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

ABSTRACT

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 multidecade FUor-type outbursts. Here, we present a near-infrared spectroscopic follow-up study of 33 high-amplitude (most with Δ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 (ΔK_s and ΔW2), faster eruptive time-scales and bluer infrared colours than the other outburst types. In addition, we identified seven post-MS variables apparently associated with deep dipping events and an eruptive star with deep aluminium monoxide absorption bands resembling those seen in the V838 Mon stellar merger.

stars: AGB and post-AGB, stars: pre-main-sequence, stars: protostar, stars: variables: T Tauri, Herbig Ae/Be, infrared: stars

1 INTRODUCTION

The mass accretion process plays a critical role on many astrophysical scales, from young stellar objects (YSOs) to supermassive black holes, with variability as a main character (see Scaringi et al. 2015). The protostellar and pre-main-sequence (MS) stellar evolution is dominated by the mass accretion (reviewed by Hartmann, Herczeg & Calvet 2016) with unstable nature on time-scales from hours to decades (see Hillenbrand & Findeisen 2015). Stellar evolutionary models (Hartmann & Kenyon 1996) have predicted multiple accretion bursts throughout the protostellar 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, Belloche & Garrod 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 magnetohydrodynamic turbulence at the inner accretion disc (Armitage, Livio & Pringle 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 magnetorotational instability (MRI; Zhu, Hartmann & Gammie 2009b; Elbakyan et al. 2021), fragmentation of massive protostellar discs (Vorobyov & Basu 2005, 2010), infalling of piled-up disc materials outside the star–disc corotation 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 time-scale (reviewed by Audard et al. 2014; Fischer et al. 2023). 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, FUor-type objects experience high-amplitude outbursts with rapid-rising (t_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 e.g. in Connelley & Reipurth 2018, and the list from Contreras Peña in preparation), 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⁵ yr on T Tauri stars (Scholz, Froebrich & Wood 2013; Hillenbrand & Findeisen 2015; Contreras Peña, Naylor & Morrell 2019). During an FUor-type event, the stellar mass accretion rate reaches an extremely high value (10⁻⁵–10⁻⁴ M_⊙ yr⁻¹), and the viscous heated inner accretion disc becomes self-luminous and overshines the stellar photosphere (Zhu et al. 2009a; Hartmann, Zhu & Calvet 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 time-scales of several hundred days and are sometimes repeatable (Herbig 1989, 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 H i 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 K_s 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 K_s time-series, including eruptive YSOs, Miras, symbiotic systems, and ‘blinking’ post-main-sequence (MS) 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 (ΔK_s > 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 subcategory 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 summarized 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 (t ≥ 5 yr), emission-line objects (under the magnetically controlled accretion mode) still predominate among the long-lasting events, with K_s 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. (2024, hereafter LSG23) identified 222 high-amplitude (ΔK_s ≥ 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 long-lasting 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 organized as follows: basic information on target selection, the design of spectroscopic observation and data reduction, and archival photometric data are presented in Section 2. Spectroscopic classifications of eruptive events and other sources are described in Section 3. We present discussions on different types of eruptive behaviours in Section 4, including statistical views of eruptive time-scales and variation amplitudes. An unusual source is described in Section 5. Finally, this paper is concluded in Section 6.

2 OBSERVATION AND DATA REDUCTION

2.1 Target selection and basic information

In the previous work of this series, 222 high-amplitude (ΔK_s ≥ 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 preparation). 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 (ΔK_s > 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 (K_s < 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₂ 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 Table 1, including coordinates, photometric amplitudes, and d_SFR as the star-forming region (SFR) distance. As defined in Contreras Peña et al. (2017b), d_SFR is the distance to the SFRs that are within 5 arcmin from each target. The list of SFRs around each target is presented in the online Supporting Information. Although, the spatially nearby SFR might be a projection effect, d_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 multitime-scale variables (MTVs; see Guo et al. 2020). Specifically, light curves in the eruptive category have ΔK_s > 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 t_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 t_var. The t_var of dippers is measured between the beginning and the ending of the dipping events. The simultaneous H − K_s colours measured near photometric maxima and minima (when available) are also presented in Table 1.

Table 1.

Basic information and photometric behaviours of our targets.

Name	RA (J2000)	Dec. (J2000)	H − Ks	H − Ks	LC type	ΔKs	tvar	αclass	dSFR	W1 − W2
	deg	deg	faint stage	bright stage		mag	d		kpc	faint stage
L222_1	175.78945	−62.35367	1.49	1.18	Eruptive	3.9	>2451	0.82	9.6	2.0
L222_4	185.22517	−62.63942	1.09	2.44	Eruptive	6.2	>2651	0.59	7.6	0.9
L222_6	193.73921	−61.04417	–	1.53	Eruptive	4.0	>2251	2.07	3.2	2.2
L222_10	201.45950	−62.79639	1.42	1.34	Eruptive	4.0	>2953	−0.58	–	–
L222_13	203.04039	−62.73005	0.83	0.80	Eruptive	4.1	>1921	–	3.0	–
L222_15	204.54729	−62.48267	2.92	0.93	Eruptive	4.8	>2951	–	2.0	–
L222_18	214.07488	−61.37306	1.97	1.47	Eruptive	4.0	>3051	−0.03	–	0.2
L222_25	230.39117	−57.88889	−0.09	1.10	Eruptive	3.3	>2951	–	3.1	–
L222_28	232.57471	−55.58203	2.79	1.48	Dipper	4.5	2548	0.09	–	1.4
L222_32	239.45987	−53.95964	–	1.97	Eruptive	4.2	>2151	1.29	2.8	2.6
L222_33	239.85950	−51.95328	2.15	1.58	Eruptive	4.2	>3601	0.39⋆	3.4	0.3
L222_37	241.77933	−49.40261	–	1.66	Eruptive	>4.3	>4151	1.36	–	3.1
VVV1636−4744	249.15808	−47.74558	2.27	2.03	Eruptive (MTV)	3.7	∼300	0.96⋆	4.5	–
VVV1640−4846	250.04900	−48.78150	2.33	1.52	Eruptive	3.2	>1627	0.96⋆	2.1	0.6
L222_59	253.43487	−43.47205	1.81	1.55	Eruptive	5.7	>2251	0.89	–	3.1
VVVv800	258.19183	−38.42350	–	1.88	Eruptive	3.2	1851	0.84	1.4	2.3
L222_73	258.65950	−38.49140	–	3.02	Eruptive	>3.9	>4151	1.39	6.2	0.7
L222_78	259.58187	−32.38142	0.33	0.46	Eruptive	4.6	>2642	–	1.1	0.4
L222_93	261.73025	−34.14661	1.37	1.18	Eruptive	4.5	>3151	–	6.5	–
L222_95	262.28616	−33.52969	1.27	0.43	Eruptive	5.3	>2751	1.41	6.6	0.6
L222_120†	265.37971	−31.43672	3.50	1.79	Eruptive	5.3	> 2463	0.82	4.5	2.5
L222_130	266.11000	−29.45197	3.91	3.62	Dipper?	>4.2	>2451	–	–	–
L222_144	266.53671	−28.68197	–	3.56	Dipper	4.2	>1450	2.08⋆	–	–
L222_145	266.58596	−28.56781	3.83	2.91	Dipper	4.7	>2551	0.85⋆	–	–
L222_148	266.64095	−29.37922	3.30	0.38	Eruptive	>4.5	>2751	1.67	9.1	2.5
L222_149	266.82162	−28.32786	–	3.98	Dipper	5.4	1157	−0.93	–	2.5
L222_154	266.95925	−28.24306	–	4.19	Dipper	4.8	>2165	2.16⋆	–	–
L222_165	267.41287	−28.44864	0.85	0.95	Eruptive	4.3	>2151	–	1.6	–
L222_167	267.60937	−28.87519	–	2.78	Eruptive	5.7	>2501	2.33	7.8	3.9
L222_168	267.66867	−28.06083	2.79	1.68	Dipper	4.4	∼400	−0.93⋆	8.3	1.4
L222_172	267.86138	−26.12528	–	3.13	Dipper?	>5.2	>4151	1.87⋆	–	–
L222_192	269.43525	−24.34244	1.07	0.85	Eruptive	4.8	>1881	0.40⋆	6.1	0.3
L222_210	271.91033	−21.76933	2.37	2.21	Eruptive	3.6	>2951	0.62⋆	5.0	1.2

Notes.^⋆Targets without W4 (22 μm) and MIPS (24 μm) detections. ^† Also known as DR4_v67 from Guo et al. (2021).

Multitime-scale variable (MTV): defined in Guo et al. (2020), MTVs have a general variation trend, but their short-term variability (i.e. times-cales < 1 yr) is comparable to the long-term variability.

The H − K_s colours are measured from the VVV photometry.

2.2 Spectroscopic observation and data reduction

Two spectroscopic campaigns were obtained by the X-shooter 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 z-band acquisition image. All targets were observed simultaneously with 0.9 arcsec slit width in optical and 0.6 arcsec 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 R = 8000 in K_s bandpass. The X-shooter 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 (Kausch et al. 2015; Smette 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⁻¹ 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 second-order polynomial function to generate the wavelength solution.¹ 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 Figs B1–B4.

Table 2.

Spectroscopic observations of YSOs.

Name	tobs	texp	Type	|$A_{K_s}$|	RV	dk
	d/m/y	s		mag	km s−1	kpc
X-shooter/VLT
L222_1	16/04/22	600	FUor	1.1	40	10.2|$^{+1.0}_{-0.8}$|
L222_4	28/04/21	960	FUor?	–	−20	6.6|$^{+0.9}_{-0.9}$|
L222_10	16/04/22	600	FUor	1.5	−47	3.7|$^{+0.5}_{-1.0}$|
L222_13	28/04/21	1200	FUor	0.7	−30	8.4|$^{+0.4}_{-1.1}$|
L222_15	03/05/21	960	FUor	0.7	−15	1.2|$^{+0.4}_{-0.7}$|
L222_18	16/04/22	960	FUor	1.7	−30	2.0|$^{+0.5}_{-0.6}$|
L222_25	15/04/22	960	FUor	0.9	−10	0.8|$^{+0.3}_{-0.6}$|
L222_33	29/04/21	960	FUor	1.9	−70	4.2|$^{+0.5}_{-0.5}$|
VVV1640-4846	15/05/22	960	FUor	2.0	−40	2.8|$^{+0.3}_{-0.2}$|
L222_73	15/05/22	960	FUor?	–	−30	3.8|$^{+0.6}_{-1.1}$|
L222_78	16/04/21	600	FUor	0.1	−10	4.2|$^{+1.2}_{-1.5}$|
L222_93	29/04/21	960	FUor	1.0	−95	7.8|$^{+0.3}_{-0.4}$|
L222_95	29/04/21	960	FUor	1.8	−50	6.9|$^{+0.4}_{-1.0}$|
L222_165	04/05/21	960	FUor	0.4	−30	–
L222_192	01/07/21	960	FUor	0.9	5	3.0|$^{+1.1}_{-2.0}$|
L222_6	16/04/22	960	Em. line	2.4	–	–
L222_28	03/05/22	960	Em. line	–	−10	0.9|$^{+0.3}_{-0.8}$|
VVV1636-4744	11/05/22	960	Em. line	–	–
L222_148	04/05/21	960	Em. line	1.4	−10	8.3|$^{+0.2}_{-0.2}$|
L222_167	29/04/21	960	Em. line	0.6	−5	7.8|$^{+0.8}_{-3.1}$|
L222_167	09/06/21	960	Em. line	0.4	−5	7.8|$^{+0.8}_{-3.1}$|
L222_210	10/07/21	960	Em. line	0.1	–	–
L222_32	29/04/21	960	Outflow	1.5	–	–
	08/06/21	960	Outflow	1.8	–	–
L222_37	07/05/22	960	Outflow	–	–
VVVv800	29/04/21	960	Outflow	1.2	–	–
L222_120	11/05/22	960	Outflow	3.4	–	–
FIRE/Magellan
L222_4	20/05/19	475	FUor	3.1	–	–
	21/05/19	1056	FUor	3.1	–	–
L222_15	22/05/19	634	FUor	0.6	–	–

Note. d_k: kinematic distance calculated using the online tool provided in Wenger et al. (2018).

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 arcsec slit was used, providing a spectral resolution of R = 6000 in K bandpass. The data were 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 K 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. 2020, 2021, for more details). Finally, we used the OH skylines from the scientific exposures as wavelength calibration.

2.3 Photometric data

2.3.1 Near-infrared photometric data

The pre-release VIRAC2β version of the VIRAC2 catalogue provides near-infrared photometry that contains dozens of K_s-band measurements and a few multicolour epochs (Smith et al., in preparation). Additional processing of a few sources was conducted to correct for saturation or spatially resolved nebulosity. The first multicolour 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 K_s light curves were resampled into 1-d bins by the median brightness after removing outlier detections (3σ in each bin).

In 2021 April and 2022 May, two epochs of J, H, K_s 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; 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 K_s brightness.

2.3.2 Archival-infrared photometric data

We obtained photometric data from Wide-field Infrared Survey Explorer (WISE) and Spizter surveys via the NASA/IPAC Infrared Science Archive (IRSA²). The WISE W1 to W4 photometry is adopted from the All-Sky data release when most sources were at the quiescent stage (Wright et al. 2010). The W1- and W2-band light curves were obtained from the multi-epoch photometry table of the AllWISE and the Near-Earth Object WISE (NEOWISE) surveys (Mainzer et al. 2014, 2014–2022). The W1 and W2 time-series were sampled into 1-d 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 I1 to I4 photometry was obtained from the Glimpse survey (Benjamin et al. 2003, prior to the VVV survey), and the 24 μm data were adopted from the Multiband Infrared Photometer Spitzer Galactic plane survey (MIPSGAL, 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. W1 and W2 bands), sources located within 8 arcsec are not well distinguished. We set a threshold of 2 arcsec 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., 2024).

Table 3.

Spectroscopic information of post-MS sources.

Name	tobs	EWCO	EWNa	[Fe/H]	Teff	TCO
	d/m/y	Å	Å		K	K
L222_59⋆	29/04/21	19.1	1.8	−0.68	3648	1200
L222_130	15/05/22	18.2	5.0	0.62	3744	3800
L222_144	15/05/22	18.0	7.0	1.00	3765	3600
L222_145	16/05/22	22.1	5.2	0.28	3325	3800
L222_149	04/05/21	24.1	5.2	0.10	3112	2800
L222_154	16/05/22	23.2	6.0	0.30	3208	2400
L222_168	23/05/22	14.9	2.1	−0.42	4097	3800
L222_172	16/05/22	17.8	5.9	0.97	3788	3600

Note. T_eff is measured using the methods from Feldmeier-Krause et al. (2017), with a typical uncertainty of 163 K. T_CO is measured from the analytical models of CO bandheads (Contreras Peña et al. 2017b), with an error bar of 200 K. ^⋆The measurements of T_eff and [Fe/H] are not reliable due to the unique spectral features of L222_59, for example, circumstellar CO.

Table 4.

Variation amplitudes and pre-outbursting colour of eruptive YSOs.

Name	Type	trise	ΔKs, rise	ΔKs, res	ΔW2	Ks − W2
		d	mag	mag	mag	mag
L222_1	FUors	609	3.2	0.7	2.1	3.6
L222_4⋆	FUors	793	6.1	1.0	>4.8	<4.6
L222_10	FUors	939	3.3	0.5	–	–
L222_13⋆	FUors	481	4.0	0.4	>2.5	<3.1
L222_15	FUors	995	4.6	1.2	>4.0	<3.4
L222_18†	FUors	510	4.0	0.6	3.8	3.3
L222_25	FUors	849	3.7	1.2	–	–
L222_33†	FUors	751	4.1	0.7	>3.5	2.8
VVV1640-4846†	FUors	747	3.0	1.0	1.9	5.2
L222_73‡	FUors	>589	>4.3	0.9	>3.4	5.3
L222_78	FUors	500	4.6	0.6	3.3	2.8
L222_93	FUors	527	4.6	1.4	>3.2	<3.1
L222_95†	FUors	1439	5.4	1.3	3.3	4.8
L222_165	FUors	705	3.9	1.6	>3.3	<2.9
L222_192†	FUors	463	4.9	0.6	2.8	2.6
VVVv322	FUors	129	2.3	0.6	–	3.3
VVVv721‡	FUors	>2140	>3.1	0.6	>2.0	4.3
DR4_v20‡	FUors	>838	>3.3	0.6	>1.1	5.8
VVVv270	Em. Line	895	3.3	0.7	1.3	5.7
VVVv631	Em. Line	1926	2.6	1.2	1.5	3.2
L222_6	Em. Line	900	3.4	1.4	0.6	6.6
L222_148⋆	Em. Line	2040	3.6	2.1	>0.7	<6.9
L222_167	Em. Line	1038	5.6	1.3	1.7	7.5
DR4_v10‡	Em. Line	>1173	>2.9	1.5	>0.8	6.1
DR4_v34	Em. Line	817	3.7	1.2	0.3	6.5
Stim1‡	Em. Line	>2126	>2.2	0.8	>1.6	3.6
Stim5‡	Em. Line	>874	>3.4	1.3	>0.7	5.1
L222_32	Outflow	1088	3.3	2.4	0.4	7.9
L222_37‡	Outflow	>808	>4.9	1.2	>2.6	5.0
L222_120	Outflow	808	4.5	2.2	1.8	4.8
L222_210⋆	Outflow	941	3.7	1.9	>1.0	<6.3
VVVv800	Outflow	800	2.0	0.8	1.6	6.2

Notes. This table contains 32 large-amplitude and long-duration outbursting YSOs discovered in the VVV survey (Contreras Peña et al. 2017b; Guo et al. 2021;,LSG23, and this paper).

^⋆Eruptive objects without clear detections from the WISE images due to blending with a spatially nearby source.

^†Sources without W2 detections during the photometric minima, but have I2-band detections from Spitzer.

^‡Eruptive events started their rising stage before the VVV survey.

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 (Classes 0 and I) have positive values as α_class > 0.3, while T Tauri stars (Classes II and 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 K_s and WISE All-Sky survey. The line-of-sight extinction could lead to a higher α_class as K_s is more sensitive to extinction than the longer wavelength bands. However, after examining the SEDs (see Appendix A), we found that including K_s 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.

3 RESULT

3.1 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 subcategories based on their spectral features, which also reveal their physical conditions during the eruptive stage.

3.1.1 FUor-type

During an 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 an 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 H-band continuum due to low gravity water vapour absorption, and CO overtone bandheads beyond 2.3 μm (see Fischer et al. 2023). 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–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 H-band continuum, to align with the spectrum of FU Orionis (provided by Connelley & Reipurth 2018), using the extinction law from Wang & Chen (2019) and A_V = 1.5 mag estimated in Zhu et al. (2007). Subsequently, we estimated the K_s-band extinction (⁠|$A_{K_s}$|⁠) for each FUor-type object by applying the minimum χ² method to compare the de-reddened spectra with the standard H-band spectrum, with |$A_{K_s}$| values ranging between 0 and 10 mag. However, we acknowledge that this is only a rough estimation of |$A_{K_s}$| as the disc spectra of FUors are not necessarily the same from one case to another. The de-reddened H- and K-bandpass spectra of FUor-type sources are presented in Fig. 1.

Figure 1.

Normalized near-infrared (1.5–2.4 μm) X-shooter/VLT spectra of 13 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 H bandpass continuum, to match the spectrum of FU Orionis (provided by Connelley & Reipurth 2018). The spectra are ordered by their H − K_s colour after de-reddening, with the reddest on top.

The X-shooter spectra of two FUors (L222_4 and L222_73) were taken during their fading stages (see Fig. 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 H-band continuum is absent or less prominent due to the low brightness in H, and the H₂ 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 K_s-band brightness was higher (see Fig. B5, left panel). After the original 6 mag outburst in K_s (the highest amplitude in this sample), L222_4 entered a rapidly fading stage that decayed 3 mag in 5 yr. 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-fitting models were confirmed through visual inspection. The grids of the model are ΔRV  = 5 km s⁻¹, smaller than the pixel size of X-shooter (12 km s⁻¹). The kinematic distance (d_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 d_k agrees with d_SFR, suggesting that the variable source is associated with nearby SFRs. When comparing d_k with d_SFR (listed in Table 1), many FUors have d_SFR beyond the 1σ error of d_k. However, d_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.

3.1.2 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 time-scale of 1000 d 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. 2023). In the near-infrared, a group of eruptive YSOs was discovered by the VVV survey with long-lasting 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 H₂ 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 Fig. 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 mass-accretion rate. One source, L222_167 was observed on two epochs 40 d apart and the continuum on the first epoch is 13 per cent brighter than the second epoch (Fig. 2). Stronger mass-accretion indicators are seen in the brighter epoch, whilst the fainter epoch has stronger indicators of the outflow (H₂ 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 H₂ lines, including a blueshifted high-velocity component around −130 km s⁻¹, a Gaussian component around −33 km s⁻¹, and a weak redshifted tail. This complex line feature indicates that winds/outflows were launched towards different directions, and sometimes at perpendicular angles.

Figure 2.

Left: K bandpass X-shooter 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 visualization. Right: normalized line profiles of H₂ S(1) lines at 2.12 μm. 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.

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. (2020, 2021). In comparison with the latest spectra, there is no significant spectroscopic variability. In the case of L222_32, two X-shooter spectra were obtained in 40 d. Despite very similar continuum levels, the H₂ S(1) line in the second epoch is 35  per cent 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 time-scale, 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 H₂ 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, Kirby-Docken & Dalgarno 1977). We applied the power-law extinction function from Wang & Chen (2019). The typical precision of the |$A_{K_s}$| measurement is 0.4 mag, with uncertainty rising from the flux measurement of the Q(3) line (σ = 8  per cent), which is veiled by the continuum and the CO overtone emission. The colour variation before and after the outburst is discussed further in -Section 4.3.

3.1.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 e.g. 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 Two Micron All Sky Survey (2MASS), for example, sources L222_144, L222_145, L222_154, L222_144, L222_168, and 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 (H − K_s > 3.5 mag) or lack J-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 ¹³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 Fig. 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 Fig. 3. Most post-MS sources exhibit redder H − K_s 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 Section 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 and 4000 K and supersolar 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 Appendix  D.

Figure 3.

H − K_s and K_s − W2 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 K_s 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 W2 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.

We highlight an eruptive source (L222_59, a.k.a. VVVv746) that exhibits rarely seen deep aluminium monoxide (AlO) absorption bands (see Fig. 4 and Fig. B4). This source experienced a ΔK_s > 5 mag outburst and reached the photometric maximum in less than three years. According to our follow-up K_s photometry, after spending nearly 1000 d on its brightness plateau, L222_59 dramatically faded nearly 3 mag in the year 2021–2022. The X-shooter spectrum presented in this paper was taken during the outbursting stage, with many absorption features, including cool ¹²CO overtone bands in the K bandpass (as fitted temperature, T = 1200 K, assuming local thermodynamic equilibrium) and AlO in the J and H bandpasses, confirmed by the ExoMol database (Patrascu, Yurchenko & Tennyson 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 H i 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 time-scale of 50 d (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.

Figure 4.

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

4 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 time-scales of the outburst, and the near to mid-infrared amplitudes, are measured to distinguish FUors from other eruptive events.

4.1 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 minimize the selection bias, we included all sources with ΔK_s ≥ 2 mag and have eruptive events lasting longer than 1000 d. 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 lightcurve 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 K_s light curves of these sources, including data from both the VVV survey and SOFI observations, are shown in Fig. C1. We note that three FUors previously identified in the VVV survey are not included in this sample, due to the absence of K_s 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., 2024), we will further discuss the occurrence rate of eruptive behaviours on Class I protostars, by applying mid-infrared SED-based catalogues (e.g. Spitzer/IRAC Candidate YSOs, Kuhn et al. 2021) and including archival K_s-band detections (e.g. 2MASS) to form a two-decade-long timeline.

4.2 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 yr (Herbig 1977). Here, we apply linear functions to describe the quiescent, rising, decaying, and post-outbursting stages of the aforementioned 32 eruptive events (see e.g. Fig. 5). Two parameters were defined to quantify the rising stage: the rising amplitude (ΔK_{s, rise}) and the rising time-scale (t_rise). Specifically, ΔK_{s, rise} is the K_s amplitude of the rising stage, measured from the linear fitting result. The t_rise is calculated as ΔK_{s, rise} divided by the measured linear slope (mag d⁻¹) of the rising light curve. The measured t_rise and ΔK_{s, rise} are presented in the right panel of Fig. 5, with numbers listed in Table 4.

Figure 5.

Left: K_s 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 time-scales 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 K_s (filled stars). Outbursts started in the pre-VVV era were labelled with lower limits. The distributions of rising time-scales 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 (p_WMW) for each distribution is labelled on the histogram.

We present the statistical differences of t_rise and ΔK_{s, rise} between FUors and non-FUors (emission-line/outflow-dominated objects). In general, we find that FUors experienced shorter t_rise compared to non-FUors. Among the examined samples, 14/18 FUors have t_rise < 850 d (mean at 747 d), whilst all non-FUors have t_rise > 800 d (mean at 1038 d). To support this finding, we employed the Wilcoxon–Mann–Whitney test (Mann & Whitney 1947), a non-parametric rank-based test. The possibility of FUor and non-FUor sharing the same distribution in t_rise is only p_WMW = 0.001. There should not be any selection bias on the t_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 time-scale less than 1 yr. Among the longest duration objects, we detect eruptive sources with t_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 ΔK_{s, rise}, although FUors are more abundant towards the highest amplitude end, as 61 per cent of FUors and 21 per cent of non-FUors have ΔK_{s, rise} > 3.8 mag. However, FUors and non-FUors are not distinguishable at the intermediate-amplitude range (2.0 < ΔK_{s, rise} ≤ 3.8 mag). The p_WMW on the ΔK_{s, rise} distribution is 0.07, larger than the one measured on t_rise. Since there is no strongly justified null hypothesis regarding the relative amplitudes of the FUors and non-FUors, we regard this p_WMW-value as fairly good evidence that the amplitude distributions are different.

We further compare the t_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 fast-eruptive sources with t_rise < 300 d, other archived FUors share similar t_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 t_rise in the optical than in the infrared, in the order of years (Cleaver, Hartmann & Bae 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 t_rise > 2140 d. 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 K_s 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, and 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 time-scales 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 light curves to measure the short-time-scale residual variability (ΔK_{s, res}), as the residual amplitude within a time window of two years in each measurement (see Fig. 6). In general, FUors have smaller ΔK_{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-time-scale 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 (ΔK_s) variation seen among non-FUors.

Figure 6.

The ratio between the short-term residual amplitude and the rising amplitude versus the rising time-scale of eruptive sources. Colour codes and lower limits are as same as in Fig. 5.

4.3 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) summarized 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 K_s to W2, 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 W2-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 ΔW2 and K_s − W2 can not be measured. The results are listed in Table 4.

The correlation between ΔW2 and ΔK_s is shown in Fig. 7. FUors exhibit systematically higher ΔW2 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 ΔW2 > 2 mag, even among those with ΔK_s > 3.5 mag. The infrared colour indices are presented in Figs 7 and 8. Lower limits are put on objects that erupted before the VVV survey. All FUors have K_s − W2 < 6 mag at the quiescent stage, whilst the majority of long-term non-FUors have a redder K_s − W2 colour. Similar trends are seen in the W1 − W2 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 H − K_s colour distribution shown in Fig. 3, where most FUors are bluer than non-FUors during the quiescent phase. As comparisons, we added some archived eruptive YSOs in Fig. 8, including three FUors and three recently confirmed eruptive objects from the NEOWISE data (with prefix SPICY, Contreras Peña et al., submitted).

Figure 7.

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

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 ΔW2 seen among red non-FUors (K_s − W2 > 6 mag). However, the origin of their substantial excess in W2 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 Fig. 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 ΔK_s than ΔW2. We also admit that there is a selection bias in our spectroscopic sample (i.e. only sources with high ΔK_s 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.

Figure 8.

Quiescent W1−W2 colour and ΔW2 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 W1−W2 and ΔW2 are presented as histograms.

In summary, we find FUors have bluer K_s − W2 and W1 − W2 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 preparation) will obtain spectroscopic confirmation of a group of extremely red candidate FUors discovered from the mid-infrared NEOWISE time-series.

5 THE PERIODIC VARIATION ON L222_148

One emission-line object, L222_148, exhibits repeated lower amplitude variation in addition to the overall long-term eruptive behaviours. We applied an analytical fit to the K_s light curves of L222_148, composed of a sinusoidal function and three linear slopes as the quiescent stage, rising stage, and the brightness plateau (see Fig. 9). A 378-d quasi-periodic signal is detected with ΔK_s ∼ 2.0 mag after applying the generalized 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 recognize 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 d 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).

Figure 9.

Upper left: K_s 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 H₂ emission lines and the Brγ line. The RV 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.

The near-infrared spectrum of L222_148, taken during the outbursting stage, shows emission features associated with the magnetospheric accretion process (e.g. Na ii and CO bandheads) and coupled with strong stellar outflow (H₂ line series). In contrast, the Brγ line is in absorption. Inverse P Cygni profiles are seen on the H₂ lines, suggesting a combination of blueshifted stellar outflow coupled with a stream of cold disc materials (see the right panels of Fig. 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₂ 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 GIs 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 time-scales 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).

6 SUMMARY

In our series of works, over 200 high-amplitude variable stars were identified from the decade-long VVV/VVVX K_s 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 X-shooter spectrograph on VLT and FIRE on the Magellan telescope. These sources are categorized into a few groups based on their spectral features and the nature of their variability. We summarize 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 H-band continuum.

• Four targets are classified as outflow-dominated with H₂ 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-MS sources with the existence of ¹³CO absorption features and extremely red near-infrared colours (five of them have H − K_s > 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 long-duration eruptive YSOs originally identified from the VVV survey. We find that in the near-infrared, FUors have faster eruptive time-scales and predominate the highest amplitude end (ΔK_s > 3.8 mag). In the mid-infrared, FUors have systematically higher W2 amplitude than other eruptive objects. It suggests that FUors can heat their inner accretion disc more efficiently.

• FUors have comparable near- and mid-infrared amplitudes. By contrast, most non-FUors have much lower amplitudes in the mid-infrared than in the near-infrared, with redder K_s − W2 and W1 − W2 colours, indicating early evolutionary stages or high line-of-sight extinction.

• In Guo et al. (2021), we found that most multiyear 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 ΔK_s > 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 time-scales and amplitudes of eruptive events in two subcategories. To further reveal the physical nature behind these eruptive events, high-cadence multiwavelength observations and numerical simulations are desired.

ACKNOWLEDGEMENTS

We thank the anonymous referee for providing valuable comments and suggestions to improve this work. ZG is supported by the ANIDFONDECYT Postdoctoral program no. 3220029. ZG and KM acknowledge support by ANID, – Millennium Science Initiative Program– NCN19_171.

ZG, PWL, and CJM acknowledge support by STFC consolidated grants ST/R00905/1, ST/M001008/1, and ST/J001333/1 and the STFC PATT-linked grant ST/L001403/1. This work has made use of the University of Hertfordshire’s high-performance computing facility (http://uhhpc.herts.ac.uk).

We gratefully acknowledge data from the ESO Public Survey program ID 179.B-2002 taken with the VISTA telescope, and products from the Cambridge Astronomical Survey Unit (CASU). This work contains data from ESO programs 105.20CJ and 109.233U.

JB and RK thank the support from the Ministry for the Economy, Development and Tourism, Programa Iniciativa Cientifica Milenio grant IC120009, awarded to the Millennium Institute of Astrophysics (MAS). DM gratefully acknowledges support from the ANID BASAL projects ACE210002 and nFB210003, from Fondecyt Project No. 1220724, and from CNPq Brasil Project 350104/2022-0. JA-G acknowledges support from Fondecyt Regular 1201490 and from ANID’s Millennium Science Initiative ICN12_009, awarded to the Millennium Institute of Astrophysics (MAS). ACG acknowledges support from INAF-GOG ’NAOMY: NIR-dark Accretion Outbursts in Massive Young stellar objects’ 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. R.K.S. acknowledges support from CNPq/Brazil through projects 308298/2022-5, 350104/2022-0 and 421034/2023-8. AA acknowledge the support through a Fellowship for National PhD students from ANID, grant no. 21212094. CCP was supported by the National Research Foundation of Korea (NRF) grant funded by the Korean government (MEST) (no. 2019R1A6A1A10073437) Support for MC 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. KM 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.

7 DATA AVAILABILITY

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

Footnotes

1

Originally proposed in the master thesis of A. Aguayo.

2

https://irsa.ipac.caltech.edu

APPENDIX A: NEAR TO MID-INFRARED SEDS

Here, we present the SED of our targets between 1–24 μm. Wherever possible, we applied contemporaneous J, H, and K_s photometry from VVV or NTT/SOFI observations, and contemporaneous mid-infrared data from ALLWISE (preferred) or Spitzer surveys. The spectral index (α) is measured as the slope of the linear fit to the quiescent SED between 2–24 μm. There are eight sources without the measurement of α, including three sources that were only detected by WISE or Spitzer during their photometric maxima, three sources that were blended during their photometric minima and two sources without any mid-infrared detections. The SEDs of these eight sources are not shown. The spectral index is not measured on sources without detections beyond 5 μm.

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; and Spitzer: red. Linear slopes are fit to the quiescent SEDs at λ > 2μm, among sources having detections beyond 5 μm.

APPENDIX B: NEAR-INFRARED SPECTRA

We present the near-infrared spectra obtained for this work. Sources are sorted into four groups, FUors, emission-line objects, outflow-dominated sources, and post-MS objects. For each spectrum, we also present the K_s light curves obtained from the VVV survey and our follow-up observations on NTT/SOFI. The date of the spectroscopic observation is specified in each subplot by the vertical lines. For the best visualization, all spectra are smoothed with a boxcar average of 5 pixels.

Figure B1.

X-shooter spectra and K_s 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.

X-shooter spectra and K_s light curves of emission-line objects. Symbols are the same as in Fig. B1.

Figure B3.

X-shooter spectra and K_s light curves of outflow dominated objects. Symbols are the same as in Fig. B1. Bottom right: line profiles of three H₂ lines from two spectral epochs of L222_32. The vertical dashed lines show the RVs of Gaussian fittings.

Figure B4.

X-shooter spectra and K_s light curves of post-MS giants. Symbols are the same as in Fig. B1. The J- and K-bandpass spectra of L222_59 are presented in the bottom two plots.

Figure B5.

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

APPENDIX C: THE NEAR-INFRARED LIGHT CURVES OF ERUPTIVE OBJECTS

We present the K_s-band light curves, from both the VVV survey and our NTT/SOFI follow-up observations, of 32 eruptive objects that exhibit long-duration outbursts (see Fig. C1). In most cases, we artificially divided the light curves into four stages: pre-outbursting/quiescent, rising, decaying, and post-outbursting stages. The post-outbursting stage is only applied to sources that exhibit further eruptive or rapid decaying stage after the original outburst/photometric plateau. As mentioned in the main text, we used linear fits to measure the amplitude and duration of each stage. We noticed that not all light curves are well presented by linear trends, especially the ones with outstanding short-time-scale variability as opposed to the overall eruptive behaviour, such as L222_6, L222_32, and L222_148.

Figure C1.

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

APPENDIX D: MEASUREMENTS OF STELLAR PARAMETERS OF POST-MS STARS

In this section, we briefly describe our methods to fit stellar parameters on the post-MS stars. First, we applied the theoretical models of the CO bandhead emission beyond 2.29 μm. This model was initially designed in Contreras Peña et al. (2017b), using the rovibrational states and coefficients from Farrenq et al. (1991). Three free parameters were initially designed in the model, including the RV, the temperature of CO gas (T_CO), and the column density of CO molecules (N_CO). To measure these parameters from the observed spectra, we generated synthetic spectra from the models, with T_CO ranging between 800 and 6000 K (100 K grids) and N_CO between 1e20 and 5e21 cm⁻² (1e20 cm⁻² grids). We first find the RV of the CO absorption following the methods described in the main text. Then, we removed the spectral continuum by applying polynomial fits and used the median flux between 2.285 and 2.290 μm to normalize the spectrum. Finally, we applied the minimum-χ² method to locate the best-fitting models among our synthetic spectra. In the last step, to reduce the noise level, we smoothed both the observed and synthetic spectra by five pixels. The measured CO temperatures are presented in Table 3 of the main text. The typical error of this fitting is ±200 K, mainly constrained by the signal-to-noise ratio of the spectra and the accuracy of the continuum fitting. Among these post-MS sources, L222_59 has a unique CO morphology that indicates a very cool effective temperature (1200 K) and high column density (8e20 cm⁻²). This cool T_CO is consistent with the AlO model presented in Fig. 4 of the main article. Notably, such a low temperature implies the feature arises in circumstellar matter rather than the stellar photosphere.

We applied a method based on the EW of the CO first overtone bands and Na i absorption features to measure the metallicity, as described in Fritz et al. (2021). The measured EWs and [Fe/H] are listed in Table 3, with a typical error bar of 10 per cent introduced by the standard deviation of the continuum flux. Notably, in the table, we changed the EW of absorption features to positive, to be consistent with standard practice. In addition, we used the method provided by Feldmeier-Krause et al. (2017), to calculate the effective temperature of the stellar atmosphere (T_eff) from the EW of CO bandhead, as

$$\begin{eqnarray} T_{\rm eff}=5677.0 \, K-106.3*EW_{\rm CO}. \end{eqnarray}$$

(D1)

The intrinsic scatter of this method is 163 K according to Feldmeier-Krause et al. (2017), which is similar to the uncertainty introduced by the measurements of EW. The measured T_eff are also presented in Table 3 of the main article. Except for one source with very cool CO, L222_59, the T_eff measured from EW_CO are reasonably similar to T_CO measured from the shape of the CO absorption features, considering the errors. T_CO is slightly lower than T_eff in 5/7 cases, as one would expect for an absorption feature. According to the parsec stellar evolutionary models (see Pastorelli et al. 2019, and references therein), sources with T_eff ranging between 3000 and 4000 K have typical intrinsic H − K_s colours between 0.15 and 0.35 mag. In comparison with the observed H − K_s colours from the VVV survey (listed in Table 1), we found that these post-MS sources are indeed severely embedded. Applying the extinction law from Wang & Chen (2019), A_V = 18.8E(H − K_s), the average extinction of these sources are A_V = 57 mag or |$A_{K_s} = 4.5$| mag.