The most variable VVV sources: eruptive protostars, dipping giants in the Nuclear Disc and others

We have performed a comprehensive search of a VISTA Variables in the Via Lactea (VVV) database of 9.5 yr light curves for variable sources with $\Delta K_s \ge 4$ mag, aiming to provide a large sample of high amplitude eruptive young stellar objects (YSOs) and detect unusual or new types of infrared variable source. We find 222 variable or transient sources in the Galactic bulge and disc, most of which are new discoveries. The sample mainly comprises novae, YSOs, microlensing events, Long Period Variable stars (LPVs) and a few rare or unclassified sources. Additionally, we report the discovery of a significant population of aperiodic late-type giant stars suffering deep extinction events, strongly clustered in the Nuclear Disc of the Milky Way. We suggest that these are metal-rich stars in which radiatively driven mass loss has been enhanced by super-solar metallicity. Among the YSOs, 32/40 appear to be undergoing episodic accretion. Long-lasting YSO eruptions have a typical rise time of $\sim$2 yr, somewhat slower than the 6-12 month timescale seen in the few historical events observed on the rise. The outburst durations are usually at least 5 yr, somewhat longer than many lower amplitude VVV events detected previously. The light curves are diverse in nature, suggesting that multiple types of disc instability may occur. Eight long-duration extinction events are seen wherein the YSO dims for a year or more, attributable to inner disc structure. One binary YSO in NGC 6530 displays periodic extinction events (P=59 days) similar to KH 15D.


