Spectroscopic confirmation of high-amplitude eruptive YSOs and dipping giants from the VVV survey

During the pre-main-sequence (pre-MS) evolution stage of a star, significant amounts of stellar mass are accreted during episodic accretion events, such as multi-decade FUor-type outbursts. Here, we present a near-infrared spectroscopic follow-up study of 33 high-amplitude (most with $\Delta K_s$>4 mag) variable sources discovered by the Vista Variables in the Via Lactea (VVV) survey. Based on the spectral features, 25 sources are classified as eruptive young stellar objects (YSOs), including 15 newly identified FUors, six with long-lasting but EXor-like bursts of magnetospheric accretion and four displaying outflow-dominated spectra. By examining the photometric behaviours of eruptive YSOs, we found most FUor-type outbursts have higher amplitudes ($\Delta K_s$ and $\Delta W2$), faster eruptive timescales and bluer infrared colours than the other outburst types. In addition, we identified seven post-main sequence variables apparently associated with deep dipping events and an eruptive star with deep AlO absorption bands resembling those seen in the V838 Mon stellar merger.


INTRODUCTION
The mass accretion process plays a critical role on many astrophysical scales, from young stellar objects (YSOs) to supermassive black ★ E-mail: zhen.guo@uv.clholes, with variability as a main character (see Scaringi et al. 2015).The protostellar and pre-MS stellar evolution is dominated by the mass accretion (reviewed by Hartmann et al. 2016) with unstable nature on timescales from hours to decades (see Hillenbrand & Findeisen 2015).Stellar evolutionary models (Hartmann & Kenyon 1996) have predicted multiple accretion bursts throughout the pro-tostellar and pre-MS stage, so-called episodic accretion, suggesting that most stellar mass is accumulated during several episodes of accretion outbursts.The episodic accretion model provides a solution to the problem of large luminosity spreads among low-mass protostars (Dunham & Vorobyov 2012), as well as the luminosity problems in these objects (Kenyon et al. 1990), though the latter issue has been partially solved by a longer duration of the Class I stage of evolution (Fischer et al. 2017).Moreover, episodic accretion events can efficiently heat the inner accretion disc, which results in the expansion of the snowlines of various molecules that leave impacts on the grain growth and planet formation (e.g. Lee et al. 2019;Jørgensen et al. 2020;Kóspál et al. 2023).
Theoretical studies have explored various physical mechanisms triggering episodic accretion on pre-MS sources, such as combinations of gravitational instability (GI) and magneto-hydrodynamic turbulence at the inner accretion disc (Armitage et al. 2001;Kratter & Lodato 2016;Bourdarot et al. 2023), self-regulated thermal instability (Bell & Lin 1994), thermal instability introduced by a massive young planet inside the accretion disc (Lodato & Clarke 2004;Clarke et al. 2005), the imbalance between gravitational and magneto-rotational instability (MRI, Zhu et al. 2009b;Elbakyan et al. 2021), fragmentation of massive protostellar discs (Vorobyov & Basu 2005, 2010), infalling of piled-up disc materials outside the star-disc co-rotation radius (D'Angelo & Spruit 2010, 2012), and flybys from stellar perturber (Cuello et al. 2019;Borchert et al. 2022).Nevertheless, there is no universal solution that responds to trigger all types of episodic accretion on YSOs, and more observational evidence is desired to differentiate these mechanisms.
Large-amplitude eruptive events, often attributed to episodic accretion, are observed on YSOs via time-domain surveys from optical to infrared.Dozens of such events have been captured among Class I and II YSOs with extensive diversity in the variation amplitude and timescale (reviewed by Audard et al. 2014;Fischer et al. 2022).Traditionally, these events are sorted into two main categories, FUor-type and EXor-type, separated by distinct observational and physical characteristics.Named after FU Orionis, FUortype objects experience high-amplitude outbursts with rapid-rising ( rise < 1000 d) and long-lasting (duration some decades) morphologies (Hartmann & Kenyon 1996;Kenyon et al. 2000).Until recently, there have been only 30 FUors that have been confirmed by both photometric eruption and spectroscopic observations (see examples in Connelley & Reipurth 2018, and the list from Contreras Peña in prep), including recent discoveries via Gaia time series (e.g.Hillenbrand et al. 2018;Szegedi-Elek et al. 2020;Hillenbrand 2021).The occurrence rate of such events is about once per 10 5 years on T Tauri stars (Scholz et al. 2013;Hillenbrand & Findeisen 2015;Contreras Peña et al. 2019).During a FUor-type event, the stellar mass accretion rate reaches an extremely high value (10 −5 to 10 −4  ⊙ yr −1 ), and the viscous heated inner accretion disc becomes self-luminous and over-shines the stellar photosphere (Zhu et al. 2009a;Hartmann et al. 2011;Liu et al. 2022).Molecular absorption features, such as water and CO bands, are seen in the near-infrared spectra of FUors, arising in the cool surface of the disc.By contrast, EXor-type outbursts have typical timescales of several hundred days and are sometimes repeatable (Herbig 1989(Herbig , 2008;;Kuhn et al. 2023), with significant diversity on a case-by-case basis (e.g.Aspin et al. 2009;Lorenzetti et al. 2012;Hodapp et al. 2020;Park et al. 2022).Spectroscopic follow-ups show that most EXors still maintained the magnetospheric accretion mode during the outburst, with emission features such as Hi recombination lines and CO bandhead emission (e.g.Lorenzetti et al. 2009).
The decade-long VISTA Variables in the Via Lactea (VVV) survey and its extension VVVX survey obtained near-infrared photometry of the inner Galactic plane and the Galactic bulge from 2010 to 2020 (Minniti et al. 2010;Saito et al. 2012;Minniti 2016).Most of the VVV images are observed through the   filter with tens of visits per year during 2010 -2015 and a lower cadence after 2015.In the past years, thousands of high-amplitude near-infrared variables have been identified from the VVV   time series, including eruptive YSOs, Miras, symbiotic systems, and "blinking" post-main-sequence giants (e.g.Contreras Peña et al. 2017a;Guo et al. 2021;Smith et al. 2021;Nikzat et al. 2022).Specifically, Contreras Peña et al. (2017a) discovered that large-amplitude eruptive events (Δ  > 1 mag, named therein with the prefix "VVVv") are an order of magnitude more commonly detected among embedded protostars than more evolved disc-bearing T Tauri systems.In Contreras Peña et al. (2017b), a new sub-category of young eruptive objects was identified as MNors, which presented mixed photometric and spectroscopic features of FUors and EXors.In Guo et al. (2021) we summarised photometric and spectroscopic behaviours of 61 highly variable YSOs discovered from the VVV survey.We found that despite all FUor-type outbursts having long duration ( ≥ 5 yr), emission line objects (under the magnetically controlled accretion mode) still predominate among the long-lasting events, with   amplitudes between 2 to 4 mag.Questions have been raised on the physical origination and the frequency of these long-lasting magnetospheric controlled accretion bursts, and their impacts on the stellar evolution history.
In this series of works, Lucas et al. (submitted,hereafter LSG23) identified 222 high-amplitude (Δ  ≥ 4 mag) variable sources from the latest VVV catalogues.In this paper, we present the near-infrared spectroscopic follow-up observation of 33 sources, most of which are eruptive YSO candidates from the aforementioned catalogue.The primary goal of this work is to provide a distinctive spectroscopically confirmed sample of high-amplitude and longlasting eruptive YSOs, to examine their accretion nature during the eruptive stage, and to link the episodic accretion behaviours to their astrophysical origins.Furthermore, we aim to provide comprehensive statistical views of these long-lasting eruptive events, by combining near-to mid-infrared light curves and spectroscopic characteristics.Highlighted sources will be further studied in follow-up works.
This paper is organised as follows: basic information on target selection, the design of spectroscopic observation and data reduction, and archival photometric data are presented in §2.Spectroscopic classifications of eruptive events and other sources are described in §3.We present discussions on different types of eruptive behaviours in §4, including statistical views of eruptive timescales and variation amplitudes.An unusual source is described in §5.Finally, this paper is concluded in §6.

Target selection and basic information
In the previous work of this series, 222 high-amplitude (Δ  ≥ 4 mag) variable or transient sources were selected from a pre-release version of the VVV Infrared Astrometric Catalogue (VIRAC2; Smith et al. 2018, and in prep).In this work, we selected 31 new targets for observation from LSG23, of which 30 were regarded as eruptive YSO candidates, based on their red near-infrared colours and a VVV light curve morphology that indicated a possible outburst, or possibly an extinction-driven dip in some cases that were ambiguous.Among them, 29 sources have high variation amplitude (Δ  > 4 mag, with suffix L222_), and the other two sources have slightly lower amplitudes but also with eruptive light curve morphologies (VVV1640-4846 and VVV1636-4744).
As discussed in LSG23, there was a cluster of sources near the Galactic centre with ambiguous light curves that were thought to be either eruptive YSOs, dipping YSOs or dipping giant stars, the latter two options involving variable extinction.The 31st source, L222_149, clearly had a dip in the VVV light curve so it was observed as an exemplar for comparison.An additional criterion was that the targets were estimated to be bright enough (  < 14.5 mag), based on photometry taken in previous years, that a useful spectrum could be obtained in a short observation.
We also observed two additional sources, VVVv800 (Contreras Peña et al. 2017b) and L222_120 (=DR4_v67 Guo et al. 2021) that were previously confirmed as probable eruptive YSOs but had been observed in a relatively faint state and showed only strong outflow signatures such as H 2 and [Fe ii], rather than the spectral features associated specifically with FUor-type or EXor-type outbursts.Unfortunately, both targets were again found to be in a faint state at the time the spectra were taken so the new spectra added little to earlier observations.
Some key information about these targets is presented in Ta-ble 1, including coordinates, photometric amplitudes, and  SFR as the star-forming region (SFR) distance.As defined in Contreras Peña et al. (2017b),  SFR is the distance to the SFRs that are within 5 ′ from each target.The list of SFRs around each target is presented in the online supplementary data.Although, the spatially nearby SFR might be a projection effect,  SFR provides us with a rough estimation of the distance to the target.Based on the VVV light curves, we classified our targets into three groups: eruptive, dipper and multi-timescale variables (MTVs; see Guo et al. 2020).Specifically, light curves in the eruptive category have Δ  > 2 mag and an overall rising trend.Plus, there is no large-amplitude decaying trend before the "outburst´´.The total duration of the variation (either eruption or dip) is described by  var .The majority of our eruptive sources are either in a post-outbursting decaying stage or still on their brightness plateau, hence we only have the lower limit of  var .The  var of dippers is measured between the beginning and the ending of the dipping events.The simultaneous  −   colours measured near photometric maxima and minima (when available) are also presented in Table 1.

Spectroscopic observation and data reduction
Two spectroscopic campaigns were obtained by the XSHOOTER spectrograph on the European Southern Observatory Very Large Telescope (ESO VLT, program ID: 105.20CJ, 109.233U,PI: Lucas) in the years 2021 and 2022.During the target acquisition, blind offsetting was applied as most targets are not visible in the -band acquisition image.All targets were observed simultaneously with 0.9 ′′ slit width in optical and 0.6 ′′ slit width in near-infrared arms.
The exposure time of each target was estimated by the expected infrared brightness.Some targets were observed multiple times due to poor weather conditions.The typical spectral resolution is  = 8000 in   -bandpass.The XSHOOTER spectra were extracted using the pipeline built on the Recipe Flexible Execution Workbench software (Reflex; Freudling et al. 2013).The telluric features were further corrected by the molecfit software (Smette et al. 2015;Kausch et al. 2015) that provides a better calibration than simply applying the telluric standard provided per night.In Guo et al. (2021), a constant 83 km s −1 shift was found between the linear wavelength solution from the pipeline and the measured wavelengths of the skylines.Here, we removed this shift by setting the wavelength_frame parameter to air in molecfit and applying a 2nd-order polynomial function to generate the wavelength solution 1 .During the final examination of the spectra, some non-astrophysical peaks, likely caused by cosmic rays, were removed by replacing them with the median value of neighbouring pixels.The accurate flux calibration could not be performed, since most targets were acquired by blind offsetting and some were observed during cloudy weather conditions.The information on spectroscopic observations is presented in Table 2.The extracted spectra are presented in Figure B1 to B4. Near-infrared spectra of two sources (L222_4 and L222_15) were observed using the Folded-port Infra-Red Echellette (FIRE) spectrograph (Simcoe et al. 2013) on the Magellan Baade Telescope in 2019 (PI: Kurtev).A uniform 0.6 ′′ slit was used, providing a spectral resolution of  = 6000 in -bandpass.The data was reduced using the firehose v2.0 pipeline designed by Gagne et al. (2015).Some breaks were seen at the joint of orders, especially between the last two orders in -bandpass, caused by the mistracking of spectra towards the edge of the detector.We fixed these breaks by assuming a constant slope at the edge of the two orders (<0.005 m wide).We first applied polynomial functions to fit the curved edges between the two orders and then replaced the curved continuum with a linear slope throughout the edge area (see Guo et al. 2020Guo et al. , 2021, for more details).Finally, we used the OH skylines from the scientific exposures as wavelength calibration.

Near-infrared photometric data
The pre-release VIRAC2 version of the VIRAC2 catalogue provides near-infrared photometry that contains dozens of   band measurements and a few multi-colour epochs (Smith et al., in prep).
Additional processing of a few sources was conducted to correct for saturation or spatially resolved nebulosity.The first multi-colour epoch was taken in 2010 when most eruptive sources were at the quiescent stage, the second coloured epoch was taken in 2015 and sources in the inner bulge benefit from additional subsequent epochs.To increase the photometric accuracy, the   light curves were resampled into 1-day bins by the median brightness after removing outlier detections (3 in each bin).
In April 2021 and May 2022, two epochs of , ,   photometric follow-ups were carried out on the ESO 3.58 m New Technology Telescope (NTT), using the imaging mode of the Son OF ISAAC (SOFI) instrument.As extensions of the VVV survey, these two epochs targeted highlighted variables that were discovered from the VVV light curves (Contreras Peña et al. 2017a;Teixeira et al. 2018, and LSG23).In total, 136 targets were observed in the 2021 campaign, and 97 were observed in 2022.The photometric data were extracted by the aperture photometry methods, similar to the routine designed in Guo et al. (2018).Photometric calibration stars were selected based on the stability of their   brightness.

Archival-infrared photometric data
We obtained photometric data from WISE and Spizter surveys via the NASA/IPAC Infrared Science Archive (IRSA2 ).The WISE 1 to 4 photometry is adopted from the All-Sky data release when most sources were at the quiescent stage (Wright et al. 2010).The 1and 2-band light curves were obtained from the multi-epoch photometry table of the ALLWISE and the NEOWISE surveys (Mainzer et al. 2014(Mainzer et al. , 2014(Mainzer et al. -2022)).The 1 and 2 time series were sampled into 1-day bins, where outliers lie more than 3 from the median of the bin were excluded, and the final output light curve is the unweighted mean value of each bin.We also visually inspected the WISE images to identify sources that have suspicious detections from surrounding nebulosity, saturated nearby sources, or spatially close companions.Similarly, Spitzer 1 to 4 photometry was obtained from the Glimpse survey (Benjamin et al. 2003, prior to the VVV survey), and the 24 m data was adopted from the MIPSGAL survey (Carey et al. 2009;Gutermuth & Heyer 2015).
Spatially close companions are commonly seen among our targets.In the low spatial resolution WISE images (e.g.1 and 2 bands), sources located within 8 ′′ are not well distinguished.We set a threshold of 2 ′′ on the spatial distance when retrieving data from archival catalogues since the blending of two or more sources affects the location of archival detections when being extracted as a point source.Such blended sources are L222_4, L222_13, L222_25, L222_78, L222_148 and L222_210.In the following analysis, readers should keep in mind that the mid-infrared brightness of these targets, especially during photometric minima, is probably overestimated.The amplitudes of these sources are lower limits due to the shallower photometric minima (see the notes on Table 4).In a special case, L222_78, we performed custom-designed PSF photometry on the binned WISE images.The details of this procedure are presented in another publication (Guo et al., submitted).
We applied the slope of the infrared spectral energy distributions (SED), namely as the spectral index  class (Lada 1987;Greene et al. 1994), to identify the evolutionary stage of each source.Protostars (Class 0 & I) have positive values as  class > 0.3 while T Tauri stars (Class II & III) have negative slopes ( class < −0.3).The transitional period, also called the "flat-spectrum" stage named after the shape of their SEDs, is defined as −0.3 ≤  class ≤ 0.3.In this work,  class was calculated using the quiescent infrared SEDs (2 − 24 m, where available).Whenever possible, we preferred the SEDs from contemporary VVV   and WISE All-Sky survey.The line-of-sight extinction could lead to a higher  class as   is more sensitive to extinction than the longer wavelength bands.However, after examining the SEDs (see appendix), we found that including   has not significantly affected the result of the linear fitting, in terms of the classification of YSO evolutionary stages.The measured  class is presented in Table 1.Sources without mid-infrared detections during their quiescent stage are not measured.In total, 22 targets are classified as embedded Class I objects, and four targets are identified as "flat-spectrum" and Class II objects.In the next section, we will further refine the number of YSOs via spectral features.In the end, 16 targets are eventually confirmed as Class I sources, and two targets are "flat-spectrum" and Class II objects.

Spectroscopic classification
Based on the infrared spectral features, the 33 variables observed in this paper are classified into YSO and post-MS groups.Specifically, the YSO group is further divided into sub-categories based on their spectral features, which also reveal their physical conditions during the eruptive stage.

FUor-type
During a FUor-type outburst, the mass accretion process, often enhanced by three orders of magnitudes, is no longer controlled by the stellar magnetic field (see Audard et al. 2014, and references therein).Instead, the circumstellar material is directly accreted onto the stellar surface, and the inner accretion disc (up to a scale of 1 AU) is heated internally by the viscous dissipation of accretion energy, with a hotter mid-plane and a cooler disc surface (see Zhu et al. 2007).Numerical models show that, during a FUor-type event, the observed surface brightness of the viscous heated accretion disc could be five times greater than the stellar photosphere, dominating the infrared spectral features (Liu et al. 2022).On the observational side, FUor-type objects are diagnosed with rich absorption features in their near-infrared spectra, such as hydrogen recombination lines, the triangular-shaped -band continuum due to low gravity water vapour absorption, and CO overtone bandheads beyond 2.3 m (see Fischer et al. 2022).In the optical domain, the spectra of FUors resemble supergiants with low gravity and an early spectral type (G to K) at shorter wavelengths but a later type (M) at longer wavelengths.The infrared spectra of FUors resemble those of young brown dwarfs with low gravity, although FUors have much higher bolometric luminosity and usually show strong Pa and He i absorption features in the 1 to 1.3 m region, attributed to a wind (see Connelley & Reipurth 2018).
In this work, 15 eruptive sources are identified as FUor-type objects.We de-reddened their spectra based on the triangular-shaped Figure 1.Normalised near-infrared (1.5 to 2.4 m) XSHOOTER/VLT spectra of thirteen new FUor-type objects originally discovered by the VVV survey.In addition, spectra of two previously identified FUors (VVVv322 and VVVv721) are also presented.All spectra were de-reddened based on the shape of their  bandpass continuum, to match the spectrum of FU Orionis (provided by Connelley & Reipurth 2018).The spectra are ordered by their  −   colour after de-reddening, with the reddest on top.-band continuum, to align with the spectrum of FU Orionis (provided by Connelley & Reipurth 2018), using the extinction law from Wang & Chen (2019) and   = 1.5 mag estimated in Zhu et al. (2007).Subsequently, we estimated the   -band extinction (   ) for each FUor-type object by applying the minimum  2 method to compare the de-reddened spectra with the standard -band spectrum, with    values ranging between 0 and 10 mag.However, we acknowledge that this is only a rough estimation of    as the disc spectra of FUors are not necessarily the same from one case to another.The de-reddened  and -bandpass spectra of FUor-type sources are presented in Figure 1.
The XSHOOTER spectra of two FUors (L222_4 and L222_73) were taken during their fading stages (see B1), and show spectral features that slightly differ from those of classical FUors.Some weak and cool CO absorption features were detected, however, the triangular-shaped -band continuum is absent or less prominent due to the low brightness in , and the H 2 lines are seen in emission.In the case of L222_4, the FIRE spectrum taken in 2019 matches a standard FUor-type, where deeper absorption features were observed when the   -band brightness was higher (see Figure B5,left panel).After the original 6 mag outburst in   (the highest amplitude in this sample), L222_4 entered a rapidly fading stage that decayed 3 mag in 5 years.The spectroscopic and photometric variability of L222_4 resembles another rapidly faded FUor (VVVv322, Contreras Peña et al. 2017b).It agrees with the physical model of FUors, as a warm circumstellar disc with a steeper vertical temperature gradient at the peak of the outburst.
The radial velocity (RV) of FUors is measured by the CO bandhead absorption feature.We applied the rovibrational CO model from Contreras Peña et al. (2017b), with free parameters including the temperature and column density of the CO gas, and the RV of the system.The best-fit models were confirmed through visual inspection.The grids of the model are ΔRV = 5 km/s, smaller than the pixel size of XSHOOTER (12 km/s).The kinematic distance ( k ) is then calculated using the "Monte Carlo simulation mode" of the online tool provided by Wenger et al. (2018), based on the Galactic coordinates and RV of each source.In most cases, either the near or far  k agrees with  SFR , suggesting that the variable source is associated with nearby star-forming regions.When comparing  k with  SFR (listed in Table 1), many FUors have  SFR beyond the 1 error of  k .However,  k could be wrong if the target is in a binary system due to the orbital motions or if it has an intrinsic velocity, but our data are not enough to confirm or reject the existence of any close companions.

Emission line and outflow-dominated objects
Some YSOs exhibiting long-lasting eruptive events have spectroscopic features disagreeing with FUors.One highlighted source is V1647 Ori, whose eruptive light curve resembles FUors.However, a local dimming event with a timescale of 1000 days was seen after the original outburst.The spectra of V1647 Ori have water vapour absorption accompanied by weak CO absorption in the first outburst but with Br and CO bandheads emission in the second outburst (see Aspin et al. 2009;Connelley & Reipurth 2018).Since then, eruptive sources with features between FUors and EXors were named MNors or V1647 Ori-type (Contreras Peña et al. 2017b;Hodapp et al. 2020;Fischer et al. 2022).In the near-infrared, a group of eruptive YSOs was discovered by the VVV survey with longlasting outbursts but EXor-like spectra (e.g.Na I doublet, Br and CO bandheads emission), which are named as emission line objects in Guo et al. (2021).Another group of young variables are classified as outflow-dominated objects since only  2 and/or [Fe II] emission lines were identified from the spectra (see Section 3 and Table 5 in Guo et al. 2021).In the VVV survey, these novel categories of embedded outbursting YSOs are more populated than FUors, even towards the longest duration end.
In this paper, six sources are identified as emission line objects with signatures of magnetospheric accretion (see Figure B3).Among them, L222_148 and L222_167 were observed during their brightness plateau, whilst the rests were observed at least 1.5 mag fainter than their photometric maxima.CO bandheads were detected (beyond 3) on L222_148 and L222_167 in emission, presenting a positive correlation between the infrared brightness and the existence of hot gas at the accretion disc, aligned with a high massaccretion rate.One source, L222_167 was observed on two epochs 40 days apart and the continuum on the first epoch is 13% brighter than the second epoch (Figure 2).Stronger mass-accretion indicators are seen in the brighter epoch, whilst the fainter epoch has stronger indicators of the outflow ( 2 lines), suggesting non-correlation between the mass-accretion rate and the strength of outflow or there is a time delay between them.Typically, emission from a jet is rarely detected among FUor-type outbursts, as the strong accretion disrupts the magnetic field of the inner disc, preventing the launching of a jet.Complex line profiles are detected on the  2 lines, including a blue-shifted high-velocity component around -130 km/s, a Gaussian component around -33 km/s, and a weak red-shifted tail.This complex line feature indicates that winds/outflows were launched towards different directions, and sometimes at perpendicular angles.
Four targets are classified as outflow-dominated objects, only having spectral indicators of stellar wind/outflows.Among them, the spectra of two YSOs (VVVv800 and L222_120, a.k.a.DR4_v67) were previously published in Guo et al. (2020Guo et al. ( , 2021)).In comparison with the latest spectra, there is no significant spectroscopic variability.In the case of L222_32, two XSHOOTER spectra were obtained in 40 days.Despite very similar continuum levels, the  2 S(1) line in the second epoch is 35% stronger in flux and 11 km/s redder in RV than the first epoch.An interpretation is that the variability of the outflow is not correlated with or delayed response to the variability of the infrared continuum.Combined with the list published in Guo et al. (2021), 12 VVV targets are classified as outflow-dominated sources.In most cases, the VVV light curves of outflow-dominated sources are irregular, with significant variations on (sub)year-long timescale, in contrast to classic outbursts (see definitions in LSG23).This could be a geometric effect as one is looking through the outflow, as part of the variability might be attributed to the winds/outflows lifting dust from the disc, hence affecting the extinction towards the source.
Wherever available, we measured the extinction using the flux ratio between  2 1-0 Q(3) (2.42 m) and 1-0 S(1) (2.12 m) lines, as it is independent of the excitation conditions (Q(3)/S(1) = 0.7; Turner et al. 1977).We applied the power-law extinction function from Wang & Chen (2019).The typical precision of the    measurement is 0.4 mag, with uncertainty rising from the flux measurement of the Q(3) line ( = 8%), which is veiled by the continuum and the CO overtone emission.The colour variation before and after the outburst is discussed further in §4.3.

Post-main-sequence objects
Many post-MS variables in the sample have VVV light curve morphologies that appear either to show a dip or an eruption: this is ambiguous in some cases (see examples in Contreras Peña et al. 2017b;Guo et al. 2021).Spectra were obtained to determine which of the two is occurring or whether these sources are YSOs or giant stars.However, most such sources with ambiguous light curves were detected in a bright state by 2MASS, e.g.sources L222_144, L222_145, L222_154, L222_144, L222_168, L222_172, which tended to indicate that we are observing the dips due to extinction rather than photometric eruptions (see LSG23).These sources often have a red colour in the near-infrared ( −   > 3.5 mag) or lack -band detections in VVV data.Mid-infrared colour-colour and colour-magnitude diagrams were applied to distinguish YSOs with post-MS sources, although no rigid boundaries could be drawn (see Lucas et al. 2017;Guo et al. 2022).Due to nucleosynthesis processes and the first dredge-up, most post-MS stars display a fairly strong 13 CO absorption bandhead at 2.345 m in our intermediate resolution spectra, whereas this feature is rarely seen in the YSO spectra presented herein.
In the LSG23 catalogue, 21 sources are classified as dipping giants based on their infrared light curves, with a majority of them projected in the Nuclear Disc of the Milky Way.We obtained spectroscopic verification of eight sources (see Figure B4), including seven dipping giants and one eruptive source (L222_59).All dipping giants have similar spectroscopic features resembling O-rich post-MS stars.The quiescent colour-colour diagrams of our targets are shown in Figure 3.Most post-MS sources exhibit redder  −   colours compared to YSOs.Furthermore, it is noteworthy that FUors consistently display smaller colour indices than non-FUor eruptive YSOs, which will be explored further in §4.3.We derived the stellar parameters for post-MS sources by fitting the CO bandhead models plus equivalent width (EW CO and EW Na ) based on methodologies derived from Feldmeier-Krause et al. ( 2017) and Fritz et al. (2021).The outcomes of this analysis are listed in Table 3.Most dipping giant sources have an effective temperature between 3000 to 4000 K and super solar metallicity (except L222_168).This high metallicity is consistent with the location in the Nuclear Disc, for which the metallicity distribution is distinctly metal-rich (Fritz et al. 2021).More details are presented in the Appendix D.
We highlight an eruptive source (L222_59, a.k.a.VVVv746) that exhibits rarely seen deep aluminium monoxide (AlO) absorption bands (see Figure 4 and B4).This source experienced a Δ  > 5 mag outburst and reached the photometric maximum in less than three years.According to our follow-up   photometry, after spending nearly 1000 days on its brightness plateau, L222_59 dramatically faded nearly 3 mag in the year 2021-2022.The XSHOOTER spectrum presented in this paper was taken during the outbursting stage, with many absorption features, including cool 12 CO overtone bands in the  bandpass (as fitted temperature,  = 1200 K, assuming LTE) and AlO in the  and -bandpasses, confirmed by the ExoMol database (Patrascu et al. 2015;Tennyson et al. 2016;Bowesman et al. 2021) to be the A-X electronic transitions of the AlO radical.Interestingly, it also has narrow Hi lines in absorption, with different RV compared to the CO bandhead.These spectral features indicate the complexity of a cool envelope and a warm central star.Hitherto, such deep AlO absorption bands have only been observed in the stellar merger candidate V838 Mon (Evans et al. 2003;Banerjee et al. 2005), which exhibits a Nova-like eruption on a timescale of 50 days (Munari et al. 2002).The hypothesis of a candidate stellar merger will be further discussed in a follow-up work, including up-to-date photometric and spectroscopic data.

DISCUSSION
In the following sections, we present analyses based on the photometric and spectroscopic behaviours of eruptive YSOs discovered in the VVV survey.Individual characters, such as the rising timescales of the outburst, and the near to mid-infrared amplitudes, are measured to distinguish FUors from other eruptive events.

Eruptive samples
We selected a sample of 32 long-lasting eruptive events discovered by the VVV survey, with near-infrared spectroscopic confirmation by our group.To minimise the selection bias, we included all sources with Δ  ≥ 2 mag and have eruptive events lasting longer than 1000 days.Among them, 18 sources are spectroscopically confirmed as FUors and the rest are either emission line  objects or outflow-dominated sources (non-FUors).This selection of confirmed long-lasting eruptive YSOs is roughly consistent with the sample discussed in LSG23, although the samples in this work include three emission line objects with light curve morphologies different from the "classic outbursts".This is to maintain enough non-FUors in our sample.However, the results of the measured distribution of several characters (between FUors and non-FUors) are not severely impacted if these three sources were not included.
According to their quiescent SEDs and estimated distance, the vast majority of these sources are embedded low-mass Class I protostars.
Our sample forms a valuable laboratory to investigate outbursting events on embedded sources, especially to explore the relationship between photometric and spectroscopic features.The   light curves of these sources, including data from both the VVV survey and SOFI observations, are shown in Figure C1.We note that three FUors previously identified in the VVV survey are not included in this sample, due to the absence of   detections during the rising stage, hence their eruptive stages are poorly described.The complexity of our eruptive sample plays a key role in our following statistical analysis.A few factors could affect the detection of long-duration outbursts, such as saturation in the VIRAC catalogue during the outbursting stage, eruption prior to the beginning of the VVV survey (the year 2010), and lower amplitude events that are yet to be confirmed by spectroscopic follow-ups.In another paper of this series (Contreras Peña et al., submitted), we will further discuss the occurrence rate of eruptive behaviours on Class I protostars, by applying mid-infrared SED-based catalogues (e.g.SPICY, Kuhn et al. 2021) and including archival   -band detections (e.g.2MASS) to form a two-decade-long timeline.

Near-infrared light curves of eruptive samples
Eruptive events on YSOs are known to have various photometric characteristics, due to different triggering mechanisms and disc configurations (Vorobyov et al. 2021).For instance, FU Ori reached the peak brightness in eight months while the rising stage of V1057 Cyg lasted 15 years (Herbig 1977).Here, we apply linear functions to describe the quiescent, rising, decaying, and post-outbursting stages of the aforementioned 32 eruptive events (see examples Figure 5).Two parameters were defined to quantify the rising stage: the rising amplitude (Δ s,rise ) and the rising timescale ( rise ).Specifically, Δ s,rise is the   amplitude of the rising stage, measured from the linear fitting result.The  rise is calculated as Δ s,rise divided by the measured linear slope (mag/d) of the rising light curve.The measured  rise and Δ s,rise are presented in the right panel of Figure 5, with numbers listed in Table 4.
We present the statistical differences of  rise and Δ s,rise between FUors and non-FUors (emission line/outflow-dominated ob-Figure 6.The ratio between the short-term residual amplitude and the rising amplitude versus the rising timescale of eruptive sources.Colour codes and lower limits are as same as in Figure 5. jects).In general, we find that FUors experienced shorter  rise compared to non-FUors.Among the examined samples, 14/18 FUors have  rise < 850 days (mean at 747 d), whilst all non-FUors have  rise > 800 days (mean at 1038 d).To support this finding, we employed the Wilcoxon-Mann-Whitney test (Mann & Whitney 1947), a nonparametric rank-based test.The possibility of FUor and non-FUor sharing the same distribution in  rise is only  WMW = 0.001.There should not be any selection bias on the  rise among our sample, as we examined almost all high-amplitude eruptive candidates, however, the sparse cadence of the VVV survey (after 2015) limits our capability to measure any timescale less than 1 year.Among the longestduration objects, we detect eruptive sources with  rise > 2000 d, including two sources that entered their outbursting stage before 2010.We find that FUors and non-FUors share a comparable range of Δ s,rise , although FUors are more abundant towards the highest amplitude end, as 61% of FUors and 21% of non-FUors have Δ s,rise > 3.8 mag.However, FUors and non-FUors are not distinguishable at the intermediate amplitude range (2.0 < Δ s,rise ≤ 3.8 mag).The  WMW on the Δ s,rise distribution is 0.07, larger than the one measured on  rise .Since there is no strongly justified null hypothesis regarding the relative amplitudes of the FUors and non-FUors, we regard this  WMW -value as fairly good evidence that the amplitude distributions are different.
We further compare the  rise between VVV and FUors discovered from optical time series, including FU Ori, V1057 Cyg, V960 Mon, V2775 Ori, Gaia17bpi, Gaia18dvy, and Gaia21bty (Kenyon et al. 2000;Connelley & Reipurth 2018;Hillenbrand et al. 2018;Szegedi-Elek et al. 2020;Siwak et al. 2023).Despite three fasteruptive sources with  rise < 300 days, other archived FUors share similar  rise with the FUors discovered in the infrared (see more discussion on this topic in LSG23).However, some theoretical works predicted that outbursts triggered by outside-in migrated instabilities should have shorter  rise in the optical than in the infrared, in the order of years (Cleaver et al. 2023).Intriguingly, our VVV FUors have two outliers, one has a rapid rise followed by a rapid decay (VVVv322), and the other (VVVv721) experienced a long rising stage as  rise > 2140 days.The long rise time of VVVv721 indicates that the instability originated further out from the accretion disc or else the viscosity parameter was smaller (see discussions in Liu et al. 2022).
Some eruptive sources exhibit long-term variability during their quiescent stage, such as L222_1, L222_10, L222_25 and VVV1640-4846.All of these sources have a gradual increase in   before the main outburst, which can be interpreted as the piling up of warm disc material.This interpretation aligns with the first stage of the MRI-GI model, in which disc material is first gathered outside the magnetic dead zone through GI (Bourdarot et al. 2023).We have seen different cooling slopes among VVV FUors, including a few rapidly decaying sources (e.g.VVVv322, L222_4, VVV1640-4846) that differ from FU Ori, and even faster than the fading stage of V1057 Cyg.Apart from VVVv322, we have not observed other FUors returning completely to the quiescent brightness or the re-establishment of magnetospheric accretion.
Short-term variations are detected in eruptive YSOs with timescales from hours to hundreds of days, in addition to their overall eruptive behaviour, on both FUors (e.g.V1037 Cyg Szabó et al. 2021) and non-FUors (e.g.V899 Mon Ninan et al. 2015).We subtracted the long-term linear trends in the lightcurves to measure the short-timescale residual variability (Δ s,res ), as the residual amplitude within a time window of two years in each measurement (see Figure 6).In general, FUors have smaller Δ s,res than non-FUors, consistent with their physical nature, as an abrupt and persistent enhancement of the mass accretion rate which efficiently heated the inner accretion disc.Short-timescale extinction and accretion variabilities are not commonly seen among FUors in our sample.On the contrary, some local dipping events are observed on non-FUors, such as the quiescent stage of L222_6, and the risings stage of Stim 5, DR4_v10, DR4_v34 and DR4_v67.These features could be attributed to extinction events (like V582 Aur Ábrahám et al. 2018) or the discontinued accretion stream.We infer that the reduction of line-of-sight extinction may in fact contribute to some high-amplitude (Δ  ) variation seen among non-FUors.

Mid-infrared amplitude and colour of eruptive samples
Many eruptive YSOs show large amplitudes across a wide range of wavelengths, from optical to mid-infrared (e.g.Kóspál et al. 2016;Lucas et al. 2020;Contreras Peña et al. 2023).Guo et al. (2021) summarised that eruptive YSOs have comparable amplitudes in the near-and mid-infrared, in contrast to the expectations from variable extinction.In this section, we will compare the contemporary amplitudes and colours of these eruptive events, from   to 2, to further investigate the physical nature.Since the VVV and NEOWISE time series have different cadences, we measured the contemporary amplitudes via linear interpolation, based on the photometric maximum identified from 2-band light curves.For sources without WISE detections during the quiescent stage, we converted the fluxes measured in Spitzer filters to WISE filters using the empirical photometric transformations from Antoniucci et al. (2014).In two cases (L222_10 and L222_25), there are no archival quiescent mid-infrared detections, hence the Δ2 and   − 2 can not be measured.The results are listed in Table 4.
The correlation between Δ2 and Δ  is shown in Figure 7. FUors exhibit systematically higher Δ2 compared to non-FUors, indicating efficient heating of the circumstellar disc during the outburst.This aligns with the spectral feature of FUors, where their optical to mid-infrared emission is primarily from a self-luminous inner accretion disc.By contrast, there are no emission line objects with Δ2 > 2 mag, even among those with Δ  > 3.5 mag.The infrared colour indices are presented in Figure 7 and Figure 8. Lower  limits are put on objects that erupted before the VVV survey.All FUors have   − 2 < 6 mag at the quiescent stage, whilst the majority of long-term non-FUors have a redder   −2 colour.Similar trends are seen in the 1 −2 colour.During the outbursting stage, most non-FUors are still redder than most FUors although the gap is smaller.This colour difference is consistent with the  −   colour distribution shown in Figure 3, where most FUors are bluer than non-FUors during the quiescent phase.As comparisons, we added some archived eruptive YSOs in Figure 8, including three FUors and three recently confirmed eruptive objects from the NEOWISE data (with prefix SPICY, Contreras Peña et al., submitted).
Due to the logarithmic nature of the magnitude measurement, the amplitude is smaller when the source is already bright in a particular band.This aligns with the low Δ2 seen among red non-FUors (  − 2 > 6 mag).However, the origin of their substantial excess in 2 is still uncertain.These objects could be heavily obscured during the quiescent stage.The high line-of-sight extinction may arise from a dusty envelope, an edge-on circumstellar disc, or even the foreground cloud.However, the disc inclination and random foreground extinction cannot account for the fact that most non-FUors exhibit a redder colour than most FUors (see distribution in Figure 8).Nevertheless, the extinction measured from in-outburst spectra between FUors and non-FUors is comparable (see Table 2).The reduction of extinction during the outburst might play a role here, where the stellar outflow clears the immediate vicinity around the star, and results in a greater Δ  than Δ2.We also admit that there is a selection bias in our spectroscopic sample (i.e.only sources with high Δ  were selected from LSG23).Alternatively, a red colour index typically indicates an earlier stellar evolution stage, suggesting that most long-lasting emission line objects could be Class I systems characterized by massive accretion discs and envelopes.
In summary, we find FUors have bluer   − 2 and 1 − 2 colours than most non-FUors during the quiescent phase.However, we can not simply draw a solid conclusion that there are no heavily embedded FUors or that all FUors are slightly more evolved than other eruptive sources, due to the limited size of samples and contribution from variable extinction.In a forthcoming work, Morris et al, (in prep) will obtain spectroscopic confirmation of a group of extremely red candidate FUors discovered from the mid-infrared NEOWISE time series.as the quiescent stage, rising stage and the brightness plateau (see Figure 9).A 378-day quasi-periodic signal is detected with Δ  ∼ 2.0 mag after applying the generalised Lomb-Scargle periodogram (Zechmeister & Kürster 2009) to the residual light curve.Since the photometric minima are poorly covered, the M-index is not accurately calculated, hence we did not recognise this target as a periodic outbursting candidate (Guo et al. 2022, see the definition of M-index herein).The hybrid photometric behaviour of L222_148 resembles the eruptive source LkH 225 South (Hillenbrand et al. 2022), with a period of 43 days and an overall 7 mag eruption in optical bands.The photometric period of L222_148 is likely attributed to the Keplerian rotation of an asymmetric structure in the circumstellar disc, which may be triggered by a secondary body in the system, similar to the 130 candidate periodically outbursting YSOs discovered in Guo et al. (2022).
The near-infrared spectrum of L222_148, taken during the outbursting stage, shows emission features associated with the magnetospheric accretion process (e.g.Naii and CO bandheads) and coupled with strong stellar outflow (H 2 line series).In contrast, the Br line is in absorption.Inverse P Cygni profiles are seen on the H 2 lines, suggesting a combination of blue-shifted stellar outflow coupled with a stream of cold disc materials (see the right panels of Figure 9).The RV of the absorption components are consistent with the CO bandhead emission (∼ −10 km/s).In a recent spectrum of this target, we do not detect any absorption features around H 2 lines, and the Br line (previously in absorption) is in emission.We infer that the changing of the spectral features results from the evolution of the inner disc structure, which might be triggered by gravitational instabilities or perturbation from a secondary object.The new spectrum will be published in a follow-up work.
This unique target, L222_148, provides us with a complex eruptive phenomenon, as a slow-rising accretion burst happens on a possible stellar binary or a star-brown dwarf system.We have detected the hot inner accretion disc, a cold material stream, and stellar wind/outflow from the near-infrared spectrum.It is an intriguing case to explore the physics corresponding to variable mass accretion processes on timescales from years to a decade.In addition, we can not rule out the possibility that this object is in a stellar binary (explains the two spectral components), in which one source exhibits periodic variation and the other source undergoes a slow-rising outburst, like the case of Z CMa (Bonnefoy et al. 2017).

SUMMARY
In our series of works, over two hundred high-amplitude variable stars were identified from the decade-long VVV/VVVX   time series (see LSG23).Among them, FUor-type outbursts are seen among the highest amplitude and longest duration events.In this work, we have presented near-infrared spectroscopic follow-up observations of 33 sources using the XSHOOTER spectrograph on VLT and FIRE on the Magellan telescope.These sources are categorised into a few groups based on their spectral features and the nature of their variability.We summarise our findings as follows.
• Fifteen new FUor-type objects are confirmed by unique spectral signatures, such as CO absorption bands beyond 2.3 m and the triangular-shaped -band continuum.
• Four targets are classified as outflow-dominated with  2 emission line series, and seven targets are identified as emission line objects with signatures of magnetospheric accretion.Most of these sources have a longer duration than the optically seen EXors.
• Eight sources, including seven dipping giants, are identified as post-main-sequence sources with the existence of 13 CO absorption features and extremely red near-infrared colours (five of them have  −   > 3.5).The rarely seen deep near-infrared AlO absorption bands are detected on an eruptive source, L222_59.
• We form a sample of 32 spectroscopically confirmed longduration eruptive YSOs originally identified from the VVV sur- vey.We find that in the near-infrared, FUors have faster eruptive timescales and predominate the highest amplitude end (Δ  > 3.8).In the mid-infrared, FUors have systematically higher 2 amplitude than other eruptive objects.It suggests that FUors can heat their inner accretion disc more efficiently.
• FUors have comparable near-infrared and mid-infrared ampli-tudes.By contrast, most non-FUors have much lower amplitudes in the mid-infrared than in the near-infrared, with redder   − 2 and 1 − 2 colours, indicating early evolutionary stages or high line-of-sight extinction.
• In Guo et al. (2021), we found that most multi-year duration eruptive events on YSOs are emission line objects.However, in this work, we notice that the distribution is reversed towards the highest amplitude end, as FUors are more abundant among the eruptive samples with Δ  > 4 mag.
In summary, this work significantly increased the total number of spectroscopically known FUor-type events.Our photometric and spectroscopic surveys provide several key characteristics of the eruptive timescales and amplitudes of eruptive events in two sub-categories.To further reveal the physical nature behind these eruptive events, high-cadence multi-wavelength observations and numerical simulations are desired.

Figure 2 .
Figure 2. Left:  bandpass XSHOOTER spectra of L222_167 with spectral features marked out individually.Two spectral epochs are presented in different colours.The flux of the second epoch is shifted by 1 unit smaller for better visualisation.Right: Normalised line profiles of H 2 S(1) lines at 2.12 .Two spectral epochs are shown by the black and blue solid lines.Three Gaussian components, presented as the dotted lines, are fit to the line profile of the second epoch.

Figure 3 .
Figure 3.  −   and   − 2 colour-colour diagram of targets observed in this work (see name tags on the plot, the suffix "L222_" is omitted), colour-coded by spectral features.Targets with classic eruptive   light curves are shown by dots and others are shown by triangles.Filled dots represent YSOs with colours measured during the pre-outbursting stage.Sources without any W2-band detections are not presented.Four FUors that only have 2 detections during the outbursting stage are presented by open circles.Typical photometric error bars (<0.05 mag) are about the same size as the data points.

Figure 4 .
Figure 4.The -bandpass spectrum of L222_59 with an analytical model of the AlO absorption bands from ExoMol (red lines; Bowesman et al. 2021).

Figure 5 .
Figure 5. Left:   light curves of four eruptive YSOs observed in this work.Linear functions are applied to fit the quiescent/pre-eruptive (blue), rising/eruptive (green), and post-outbursting (red) stages.Right: The rising timescales and amplitudes of large amplitude and long-lasting eruptive YSOs.Archival FUors are shown as comparisons, among which the amplitudes of three FUors were measured in optical (open stars) and the rest were measured in   (filled stars).Outbursts started in the pre-VVV era were labelled with lower limits.The distributions of rising timescales and amplitudes of VVV sources (Archival sources are not included) are presented in the top and right panels.The possibility from the Wilcoxon-Mann-Whitney test (    ) for each distribution is labelled on the histogram.

Figure 7 .
Figure 7. Left: Contemporary near-infrared (  ) and mid-infrared (2) amplitudes of long-term eruptive YSOs.Sources without quiescent detections are marked as lower limits.Photometric uncertainties are smaller than the symbols.Right:   − 2 colour versus the mid-infrared amplitudes.The quiescent and in-outburst colours are presented by filled and open dots, respectively.

Figure 8 .
Figure 8. Quiescent 1-2 colour and Δ2 of eruptive YSOs.Sources without WISE detections during their quiescent stage are shown with lower limits.Archived eruptive YSOs are presented as open dots and triangles.The typical uncertainties on the colour and amplitude are 0.15 mag.Distributions of 1-2 and Δ2 are presented as histograms.

Figure 9 .
Figure 9. Upper Left:   light curve of L222_148, colour coded by the observation time.Overall linear trends are shown by the dashed lines.Middle Left: residual light curves after removing the linear trends.The locations of periodic photometric maxima are marked by dotted vertical lines.Lower Left: phase-folded residual light curves.Right: profiles of  2 emission lines and the Br line.The radial velocity of CO bandhead emission is marked by the black dashed line.The central wavelengths of the emission/absorption features are shown by the dotted lines.

Figure A1 .
Figure A1.Near to mid-infrared SEDs of 28 out of 33 sources observed in this work.Data obtained during the photometric minima are shown by dots while data obtained during photometric maxima are presented by triangles.The data point is colour-coded by the origin, as VVV or NTT/SOFI: blue; WISE: brown; Spitzer: red.Linear slopes are fit to the quiescent SEDs at  > 2m, among sources having detections beyond 5 m.MNRAS 000, 000-000 (2023)

Figure B1 .
Figure B1.XSHOOTER spectra and   light curves of FUor-type objects confirmed in this paper.The detections from NTT/SOFI are marked in red.The observation time of each spectrum is presented as the dashed line in the lower panels.

Figure B2 .
Figure B2.XSHOOTER spectra and   light curves of emission line objects.Symbols are the same as in Figure B1.

Figure B3 .Figure B4 .
Figure B3.XSHOOTER spectra and   light curves of outflow dominated objects.Symbols are the same as in Figure B1.Bottom Right: Line profiles of three  2 lines from two spectral epochs of L222_32.The vertical dashed lines show the radial velocities of Gaussian fittings.

Figure B5 .
Figure B5.Magellan/FIRE spectra of two FUors taken in 2019.Symbols are the same as in Figure B1.These two spectra were originally published by R. Kurtev as a poster at the Crete III conference in 2019 without permanent records.

Figure C1 .
Figure C1.  light curves of long-duration eruptive YSOs discovered from the VVV survey.Linear functions are applied to fit different stages on the light curve.The eruptive stage, both the light curve and the best fit, is highlighted in green colour.

Table 1 .
Basic information and photometric behaviours of our targets

Table 2 .
Spectroscopic observations of YSOs : kinematic distance calculated using the online tool provided inWenger et al. (2018).

Table 3 .
Peña et al. 2017b)rmation of post-MS sources eff is measured using the methods from Feldmeier-Krause et al. (2017), with a typical uncertainty of 163 K.  CO is measured from the analytical models of CO bandheads (ContrerasPeña et al. 2017b), with an error bar of 200 K. ★ : the measurements of  eff and [Fe/H] are not reliable due to the unique spectral features of L222_59, e.g.circumstellar CO.

Table 4 .
Variation amplitudes and pre-outbursting colour of eruptive YSOs Name Type  rise Δ s,rise Δ s,res Δ2   − 2 Eruptive objects without clear detections from the WISE images due to blending with a spatially nearby source.† : Sources without 2 detections during the photometric minima, but have I2-band detections from Spitzer.‡ : Eruptive events started their rising stage before the VVV survey.
WISE and Spitzer data underlying this article are publicly available at the IRSA server https: //irsa.ipac.caltech.edu/Missions/wise.html,https: //irsa.ipac.caltech.edu/Missions/spitzer.html,and https://irsa.ipac.caltech.edu/Missions/2mass.html.The VVV and VVVX data are publicly available at the ESO archive http://archive.eso.org/cms.html.The VIRAC2 version of the VVV/VVVX light curves has not yet been publicly released but is available on request to the first author.Raw spectra are available on the ESO archive service when released to the public.Reduced spectra are provided at http://star.herts.ac.uk/~pwl/Lucas/GuoZ/VVVspec/.and PRIN 2022 20228JPA3A -PATH.SY and JT thank the support of the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme through Advance Grant number 883830 and of STFC under the project ST/R000476/1.AA acknowledge the support through a Fellowship for National Ph.D. students from ANID, grant number 21212094.CCP was supported by the National Research Foundation of Korea (NRF) grant funded by the Korean government (MEST) (No. 2019R1A6A1A10073437) Support for M.C. is provided by ANID's Millennium Science Initiative through grant ICN12_009, awarded to the Millennium Institute of Astrophysics (MAS); by ANID/FONDECYT Regular grant 1231637; and by ANID's Basal grant FB210003.K. M. is funded by the European Union (ERC, WANDA, 101039452).Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Research Council Executive Agency.Neither the European Union nor the granting authority can be held responsible for them.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.