INTRODUCTION
Time domain astrophysics has led to the discovery of several new types of variable star in recent years.Mid-infrared studies have found ★ E-mail: p.w.lucas@herts.ac.uk (PWL) transients known as "eSPecially Red Intermediate Luminosity Transient Events (SPRITES)" in nearby galaxies, having luminosities between classical novae and supernovae (Kasliwal et al. 2017).In optical studies of variable stars, the abundant OGLE Small Amplitude Red Giants (OSARGs, Soszynski et al. 2004;Soszyński et al. 2013) and the rarer Blue Large Amplitude Pulsators (BLAPs, Pietrukowicz et al. 2017) have been found by the Optical Gravitational Lensing Experiment (Udalski 2003).Circumstellar matter is the cause of other new types of variable star.For example, ASAS J140748−3945.7 (J1407) hosts a transiting giant ring system bound to the substellar companion of a pre-main sequence star (Mamajek et al. 2012).Unpredictable dips and slow variations in brightness are also seen in main sequence stars and giant stars (e.g.Boyajian et al. 2016Boyajian et al. , 2018;;Schmidt 2019;Rodriguez et al. 2016;Smith et al. 2021), attributed to occultation by irregular distributions of circumstellar matter in the main sequence stars or a circumsecondary disc in the case of the giants.
The VISTA Variables in the Via Lactea (VVV) survey (Minniti et al. 2010) is the first major near infrared variability survey of the Milky Way.In view of the above list of discoveries, we might reasonably hope that it will discover new types of infrared variable source.Due to the Wien displacement law, it is perhaps likely that the physical mechanisms behind any such infrared discoveries would involve circumstellar matter.The most numerous high amplitude variable stars discovered by VVV thus far are YSOs with circumstellar accretion discs, see Contreras Peña et al. (2017a,b), hereafter CP17a and CP17b and Medina et al. (2018); Teixeira et al. (2018); Guo et al. (2020Guo et al. ( , 2021Guo et al. ( , 2022)).These generally have amplitudes Δ  ≤ 4 mag in these previous works (measured as  , −  , ) and a diverse range of light curve types.Short term variability on timescales of a few days is common at 1 < Δ  < 1.5.At higher amplitudes, longer timescale variability predominates: mostly cases of eruptive variability but also cases of variable extinction, sometimes in the same star (see e.g.Kóspál et al. 2011).CP17a found that eruptive variability is dominated by optically faint class I and flat spectrum YSOs, with very few examples of class II YSOs being observed.Variability in embedded YSOs is more easily studied in the near infrared, the mid-infrared (e.g. Park et al. 2021;Contreras Peña et al. 2023) or the submillimetre waveband (Lee et al. 2021).
Most past studies of YSO outbursts have focussed on the most dramatic events, either long duration events (FUors) with optical amplitudes up to 6 mag or more (e.g.FU Ori, V1057 Cyg, V1515 Cyg, V346 Nor and PTF14jg Herbig 1977;Graham & Frogel 1985;Kóspál et al. 2020;Hillenbrand et al. 2019) or shorter duration events (EXors) with amplitudes up to 5 mag such as EX Lup and V1118 Ori (Herbig 2008).In recent years, the Gaia Alerts data (Kostrzewa-Rutkowska et al. 2020;Hodgkin et al. 2021) have enabled more systematic searches that have discovered two new FUors, Gaia 17bpi (Hillenbrand et al. 2018) and Gaia 18dvy (Szegedi-Elek et al. 2020) amongst other things, complementing the more general Gaia classification of YSOs and other variable stars (Eyer et al. 2023;Rimoldini et al. 2023;Marton et al. 2023).VVV and its extension, VVVX (vvvsurvey.org)offer the chance to establish a large sample of events that characterises the outburst onset better than in most previous (often accidental) discoveries and lacks any significant bias against lower amplitude events.The great majority of outbursts have lower amplitude, especially in the near infrared.E.g.CP17a listed 106 candidates with 1 < Δ  < 4, with mean Δ  = 1.72.
Outburst durations do not always fit neatly into the original FUor (multi-decade) and EXor (8-16 month) categories (Herbig 1977(Herbig , 2008)).Spectroscopically confirmed VVV outbursts (CP17b) appeared to have typical durations of 1-5 yr (though sometimes longer) and the spectroscopic characteristics of the two traditional outburst classes were sometimes mixed or swapped.This led to the tentative suggestion that a new "MNor" outburst category might be needed to describe mixed outbursts (CP17b).More recently Guo et al. (2021) showed that within the broad category of mixed outbursts there is a distinct class of long duration eruptions with EXor-like spectra, apparently with the accretion luminosity of FUors.In fact these long lasting but EXor-like outbursts were more numerous than FUors amongst the long duration events in that work.
In this work, we perform a thorough search of the 562 deg 2 VVV area for very high amplitude (Δ  ≥ 4) variable stars and transients.This threshold was chosen because previous VVV-based searches had found very few eruptive YSOs with such high amplitudes.A complete search for all types of variable sources is more difficult at lower amplitudes, if all false positives are to be removed, due to the larger number of candidates.While our main interest was eruptive YSOs, unexpectedly we also discovered a population of aperiodic dipping giants in the Nuclear Disc of the Milky Way.Consequently, this paper focuses on these two topics.Supporting spectra of YSOs and giants are found in a companion paper (Guo et al. submitted, hereafter GLK).Separately, in Contreras Peña et al.(submitted) we present results of a VVV-based search for eruptive YSOs with Δ  >2, within published catalogues of YSO candidates.
VVV observed much of the Galactic bulge and the adjacent Galactic disc from 2010-2015 and this study benefits from the ongoing VVVX extension in the same area, though we have not included the new areas of the Galactic plane area monitored by VVVX since 2016.This new search benefits not only from longer duration light curves and larger area than the 119 deg 2 search of CP17a but also from the new VVV Infrared Astrometric Catalogue version 2 (VIRAC2) profile-fitting photometry database (Smith et al. 2018, Smith et al., in prep.) which is deeper and more reliable than the standard VISTA data pipeline in crowded fields.VIRAC2 light curves have reduced noise due to a carefully optimised relative photometric calibration and the absolute photometric calibration is also improved, particularly in crowded inner-bulge fields (Smith et al., in prep.).
In §2 and Appendix A we describe the data and the variable star searches.§3 summarises our initial classifications of the variable sources and their distribution on the sky.In §4 we examine the YSOs, with emphasis on the timescale of the rise to maximum light and the form of the light curve in FUor-like events, and in section §5 we discuss the aperiodic dipping giants.The dipping YSOs are presented fully in Appendix E. Similarly, a few noteworthy findings among other types of variable star (LPVs, novae and unusual sources) are summarised in §6 and discussed more fully in appendices F, G and I. Microlensing events are discussed only briefly in Appendix H. Our conclusions are given in §7.

Observations
The VVV and VVVX surveys employed the VISTA telescope (Sutherland et al. 2015) at Cerro Paranal observatory and the VIR-CAM camera (Dalton et al. 2006), which houses 16 HgCdTe 2048 2 Raytheon infrared array detectors.From 2010-2015, VVV surveyed a region of the Galactic bulge of at Galactic coordinates |ℓ| < 10 • , −10 • <  < 5 • and the adjacent southern Galactic disc at 295 The region is divided into filled 1.5×1.1 deg 2 tiles, each of which comprises six unfilled pawprint stack images.The tiling pattern is such that most parts of the tile are observed by at least two pawprint stack images within a few minutes, with the two or more images of each star falling on different arrays or widely separated parts of the same array.VVV typically has between 50 and 80 epochs of observation for each tile in   , with additional observations in the ,  ,  and  filters at the beginning and the end of the survey.The   observing blocks in the VISTA observing queue each contained a group of four adjacent, slightly overlapping tiles that could be observed in about an hour.From 2016-2022, VVVX extended the survey to a larger area, whilst continuing to take a few more epochs of   data in the original VVV area, supplemented by some deeper  and  observations in the inner bulge.Ultimately, most sources in the VVV area have light curves comprising ∼200   pawprint stack-based measurements over the 2010-2019 interval.
The VISTA   pawprint stack images have a typical spatial resolution of 0.75 ′′ but data from images with much poorer seeing are retained in the VIRAC2 light curves.The sensitivity limit is   ≈ 17.5 (Vega system) in relatively uncrowded fields and the saturation limit is in the range   = 11-12, with some variation depending on the seeing, the array detector, the brightness of the background and the transparency.Brighter sources are often included in the database, albeit with relatively few epochs of data, taken in poor seeing conditions.

VIRAC2
Our searches for highly variable stars and transients mainly used two preliminary versions of the VIRAC2 database, VIRAC2- and VIRAC2-, which contain time series profile fitting photometry taken over 2010-2018 and 2010-2019, respectively.The light curves presented herein are almost all from VIRAC2-, supplemented by additional saturation-corrected photometry and photometry on stacked images in cases where this was needed for bright or very faint sources lying outside the single-epoch dynamic range of the survey.VIRAC2- light curves are presented only for a few poorly sampled transients that were were not included in VIRAC2- because they were detected in fewer than ten pawprint stack images.The dophot software (Schechter et al. 1993;Alonso-García et al. 2012) was employed for the photometry on each pawprint stack image.Photometric and astrometric calibration were performed as part of the production of VIRAC2 (Smith et al., in prep.), see the summary in Appendix A.

VVV 4th Data Release
An earlier search employed the SQL database of VVV tile-based aperture photometry available in the 4th VVV Data Release (VVV DR4), hosted at the VISTA Science Archive (Cross et al. 2012, vsa.roe.ac.uk), supplemented by further analysis using the VVV pawprint stack-based aperture photometry for each candidate variable source (see Appendix A).VVV DR4 comprises aperture photometry taken in the 2010-2013 time period, processed with version 1.3 of the casutools VISTA data reduction pipeline (González-Fernández et al. 2018).The 20 sources having Δ  > 4 mag in the VVV DR4 light curves were all recovered by the later VIRAC2-based searches.However, the VVV DR4 search was extended to lower amplitudes (Δ  > 3 mag), yielding a somewhat larger number of discoveries (see §3.1).

Search methods
Details of how the above databases were searched for bona fide large amplitude variable stars and transients are given in Appendix A.

OVERVIEW OF RESULTS
In Table 1a we present the list of 222 bona fide variable stars and transients with Δ  > 4 mag discovered by our VIRAC2-based searches.Light curves for all these sources are provided in catalogue and PDF forms in the Supplementary Online Information.Of these, 162 (73 per cent) are VVV discoveries, ∼14 per cent of which were previously reported in VVV-based searches for novae (Saito et al. 2015(Saito et al. , 2016a;;Montenegro et al. 2015;Gutierrez et al. 2016), highly variable stars of all types (CP17a, Medina et al. 2018) and periodic variable stars (Guo et al. 2022) 1, 2 .Table 1a gives a running number for each source and any other name by which the source is known, the VIRAC2 J2000 equatorial coordinates of each source, a measurement of the amplitude of variability (Δ  ), the median   magnitude, an initial classification of the type of astrophysical source, a description of the light curve, a comment, and the time baseline in years over which the source was detected in the VVV/VVVX surveys.In Table 1a, Δ  is computed using flux measurements that are averaged within 1 day time bins 3 to reduce measurement error, after rejecting any measurements with discrepant VIRAC2 astrometry or a poor dophot fit to a stellar profile (see Appendix A) and removing any measurements that differ by more than 3 from the mean in each bin.A small number of sources have Δ  just under 4 mag as measured in the binned light curves (see Appendix B for details) since the sparse sampling and finite sensitivity of VVV would tend to cause the true amplitude of most sources to be slightly under-estimated.Lower limits on Δ  are given for sources that dropped out in part of the survey; these are based on upper limits on source flux in images constructed by stacking images in time groups, grouped by the calendar year of observation.
classified "VIVACE" list of 1.4 million periodic VVV variable star candidates (Molnar et al. 2022), classified as LPVs in almost all cases.However, 72 per cent of the 127 (i.e.92 sources) are not in fact periodic variable stars.This suggests that the VIVACE classification procedures had some difficulty with aperiodic variable sources, which were not represented in the training sets used in that work.The authors of that work noted that some degree of contamination by aperiodic YSOs was expected in their catalogue; it now appears that other types of aperiodic variable source such as the novae and microlensing events reported herein can also appear as contaminants, though YSOs are likely to be the most numerous type. 2 For completeness, we note that 133 of the 222 variable sources are included in a list of 45 million VVV variable star candidates published by (Ferreira Lopes et al. 2020), based on the aperture photometry in VVV DR4 (2010 to 2013 time series).Bona fide variable stars are not clearly distinguished from the far more numerous false positives in that database, nor given astrophysical classifications, but it provides a basis for searches, like VIRAC2. 3 For source 185, Δ  was computed with the unbinned light curve due to large intra-night variability of this brief transient.The initial classification of each source is based on an assessment of all available data, beginning with the light curve morphology, near infrared and mid infrared colours (the latter from the Spitzer Galactic Legacy Infrared Mid-plane Surveys, GLIMPSE I-II , Benjamin et al. 2003;Churchwell et al. 2009, and 2017), but we treat the results only as one indicator amongst several when assigning initial classifications.
A period search was performed with elements of a software pipeline developed to search the entire VIRAC2- database for periodic sources (Miller et al., in prep).For sources in the present work, the pipelien uses three period search methods after first removing any overall linear trend: Phase Dispersion Minimisation (Stellingwerf 1978, specifically the PDM2 code), Lomb Scargle (Scargle 1982), Conditional Entropy (Graham et al. 2013) and selects the one that best fits the data.For this work, periods were confirmed by visual inspection of the light curves and independently confirmed by a novel machine learning procedure that computes a false alarm probability (Miller et al., submitted) after being trained to distinguish between aperiodic and periodic light curves.Periods of LPVs are given in Appendix F and the one other periodic source is discussed in Appendix E.
Our initial classifications break down as follows: 70 CVs, 40 YSOs, 38 microlensing events, 35 LPVs, 21 sources that we call "dipping giants", 10 sources that we classify as "unusual" and eight transients denoted "sparse", where there are too few detections over time to attempt an astrophysical classification.We discuss each of these classes below in turn.To note, there is relatively little useful optical time series data to aid classification of the new discoveries, from Gaia (Eyer et al. 2023) or other public optical time domain surveys, owing to the red nature of the sample.Most of the 50 Gaia DR3 matches are novae (28 sources, previously published in all but one case), microlenses (11 sources, none of which have multi-epoch photometry in Gaia DR3), or YSOs (7 sources, among which only source 78 has a Gaia DR3 light curve, presented in a single object study, Guo et al., submitted).Additionally, source 192 had a Gaia Alert, Gaia 17bzk, though only the fading portion of this eruptive YSO event was sampled and the star is not in Gaia DR3.
Figure 1 shows the spatial distribution of all the sources in Galactic coordinates, with the VVV survey area indicated.The YSOs are fairly uniformly spread along the mid-plane of the inner Galactic disc.Most YSOs (35/40) lie at Galactic latitudes || < 1 • , as expected for a young population with a small Galactic scale height.By contrast, CVs are more concentrated in the VVV bulge region, especially the inner bulge.Microlenses are even more strongly concentrated in the bulge, consistent with past optical surveys such as OGLE (Wyrzykowski et al. 2015) and prior work using VVV (Navarro et al. 2018(Navarro et al. , 2020a,b),b).
The dipping giants stand out as a very tight cluster of 21 points near the Galactic centre, most of them projected in the Central Molecular Zone and the Nuclear Disc of the Milky Way that lie between longitudes 358 • and 2 • .At first, there were three possible interpretations of these sources: YSOs with long-lasting dips like AA Tau (Bouvier et al. 2013;Covey et al. 2021), giant stars undergoing unusually strong dimming by circumstellar dust, or eruptive YSOs (in cases where the VVV light curve morphology was ambiguous).However, spectroscopy of seven of these sources (see GLK) and use of photometry from earlier surveys (see §5) helped to distinguish giant stars from YSOs and resolve the nature of ambiguous light curves: such sources are assumed to be undergoing a dip in VVV if they were bright in earlier surveys.In some cases there is supporting evidence of variable reddening from the VVV multi-filter observations.Moreover, the recognition that the Central Molecular Zone has a very high surface density of molecular cloud cores and clumps (e.g.Parsons et al. 2018) helped us to remove YSO status for some sources lacking spectroscopic data that are projected adjacent to one or more cloud cores.
Given the spectroscopic evidence for evolved status of seven sources we opt for a giant star classification for dipping sources near the Galactic centre that presently lack spectra in all but one case (source 134) where YSO status is supported by a previous study, see §4.3.In total we classify eight sources with long-lassting dips as YSOs but only source 134 is projected in the Nuclear Disc.The dipping giants are discussed further in §5 and in GLK.
The LPVs are spread along the Galactic mid-plane but with a greater concentration towards the inner Galaxy than the YSOs.This is consistent the Akari-based study of Ishihara et al. (2011), which was sensitive to LPVs at all Galactocentric radii and quantified the high concentration of O-rich LPVs in the inner Galaxy, where the more uniformly spread C-rich LPVs are greatly outnumbered.
Interestingly, the eight sparsely sampled transients (see lower right panel of Figure 1) are spread along the mid-plane (7/8 lie at || < 1 • ) in a manner that appears different to the bulge-centric CV and microlens spatial distributions.This suggests that not all of these sources are poorly sampled examples of those two classes of variable sources.Finally, the 10 "unusual" sources are a very diverse group (see Appendix I) so their distribution has little meaning.
Table 1b lists two eruptive YSOs with amplitudes in the 3 to 4 mag range that were removed from Table 1a at a fairly late stage, after undertaking spectroscopic follow up (see GLK), when it was recognised that the amplitudes of the light curves were significantly below the 4 mag threshold after 1-day binning with 3 rejection of outliers.In the case of VVV J164011.76-484653.4, hereafter VVV 1640-4846, the YSO is embedded in a compact nebula so the VIRAC2- photometry were replaced with aperture photometry (aperture diameter 1.414 ′′ ) using a carefully selected background annulus.

DR4 sources with Δ𝐾 𝑠 > 3
Our search of the public DR4 release for sources with Δ  > 3 mag (see §2.3) yielded 105 sources.These 105 sources are listed in Appendix C. They are denoted with the prefix "DR4_v", consistent with Guo et al. (2021).Of the 105, the 20 with Δ  > 4 in the DR4 aperture photometry data were all recovered in the VIRAC2-based searches.A further 15 have Δ  > 4 in the deeper VIRAC2 dophotbased data so 35 of the DR4 sources are also listed in Table 1a.We very briefly discuss the 70 DR4 sources with   = 3 to 4 mag in order to provide some information about the types of sources seen, since a general Δ  > 3 mag search has not yet been done with VIRAC2.Our initial classifications break down as follows: 23 LPVs, 22 YSOs, 1 CV (nova candidate VVV-NOV-012, Saito et al. 2016a), 15 microlenses, five listed as unusual and four in the sparse category.We see that, as expected, there are larger proportions of YSOs and LPVs and a much smaller proportion of CVs with measured amplitudes in this range than in the VIRAC2 Δ  > 4 sample.Spectroscopic follow up and analysis of 15 eruptive YSO candidates from this DR4 set was undertaken in 2019 (see Guo et al. 2021).Three eruptive YSO candidates in the set were noted earlier in CP17a and followed up in CP17b.

YOUNG STELLAR OBJECTS
The 40 YSOs in Table 1a, comprise 32 eruptive YSOs and eight systems with dips that appear to be due to variable extinction.Our main focus in this work, and in GLK, is the eruptive YSOs.The 32 eruptive systems have some variety in light curve morphology but 20 (63%) can be broadly described as "classic" eruptive YSO events, see Figure 2.These have a reasonably smooth rise from a faint state to a bright state that then lasts for at least several years, i.e. continuing after the end of the VVV survey, though some fade by > 1 mag 56000 57000 58000 Modified Julian Date during that time 4 .Nine eruptive systems have more irregular and diverse light curves, see § 4.2.
A further three eruptive systems may well have been of the classic type: these are cases where the outburst reached maximum in 2010, or earlier, so there was little or no sampling of the rise by VVV and the high amplitude is due to fading of the outburst by several mag within 4 Source 165 had a smaller burst two years before the main event and source 15 shows slightly more complicated behaviour after the initial rising stage but we include these in the "classic outburst" group despite these differences.
the 2010 to 2019 time period.Such a rapid rate of decay is faster than the three historical sources often called the "classical FUors" (FU Ori, V1057 Cyg and V1515 Cyg, see e.g.Hartmann & Kenyon 1996) but within the range of the 20 classic VVV events where the rise was observed, as shown by recent photometry in GLK.These three systems are source 19 (= WISEA J142238.82-611553.7, Lucas et al. 2020, noted   1b, comprise one classic eruptive system, VVV 1640-4846, and one irregular system, VVV J163637.94-474444.1 (hereafter VVV 1636-4744).The latter is one of the few irregular systems having spectroscopic follow up (GLK) and our analysis of classic eruptions includes several other spectroscopically confirmed lower amplitude events so these two are discussed along with the rest.

Classic outbursts
A key finding of this work is that high amplitude outbursts with a long duration typically have a slow rise over ∼2 yr or more, contrasting with the classical FUor systems FU Ori and V1057 Cyg where the rise time was only about 6 months (FU Ori, see Kenyon et al. 2000) or 13 months (V1057 Cyg Herbig 1977).By contrast, the third classical FUor, V1515 Cyg, had a uniquely slow rise time of at least 15 years, perhaps longer (see Herbig 1977), considerably slower than any in this work or any other published FUor.In Figure 3 (left panel) we plot the rise time of the outbursts against amplitude for the 20 classic outbursts in Table 1a, and a further 12 spectroscopically confirmed VVV discoveries with classic morphology and amplitudes in the 2 to 4 mag range, drawn from Table 1b, Appendix C and CP17b, see Guo et al. (2021Guo et al. ( , 2020) ) and CP17b for spectra and VIRAC2 light curves of the Appendix C sources and earlier discoveries.Here we measure the rise time from the beginning of the outburst until the bright plateau or turnover is reached, neglecting any short timescale/low amplitude scatter in the faint state or the bright state.Due to the sparse VVV sampling there is some uncertainty in the rise times but the error bars encompass the plausible range of values based on visual inspection.The numerical data are given in Table 2, along with a fitted rise timescale,    , see below.
We see that the typical rise time is 600 to 1000 d, though in some cases it is ∼2000 d or more and in a smaller number of cases it is only ∼200 to 400 d.
We inspected the literature on long duration outbursts, for which a helpful list was compiled in Connelley & Reipurth (2018).While most past eruptions have little photometry during the rise to maximum, useful (optical) data exist for V2493 Cyg (=HBC722), V2775 Ori (=CTF93 216-2), V582 Aur and V960 Mon, in addition to the three classical FUors mentioned above.The rise times from quiescence to the bright state are: 2 months (V2493 Cyg), < 18 months  The "Rise Time" is based on visual inspection, measuring the total time from the start of an outburst until photometric maximum, discounting short-timescale scatter.   is a fitted rise time that encompasses the bulk of the rise in brightness but excludes the asymptotic tail that occurs at the beginning, and sometimes at the end also.   is only given if a useful fit was possible It is possible that the slower rise times seen in the near infrared reflect a difference between optical and infrared measurements, if the accretion burst is detected at an earlier stage in the infrared, at lower disc temperatures, and only later becomes obvious in the optical as the temperature rises further.For example, this occurs naturally in disc models where FUor outbursts are due to sudden triggering of magneto-rotational instability in the dead zone and the outburst then propagates inward (Cleaver et al. 2023).However, for the more recent FUors Gaia 17bpi and Gaia 18dvy (Hillenbrand et al. 2018;Szegedi-Elek et al. 2020), inspection of the Gaia Alerts data shows that these had total optical rise times of ∼900-1200 d (Gaia 17bpi) and almost 700 d (Gaia 18dvy), in the same range as the near infrared data for the VVV sample.
In view of these two recent Gaia discoveries, we cannot be certain that the slow infrared rise times seen in the VVV sample are different to those of optically bright, optically measured events.The slightly faster rise times of the six historical outbursts mentioned above may simply reflect small number statistics rather than a real distinction between optical and infrared measurements, or a difference between optically bright events and embedded systems.
In this context it is perhaps significant that there appears to have been a slower rise in the mid-infrared than the optical in Gaia 17bpi and Gaia 18dvy (Hillenbrand et al. 2018;Szegedi-Elek et al. 2020), interpreted as evidence for inward propagation of the disc instability.At present, we tend to the view that either (i) there is a genuine difference in optical and infrared rise timescales in individual eruptive YSOs due to the different temperature sensitivities of the wavebands, or (ii) the faster optical timescales in most historical outbursts are a fluke of small number statistics.Whilst the embedded nature of the VVV systems might also be significant, it is unclear whether these sources are truly at an earlier evolutionary stage than optically bright FUors.Historical FUors are usually associated with large amounts of circumstellar matter, seen in the far infrared and mm wavebands (Green et al. 2013;Fehér et al. 2017) so their less embedded nature may simply be due to the lower foreground interstellar extinction towards nearby star forming regions or better clearing of circumstellar material along the individual lines of sight by the stellar wind.
The spectral index, , measured at a time prior to the outburst for 15 systems is given in Table 2 and Figure 4.This was measured using either 3-22 m WISE photometry from 2010 (preferred if available and pre-outburst) or 3-24 m Spitzer photometry otherwise.5(VVV   photometry were excluded in order to reduce the potentially significant effect of foreground extinction in the Galactic plane, see CP17b).We also include the index for a sixteenth system, source 37, which appears to have been in the early stages of the rise at the time of the measurement in 2010.The mean spectral index is  = 0.16, or 0.13 if source 37 is excluded, and (14/16) systems fall in the Class I category ( > 0.3) or Flat Spectrum category (-0.3 <  < 0.3), see Greene et al. (1994).Most historical FUors lack a measurement of the spectral index in the quiescent state, preventing a direct comparison between optically bright and embedded systems.These data are therefore a useful indicator of the evolutionary stage of the progenitors of classic outbursts, for cases where the system is not too deeply veiled to be detected in the near infrared.The systems lacking a spectral index due to non-detection in the [24] and 4 passbands were in most cases veiled from view by strong nebulosity gradients or blending, so their spectral energy distributions (SEDs) cannot be assumed to be bluer.
The morphology of the   light curves in the rising phase of classic outbursts can be described as a rise that begins slowly, accelerates to a faster, approximately constant rate of increase and then changes slope abruptly as the bright state is reached.There is some variation here: a minority of light curves show a gradual decline in gradient as the rise approaches the bright state, rather than an abrupt change.An example of this "S-shaped curve of growth" (to borrow a term from economics) is found in source 10, which rises asymptotically towards maximum light over several years (making it hard to define when the bright state is reached).Other examples are source 165 and perhaps sources 6, 95 and VVV 1640-4846, though the sampling is too sparse for certainty in these last three cases.Gaia 18dvy has a sufficiently well-sampled Gaia light curve to inform this topic: the outburst has the characteristic slow start over ∼400 d, followed by a fast linear rise to within ∼0.5 mag of peak brightness over ∼200 d, but the final stage of the rise then proceeds at a more gradual rate over ∼ 100 d.Given the sparse VVV sampling, a similar brightening profile to Gaia 18dvy is certainly possible in many cases.Finally, we note that source 167 lacks the slow start and shows an approximately linear rise from quiescence to the bright state.
In Figure 3 (right panel) we plot the rise time against total duration, though almost all of these outbursts were ongoing at the end of the VIRAC2 time series so the durations are usually lower limits.We see that the outbursts typically have a total duration of at least several years, though the decay rates vary considerably between systems (discussed further in GLK).In terms of duration then, these high amplitude outbursts appear to match the expectation for FUor events.In GLK we confirm spectroscopically that most of the highest amplitude classic events are indeed FUors, the remainder being eruptive YSOs of other types.Only sources 3 and 41 presently lack spectra.The companion paper by Contreras Peña et al. (submitted) includes a small subset of the high amplitude sample presented here, likewise noting their long outburst durations in comparison with lower amplitude events, which have a spread of durations (as seen in Guo et al. 2021 also).
For 18 sources in Table 2 where the full rise was adequately sampled from quiescence, we fitted the source magnitude, (), during the rising portion of the outbursts with the following formalism: where   is the quiescent magnitude,  is the amplitude of the outburst,  0 is the time at which the source has brightened by 0.5 and  is the e-folding timescale.With this form, the bright state is reached at  =  0 + 2, at which point, () =   − .To aid the fitting of  and  0 , we fixed   and  and discarded data taken after the end of the rise and data taken more than 1-2 yr before the rise began (which otherwise had a tendency to dominate the fit due to denser VVV sampling in 2012-2013 than later years).Levenberg Marquardt fitting was performed with the Python lmfit package and the piecewise routine in Numpy.Sources that were better described by an S-shaped curve of growth (sources 10, 165, 6 and 95) were fitted using equation 1a for all times, up to the end of the VIRAC2- time series.Equation 1a is simply a hyperbolic tangent, scaled by the outburst amplitude and shifted appropriately.Examples of these fits are shown in Figure 2.
Figure 5 shows the distribution of rise times resulting from these fits, plotted as a histogram of    = 4.Using the form in equations 1a and 1b, the time interval    , from  0 − 2 <  <  0 + 2, encompasses most of the rise in flux, specifically 0.88.When using equation 1a only, for an S-shaped rise, this time interval corresponds to 0.76, whereas for a linear rise,    corresponds to the full rise, .We see that the timescale represented by    is typically between 400 and 1000 days, with the peak at 600 to 800 days.The distribution and this typical timescale are broadly consistent with the 600-1000 day typical rise time found by visual inspection, given that (i) the full rise does take longer than    and (ii) three of the slowest rises plotted in Figure 3 are lower limits and therefore not fitted.Rise times longer than 1000 days are more significantly under-represented by    but since this parameter encompasses most of the change in brightness, aside from the slowest-changing portions of the S-shape, it may be of use to theorists.GLK use piecewise linear fits to describe the full light curves and come to similar conclusions about the rise timescale.
The faint pre-outburst magnitudes shown in Figure 2 allow us to infer the relevant range of stellar masses.The median   = 16.5 and source distances given in GLK for classic outbursts are usually between 2 and 8 kpc according to two separate methods (see Tables 1 and 2 of that work).Absolute magnitudes then lie mostly in the range 2 <    < 5. Comparison with YSOs with well measured SEDs, similar spectral indices (−0.3 <  < 0.6) and absolute magnitudes in the "Cores to Disks" sample of Evans et al. (2009) suggests that extinction-corrected bolometric luminosities are probably in the range 0.2 <  bol / ⊙ < 20.This corresponds to low mass YSOs: if the quiescent luminosity is dominated by the star rather than accretion, the range of masses is approximately 0.2 < / ⊙ < 3, based on the 0.5 to 1 Myr isochrones of Baraffe et al. (2015) andFeiden (2016).If accretion dominates then the masses would be lower.
The effect of FUor outbursts on source SEDs has been modelled in detail by MacFarlane et al. ( 2019), Rodriguez & Hillenbrand (2022), Hillenbrand & Rodriguez (2022) and Liu et al. (2022).The first of these studies performed a full hydrodynamic simulation and radiative transfer modelling, focussing on class 0 YSOs, while the latter three studies explore a larger parameter space and focus on optical to mid-infrared observations.These works show that the outburst amplitude behaviour depends on many factors, such as the properties of the central star, the temperature structure of the disk, and the envelope density structure.Hillenbrand et al. (2022) find a dependence of the amplitude of the outburst with wavelength, with optical amplitudes being typically larger than infrared ones, in agreement with observations (Hillenbrand et al. 2018;Szegedi-Elek et al. 2020).Liu et al. (2022) model the full range of stellar masses mentioned above, and quantify the changes in 1 and 1 − 2, finding a smaller difference between the infrared and optical amplitudes.The authors attribute this to the assumption of a dipole magnetic field, which leads to a fainter quiescent disk and larger mid-infrared amplitudes.Nevertheless the amplitudes in their models depend on the maximum accretion rate reached during the outburst and the mass of the central star.From Figure 2, only sources 1, 6, 10, 18 and 37 have well-sampled WISE/unTimely light curves (Meisner et al. 2023) and these typically have Δ1 ≥ 3 mag (not shown).Inspection of Figure 16 in Liu et al. (2022) then suggests at least a 1.5-2 dex increase in accretion rate.
Unfortunately, the bolometric luminosity in outburst is not very well constrained for the VVV sources by the available VVV and WISE photometry at  ≤ 4.6 m and the somewhat uncertain individual distances.Moreover, most of these sources have pre-outburst 1 − 2 colours that are considerably redder than predicted and change (become bluer) by a larger amount than predicted, a behaviour also seen in the embedded FUor SPICY 97855 (Contreras Peña et al. 2023).Liu et al. (2022) note that the 1 − 2 colours of FUors are sometimes considerably redder than predicted, which they attribute to various effects of the circumstellar envelope that were not included in the model.High and variable extinction alone cannot resolve these issues in the VVV sources, given the measured  −   colours (see GLK).Some recently discovered embedded FUors actually become redder when brighter during outburst (listed in Contreras Peña et   2023).In view of these difficulties, we do not attempt to quantify the accretion rate in outburst or analyse the systems fully at this stage.

Irregular Eruptive Behaviour
In Figure 6 we show the light curves of nine eruptive variable YSOs that display irregular outbursts and three (in the bottom row) that show fading behaviour within VVV because the start of the outburst was missed.Eight of the nine irregular outbursts in Figure 6 are drawn from Table 1a and one slightly lower amplitude system, VVV 1636-4744, is drawn from Table 1b.A tenth irregular outburst, in source 148, is not shown here.That system displays strong periodic oscillations atop a rising trend, somewhat resembling the recent outburst of LkH 225 South, also known as V1318 Cyg South (Hillenbrand et al. 2022;Magakian et al. 2019).Source 148 is discussed further in GLK, where a spectrum confirming its nature as an eruptive YSO is presented.
The light curves of the six sources plotted in the upper two rows of Figure 6 (sources 32, VVV 1636-4744, 55, 72, 77 and 164) have some visual resemblance to each other.There are signs of two or more photometric outbursts, or at least two times when each source is relatively bright (though this is unclear in source 77) and significant variability on short timescales (< 100 d).None of these six sources have a detection in 2MASS PSC, indicating that they were relatively faint at the time of that survey.This supports the view that the variability is mainly due to accretion driven luminosity increases from a faint state rather than a reduction in luminosity due to variable extinction.There is no useful data on the near infrared colour changes of these sources (all of which have red 1 to 5 m SEDs): in most cases VVV detected them only in the  and   filters, with only one  band detection and non-detections at other times when they were fainter.Next we give some detail on individual sources, most of which are not discussed in GLK, before analysing the timescales of variability collectively.

Individual sources
VVV 1636-4744 and source 72.We consider the eruptive status of these two stars to be doubtful.In VVV 1636-4744 this is because the 2MASS images shows signs of an uncatalogued, heavily blended detection at   ≈ 14-14.5, slightly brighter than the mean brightness level for this source, based on visual comparison with adjacent sources.The source has a spectrum in GLK, observed in an intermediate state of brightness with a modest signal to noise ratio.It shows fairly strong Br emission and a red continuum, consistent with a veiled YSO undergoing magnetically controlled accretion, but no CO emission.Source 72 simply has an ambiguous VVV   light curve morphology.There is no spectrum available and little additional information in other bandpasses.
Sources 32, 55, 77 and 164.Source 32 has an "outflow dominated" spectrum (see GLK) 6 , with exceptionally strong H 2 emission lines.Sources 55, 77 and 164 lack spectra but the WISE unTimely and NE-OWISE data (Meisner et al. 2023;Mainzer et al. 2014) are useful for the latter two sources and for source 32, see Appendix D. Sources 32 and 77 become bluer when they are brighter, with magnitude changes approximately in the ratio Δ2/Δ  = 0.25 and 0.52, respectively (with some uncertainty since the measurements are not quite contemporaneous).Source 164 becomes redder when it is brighter, with Δ2/Δ  ≈ 1.7.According to Wang & Chen (2019) the expected ratio in the case of variable extinction is Δ2/Δ  = 0.33, based on measurements of Cepheids with high extinction (Chen et al. 2018).Only source 32 has colour variations close to this ratio and the exceptionally strong H 2 emission indicates the presence of a strong molecular outflow or disc wind, implying a high time-averaged accretion rate.Source 55 is blended with brighter neighbours in the ∼6 ′′ WISE beam (as are VVV 1636-4744 and source 72) but in this case the light curve has a convincing eruptive morphology.Hence the available data tend to support an accretion driven origin for the variability of these four systems.
Sources 42 and 48.These two stars display only relatively brief outbursts.In VVV the decline timescales were ∼100 d but the rise and the photometric maximum were not sampled.For source 42, WISE provides data for the second outburst sampled by VVV, at MJD=57925 (see Appendix D) indicating a duration of 1 yr.There is evidence for earlier outbursts of source 42: the 2MASS survey measured   = 13.52 in 1999 and the Spitzer GLIMPSE survey in 2003 measured [4.5] = 9.58 at 4.5 m, similar to the WISE/unTimely measurement 2 = 9.48 at 4.6 m during the outburst in 2017.During the decline of the first outburst sampled by VVV there was a brief ∼0.5 mag re-brightening (see inset in Figure 6).In source 48, four outbursts are seen in the WISE unTimely data, each observed at a single unTimely epoch (two with ∼2 mag amplitude in 2 and two with ∼1 mag amplitude).This indicates that the durations were of order 6 months.Source 55 may also have outburst durations of only 6 to 12 months but the sparse sampling makes this uncertain.Among previous VVV discoveries, only VVVv118 shows somewhat similar behaviour (CP17b, Guo et al. 2020) though in that case the duration of the outbursts was only ∼50 d.
Source 120 (=DR4_v67).This star shows rapid variability, in addition to larger inter-year variability, resembling that seen in the "multiple timescale variables" (MTVs) reported in Guo et al. (2020Guo et al. ( , 2021)).Spectra shown in Guo et al. (2021) and GSK are outflow dominated, with unusually strong H 2 emission.Comparison of the VVV and WISE light curves, shown in Appendix D, shows bluer-when-brighter behaviour with ratios Δ2/Δ  = 0.48 and Δ1/Δ  = 0.75, somewhat larger than would be expected for variable extinction by interstellar dust.

Timescales
Figure 7 illustrates the variability of the eruptive YSOs as a function of time baseline, comparing 21 classic outbursts (grey dashed lines) with the six irregular eruptive systems in the upper two panels of Figure 6 (solid blue lines) and sources 42 and 48 (solid red lines).These plots are constructed by computing the r.m.s.variability (upper panel) or maximum variability (lower panel) over all pairs of   points,   ,   , in each light curve and binning the results in 0.5 dex bins in the logarithm of time baseline.
The r.m.s.variability was computed using the unbinned pawprint light curves and the uncertainty in each value of (  −   ) 2 was subtracted from the mean square variability in order to debias the result.The r.m.s.variability of source magnitudes is the square of the structure function as defined in Lakeland & Naylor (2022), though some authors have preferred to use flux differences, see e.g. the YSO optical variability studies of Sergison et al. (2020), Venuti et al. (2021).An approach based on magnitudes, i.e. flux ratios, is more suitable than flux differences for large amplitude variability, as discussed in Lakeland & Naylor (2022).The maximum variability vs. time baseline (also called the "accumulation function", see Guo et al. 2020) is a useful alternative for the sparsely sampled VVV light curves.This was computed using the 1-day binned light curves plotted in this work, in order to reduce the effect of individual 2-3 outliers.
We see that the six irregular systems and the two fast eruptive systems tend to have more variability than classic outbursts on timescales up to 100 d.The irregular outbursts generally show maximum brightness changes of 1 mag or more on a 100 d timescale (Figure 7, lower panel).The variability of the classic outbursts rises sharply on baselines from 100-1000 d.Source 42 displays the largest r.m.s.variability on 2 d to 300 d time baselines, though this is partly due to the good sampling of the first outburst seen by VVV.Source 93 stands out slightly from the other classic outbursts over 30 d to 100 d time baselines due to the unusually fast rise time of the event (see Figure 2).
YSOs having irregular high amplitude accretion-driven behaviour clearly require further study and spectroscopic follow up, as do the MTVs reported in Guo et al. ( 2020 clearer.There appear to be two main types of these irregular light curve: those with repeating outbursts that last 1 yr or less (sources 42, 48, VVVv118 and perhaps source 55) and those those with both intra-year and inter-year variability.

YSO Extinction Events
In eight YSO candidates there is good evidence that the variability is due to changing extinction.Lucas et al. 2008) finds that 18/21 sources were detected in a bright state by one of the two surveys some years prior to the VVV survey (i.e.close to or brighter than the brightest VVV   measurements).This helps to verify that we are indeed observing dipping events rather than photometric outbursts.N.B. non-detection in 2MASS PSC is sometimes due to source confusion in this region, which also increases scatter in the photometry of detected sources.2MASS PSC typically reports magnitudes slightly brighter than either UKIDSS or the brightest VVV measurement for the same source in this set, most likely due to blending.Only 16/21 sources lie in the UKIDSS footprint but 14/16 were in a bright state in 2006 or 2007.
GLK present near infrared spectra of seven of the dipping giants (sources 130, 144, 145, 149, 154, 168, 172), six of which were taken during the bright state after the dip had ended.All of these are latetype giant spectra suffering high infrared extinction, such that high quality data were obtained only at  > 2 m.The initial analysis in that work finds effective temperatures from 3100 K<     < 4100 K and suggests that these are O-rich stars: there was no sign of C-rich features such as the 1.77 m C 2 absorption band, with the caveat that the spectra are of poorer quality in the  bandpass.The high foreground extinction towards the 21 sources is further evidenced by their very red ( −   ) colours.In Figure 10 blue points are plotted for 20 sources, corresponding to data taken near maximum brightness within the VIRAC2 time series.The 11 red points correspond to data taken at the faintest available multi-filter detection in each light curve, the remaining sources lacking a second detection in .A reddening vector is also plotted, based on the Wang   Using the Wang & Chen (2019) extinction law, the extinction towards these 21 sources outside the dips is typically 4 to 5 mag in   , based on observed colours 3 < -  < 4 (see Figure 10) and intrinsic -  ≈ 0.2 to 0.4 for post-main sequence stars at these temperatures, according to the PARSEC-COLIBRI isochrones (Chen et al. 2015;Marigo et al. 2017).This is supported by fitting the dereddened spectra to comparison giant stars from the Xshooter Spectral Library (Verro et al. 2022), showing no obvious contribution of emission by circumstellar dust to the red continuum (to be discussed further in a future work).Observed magnitudes are   ≈ 11 to 13 mag outside the dips.After correcting for extinction and adopting an 8 kpc distance for the Nuclear Disc, the distance modulus is 14.5 and the absolute magnitudes are in the approximate range −8 <    < −6.From inspection of the PARSEC-COLIBRI isochrones, the range of    and     is consistent with an AGB interpretation.
Mid-IR colours from Spitzer/GLIMPSE-II and GLIMPSE-3D, where available, are plotted in Figure 11.These were measured prior to the dips observed by VVV.The very high foreground extinction prevents us from distinguishing O-rich and C-rich AGB stars via the colours but we see that only a minority (6/19) of the sources in the left panel have [3.6]− [4.5]> 1.The very red colours of this smaller subset suggest the presence of hot dust close to the photosphere.However, no source has [4.5] − [8.0]> 1 (the threshold for extremely dusty AGB stars and other intrinsically red sources listed in Robitaille et al. 2008).None are detected in the Spitzer MIPSGAL survey or the far-infrared Herschel Hi-Gal survey (Molinari et al. 2010;Elia et al. 2017), and there is only one detection at 22 m in WISE (source 91, 4=1.2).This indicates that there is no large mass of cooler dust in almost all cases.There is little other useful data from the WISE satellite and the associated unTimely light curves, owing to severe source confusion.
The VIRAC2 proper motions in Galactic coordinates have a cluster around   = −7 mas/yr,   = −3 mas/yr (not shown).The  and  components of Galactic velocity can be closely approximated by adopting the Galactic centre sky coordinates for these sources so that the velocity components could be determined simply from the proper motions, assuming  = 8 kpc.The  component (parallel to Galactic rotation in the solar neighbourhood) is then usually between 200 and 300 km/s, suggesting locations on the near side of the Nuclear Disc and consistent with typical Nuclear Disc stellar motions measured by VIRAC2 (Sormani et al. 2022) for the Nuclear Disc stars in the spectroscopic sample of Fritz et al. (2021).However, some sources have highly discrepant velocities.On investigation, this appears to be due to a rare type of astrometric error that affects only highly variable sources that are in very crowded fields, or otherwise blended with adjacent stars.A comparison with Hubble Space Telescope astrometry has confirmed that the uncertainties on VIRAC2 proper motion and parallax measurements are usually correct even in crowded fields (Luna et al. 2023), partly because the error calculations include a comparison with Gaia DR3 rather than internal errors alone (Smith et al., in prep).However, the effect of large changes in source flux is to change the number of stars in a blend that are detected by the dophot software, which can cause a systematic shift in source coordinates in sources like dipping giants that are faint for a significant fraction of the light curve.
Variability due to changing extinction by circumstellar matter is often seen in pulsating dusty Mira variables.In these stars, mass loss is attributed to radiative pressure on dust grains that form in the outer layers of the stellar atmosphere (and above the photosphere) and drag the gas along with them.Pulsations and convection are thought to assist the process by providing an initial upward motion and causing shocks to form, creating high density regions that enhance dust formation (Höfner & Olofsson 2018).This process is inherently three-dimensional, a fact confirmed by recent simulations (Freytag & Höfner 2023) and several high resolution observations at infrared and radio wavelengths (see e.g.Adam & Ohnaka 2019;Matthews 2023).Consequently, ejection of "dust puffs" can be directional, providing an explanation for the deep aperiodic dips occasionally seen in Mira variables that lack the very thick dust shell and the corresponding very red colours seen in OH/IR stars (the reddest Mira variables).E.g.Whitelock et al. (2006) note that such dips are seen in C-rich Miras, appearing observationally similar to the R Cor Bor events seen in H-deficient C stars, where random fading events of up to 7 mag in the visible (∼0.7 mag in   ) are attributed to ejection of puffs of C-rich dust at random times and in random directions (Clayton 1996(Clayton , 2012)).They went on to note that dust puffs had not been seen in solitary O-rich Miras but occur commonly in O-rich symbiotic Miras (Whitelock 1987), thought to be caused by interaction between the binary components.While the spectra presented in GLK do not show the emission lines that would be associated with accretion on to a compact companion, interaction with a main sequence companion is certainly a possibility to consider.
In Appendix F we highlight two Mira variables with exceptionally deep dips, also projected in the Nuclear Disc.Deep dips in infrared light curves due to circumstellar matter in non-pulsating giant stars have been seen much more rarely hitherto.Outside the nuclear disc, perhaps the best examples are found in Whitelock et al. (2006), where the carbon-rich star IRAS 09164-5349 was discussed.That system showed a slow fading trend of at least 2 mag in  after previously showing only very low amplitude semi-regular variations.Two other carbon-rich giants, IRAS 10136-5743 and IRAS 16406-1406, showed dips of 1 to 1.7 mag in  and no obvious pulsations, though the light curves were sparsely sampled.In a single 11.5×11.5 arcmin field close to the Galactic centre, aperiodic dips with ∼1 mag depth in   were reported in several apparently non-pulsating giant stars in the VVV-based study of Molina et al. (2019), along with two cases of ∼2.5 mag dips (the stars NV232 and NV261).One of the stars with a ∼1 mag dip, NV062, was spectroscopically classified as M4-5 III type (though the 2.0-2.4 m spectrum did not unambiguously distinguish O-rich and C-rich giants).The prior study of the Galactic centre region by Matsunaga et al. (2009) listed a large number of variable giant stars for which a period could not be determined, though these were not discussed.In fact, two of the 21 dipping giants in our sample, sources 137 and 146, are listed as variable stars with no known period in Table 5 of Matsunaga et al. (2009).These are presently listed in SIMBAD as [MKN2009] 37 and [MKN2009] 1313, respectively, arguably misclassified there as LPVs.
Interestingly, [MKN2009] 37 shows a 1.3 mag rising trend from 2002 to 2008 in the published light curve, at which time it reached the flat plateau at   = 11.5 seen in the VIRAC2 light curve from 2010 to 2014, prior to the deep dip that followed (see source 137 in Figure 9).In this star then, we have evidence for an earlier dip, where the minimum occurred at least 17 years prior to the minimum of the dip seen by VVV.By contrast, [MKN2009] 1313 showed only small variations in   in the 2002 to 2008 interval.
The key new feature of the present sample is the clustering within the Nuclear Disc.The metallicity distribution in this region is largely in the super-solar range 0 < [Fe/H] < 1 (Fritz et al. 2021), exceeding even the nuclear bulge of the Milky Way (Queiroz et al. 2021) which has a broader distribution.In view of this difference from other parts of the Milky Way, it seems plausible that high metallicity is aiding the mass loss e.g. by strongly enhancing production of refractory dust species.Increased dust opacity in the photosphere might also lead to more energetic convection.The analysis of the seven dipping giant spectra in GLK, based on CO and NaI indices, appears to favour high metallicity.Additional observations, searches for less extreme dipping events in the VIRAC2 database, and 3D modelling of mass loss in metal-rich atmospheres will be needed to better understand this population and distinguish the various possibilities.
An explanation related to high metallicity would indicate that these stars are O-rich, since dust production in O-rich stars is metallicity dependent, whereas dust production in C-rich stars is due to nuclear burning and dredge-ups from the interior, rather than the initial metallicity of the star (see e.g.Bladh et al. 2019).
To finish this topic, we should mention the possibility of periodic dips due to eclipses by either (i) a circumsecondary disc around a very low mass companion on a wide orbit, or (ii) matter spread around the orbit of such a companion.The high amplitude samples presented herein (Table 1a and Appendix C) include three late type giant stars displaying a deep dip attributed to a circumsecondary disc: VVV-WIT-08 (=source 199 in Table 1a), VVV-WIT-10 (=DR4_v16 in Appendix C) and VVV-WIT-11 (=source 101 in Table 1a), see Smith et al. (2021).Two of those three sources show a single symmetric dip, whilst the third (VVV-WIT-11) was undetected for a year during the dip and is regarded as a candidate of that type.None of the three are located in the Nuclear Disc.The light curves of the 21 dipping giants discussed above differ from those three sources in that the dips are not symmetric and in most cases there is significant variability outside the deep dip.Consequently, an eclipse by a circumsecondary disc can be ruled out.We also note that periodic dips (a "long secondary period") are seen in up to a third of LPVs (Wood et al. 1999;Soszyński et al. 2013;Pawlak 2023), these having a longer timescale than the pulsation period.These events are due to either option (i) or option (ii) above (Wood et al. 1999;Soszyński et al. 2021), the latter option being the favoured explanation for asymmetric and longer lasting dips.Option (ii) is a possible explanation for 21 dipping giant systems.However, the range of periods seen in Miras with a long secondary period is usually only 200 d to 1600 d (see Figure 2 of Soszyński et al. 2021), whereas the 21 dipping giant systems would have to repeat on a longer timescale, if they are periodic.A few Mira variables are known with very long secondary periods, see Appendix F so this explanation cannot be ruled out entirely.However, long secondary periods are found in Miras all over the Milky Way so it seems more likely that the explanation is related to the high metallicity of the Nuclear Disc, with no obvious requirement to invoke a companion object as yet.

Long Period Variables
The 35 sources classified as LPVs are dusty Mira variables with extremely high amplitudes, ∼5 mag from peak to trough in   in two cases.Separately, sources 163 and 147 are very interesting in the context of the Dipping Giants discussed above, since both show a very deep, apparently aperiodic dip and both are projected in the Nuclear Disc Other LPVs with dips that are slightly less deep are also of interest in this context.See Appendix F for the light curves and further discussion.

Novae and other transients
In Table 1a, 70 sources are classified as CVs.Of these, 45 have been previously classified as classical novae, recurrent novae or nova candidates and one other, source 219 (aka ASASSN-17fm) is a dwarf nova.Individual references and identifications are given in Table 1a.In Appendix G we briefly discuss four of the 24 new discoveries with unusual light curve features, such as the oscillations and cusps occasionally seen at optical wavelengths (Strope et al. 2010) but to our knowledge not previously seen in the infrared.

Unusual sources
Of the 222 sources in Table 1a, there are 10 that have fairly well sampled light curves but do not fit into the categories discussed hitherto.Two of these were known as unusual sources independently of the VVV survey (V4334 Sgr (= Sakurai's star) and the X-ray source MAXI J1348-630).The nature of two others (VVV-WIT-08 and VVV-WIT-10) is thought to be at least partly understood (Smith et al. 2021) in terms of an eclipse of a giant star by a circumsecondary disc.The nature of the remaining six is unclear at present.Each of the 10 is discussed briefly in Appendix I.

SUMMARY AND CONCLUSIONS
A thorough search of the VVV/VIRAC2 survey database covering 562 deg 2 of the Galactic bulge and adjacent Galactic disc has yielded a sample of 222 variable stars and transients with Δ  ≥ 4 mag, most of which are new discoveries.Novae and nova candidates make up the largest proportion but YSOs are slightly more numerous in the disc fields at 295 • <  < 350 • (see Figure 1).Whilst detection of eruptive YSOs was the principal motivation for this work, this was the first panoramic infrared search for variable sources of all types to cover the most obscured regions of the Milky Way, so it is not entirely surprising that it resulted in detection of an unexpected new population: the aperiodic late-type dipping giant stars in the Nuclear Disc.This discovery comes soon after the recent detection of the population of post-Red Giant Branch stars (Kamath et al. 2016;Oudmaĳer et al. 2022), thought to be products of post-common envelope evolution.
Our conclusions regarding YSOs are as follows.
• Of the 40 sources classified as YSOs, 32 are eruptive variables presumed to be undergoing an episodic accretion event.We describe most of these as "classic" outbursts with a single long-lasting event, the majority of which are spectroscopically confirmed to have FUortype or EXor-type spectra in the companion paper by GLK.
• The rise time of these long-lasting accretion-driven outbursts as measured in   is most often between 600 and 1000 d, notably slower than the 6 to 12 month timescale reported in the optical for the small number of historical FUor events where the timescale is known.The typical outburst duration is at least 5 yr, rather longer than many lower amplitude events (see the companion paper by Contreras Peña et al. submitted for a sample with a broad range of amplitudes.)There is considerable diversity in the light curves of the classic outbursts, in terms of rise time, rate of decline postmaximum, intra-year variability and overall morphology.The spread of rise times might be explained by the onset of the accretion disc instability occurring at different radial locations, or different disc viscosity parameters (Liu et al. 2022).However, the overall diversity may indicate that more than one disc instability mechanism can trigger long-lasting outbursts.
• The form of the classic outburst light curves during the rise is usually described by a slow start that then accelerates to a rapid linear rise in brightness that stops rather quickly when the bright state is reached.A minority of classic outbursts have a more symmetric "S-shaped" rise, with a slow start and a slow finish.
• A small proportion of the eruptive YSOs have irregular light curves, showing strong variability on both intra-year and inter year timescales and typically at least two clear maxima within the 9.5 yr VIRAC2 time series.These somewhat resemble the Multiple Timescale Variables detected in Guo et al. (2020Guo et al. ( , 2021)).Most of these sources presently lack spectroscopic follow up.
• Periodic outbursts are not seen at these high amplitudes, despite making up a significant proportion of outbursts up to 4 mag in amplitude (Guo et al. 2022).However, one YSO, source 148, shows periodic 2 mag outbursts atop a 3 mag classic outburst light curve (discussed further in GLK).
• Eight YSOs with deep, long-duration extinction events (dippers) are seen wherein the source is dimmer and usually redder for a year or more, attributable to inner disc structure (e.g.Bouvier et al. 2013;Rice et al. 2015).
• Notable amongst the dipping YSOs is source 207, a candidate KH 15D-like periodic system in the Lagoon Nebula Cluster (NGC 6530) that showed two almost equal peaks in flux within each 59 d period, attributed to alternating extinction of the two components of a binary YSO by a circumbinary disc (Chiang & Murray-Clay 2004;Arulanantham et al. 2016).
The aperiodic dipping giants will be an interesting topic for future study.The main things to note concerning these are as follows.
• These sources are strongly clustered in the Nuclear Disc, evidenced by the projection of 18/21 dipping giant candidates in that region and their high infrared extinction even outside the dipping event, Of the 21 candidates, 7 are spectroscopically confirmed as late type giants in GLK.This informed the classification of the 21 candidates, indicating that the remaining 14 that lack spectra at present are unlikely to be YSOs.
• The light curves of the dips are diverse but they are generally asymmetric (hence not caused by a circumsecondary disc) and last for at least a year, more often a few years or longer.Their (−  ) colours approximately follow the interstellar reddening law, indicating that the variable extinction is due to small dust grains.Some variability is almost always seen outside the main dip.Only a minority of the dipping giants have red mid-infrared colours in the Spitzer/GLIMPSE-II survey, indicating that there is no pre-existing large mass of warm dust in most cases.This tends to suggest a scenario of dust puffs along the line of sight.
• In view of their location and the uniquely high metallicity (within the Milky Way) of many stars in the Nuclear Disc, we suggest that the aperiodic dips are caused in some way by super-solar metallicity, e.g. by an increase the production of highly refractory dust species.The initial spectroscopic analysis presented in GLK appears to support this, and the spectra, effective temperatures and luminosities appear to be consistent with an AGB population that is probably O-rich.The radiatively driven mass loss in AGB stars proceeds via radiative pressure on dust but is usually linked to pulsation (Höfner & Olofsson 2018).However, it depends on luminosity, chemical composition and the 3-D process of dust production (Dell'Agli et al. 2015), which is influenced by convection and shocks (Freytag & Höfner 2023).Further observational and theoretical work are clearly needed to understand this population, including a search for lower amplitude examples of the phenomenon.
Other discoveries to note in the sample are as follows.
• The set of 35 LPVs includes several dusty Mira variables with deep, long-lasting dips or fading trends.All three of the LPVs that are projected in the Nuclear Disc display this behaviour, two of these having the deepest dips in the set (3.5 to 4 mag in   ) which far exceed the amplitude of the pulsations in these stars.
• The set of 70 novae and nova candidates includes 24 new candidates.Some of the light curves display unusual features within the overall decline, such as a cusp or oscillations, to our knowledge previously reported only in the optical waveband (Strope et al. 2010).
• We identify 10 unusual variable sources, several of which are not yet understood and warrant further investigation.

Data Availability
The VVV and VVVx images are publicly available in the VISTA Science Archive (VSA) (vsa.roe.ac.uk) and the ESO archive http://archive.eso.org/cms.html.VSA facilitates generation of cut-out images around a specified coordinate so we have not included images in this work.The final version of the VIRAC2 database, representing a minor improvement on VIRAC2-, is currently being prepared for publication and upload to ESO archive, expected to occur within the next few months.UKIDSS data are available at the WFCAM Science Archive wsa.roe.ac.uk.The WISE, Spitzer and 2MASS datasets are publicly available, see the

APPENDIX A: SEARCHES A1 VIRAC2−𝛼 and VIRAC2-𝛽 differences
As noted in the main text VIRAC2− light curves have data from 2010-2018, whereas VIRAC2- light curves have data from 2010-2019.The absolute photometric calibration of VIRAC2- is anchored to the 2MASS Catalogue of Point Sources (PSC) (Skrutskie et al. 2006) in a region of the lower bulge having low extinction and negligible source confusion, thereby circumventing absolute calibration issues in crowded star fields of the inner bulge (Hajdu et al. 2020).This calibration was propagated across the VVV survey area via the field overlaps.VIRAC2- was calibrated using the VVV dophot-based catalogue of Alonso-García et al. ( 2018), which provides ZYJHKs photometry averaged over two epochs, ultimately calibrated to 2MASS PSC on a per-field basis.The astrometric calibration of both versions of VIRAC2 is anchored to the absolute reference frame of the Gaia second data release (DR2).
VIRAC2- benefits from a more sophisticated detection-matching scheme for each source than VIRAC2-, allowing for proper motion and parallactic motion (Smith et al., in prep).This led to a modest but useful increase in the number of matched detections for a substantial minority of sources, within the 2010-2018 time interval.However, VIRAC2- contains a small fraction of cases where two adjacent stars were incorrectly matched as a single detection with high proper motion and/or large parallax.These occasional mismatches required some additional effort to remove false positive detections whilst making the searches as complete as possible.

A2 VIRAC2 searches
The VIRAC2 databases contain a number of light curve parameters that can help to distinguish bona fide variable stars and transients from false positives.The most productive searches were based on a combination of the well-known Stetson  index (Welch & Stetson 1993) and the slightly less well-known von Neumann  parameter (von Neumann 1941).The Stetson  index, when applied to times series measurements in a single filter, has a large value if pairs of magnitude measurements taken very close together in time differ from the light curve mean in the same direction, by amounts that are large in comparison to the measurement error.Defining  , ,  , and  , ,  , as the magnitude values and uncertainties in the th pair of measurements in the light curve, and  as the number of pairs, the Stetson  index is computed as: where , and m is the mean magnitude for the source.In VIRAC2, pairs of measurements are required to be separated by < 1 hour in order to contribute to the sum, though the interval is more typically a few minutes for consecutive pawprint stack images within a six point tile.If there were more than two measurements taken within an hour, pairs of measurements were formed by first selecting the two that were closest together in time, then forming any additional pairs by successively selecting the two remaining measurements that are closest in time.Most locations within the VVV survey area have two or more pawprint stack observations taken within an hour during the course of a six point tile observation (or during observation of an adjacent, slightly overlapping tile).However, a few highly variable sources found in our searches have relatively low Stetson  values (see Figure A1) because the number of suitable pairs of observations in the light curves is small, e.g.due to location at the edges of the group of four tiles that made up the observing block.
The von Neumann  parameter is conceptually somewhat similar to the Stetson index in that it has an outlying (small) value if the consecutive values in a light curve (whether close together in time or not) typically differ by less than the variance of the whole light curve.The  parameter is defined as: where   are the individual magnitudes,  is the number of measurements in the light curve and   is the standard deviation of all the magnitudes in the light curve.Unlike the Stetson index, the  parameter is robust against the problem of sources that have only a single pawprint stack image at each observing epoch so it helps to pick out bona fide variable stars in the lower left quadrant of Figure A1.However, the Stetson index performs better than the  parameter for detection of poorly sampled transient events, e.g.some microlensing events, which tend to lie toward the upper right of Figure A1 (right panel).7A simplistic search of the VIRAC2 databases for highly variable sources using only the value of Δ  inevitably finds numerous falsepositives where Δ  is inflated by bad data at one or more epochs.The databases also include numerous spurious sources located in the wings of the image profiles of very bright stars, though VIRAC2- mitigates this somewhat by including only sources having detections in at least 10 pawprint stacks.(This problem is partly due to the choice of dophot parameters in VIRAC2 because it was decided to opt for completeness in the database rather than exclude sources within the seeing halo of saturated stars).These spurious sources (including those relating to diffraction spikes and ghost images) are usually very obvious on inspection of the images.For example, ghost images of bright stars where the image falls just off the edge of an array typically appear as a streak or "jet" (see section 13.4 of Sutherland et al. 2015), though occasionally they resemble a point source.We also noticed that bright Long Period Variables (LPVs, which are often saturated in VVV), were a common source of false positives because dophot detected time-variable spurious sources in the image halo, the diffraction spikes or ghost images.The light curves of these false positives have the appearance of a very noisy sinusoid, the detections associated with each datum often having large values of the dophot ℎ parameter (indicating that the image profile was not well fit by a point source) and large values of the VIRAC2 __ℎ parameter, which quantifies the significance of the residual to the VIRAC2 five parameter astrometric fit as an astrometric  2 .
A selection of ∼122 000 VIRAC2- false positive candidates having Stetson >15, Δ  > 4 mag and detections spanning at least a three year time baseline is plotted in the left panel of Figure A1.The colour coding illustrates the fact that most false positive sources are detected in only a small fraction of the   pawprint stack images covering their location.It is also obvious that bona fide variable sources are concentrated at the upper left of the panels in Figure A1, typically with Stetson >500 and <0.5.Since it is impractical to visually inspect time series images of very large numbers of sources, the Stetson ,  and detection fraction parameters, supplemented by a few others, were combined in various ways to slice the large pool of candidates in all promising regions of the parameter space and yield an essentially complete list of all bona fide variable stars and transient sources in VIRAC2 with Δ  > 4 mag.
The precise details of our numerous searches of VIRAC2- and then VIRAC2- are given in the following sub-sections.Any VVV sources that were missed are either very poorly sampled or have variability that fell outside the dynamic range of the dophot-based photometry in the database.A large proportion of VVV transients are saturated during part of their light curves, requiring additional processing for this work, so a search using saturation-corrected photometry would probably yield some additional discoveries.E.g., the bright transient VVV-WIT-05, discovered by Saito et al. (2016b) at   ∼ 4, is not included in VIRAC2-.
Causes of false positives that were encountered during our searches were, in decreasing order of frequency: (i) bright stars, including bright LPVs (see above) (ii) stellar blends We attribute the last category to cosmic rays or bad pixels.Regarding category (iv) we retained sources in Table 1a if their amplitude in the one-day binned light curves remained close to or larger than 4 mag, computed after rejection of bad data (see Appendix B and 3 rejection of outliers during the binning.The logic here is that VVV is sparsely sampled and will therefore typically not measure the full amplitude of a variable source, so retaining borderline cases helps to give a more complete view.Here "close to" means Δ  > 3.6 mag in the case of transients and stars that dropped below the detection threshold when faint whilst in other cases we required that the amplitude can exceed 4 mag within the 1 errors of the photometry. The asteroid+star blends (category iii) occur in pawprint stack images taken at a single epoch in a light curve, corresponding to the spatial coincidence of a fairly bright main belt asteroid and a fainter star in the Galactic bulge, typically within a few degrees of the intersection of the ecliptic and the Galactic equator.These blends were often not detectable via inspection of the images because the spatial match was within 0.1 ′′ .They were instead identified using the Minor Planet Checker software available at the Minor Planet Centre of the International Astronomical Union.
Fortunately, the precision astrometry of VIRAC2 greatly facilitates identification of false positive candidates in all the above categories because any type of image defect, blend or spurious detection due to a bright star is often associated with large values of the __ℎ parameter, corresponding to astrometric shifts large in comparison to the 5-30 milli-arcsecond precision of the individual detections.The VIRAC2 data release paper (Smith et al., in prep) gives examples of astrometric elimination of false positives by colour coding each datum in VIRAC2 light curves according to the __ℎ value.

A3 VIRAC2-𝛼 searches
As described in Appendix A2, the most important search parameters were the Stetson  and von Neumann  parameters, in combination with a parameter approximating to the fraction of   images of a given sky coordinate in which a source is detected.In VIRAC2- the latter parameter is called pp2frac.To explain this term, we first need to define a "pawprint set" as a set of VVV+VVVX pawprint stack images that have the same telescope pointing coordinates (within 30 ′′ ).Then we recall that each location in a VVV tile is normally covered by at least two pawprint stack images so the great majority of sources have photometry in two or more pawprint sets.For each VIRAC2- source, the pawprint sets having detections were ranked according to the number of detections.(Commonly there are detections in only two pawprint sets, with a similar count in each).pp2frac is defined as the fraction of images in the second-ranked pawprint set for each source where there is a detection.For most sources this rather convoluted parameter is similar to the detection fraction in the full sequence of pawprint stack images covering the location on sky.The reader may be relieved to know that VIRAC2- enabled us to use just such a simple detection fraction in later searches.
Another VIRAC2- parameter employed in our searches is epoch_baseline, the time in years between the first detection and the last detection of each source.This parameter was helpful in searches for transients having low values of pp2frac.
Below we list the parameters employed in of the VIRAC2- searches in turn.All searches required Δ  > 4 mag in addition to the other cuts listed.Except where otherwise stated, the initial choice of parameter values was based on inspection of the parameter distributions with the aim of removing the bulk of false-positive candidates, some of which are illustrated in Figure 1 (left panel).
(i) Search 1a. 2   > 0.2 AND  > 15 AND  > 2 AND median   > 11.25 AND ( > 1000 OR  < 0.5).To remind,  is the number of pairs of contemporaneous pawprint stack images that have detections and therefore contribute to the calculation of Stetson .The cut on median   magnitude was designed to limit the number of false positive candidates caused by saturation.This was the most successful search, yielding 248 candidates of which 175 were confirmed to be real by visual inspection of the images.Sources found in this search tend to have fairly well-sampled light curves.They lie towards the upper half or the left side of Figure 1 (right panel) and are colour-coded in any hue except red.Eight of the false positive candidates were real variable sources found to have amplitudes Δ  < 4 mag after removal of a few outlying data points or, in one unusual case, the use of a separate photometric procedure (see Appendix B).
Inspection of the parameter distribution of bona fide variable sources found by this search found that the subset with  < 0.5 AND ( > 1000 OR  2   > 0.35) included 167/175 bona fide variable sources and only eight false positives.This refined set of parameters that describe the richest region of parameter space helped to define Searches 4a and 5a, see below, where the cuts on  2   > 0.2 and median   > 11.25 were removed.
The threshold median   value of 11.25 was estimated from results of the prior VVV DR4 search (see Appendix B) using aperture photometry data.Later investigation of the parameter distribution in VIRAC2- showed that a stricter threshold at median   = 12.0 would have more reliably removed false positive candidates caused by saturation.
(ii) Search 2a. 2   < 0.2 AND  ℎ_ < 1 yr AND  > 15 AND  > 2 AND median   > 11.25 AND ( > 1000 OR  < 0.5).This search was one of a few that attempted to detect transient or poorly sampled sources, which are often more difficult to pick out from the larger number of false positive candidates.Here the  2   parameter is relaxed, potentially adding a very large number of false positives detected at only a few epochs, but the  ℎ_ < 1 yr cut served to limit the number to a manageable level.This search for relatively brief transients yielded 44 candidates, seven of which were found to be real.
(iii) Search 3a. 2   < 0.2 AND  ℎ_ = 1 to 3 yr AND  > 15 AND  > 2 AND median   > 11.25 AND ( > 500 OR  < 0.5).This search for transients detected over a longer time baseline yielded 235 candidates, 11 of which were found to be real variable sources.Of these 11 sources, 10 have  < 0.5.
(iv) Search 4a. 2   < 0.2 AND  ℎ_ > 3 yr AND  > 1000 AND  < 0.5 AND median   > 11.25 AND  > 2. Most real variables detected over a time baseline longer than 3 yr have  2   > 0.2 and would be found in Search 1 so this search aimed to catch any unusual cases, aided by the requirement for both a high  value and a low  value that was informed by the rich region of parameter space identified in Search 1.This search yielded 41 candidates (after removing 2 duplicates in the database with almost identical coordinates) and seven of them were found to be real.Some of these are transients or eruptive YSOs that erupted during the later stages of the VVV/VVVX surveys, when the number of observations each year was smaller than in earlier years.
(v) Search 5a. > 15 AND  > 2 AND  < 0.5 AND ( > 1000 OR  2   > 0.35) AND median   < 11.25.This search reversed the mean magnitude cut used in Searches 1 to 4 with the aim of finding any real variable sources that are likely to be saturated in part of their light curves by exploring only the rich region of (, ,  2  ) parameter space identified in Search 1a.This effort was motivated by the discovery of the extremely high amplitude eruptive variable YSO WISEA J142238.82-611553.7 (Lucas et al. 2020) using time series data from the Wide-field Infrared Survey Explorer ( ) and Spitzer satellites.We noted that the source had been missed by Search 1a because it had median   = 11.24.Search 5 yielded 22 candidates and three of them were found to be real, one of which is WISEA J142238.82-611553.7.
(vi) Search 6a. ≤ 2 AND  ℎ > 4 AND  < 0.2 OR ( < 0.3 AND  ℎ > 10).Here  ℎ is the number of pawprint stack images having a detection in   .This search was designed to recover any well sampled variable sources that lack a useful Stetson  index value due to a location at the outer boundary of the VVV survey area or on the edge of one of the groups of four tiles that make up a typical VVV/VVVX observing block.(Some variable sources at the edge of a group of tiles were recovered in earlier searches employing the Stetson index because adjacent tile groups were sometimes observed consecutively, depending on the vagaries of the VISTA observing queue.)The above selection yielded 73 candidates, of which three were confirmed as real.

A4 VIRAC2-𝛽 searches
Our searches of VIRAC2- were mostly a repeat of the VIRAC2- searches for sources with Δ  > 4 mag, with some necessary adjustments.These searches aimed to benefit from the additional year of data and the more sophisticated source matching of the newer database.A feature of VIRAC2- is that some light curves contain values with very bright or very faint magnitudes, well outside the range of meaningful VVV dophot measurements, even for saturated stars.These sources typically have  > 1.99.All the searches listed below therefore included cuts against sources where the brightest magnitude was   < 9 or the faintest magnitude was   > 85 or  > 1.99, in order to reduce the number of candidates requiring visual inspection.Cuts against sources with a very large absolute value of parallax (|| > 500 mas) or proper motion ( > 500 mas yr −1 ) were implemented in all searches; a check for bona fide very high proper motion stars not already detected in VIRAC1 had already shown that there are none in the database.An additional limitation of VIRAC2- is that it contains a large number of false positive candidates, caused by mis-matches, in a region of 8 contiguous VVV bulge tiles (∼12 deg 2 ) at 1.6 • <  < 7.5 • , −3.6 • <  < −1.5 • where the survey has a much higher sampling cadence than elsewhere.The combination of high cadence and high source density in this region is presumed to be responsible.It was necessary to exclude this area from all the VIRAC2- searches and rely on VIRAC2-.VIRAC2- enables calculation of the fraction of   images of a given sky coordinate in which a source is detected via the ratio of the ks_n_detections and ks_n_observations database parameters.We denote this ratio as   and it is used in VIRAC2- searches in place of the  2   parameter.
(i) Search 1b.  > 0.2 AND  > 7.5 AND  > 2 AND median   > 11.25 AND ( > 100 OR  < 0.5).Here an important change from Search 1a was the reduction of the thresholds in Stetson  from 15 and 1000 to 7.5 and 100, respectively.This was required due to the typically lower values of this index in VIRAC2- than VIRAC2-, owing to the larger photometric errors in the newer database that include the calibration uncertainty.Stetson  is usually a factor of ∼2 lower in VIRAC2- than VIRAC2- but in some cases the difference is larger.This search yielded 613 candidates not selected in earlier searches, of which 17 were found to be real.
(ii) Search 2b.  < 0.2 AND  ℎ_ < 1 yr AND  > 15 AND  > 2 AND median   > 11.25 AND ( > 100 OR  < 0.5).This search for transient sources yielded only one new candidate, which did not pass visual inspection.The small number of candidates is likely to be due to the requirement of a minimum of 10 detections in papwrint stack images for a source to be included in VIRAC2-.
(iv) Search 4b.  < 0.2 AND  ℎ_ > 3 yr AND  > 100 AND  < 0.5 AND median   > 11.25 AND  > 2. This search yielded 78 new candidates, none of which were found to be real variable sources.
(v) Search 5b. > 7.5 AND  > 2 AND  < 0.5 AND ( > 100 OR   > 0.35) AND median   < 11.25.This search yielded two new candidates, neither of which were found to be real variable sources.
Whilst most of the above VIRAC2- searches were fruitless, the 17 sources found in search 1b included a few sources that are of very rare or unique nature within the full VIRAC2 sample of 222 sources presented in this work so the effort can be deemed to be justified.These are discussed in Appendix I. Examples include a born-again giant undergoing a late thermal pulse (source 175 = V4334 Sgr, aka Sakurai's star) and a black hole X-ray binary candidate (source 16 = MAXI J1348-630).In addition, one of the "false positive" candidates from search 1b was also of astrophysical interest and helped to motivate further work, though its 3.5 mag amplitude fell below the Δ  = 4 mag threshold for inclusion in the list published here, after removal of a bad datum.This source, recovered separately by Guo et al. (2022) as VVV_PB_41, was an early example of a periodically outbursting YSO candidate and the amplitude is among the highest of many such objects presented in that work.

A5 VVV DR4 searches
Candidate variable stars in VVV DR4 (comprising 2010 to 2013 VVV data) were selected in two stages.First, a search of public SQL database at the VSA at Edinburgh was done for all sources with an amplitude exceeding a 4 mag threshold (later reduced to 3 mag) with the additional criteria that in each light curve all points lie in the range 0 <   < 18.These cuts removed a large number of spurious faint detections at the noise level of the survey, whilst also removing non-detections, reported as negative magnitudes in the VSA database.Non-detections appear to be numerous in the database so it was necessary to exclude them.A final criterion was that there were at least 5 good quality   observations of the field, i.e. the parameter ksnGoodObs ≥ 5.
The first stage provided 460 candidates with Δ  > 4, and later 3343 candidates with Δ  > 3. Since VVV DR4 was based on the pipeline aperture photometry from the VISTA tiles, i.e. images created by coadding the six pawprint stack images taken for each field at each epoch, it was not possible to compute a Stetson index with the DR4 data.Therefore a local database was created at Hertfordshire by co-author Smith, using only the pipeline aperture photometry performed on the pawprint stack images by the CASU team at Cambridge.The database contained matched photometry for every source over the 2010 to 2015 time period.With this database prepared, the second stage of the search was simply to compute the Stetson  index for each candidate, discard those with  < 15 and then visually inspect the much more manageable number of images remaining.The  = 15 threshold was based on inspection of the data distribution and the images for a subset of candidates.Stetson  performed similarly well, but Stetson  produced a marginally cleaner selection for this data set.For candidates observed to be bright on a single date and otherwise faint and non-variable, a final step was to check for asteroids in the same way as for VIRAC2 candidates.This step removed several false positives of that type.

APPENDIX B: ADDITIONAL PROCESSING AND DIGITAL FORMAT LIGHT CURVES
Light curves are provided for all sources in two tar sets.As noted in the main text, these are usually based on VIRAC2- data, with further processing, except in three cases, sources 60, 125 and 185, where VIRAC2- was used because a source was missing from VIRAC2-.In addition, the final publication version of VIRAC2 became available in the final stages of this project so it was used for sources 103 and 156, which are blended sources for which VIRAC2- and VIRAC2- light curves were inadequate.In each case it was necessary to amalgamate measurements from light curves for two separate sources, aided by the per-epoch astrometry supplied in the final version of VIRAC2.Due to the substantial work it would entail for very little benefit, the publication version was not used for other sources.
The first tar set contains the VIRCAM pawprint-based light curves in the  , , ,  and  filters.The second tar set contains the binned version of the light curves that are used in the figures in this work, using 1 day bins and giving inverse variance weighted mean magnitudes.The binned version includes the most sensitive upper limit (see below), for sources that dropped out in one or more calendar years.A PDF file of all the light curves is also included in the Supplementary Information, using the binned light curves except for source 185, which has strong intra-night variability.The plots in the PDF do not include the upper limits since these tend to compress the time axis to an undesirable extent.
Both tar sets benefit from additional processing as follows.All VIRAC2 light curves were inspected for signs of saturation.Bright detections having values of the dophot stellar profile parameter ℎ > 5 or the VIRAC2 astrometric parameter __ℎ > 30 were automatically saturation corrected.An additional check was performed using the ratio of fluxes in concentric apertures of radius 1.0 ′′ and 1.414 ′′ and detection with unusually low values of this ratio (compared to other detections in the same image) were also assumed to be saturated.Where saturation occurred, we replaced the datum with a saturation-corrected value computed with the fit-sio_cat_list FORTRAN script available at casu.ast.cam.ac.uk, which estimates source magnitudes using the flux in ring-shaped apertures outside the saturated core of the image profile.Images were visually inspected to select suitable apertures for this and ensure that the fluxes were not contaminated by adjacent stars.All flux measurements having ℎ > 5 or __ℎ > 30 were set to zero in the light curves and not used further.
A few sources are spatially resolved, e.g. the eruptive YSO VVV 1640-4846 displays compact nebulosity whilst in outburst.Aperture photometry in a 1.414 ′′ diameter aperture was used in such cases since the VIRAC2 dophot-based light curves display large scatter.
Many sources such as transients, dipping giants and dipping YSOs faded below the single VIRCAM pawprint-based detection limit of VIRAC2 in some calendar years.In an attempt to increase the dynamic range of this work, we stacked cut-out images for such sources in each calendar year where there was no detection or only a single pawprint detection (in the latter case the aim is to provide confidence in the photometry, given that there are usually two or more pawprint images at each epoch.)Aperture photometry was then performed on the image and any detections were added to the binned version of the light curve.In practice, this process led to only a small number of detections since (i) the annual stacks are not significantly deeper than single pawprint images in the more crowded VVV star fields (such as bulge fields that contain most of the novae and dipping giants) and (ii) the annual stacks are often only 1 or 2 mag deeper than single pawprints even in uncrowded fields.For sources with non-detections in a calendar year the variability amplitudes given in the 8th column of Table 1a, Δ  , is given as a lower limit, e.g.Δ  > 4.16 for source 7.This is based on either (i) the faintest detection in the binned light curve, or, (ii) approximate upper limits computed using the annual stacks by generating artificial stars with a Gaussian image profile at the source's location and performing forced aperture photometry there.Artificial sources were counted as detected if the measured source magnitude agreed with the injected magnitude to within 0.5 mag.

APPENDIX E: VVV AND WISE VARIABILITY OF DIPPING YSOS
In eight YSO candidates it is likely that the variability is due to changing extinction.With one exception, the light curves plotted in Figure E1 show at least one long duration extinction event lasting a year or more.In one case the source remains in a faint state at the end of the VVV time series.In terms of their high infrared extinction, most of these eight systems can be classified as extreme examples of the aperiodic dipper phenomenon, a term used by Morales-Calderón et al. (2011) and Rice et al. (2015) to distinguish aperiodic dippers from AA Tau-like dippers that show display periodic dips (typically  < 15 d) attributed to a magnetically induced warp in the accretion disc (Bouvier et al. 2007) and sometimes long duration aperiodic dips also.Source 207 is different, being periodic.It is discussed in more detail below.
In addition to the dips in the light curves, in most cases there is further evidence that the variability is due to extinction, based on inspection of data from 2MASS, the Deep Near Infrared Sky Survey (DENIS, Epchtein et al. 1999), VVV, and to a lesser extent the WISE survey.The evidence is summarised in Table E1, see also Figure E2.Seven of the nine dippers were detected in   by 2MASS or DENIS between 1996 and 2000, all at a flux level similar to or slightly brighter than the brightest measurements in the VVV time series.8 .This supports the view that VVV observed dips below the typical brightness level rather than a gap between accretion-driven eruptions of the type seen in V346 Nor (Kóspál et al. 2020).
Most of the dipping YSOs also become redder when fainter in the   vs ( −   ) colour magnitude diagram in a manner that approximately follows the Wang & Chen (2019) reddening law for interstellar extinction, see Figure E2.(Some allowance should be made for departures from an idealised reddening behaviour due to changes in the accretion disc over timescales of years, see e.g.Covey et al. (2021).) Figure E2 illustrates the colour changes for sources 28, 55, 136, 207, 208 and 212, the six sources having colour measurements on two different dates.They are plotted with blue points in the brighter state and red points in the fainter state, each pair being joined by a dashed line.The exception is source 207, the KH15Dlike candidate, which shows no significant colour change despite a change in brightness of almost 4 magnitudes.Only source 207 and source 55 have detections in ,  and   that are contemporaneous (measured on the same night).Source 207 actually becomes bluer in the faint state ( −  = 0.71, compared to  −  = 0.94 in the bright state), which can perhaps be attributed to scattered light.
Finally, for one of the eight dippers, source 28, the WISE/unTimely time series data show large changes in 2, 1 and VVV   magnitudes that are consistent with variable reddening based on the Wang & Chen (2019) extinction law, see Figure E3.Supporting evidence from the unTimely or NEOWISE data is found in a further two systems but these cases are less clear, either because the flux changes sampled by WISE are small, or because the WISE and VVV light curves are not contemporaneous and have subtly different trends over time (not shown).In the other five systems the sources are not genuinely detected by WISE due to blending, or at least not detected when they are faint.Careful inspection of the WISE, GLIMPSE and VVV images, and any astrometric shifts in WISE detections over time, was required to guard against being misled by blends in the WISE data.
Source 207.This is a periodically dipping system (period  = 59.35d) that appears to resemble KH 15D (=V582 Aur, Chiang & Murray-Clay 2004;Herbst et al. 2002), a system wherein the two bright components of a binary YSO are extinguished by a circumbinary ring for much of their orbits but each becomes visible in alternation at certain phases.A similar type of behaviour has been seen at lower amplitude in WL4 and YLW16A in the  Ophiuchi dark cloud (Plavchan et al. 2008(Plavchan et al. , 2013, though these two are triple systems) as well as CHS 7797, ONCvar 149 and ONCvar 1226 in the Orion Nebular Cluster (Rodríguez-Ledesma et al. 2012, 2013;Rice et al. 2015) and in the recent candidates Bernhard-1 and Bernhard-2 (Zhu et al. 2022).Source 207 is currently mis-classified in the SIMBAD database as a "Cataclysmic Binary" (designated OGLE BLG-DN-652) due to being listed as a photometric dwarf nova candidate, within a list of 1091 candidates (Mróz et al. 2015) drawn from the OGLE survey (most of which are indeed high quality dwarf nova candidates).In this case though, the light curve does not resemble that of a dwarf nova (see Figure E1, bottom left and bottom centre panels) and the system is projected in the Lagoon Nebula Cluster (NGC 6530), with cluster membership confirmed by Gaia parallax and proper motion (and VIRAC2 proper motion) consistent with that of the cluster, at distance  ≈ 1.3 kpc. 9  The two components of the binary have similar but not identical luminosities, see Figure E1, such that the "outburst frequency" given 9 The Gaia DR3 parallax and proper motion of source 207 are,  = 0.71 ± 0.13 mas,    = 1.95 ± 0.09 mas/yr,   = −2.05± 0.06 mas/yr.The Gaia DR3 cluster parameters, are  = 0.76 mas,    = 1.28 mas/yr,   = −2.06mas/yr, with standard deviations of 0.66 mas/yr and 0.70 mas/yr for    and   and a systematic uncertainty on  of order 0.02 mas.by Mróz et al. (2015) is double the true frequency of the periodic variability.The system is also listed with half the true period in separate VIRAC2--based searches for periodic systems by Molnar et al. (2022) and Guo et al. (2022).(The latter study noted it as a dipping YSO but analysis was deferred to this study because the system had already been classified as KH 15D-like in the initial stages of this work).The peak near phase 0.6 in the phase-folded plot is narrower as well as shorter, but the colour-coding by Modified Julian Date shows that it broadened over the course of the survey.The few red-coloured points, corresponding to the sparse sampling at the end of the time series, indicate that the phasing of the variability continued to change and the dips may have become shallower.In KH 15D changes of this sort are explained as precession of the warp in the circumbinary ring, see e.g.Arulanantham et al. (2016), so this is consistent with our interpretation of the system.In VVV, source 207 is not very red, with ( −  ≈ 0.8), ( −   ≈ 0.45) and its variability is almost wavelength independent in the five contemporaneous epochs of   photometry, suggesting either grey extinction by large dust grains or a contribution by scattered light in the faint state, as is seen in KH 15D.The system will be analysed in more detail in a future publication, with the aid of follow up spectroscopy.

APPENDIX F: LONG PERIOD VARIABLES
The 35 sources classified as LPVs are dusty Mira variables with unusually high amplitudes.These represent the upper end of the amplitude distribution in such sources and consequently given the period vs. amplitude relation (see e.g.Nikzat et al. 2022) they also tend to have unusually long periods.The periods are given in Table F1.Given the inner Galaxy locations (see Figure 1), we expect that most of the 35 sources are O-rich.Some (10/35) are known OH/IR stars and most (24/35) have previously been identified as LPVs, either in VVV-based variable star searches (e.g CP17a), or earlier work (see references in Table 1a).
A notable example, shown in Figure F1), is source 107 (=OH 355.815-0.226 = IRAS 17328-3234) which has a peak to trough amplitude of ∼5 mag in   (Δ  ≈ 5.6 mag in total, due to a slow fading trend) and a period  ≈ 2000 d.Source 2 (=OH 259.8-1.8 = IRAS 11438-6330) also has a ∼5 mag amplitude but a shorter period and no clear long term trend.These are the highest   amplitudes of pulsation that we are aware of in Mira variables.OH/IR stars often have large amplitude pulsations and they are associated with high mass loss rates (see e.g.Gaylard et al. 1989;Höfner & Olofsson 2018).Source 2 is listed as d001-79 in a previous VVV-based study  228 (=V5084 Sgr).In that work (with a four year time series in ), a perfectly regular period of 570 d was found and a 1.5 mag peak to trough amplitude was reported, although inspection of the published light curve shows a value closer to 1 mag, consistent with the pulsations seen in Figure F1).Without this prior work, the VVV light curve of source 147 might be interpreted as showing signs of a long secondary period with a 4 mag peak to trough amplitude and a possible period  ∼ 2100 d.However, this seems unlikely in view of the previous lack of any trend over fours years, coupled with the fact that the flux in the VVV light curve remains near maximum after recovering from the dip.A 2100 d period is also a bit longer than the 200 to 1600 d range usually seen in such systems (Soszyński et al. 2021).Glass et al. (2001) noted that many of the Miras are "quite irregular" but did not discuss aperiodic trends further.
Source 163 shows a long term fading trend after 2012, with sinusoidal oscillations superimposed upon it (period  ≈1 yr, peak to trough amplitude ∼1 mag), see  The periodic variations and the overall decline are fairly well sampled by NEOWISE in 1 and 2 (better than the unTimely data in this case), with similar behaviour to that seen in VVV and no significant blending issues.Within the individual cycles, the magnitude changes between MJD = 58000 and MJD = 59000 are very similar in   , 1 and 2.This is consistent with Mira-like pulsations: Iwanek et al. (2021) found that the amplitude of such pulsations declines with increasing wavelength from the optical to the midinfrared but there was no significant change between the   , 1 and 2 passbands.(There is no useful WISE light curve for source 147: the VVV   images show it is heavily blended in both WISE and Spitzer/GLIMPSE).
Considering the longer term fading, we can compare the flux levels in similar (faint) parts of the periodic cycle in 2014 (MJD ≈ 56700) and 2018 (MJD ≈ 58350).Source 163 faded by 2.8 mag, 2.1 mag and 1.3 mag in   , 1 and 2.The change in   is smaller than expected from the 1: 0.52: 0.33 ratio predicted by the Wang & Chen (2019) extinction law but the Δ1 to Δ2 ratio is as expected.This could be explained by a model where the mid-infrared fluxes decline as extinction rises but the   flux declines more slowly due to a contribution by scattered light.VVV provides only two Dusty Mira variables often show aperiodic variability on a range of time scales (Iwanek et al. 2021), slow trends typically being attributed to mass loss and variable extinction.Whilst sources 163 and 147 could simply be seen as extreme examples of the variable extinction seen in many dusty Miras in the Milky Way, their locations in the Nuclear Disc mean that they could plausibly have super-solar metallicity, as proposed for the 21 dipping giants discussed earlier, which may enhance dust production.Such aperiodic events have been observed, albeit with considerably lower amplitude in   than source 163, in the Mira variables R For (Whitelock et al. 1997), II Lup (Feast et al. 2003) and IRAS 17103-0559 (Jiménez-Esteban et al. 2006).We should consider whether a long secondary period can explain the event in source 163.Such periods are quite common in LPVs, though less common in the bulge (Soszyński et al. 2013) and they have been convincingly attributed to eclipses by either a circumsecondary disc around a very low mass companion or by matter spread around the orbit of such a companion (Soszyński et al. 2021), see also Wood et al. (1999).It is proposed by (Soszyński et al. 2021) that the companion is a former giant planet that has gained mass from the stellar wind of the pulsating Mira variable and become a brown dwarf.A secondary period is not yet established in source 163 and it would have to be rather longer than the 200 d to 1600 d periods usually observed (see Figure 2 of Soszyński et al. 2021).However, a few rare cases of very long secondary periods are known, such as the C-rich Mira R Lep (Whitelock et al. 1997) and V Hya (Knapp et al. 1999), the latter recently revealed to have an expanding spiral circumstellar disc (Sahai et al. 2022).
Source 163 is projected in a bright H ii region seen in the WISE mid-infrared images, apparently corresponding to the giant H ii region G000.5-00.7 identified by Kuchar & Clark (1997) and Conti & Crowther (2004).It also has a large infrared excess and the Spitzer colour [8] − [24] = 3.3 is more typical of YSOs than dusty AGB stars (Robitaille et al. 2008).For these reasons it was initially thought that source 163 might well be a dipping YSO of the AA Tau type.However, those systems usually have much shorter periods ( < 15d, Rice et al. 2015), with the longer-term dip in AA Tau and some other systems attributed to occultation by non-axisymmetric matter in a much wider orbit (Bouvier et al. 2013).The VIRAC2 proper motion helps to reject the possibility of a dipping YSO: the values are    = −5.04 ± 0.42 mas/yr and   = −7.41± 0.46 mas/yr, which corresponds to   = −8.03± 0.43 mas/yr and   = −3.97± 0.45 mas/yr in Galactic coordinates.For a source projected very close to the Galactic centre at heliocentric distance /kpc, these relatively large proper motions correspond to tangential velocity components,  = 38 km/s and  = −19 km/s, parallel and perpendicular to the Galactic disc, respectively.This is incompatible with a YSO in the Galactic disc, unless it were located within ∼1 kpc of the sun.That is highly unlikely, given that no such nearby star formation region is known on the line of sight to the Galactic centre.

APPENDIX G: CATACLYSMIC VARIABLES AND OTHER TRANSIENTS
The great majority of these 70 sources are expected to be classical novae.In fact, 43 of them have been previously classified as classical novae, recurrent novae (i.e.novae for which more than one explosion has been observed) or nova candidates.These 43 were discovered either in the VVV-based searches listed in § 3, by optical surveys, e.g.OGLE in 12 cases, the All-Sky Survey for Supernovae (ASAS, Shappee et al. 2014) in one case, the Gaia Photometric Science Alerts system (Hodgkin et al. 2021) in one case, by amateur astronomers, or, by the mid-IR WISE satellite in three cases (e.g.Zuckerman et al. 2023).(Individual references and identifications are given in Table 1a).Some of the remaining 27 nova candidates are poorly sampled so their classification is tentative.While most of these 27 light curves display a typical nova-like monotonic fade, there are a few cases where a smooth-topped secondary peak attributable to dust condensation (Gehrz 1988(Gehrz , 1999;;Hachisu & Kato 2006) is observed, with varying degrees of confidence.Source 84 shows such a peak and there are hints of one in sources 45 and 222, though source 45 has noisy data due to saturation and source 222 is poorly sampled.Sources 20 and 109 show a modest rise from the initial detection up to the maximum flux: we have chosen to classify these two as nova candidates where the initial outburst maximum was missed.Sources 160 and 183 show a similar behaviour but with noteworthy differences (see below).
Four of the 29 new nova candidates have rather unusual light curves, discussed briefly below and illustrated in Figure G1.
(i) Source 88 has a sharp secondary peak in its light curve 680 days after the first observation (see Figure G1), rising almost 2 mag above the overall declining trend.This feature, distinctly not a rounded secondary dust peak, lasted for at least 100 days and it resembles the "cusp" seen in the rare C-class of novae defined by Strope et al. (2010) in optical light curves.That work stated that only three examples of such novae were known at the time (∼1 per cent of novae) and the origin of the cusp features is not yet clearly determined, though two theories are noted in that work.A more recent example is found in Wee et al. (2020).
(ii) Source 97 was first detected on 2011 August 30.The flux declined rapidly but quasi-periodic oscillations within the overall decline were seen in the 2012 observing season (data at Modified Julian Dates, MJD = 56086 to 56198).The period was ∼45 days.In subsequent years the photometry began to suffer from added noise due to blending with adjacent stars so it is difficult to determine whether the oscillations continued.Oscillatory behaviour in optical light curves is described as the O-class in Strope et al. (2010) but only two examples were mentioned.We have not seen this behaviour in other infrared light curves of novae in the literature, such as the numerous examples available in the SMARTS database (Walter et al. 2012).
(iii) Source 160 was first detected as a bright star (  ≈ 11.4 on 2011 August 4) but it then brightened by at least 2 mag over the following 70 days until observations were interrupted by the end of the observing season.From February 2012 the light curve was initially flat at   ≈ 10.8 for about 4 months but it then began to decline steadily, fading from view after 2014.This source is projected close to the Galactic centre ( = 0.74,  = −0.36)and detections in the  filter on two dates in August 2011 allow us to calculate a colour  −   ≈ 7.2.If we employ the extinction law of Cardelli et al. (1989) and assume an unreddened colour  −   = 0 for a hot source, this implies   ≈ 27 mag, consistent with a location in the Nuclear Disc of the Milky Way and the CMZ.The initial brightening behaviour can probably be accommodated within the framework of a classical nova with a secondary dust peak, from which we infer that the unobserved primary peak occurred between the last pre-outburst observation on 2010 September 12 and the initial detection almost a year later.A concern is that the proposed secondary peak has a rather long duration for a nova, given the event took 400 days to fade below the level of the initial detection.Several examples of novae with a secondary peak are shown in the SMARTS database (e.g.N Mus 2018) but in every case the peaks have a shorter duration than is seen in this source, by factor of two or more.However, the timescales over which nova shells expand and fade have a very wide range, determined by the energetics of the explosion (see e.g.Strope et al. 2010), so our initial classification is that this was a slow nova event corresponding to a low energy outburst.
(iv) Source 183 rose in brightness by ∼2 mag over ∼ 100 days up to the time of maximum flux, suggesting that we are seeing a secondary dust peak.While this is not especially remarkable, there is a also a dip in the light curve just before maximum flux, having a duration of approximately two to four weeks, though measured on only two nights in the relevant interval.This might be caused by high optical depth in the ejecta.

APPENDIX H: MICROLENSING EVENTS
We classify 38 sources from Table 1a and a further 15 sources from Table C1 as microlensing events, with varying degrees of confidence owing to the sparse sampling of the VVV survey.The classifications are based on visual inspection of the light curves for the characteristic morphology and a timescale of a few weeks.Almost all events appear to be single star lenses, as opposed to binary lens events.Three of the events (sources 180, 189 and 214) were previously reported by the OGLE early warning system (Udalski 2003) and two additional events with unusually good sampling, sources 215 and DR4_v51 from Table C1 , were verified by co-author Navarro using a microlens event-fitting code.This gives us some confidence  This paper has been typeset from a T E X/L A T E X file prepared by the author.

Figure 1 .
Figure 1.Spatial distribution of variable sources in the Δ  > 4 sample.The vertical scale is stretched for clarity.The "dipping giants" show a strong clustering towards the Nuclear Disc of the Milky Way.VVV survey boundaries are indicated with black lines.

Figure 3 .
Figure 3. Rise time plotted against amplitude (left panel) and outburst duration (right panel) for 32 classic events drawn from this work and earlier VVV discoveries.Typical rise times are 600 to 1000 days, though there are faster rises and slower rises.Durations are usually lower limits since the events continued beyond 2019.The shaded region is empty because the duration includes the rise time.

Figure 4 .
Figure 4. Distribution of spectral indices for 16 classic eruptive variables.

Figure 5 .
Figure 5.The distribution of the fitted rise timescale for the classic outbursts.The quantity    = 4 encompasses most of the rise (see main text) so it gives a similar but slightly shorter timescale than the full rise timescales plotted in Figure3, which were estimated from visual inspection.

Figure 6 .
Figure 6.Light curves of the nine irregular outbursting YSOs (upper three rows) and three others (bottom row) that mainly show a fading trend in VVV.An inset is included for source 42, to better illustrate the rapid variation of the first outburst of that source.

Figure 7 .
Figure7illustrates the variability of the eruptive YSOs as a function of time baseline, comparing 21 classic outbursts (grey dashed lines) with the six irregular eruptive systems in the upper two panels of Figure6(solid blue lines) and sources 42 and 48 (solid red lines).These plots are constructed by computing the r.m.s.variability (upper panel) or maximum variability (lower panel) over all pairs of   points,   ,   , in each light curve and binning the results in 0.5 dex bins in the logarithm of time baseline.The r.m.s.variability was computed using the unbinned pawprint light curves and the uncertainty in each value of (  −   ) 2 was subtracted from the mean square variability in order to debias the result.The r.m.s.variability of source magnitudes is the square of the structure function as defined inLakeland & Naylor (2022), though some authors have preferred to use flux differences, see e.g. the YSO optical variability studies ofSergison et al. (2020),Venuti et al. (2021).An approach based on magnitudes, i.e. flux ratios, is more suitable than flux differences for large amplitude variability, as discussed inLakeland & Naylor (2022).The maximum variability vs. time baseline (also called the "accumulation function", seeGuo et al. 2020) is a useful alternative for the sparsely sampled VVV light curves.This was computed using the 1-day binned light curves plotted in this work, in order to reduce the effect of individual 2-3 outliers.We see that the six irregular systems and the two fast eruptive systems tend to have more variability than classic outbursts on timescales up to 100 d.The irregular outbursts generally show maximum brightness changes of 1 mag or more on a 100 d timescale (Figure7, lower panel).The variability of the classic outbursts rises sharply on baselines from 100-1000 d.Source 42 displays the largest r.m.s.variability on 2 d to 300 d time baselines, though this is partly due to the good sampling of the first outburst seen by VVV.Source 93 stands out slightly from the other classic outbursts over 30 d to 100 d time baselines due to the unusually fast rise time of the event (see Figure2).YSOs having irregular high amplitude accretion-driven behaviour clearly require further study and spectroscopic follow up, as do the MTVs reported inGuo et al. (2020) andGuo et al. (2021).The connection to classical EXors(Giannini et al. 2022) may then become

Figure 9 .
Figure 9.Light curves of the dipping giants

Figure 9 -
Figure 9 -continued.Light curves of the dipping giants.

Figure 10 .
Figure 10.  vs.  −   colour-magnitude diagram for 20 of the dipping giants.Blue points correspond to the brightest available epoch of multi-filter VVV data.Red points show the faintest available VVV epoch, plotted for 11 sources.Elements of a pair are connected by dashed lines.The changes are roughly parallel to the extinction vector.

Figure A1 .
Figure A1.Stetson  vs. von Neumann  plots for candidate variable sources (small red + symbols) and bona fide variables (filled circles).The points are colour coded by a parameter that represents (approximately) the fraction of images of the field in which the source was detected.The plotted data distribution drives most of the logic of our search methods.The left panel includes numerous false positive candidates, whereas the right panel shows more clearly the distribution of the bona fide variables in our sample VIRAC2 sample of 222 Δ > 4 mag sources.Data are taken from VIRAC2- except for 17 sources discovered using VIRAC2- which are highlighted with a ring around the filled circle.These 17 tend to have slightly lower Stetson  values, see Appendix A.
The cluster parameters are based on median values from the set of cluster members in Prisinzano et al. (2019) having 5 parallaxes and reduced unit weight error,  < 1.4.

Figure D1 .
Figure D1.VVV and WISE light curves of six irregular eruptive YSOs.

Figure E1 .
Figure E1.Light curves of the eight dipping YSOs.The variation in these systems is thought to be due to variable extinction by circumstellar matter.Source 207 appears to be a KH15D-like system with periodic dips due to a circumbinary disc.The phase-folded light curve (lower middle panel) is colour-coded by Modified Julian Date.

Figure E2 .
Figure E2.  vs.  −   colour-magnitude diagram for six of the dipping YSOs, labelled with their source numbers.Blue and red points correspond to the brightest and faintest available epochs of multi-filter data, respectively.Elements of a pair are connected by dashed lines.The changes are roughly parallel to the extinction vector in most cases.One source is clearly different: source 207 (the candidate KH15D-like object in the Lagoon Nebular Cluster) shows negligible colour change, see text.
Figure F2.2MASS measured   = 11.39 in 1998, similar to VVV measurements before the decline began.The VIRAC2- light curve is moderately well fit by the sum of sine wave with  = 338±2 d and a quadratic function (plotted as a dashed curve), until the fading trend appears to bottom out near MJD = 58500.The system is classified in the VIRAC2-based VI-VACE catalogue (Molnar et al. 2022) as an LPV with  = 692.8, approximately double the true period, after fitting and subtracting a quadratic long term trend (which was standard procedure in VI-

Figure F2 .
Figure F2.VVV and WISE light curves of source 163.These appears to be a combination of pulsation and a slower, extinction-driven dip.The dashed line shows a simple model fit, see main text.The background grid aids comparison of the amplitude of pulsation in the three filters.

Figure I1 .
Figure I1.VIRAC2 light curves of six of the unusual variable stars that are poorly understood at present, see main text.

Table 1a .
Table of 222 variable sources with Δ  ≥ 4 mag Eruptive YSO, 6 year slow rise from   ≈16 to   ≈12.Only the first 10 rows of the table are shown here, split into two sections.The full table is available in the online supplementary information.J2000 coordinates are calculated at epoch 2015.0 for all sources in VIRAC2− but at epoch 2014.0 for sources 60, 103, 125 and 185, which were taken from VIRAC2−  due to not being properly recovered in VIRAC2−.
Collaboration et al. 2023urvey Explorer, WISE,Wright et al. 2010).We then considered any past studies (assisted by SIMBAD,Wenger et al. 2000, Vizier and the  International Variable Star Index, VSX), follow up spectroscopy (for many of the YSO candidates, see GLK), source location and the parallax and proper motion from Gaia Data Release 3 (DR3, GaiaCollaboration et al. 2023) and VIRAC2.Source location can help to identify YSOs, if the star is projected in a star formation region, but this method is less useful on Galactic disc sight lines within 10 • of the Galactic centre, where indicators of star formation are ubiquitous in SIMBAD and in the WISE mid-infrared images.To test for location in a star forming region, we employed the method as described in CP17a andLucas et al. (

Table 2 .
Rise times and quiescent-state spectral indices of 32 classic outbursts   amplitudes in column 2 exclude short timescale variability so they are a little less than the Δ  values in Table1a. al.
Lawrence et al. 2007;2;Fritz et al. 2021) Nuclear Disc.A 19th case is plotted at  = 3.1 • and two others lie a little further from the Galactic centre, at (,)=(350.7 • ,-1.2 • ) and (8.2 • ,-0.2 • ), respectively.As discussed in §3, we classify 21 sources as dipping giants.Of these, 18/21 are projected in the Nuclear Disc of the Milky Way, see Figure8.This structure is roughly coincident with the Central Molecular Zone, with radius 220 pc and scale height ∼45 pc, corresponding to Galactic longitudes 358.45 • <  < 1.55 • , and scale height ∼0.3 •(Launhardt et al. 2002;Fritz et al. 2021).The   light curves of the 21 sources are shown in Figure9.They are very diverse but all can be described as a single, apparently aperiodic flux dip with a duration of at least a year.In fact most of the dips either started before or ended after the 2010 to 2019 VIRAC2 time baseline.Only 6/21 dips can be said to be fully contained within the time series.A cross match to 2MASS PSC and the United Kingdom Infrared Deep Sky Survey Galactic Plane Survey (UKIDSSLawrence et al. 2007; Herbst et al. 2002;Chiang & Murray-Clay 2004) Appendix E. We highlight one case of special interest: source 207 appears to resemble KH 15D (=V582 Aur,Herbst et al. 2002;Chiang & Murray-Clay 2004), a system wherein the two bright components of a binary YSO are extinguished by a circumbinary ring for much of their orbits but each becomes visible in alternation at certain phases.5 DIPPING GIANTS IN THE NUCLEAR DISC & Chen (2019) extinction

Table E1 .
Supporting evidence for an extinction-driven dip Multi-filter light curves are shown in Fig.D1, for six irregular eruptive YSOs having useful data in the WISE/unTimely database.Source 42 shows the ∼1 yr duration of the second outburst near MJD=58000 more clearly in the WISE 2 data than in VVV.Source 48 shows four brief, single epoch outbursts of 1 to 2 mag amplitude in the WISE data, suggesting a duration of order ∼6 months in each case.

Table F1 .
Periods of the LPVs