Determining the recurrence timescale of long-lasting YSO outbursts

TL;DR

This study estimates the recurrence timescale of large accretion outbursts in young stellar objects, revealing they occur roughly every 112,000 years, with fewer events in the class II stage compared to class I.

Contribution

First robust measurement of outburst rate in class II YSOs, showing large accretion events are less frequent than in earlier stages.

Findings

Outburst recurrence timescale is approximately 112 kyr.

Large accretion events are about 10 times less frequent in class II than in class I YSOs.

Identified 3 previously unknown long-lasting accretion events.

Abstract

We have determined the rate of large accretion events in class I and II young stellar objects (YSOs) by comparing the all-sky digitised photographic plate surveys provided by SuperCOSMOS with the latest data release from Gaia (DR2). The long mean baseline of 55 years along with a large sample of class II YSOs ($\simeq$15,000) allows us to study approximately 1 million YSO-years. We find 139 objects with $\Delta R\geq1$~mag, most of which are found at amplitudes between 1 and 3 mag. The majority of YSOs in this group show irregular variability or long-lasting fading events, which is best explained as hot spots due to accretion or by variable extinction. There is a tail of YSOs at $\Delta R\geq3$~mag and they seem to represent a different population. Surprisingly many objects in this group show high-amplitude irregular variability over timescales shorter than 10 years, in contrast with the view that high-amplitude objects always have long outbursts. However, we find 6 objects that are consistent with undergoing large, long lasting accretion events, 3 of them previously unknown. This yields an outburst recurrence timescale of 112 kyr, with a 68\% confidence interval [74 to 180] kyr. This represents the first robust determination of the outburst rate in class II YSOs and shows that YSOs in their planet-forming stage do in fact undergo large accretion events, and with timescales of $\simeq$100,000 years. In addition, we find that outbursts in the class II stage are $\simeq$10 times less frequent than during the class I stage.

Tables

Table 1. Table 1: Example of the photometry obtained for each high-amplitude variable. In the tables we present the modified Julian date of the observations (estimated from the Julian date as MJD=JD − - 2400000.5), the R 𝑅 R or r 𝑟 r -band magnitudes and the survey from where they are obtained. In the last column we present a note depending on whether the magnitudes are an upper limit or if these result from transforming between different filters (marked as approximate in the table). The full table is available online.

ID	MJD	Mag	Filter	Survey	note
V4	34270.404	21.0	R	SSS(R1)	upper limit
V4	47880.208	21.0	R	SSS(R2)	upper limit
V4	49253.445	18.5	R	SSS(B)	approximate
V4	49654.347	18.5	R	SSS(I)	approximate
V4	52953.981	17.46	r	IPHAS	–
V4	56227.446	16.79	r	Pan-STARRS	–
V4	56537.365	16.49	R	PTF	–
V4	57204.000	16.5	R	Gaia	approximate

Table 2. Table 2: The 139 high-amplitude variable stars detected in our analysis. The full table is available online.

ID	$\alpha$	$\delta$	SIMBAD ID	SIMBAD class	$\Delta R$	LC class	Final class	Reference	Known?	SSS?	YSO class
V1	20:58:53.72	$+$44:15:28.3	EM* LkHA 190	FUOr	5.2	Outburst	Long-lasting (class I)	This work (Section 5)	Y	Y	Marton
V2	20:58:17.03	$+$43:53:43.3	V* V2493 Cyg	FUOr	5.0	Outburst	Long-lasting (class II)	This work (Section 5)	Y	Y	Marton
V3	20:54:07.39	$+$41:34:58.1	V* V1219 Cyg	Orion V*	4.7	Repetitive fadings	Deep fades	Stecklum & Linz (2013)	Y	Y	Marton
V4	02:33:53.40	$+$61:56:50.1	2MASS J02335340$+$6156501	YSO	4.5	Outburst	Long-lasting (class II)	This work (Section 5)	N	N	2M_AAA
V5†	05:28:54.06	$-$06:06:06.3	V* RX Ori	Orion V*	4.5	Periodic	Non-YSO	Lee & Chen (2007)	Y	Y	WISE_slt04
V6	21:43:00.01	$+$66:11:27.9	V* V350 Cep	Orion V*	4.5	Outburst	Long-lasting (class I)	This work (Section 5)	Y	Y	Marton
V7	16:11:46.00	$-$25:32:00.8	V* V931 Sco	Orion V*	4.4	Irregular	Periodic (P=215 d)	Pojmanski (2002)	Y	Y	Marton
V8	05:10:11.00	$-$03:28:26.2	2MASS J05101100$-$0328262	YSO	4.0	Irregular	Short-term outburst	Sicilia-Aguilar et al. (2017)	Y	Y	WISE_slt04
V9	08:41:06.76	$-$40:52:17.4	2MASS J08410676$-$4052174	Em*	4.0	Outburst	Long-lasting (class II)	This work (Section 5)	Y	N	Marton
V10	22:53:33.26	$+$62:32:23.6	V* V733 Cep	FUOr	4.0	Outburst	Long-lasting (class II)	This work (Section 5)	Y	Y	Marton
V11	05:37:00.10	$-$06:33:27.3	V* BE Ori	Em*	4.0	Repetitive outbursts?	–	–	Y	Y	Megeath_alpha
V12†	16:00:07.42	$-$41:49:48.4	2MASS J16000742$-$4149484	Candidate YSO	3.8	Irregular	Non-YSO	Frasca et al. (2017)	N	N	C2D
$†$ Objects not included in any subsequent analysis.

Table 3. Table 3: Final sample of YSO outbursts.

ID	Name†	$\alpha$	Envelope?
V1	V1057 Cyg	$-0.1$	Yes (Fehér et al., 2017)
V2	HBC 722	$-0.44$	No (Fehér et al., 2017)
V4	2M J0233$+$6156	$-0.4$	?
V6	V350 Cep	$-0.17$	Yes (Muzerolle et al., 2004)
V9	2M J0841$-$4052	$-0.23$	?
V10	V733 Cep	$-0.78$	No (Fehér et al., 2017)
V31	V960 Mon	$-0.19$	Yes (Kóspál et al., 2015)
V36	V582 Aur	$-0.33$	No (Ábrahám et al., 2018)
V51	WRAY 15-488	$-0.27$	?
$†$ We do not present the SIMBAD ID of the variables stars from
Table 2, but instead the column shows the most common name
found in the literature.

Table 4. Table 4: List of surveys that are part of the analysis of high-amplitude variables found in our study.

Vizier ID	Title	Reference
V/73A/catalog	Third Catalog of Emission-Line Stars of the Orion Population	Herbig & Bell (1988)
III/177/catalog	Observations of Emission-line Stars in the Orion Region I. The Kiso Area A-0904	Wiramihardja et al. (1989)
B/denis/denis	DENIS: A Deep Near-Infrared Survey of the southern sky	Epchtein et al. (1994)
J/AJ/118/1043/tables	A photometric catalog of Herbig Ae/Be stars…UX Orionis stars.	Herbst & Shevchenko (1999)
J/AJ/119/3026/table1	Circumstellar disc candidates identified in the Orion nebula cluster flanking fields.	Rebull et al. (2000)
J/AJ/121/3160/table4	Near-infrared photometric variability of stars toward the Orion A molecular cloud.	Carpenter et al. (2001)
J/AJ/124/1001/table5	Near-infrared photometric variability of stars toward the Chamaeleon I molecular cloud.	Carpenter et al. (2002)
J/AJ/123/387/phot	A photometric study… $\lambda$ Orionis star-forming region.	Dolan & Mathieu (2002)
J/AJ/123/304/table1	The young cluster IC 5146	Herbig & Dahm (2002)
II/264/var	The All Sky Automated Survey. Catalog of Variable Stars	Pojmanski (2002)
2mass-psc	2MASS All-Sky Catalog of Point Sources	Cutri et al. (2003)
J/AJ/125/1537/table2	Deep near-infrared observations… in Orion Molecular Clouds 2 and 3	Tsujimoto et al. (2003)
J/AN/325/705/table2	Optical and infrared photometry… $\sigma$ Ori cluster	Béjar et al. (2004)
J/A+A/417/557/table4	A Rotational and Variability Study of a Large Sample of PMS Stars in NGC 2264	Lamm et al. (2004)
J/AJ/128/2316/table4	The low-mass population of Orion OB1b. I. The $\sigma$ Ori cluster	Sherry et al. (2004)
J/AJ/128/1684/table1	The initial mass function and young brown dwarf candidates in NGC 2264. I.	Sung et al. (2004)
J/AJ/129/829/table1	The T Tauri star population of the young cluster NGC 2264.	Dahm & Simon (2005)
J/MNRAS/356/89/catalog	Membership, binarity and accretion… $\sigma$ Ori cluster	Kenyon et al. (2005)
II/271A/patch2	TASS Mark IV Photometric Survey of the Northern Sky	Droege et al. (2006)
J/ApJ/662/1067/sigori	A Spitzer Space Telescope study of disks in the young $\sigma$ Ori cluster	Hernández et al. (2007)
J/MNRAS/375/1220/ngc2264	Empirical isochrones and relative ages for young stars… NGC2264	Mayne et al. (2007)
J/MNRAS/375/1220/sigori	Empirical isochrones and relative ages for young stars… $\sigma$ Ori	Mayne et al. (2007)
II/319/gcs9	UKIRT Infrared Deep Sky Survey (UKIDSS) Release 9	Lawrence et al. (2007)
II/316/gps6	UKIRT Infrared Deep Sky Survey (UKIDSS) Galactic Plane Survey (GPS) Release 6	Lucas et al. (2008)
J/ApJS/183/261/table2	A multi-color optical survey of the Orion nebula cluster. I. The catalog.	Da Rio et al. (2009)
J/ApJS/191/389/table2	Precision photometric monitoring of very low mass $\sigma$ Orionis cluster members	Cody & Hillenbrand (2010)
J/ApJ/726/18/stars	Disk evolution in W5: intermediate-mass stars at 2-5 Myr	Koenig & Allen (2011)
J/ApJ/733/50/table2	YSOVAR: …mid-infrared photometric monitoring of the ONC	Morales-Calderón et al. (2011)
J/A+A/525/A47/tablec2	U-band study of the accretion properties in the $\sigma$ Ori star-forming region	Rigliaco et al. (2011)
J/AZh/88/143/oricat	Proper motions of stars in the region of the Great Nebula in Orion	Vereshchagin & Chupina (2011)
J/A+A/548/A79/omc-var	INTEGRAL-OMC optically variable sources	Alfonso-Garzón et al. (2012)
J/ApJ/751/22/table3	A significant population of candidate new members of the $\rho$ Ophiuchi cluster	Barsony et al. (2012)
J/ApJ/752/59/members	The low-mass stellar population in L1641…stellar initial mass function	Hsu et al. (2012)
II/313/table3	The Palomar Transient Factory photometric catalog 1.0	Ofek et al. (2012)
J/ApJ/754/30/youth	New isolated planetary-mass objects…$\sigma$ Orionis cluster.	Peña Ramírez et al. (2012)
J/MNRAS/434/806/catalog	Pre-main-sequence isochrones. II. Revising star and planet formation time-scales.	Bell et al. (2013)
J/ApJS/209/28/table2	The MYStIX wide-field near-infrared data: optimal photometry in crowded fields	King et al. (2013)
II/321/iphas2	The Second Data Release of the INT Photometric H-Alpha Survey…	Barentsen et al. (2014)
J/A+A/564/A29/catalog	Orion Revisited: II. The foreground population to Orion A	Bouy et al. (2014)
J/ApJ/794/36/table2	A spectroscopic census in young stellar regions: the $\sigma$ Orionis cluster.	Hernández et al. (2014)
J/other/RAA/14.1264/table1	V350 Cep UBVRI long-term photometry	Ibryamov et al. (2014b)
I/327/cmc15	The Carlsberg Meridian Catalogue, final data release (MC15)	Muiños & Evans (2014)
V/147/sdss12	The SDSS Photometric Catalogue, Data Release 12 (DR12)	Alam et al. (2015)
J/AJ/149/200/table4	An optical survey of the partially embedded young cluster in NGC 7129	Dahm & Hillenbrand (2015)
II/336/apass9	AAVSO Photometric All Sky Survey (APASS) DR9	Henden et al. (2015)
J/ApJS/220/17/table2	Wide-field infrared polarimetry of the $\rho$ Ophiuchi cloud core	Kwon et al. (2015)
J/AJ/150/132/variables	Near-infrared variability in the Orion Nebula Cluster	Rice et al. (2015)
J/A+A/581/A140/tablea2	The VISTA Orion mini-survey: star formation in the Lynds 1630 North cloud.	Spezzi et al. (2015)
II/349/ps1	Panoramic Survey Telescope and Rapid Response System (Pan-STARRS) DR1	Chambers et al. (2016)
II/341/vphasp	The VST Photometric Halpha Survey… (VPHAS+) DR2.	Drew et al. (2016)
J/A+A/587/A153/science	VISION - Vienna survey in Orion. I. VISTA Orion A Survey.	Meingast et al. (2016)
II/348/vvv2	VISTA Variable in the Via Lactea Survey DR2	Minniti et al. (2017)
not in Vizier	Northen Sky Variability Survey (NSVS)	Woźniak et al. (2004)
not in Vizier	The Catalina Real-Time Transient Survey (CRTS)	Djorgovski et al. (2011)
not in Vizier	Panoramic Survey Telescope and Rapid Response System (Pan-STARRS) DR2	Chambers et al. (2016)
not in Vizier	The HOYS-CAPS Citizen Science Project	Froebrich et al. (2018)
not in Vizier	The American Association of Variable Stars Observers (AAVSO)	Kafka et al. (2018)
not in Vizier	SkyMapper Southern Sky Survey	Wolf et al. (2018)

Equations

$${\Delta }_{\mathrm{YSO}}=(\mathbf{x}-\mathbf{u}{)}^T{\mathbf{\Sigma }}^{-1}(\mathbf{x}-\mathbf{u})-{\chi }^2(5\sigma {)}_{k=4},$$

$${\Delta }_{\mathrm{YSO}}>0,$$

$$P(k|R)=\frac{(Rt{)}^k}{k!}{e}^{-Rt}.$$

$$p(\tau |k)=\frac{p(\tau )P(k|\tau )}{P(k)}=\frac{p(\tau )}{P(k)}\frac{(Nt/\tau {)}^k}{k!}{e}^{-Nt/\tau }.$$

$$P(k)={\int }_{\tau =0}^{\infty }\frac{(Nt/\tau {)}^k{e}^{-Nt/\tau }}{k!}p(\tau )\mathrm{d}\tau =\frac{Nt}{k(k-1)}p(\tau ),$$

$$p(\tau |k)=\frac{(Nt{)}^{k-1}}{(k-2)!{\tau }^k}{e}^{-Nt/\tau }.$$

$$\tau =\frac{Nt}{k}.$$

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

$$\tau =\frac{Nt}{(k+1)}.$$

Full text

††pubyear: 2019††pagerange: Determining the recurrence timescale of long-lasting YSO outbursts–C

1 Introduction

The crucial period for planet formation is between $\simeq$1 Myr (when the first planetesimals in our own solar system were formed Pfalzner et al., 2015), and $\simeq$10 Myr when the protoplanetary discs around most stars have dissipated (e.g. Bell et al., 2013). This corresponds exactly to the Class II phase of young stellar evolution (see e.g. Dunham et al., 2014, for a definition of the stages of young stellar evolution). Therefore properties of the disc, such as surface density or temperature, play a key role in the formation and evolution of protoplanets. These properties determine, for example, the outer boundary condition for the gas accretion rate of the protoplanets, enter into migration rates or determine the location of the snowline, the latter having an impact on the surface density of solids (see e.g. Mordasini et al., 2015).

The aforementioned disc characteristics depend on the accretion rate from the disc onto the central star. Theoretical models show that accretion onto the central star is unlikely to be a steady process (see e.g. Vorobyov & Basu, 2015), with observational support arising from the analysis of knots along jets from young stars (e.g. Makin & Froebrich, 2018) as well as outbursts in YSOs (e.g. Hartmann & Kenyon, 1996). However, planet formation models fail to include this effect and generally assume that the accretion rate decreases steadily with time (e.g. Mulders et al., 2015), with some models including a dependence with the mass of the central star (Kennedy & Kenyon, 2008).

The changes in the accretion rate through circumstellar discs (if they occur) could have a profound effect on the emerging planetary systems. For example Hubbard (2017b) argues that large accretion events can help to solve the so-called meter barrier problem (Boley et al., 2014) and allow the in-situ formation of rocky planets even at distances of $\simeq$1.5 AU. It also helps to explain the small size of Mars (Chambers, 2014). These events can also promote the formation of Jupiter and Saturn-like gas giants (Hubbard, 2017a). Finally, the increase of the central luminosity due to changes of the accretion rate (as small as a factor of two) will drive the snow lines for various ices (used to predict the composition of planets) towards larger radii (Garaud & Lin, 2007). This effect has already been observed in the known eruptive young stellar object (YSO) V883 Ori, where the water snow line is found at a distance of 42 AU from the central star, far beyond the expected location for a 1 M*⊙* star (Cieza et al., 2016). The observations of calcium-aluminium rich inclusions in chondrites (Wurm & Haack, 2009) and the depletion of lithophile elements in Earth (Hubbard & Ebel, 2014) could be evidence for large eruptions in our own solar system.

Variable accretion in YSOs can be caused by a variety of different physical mechanisms and there may exist a continuum of outbursting behaviour with a range of amplitudes and timescales (Cody et al., 2017). Stochastic accretion bursts in an instability-driven accretion regime lead to changes in the accretion rate up to 700$\%$ over 1-10 days (e.g. Venuti et al., 2014; Cody et al., 2017). Large changes of accretion rate, increasing from typical rates of 10*-7* M*⊙* yr*-1* up to 10*-4* M*⊙* yr*-1*, with outburst durations of weeks to 100 years, which can increase the central luminosity to 100L*⊙* (e.g. Hartmann & Kenyon, 1996), are observed in eruptive YSOs (the FUors, EXors and MNors, e.g. Audard et al., 2014; Contreras Peña et al., 2017b). Here, some kind of disc instability, such as gravitational instabilities Vorobyov & Basu, 2005 or planet-induced thermal instabilities (e.g. Lodato & Clarke, 2004), are suggested as explanations of the large variability.

However, many aspects of the so-called eruptive YSOs are still uncertain. The frequency and amplitude of the outbursts is not well constrained (e.g. Scholz et al., 2013). In addition, there is controversy as to whether the very largest outbursts are associated with the Class II planet building phase at all, or are just limited to the pre-planet-forming (Class 0/I) phase (Sandell & Weintraub, 2001, c.f., Miller et al., 2011, with).

In this paper we present the results of the study we conducted using Gaia and the Schmidt photographic plate surveys for a large sample of YSOs. It is our aim to answer the questions of outburst amplitudes and frequencies during the planet-forming stage of young stellar evolution. This paper is divided as follows: in Section 2 we present the data used to maximise the time baseline and the sample size as well as describing how the sample was classified as class II YSOs. In Section 3 we present an outline of our method to select high-amplitude variability and discuss the issues we found when applying the method to our sample. Section 4 presents a discussion on the distribution of amplitudes for the variability of class II YSOs. In Section 5 we show the steps taken to select the YSOs showing long-term outbursts where variability is most likely driven by dramatic episodes of enhanced accretion. From the latter sample, in Section 6 we determine the outburst rate during the class II stage and compare it to previous theoretical and observational estimates. In section 7 we present a number of caveats that could affect our estimate of the outburst rate, whilst in Section 8 we discuss the implications of our result on the likely physical mechanisms driving the outburst during the class II stage. Finally, Section 9 shows a summary of our results.

2 Data

To establish the recurrence rate of accretion related outbursts during the planet-forming stage, it is important to maximise both the time baseline and the number of YSOs surveyed (see e.g. Hillenbrand & Findeisen, 2015). In this paper we achieve both by comparing the magnitudes of a large sample of class II YSOs from the photographic atlases provided by the SuperCOSMOS Sky Survey (hereafter SSS, Hambly et al., 2001) with the latest Gaia data release.

We note that the current classification system of eruptive YSOs is a largely phenomenological one based on photometric and spectroscopic characteristics, with long-term outbursts falling into the FU Orionis (FUor) class and short-term, repetitive outbursts classified as EXors. However, the more recent data show characteristics that have been difficult to classify into the original sub-classes (see e.g. Contreras Peña et al., 2017b). Given this, during the remainder of the paper we will not ascribe a class to the dramatic accretion events that we are searching for, but we will simply use a physical classification and refer to them as high-amplitude, long-term accretion events.

2.1 Time baseline

The SSS has digitised the entire sky in three colours ($B,R$ and $I$) whilst providing a second epoch in $R$, by scanning Schmidt photographic plates (Hambly et al., 2001). The observations took place during the second half of the twentieth century, reaching a depth of $R\simeq 20$ mag. The photographic surveys included in SuperCOSMOS are the SERC-J/EJ (Cannon, 1984), SERC-ER (Cannon, 1984), AAO-R (Morgan et al., 1992), ESO-R (West, 1984), POSS-I E (Minkowski & Abell, 1963) and POSS-II B, R and I (Reid et al., 1991) surveys. For a more detailed description see table 1 of Hambly et al. (2001).

The Gaia mission (Gaia Collaboration et al., 2018) has provided the first all-sky photometric survey that is able to match the depth of the photographic Palomar and ESO/SERC surveys. This provides the chance of studying the long term variability of approximately 1 billion stars with a minimum and maximum time baselines of 15 and 66 years respectively.

2.2 The sample

2.2.1 Initial selection

To obtain the largest possible sample, we assembled a very inclusive sample of YSOs, and then identified the class II YSOs using their spectral energy distribution (SED).

The first step consisted of searching the SIMBAD database (Wenger et al., 2000) for Galactic sources with classifications that make them likely to be a young stellar object. We included objects with classification as pre-main sequence stars (pr*), pre-main sequence star candidate (pr?), emission lines star (em*), variable star of Orion type (or*), Herbig Ae/Be star (Ae*), Herbig Ae/Be star candidate (Ae?), Herbig-Haro object (HH), FU Orionis star (FU*), T Tauri-type star (TT*), T Taur-type star candidate (TT?), young stellar object (YO) and young stellar object candidate (Y?). The search yielded 77873 objects.

The SIMBAD classification does not indicate the likely evolutionary stage of the objects in our sample. The observed SED of young stellar objects between 2 $<\lambda<$ 20 $\mu$m is generally used to determine the likely evolutionary stage of the system (e.g. Lada, 1987; Greene et al., 1994; Dunham et al., 2014). Objects where the infalling envelope has dissipated, but the star is still accreting material from a circumstellar disc, are known as class II YSOs, and it is during this stage where planet formation is believed to occur (see Section 1). Therefore we had to determine the fraction of our sample that show colours and SEDs consistent with those of class II YSOs.

2.2.2 Classification

We crossmatched our sample (r$<$1″) with several studies that provide a class for YSOs in known star formation regions111We note that we did not attempt to search for every existing catalogue that provides a YSO classification. However, using these catalogues we were still able to create a large sample of class II YSOs. (Evans et al., 2003; Chavarría et al., 2008; Connelley et al., 2008; Gutermuth et al., 2008; Koenig et al., 2008; Evans et al., 2009; Gutermuth et al., 2009; Kirk et al., 2009; Billot et al., 2010; Rebull et al., 2010, 2011; Rivera-Ingraham et al., 2011; Allen et al., 2012; Megeath et al., 2012; Marton et al., 2016). All of these provide a classification based on the near- to mid-IR SEDs of the YSOs. We found that most of the objects with a crossmatch were contained within the studies of Marton et al. (2016), Megeath et al. (2012), Evans et al. (2003) and Gutermuth et al. (2009). The Marton et al. study is based on the all-sky survey WISE and as such contains the largest number of sources. Given this, for better consistency when trying to classify our sample we gave the Marton et al. (2016) study preference over Megeath et al., Evans et al. and Gutermuth et al. as these are based on specific areas of star formation. The classification of YSOs in our sample is explained below.

Marton et al. (2016). Through the use of a supervised learning algorithm Marton et al. (2016) classify objects as YSOs based on 2MASS $JHK_{\rm s}$ and WISE (Wright et al., 2010) photometry. However, the study does not provide a direct classification as class II YSOs, and instead yields a classification either as class I/II or class III YSOs. To identify the class II YSOs, we made use of the classification criteria of Koenig & Leisawitz (2014) which are based on WISE bands $W1$ (3.4 $\mu$m), $W2$ (4.6 $\mu$m) and $W3$ (12 $\mu$m) bands (see Fig. 1). The latter yielded 3253 class II YSOs. We note that Koenig & Leisawitz (2014) define different criteria to discard possible contaminant sources (AGNs, resolved PAH emission, among others) before attempting to classify sources in their sample as YSOs, however pre-selection using the Marton et al. catalogue has ensured our sample comprises true YSOs and as such we do not believe the criteria select contaminant objects.

Megeath et al. (2012). This study is based on 2MASS $JHK_{\rm s}$ (Skrutskie et al., 2006), IRAC (Fazio et al., 2004) and MIPS (Rieke et al., 2004) Spitzer photometry of the Orion A and B molecular clouds. The authors provide a value for the slope of the SED at the IRAC wavelengths (spectral index, $\alpha_{\mathrm{IRAC}}$). From the latter, 1671 YSOs show values -1.6 $<\alpha<$ -0.3 (see Fig. 1); following Greene et al. (1994) these were classified as class II YSOs. From those objects with no value of $\alpha$, 605 have reliable $K_{s}$, I1 (3.6 $\mu$m) and I2 (4.5 $\mu$m) photometry as well as values for the $K$-band extinction, $A_{Ks}$. We were thus able to classify them as class II YSOs using the criteria established in Appendix A of Gutermuth et al. (2009) (see Fig. 1). For the remaining objects which cannot be classified through these two methods, we found 86 objects that are classified as disc objects in Megeath et al. (2012), thus we assumed that these were very likely class II YSOs.

Gutermuth et al. (2009). The authors provide 2MASS $JHK_{\rm s}$, IRAC and MIPS Spitzer photometry for 36 young, nearby star-forming clusters, where they also provide an evolutionary class for YSOs in these areas (see their Appendix A). Based on this classification, we found 1060 class II YSOs that were added to our list.

Evans et al. (2003). The “Cores to discs” (c2d) Spitzer legacy project provides 2MASS $JHK_{\rm s}$, IRAC and MIPS Spitzer photometry for five large and nearby molecular clouds (see e.g Evans et al., 2009). The project also provides the slope of the SED, $\alpha$, using the region between 2 and 24 $\mu$m. Using the latter information we added 994 class II YSOs to our list.

For the objects in the SIMBAD list that were not found in any of the catalogues mentioned above, we established a likely class through the following steps.

WISE. We used the ALLWISE data release (Cutri & et al., 2013). We found 4883 objects that were detected in the $W1,W2$ and $W3$ bands with uncertainties below 0.4222We note that Koenig & Leisawitz (2014) require $\sigma<0.2$ in al three bands, however we used a less strict criterion in order to include more objects in our classification. mag in all three bands and that can be flagged as class II YSOs (see Fig. 1) using the criteria of Koenig & Leisawitz (2014).

2MASS. 1536 objects have reliable 2MASS photometry (quality flag “AAA”) and can be flagged as class II YSOs based on their location in the $J-H$ vs $H-K_{\rm s}$ colour-colour diagram (see Figure 1), i.e. they show larger $H-K_{\rm s}$ colours than main-sequence stars and have colours that are consistent with those of classical T Tauri stars (CTTS) as defined by the CTTS locus of Meyer et al. (1997).

MYStIX. Many YSOs found in the Massive Young Star-Forming Complex Study in Infrared and X-ray (MYStIX Feigelson et al., 2013) project are not part of the SIMBAD database. We found that 1400 of these were classified as being in the planet-forming stage, thus they were also added to our sample.

In summary we were able to compile a list of 15404 class II YSOs.

3 Gaia vs SuperCOSMOS search

3.1 Outline method

To search for high-amplitude variables in our sample we compared their Gaia magnitudes with those found in each individual SuperCOSMOS band. The SSS data were obtained from the merged source catalogue of SSS via a structured data language (SQL) query of the SuperCOSMOS Science Archive (SSA333http://ssa.roe.ac.uk/), whilst Gaia magnitudes were obtained from the second data release (DR2) catalogue downloaded to a local machine. The method was as follows.

For each YSO we searched for a Gaia and SSS counterpart, using a separation of 1.6″and 2″respectively. If the YSO was not found in Gaia, then we set the magnitude of the object to be $G=21$ to account for the possibility that the star was fainter than the Gaia limit at the time of observations. Whenever the YSO was not detected in an SSS band, we set the magnitude of the YSO to be that of the faintest SSS source with a stellar profile (class 2) in a $2^{\circ}\times 2^{\circ}$ box around the YSO in the corresponding band. 2. 2.

We crossmatched Gaia and SSS sources found in a $2^{\circ}\times 2^{\circ}$ box around each YSO, where SSS detections with stellar classification were selected. This allowed us to determine the relation between the Gaia $G$ broadband magnitude and the $B$, $R1$, $R2$ and $I$ magnitudes obtained from SSS. 3. 3.

For each SSS band (whose corresponding magnitude we will call $M_{\rm sss}$), we determined $G-M_{\rm sss}$, for all objects with Gaia magnitudes that were found within $|G-G_{\rm YSO}|\leq 0.2$, where $G_{\rm YSO}$ is the Gaia magnitude of the YSO being analysed. The mean value of the $G-M_{\rm sss}$ distribution gave us the expected difference between Gaia and SSS for (what we assume are largely) non-variable stellar sources in the region around the YSO. A large difference between the observed $G-M_{\rm sss}$ colour of the YSO vs the expected $G-M_{\rm sss}$ colour estimated from non-variable sources, gave us an indication as to whether the YSO had shown any variability between the photographic surveys and Gaia DR2.

We quantified the variability of the YSO by treating the $G-M_{\rm sss}$ distribution as a 4-dimensional multivariate normal distribution. We define an statistic for the variability of the source, $\Delta_{\mathrm{YSO}}$, as

[TABLE]

with $\boldsymbol{x}$ the vector of the observed $G-M_{\rm sss}$ colours for the YSO, $\boldsymbol{u}$ the vector of expected $G-M_{\rm sss}$ colours estimated for non-variable stars, $\Sigma$ is the covariance matrix determined from the expected $G-M_{\rm sss}$ colours, and $\chi^{2}(5\sigma)_{k=4}$ is the value of the $\chi^{2}$ distribution with 4 degrees of freedom, and for a confidence interval of 99.99996$\%$ or 5$\sigma$.

To illustrate the type of objects we were looking for, we analysed the known eruptive YSOs that were part of the Gaia alerts sample, V2492 Cyg (e.g. Hillenbrand et al., 2013), V350 Cep (Magakian et al., 1999) and ASASSN13-db (Sicilia-Aguilar et al., 2017). In Figure 3 we show the example of eruptive class I YSO V2492 Cyg. The YSO is not detected in the first epoch of $R$ (observations on September 1952), is faint at the second epoch (observations on July 1989) with $R=19.8$ mag and is bright during Gaia observations ($G\simeq 15$ mag). Figure 2 shows a 2-dimensional representation of our method applied to V2492 Cyg, where it can be seen that the observed $G-R1$ vs $G-R2$ colours of the object are located well beyond the 5$\sigma$ confidence intervals. We find similar results for V350 Cep and ASASSN13-db.

3.2 Finding a suitable $\Delta_{\mathrm{YSO}}$

We began by classifying the YSO as a variable star candidate if it satisfied the following condition

[TABLE]

which ensures the selection of high-amplitude variable stars.

To understand the types of variability that we would detect with this method, we analysed a subsample of 5000 YSOs444This sample is selected from some of the categories described at the beginning of Section 2. We note that these do not include the FU Orionis star class. as well as a sample of known variable YSOs taken from the Gaia Photometric Science Alerts555http://gsaweb.ast.cam.ac.uk/alerts (GSA Wyrzykowski et al., 2012) which include a number of objects that are known eruptive YSOs. We note that in both samples the evolutionary class of the individual YSOs was not taken into account.

This analysis yielded 24 bonafide high-amplitude ($\Delta$R$\geq 1.5$ mag) variable stars from the SIMBAD sub-sample and successfully retrieved the Gaia YSO alerts showing high-amplitude variability with $\Delta G>2$ mag.

We found that most variable stars from the SIMBAD and Gaia alerts sub-samples show low-amplitude, recurrent aperiodic variability. A literature search revealed that the variability in these low-amplitude objects is caused by physical processes other than variable accretion. The low-amplitude objects showed lower values of $\Delta_{\mathrm{YSO}}$, whilst a few of these had positive values of $\Delta_{\mathrm{YSO}}$, but with individual values of $G-M_{\rm SSS}$, when analysed separately as a univariate distribution, that were found to be below $5\sigma$ in each band.

In the cases where we suspected variability was due to large accretion events, and in the known eruptive YSOs from the Gaia alerts, we found that $\Delta_{\mathrm{YSO}}$ showed a value much larger than 0. In addition the inspection of individual values of $G-M_{\rm sss}$, when analysed separately as a univariate distribution, were found to be above $5\sigma$ in at least one band.

Given the results from this early analysis, when applying our method to the class II YSO sample, we defined a new set of conditions for an object to be included as a candidate variable star. This was done to reduce the number of objects to inspect (see below) and to try and include only the most extreme cases of variability, and thus more likely related to large and long-lasting changes in the accretion rate. The revised conditions were as follows.

$\Delta_{\mathrm{YSO}}\geq 70$ and the individual values of $G-M_{\rm sss}$, when analysed separately as a univariate distribution, to be above $5\sigma$ in at least one band. 2. 2.

$\Delta_{\mathrm{YSO}}\geq 100$ and the individual values of $G-M_{\rm sss}$, when analysed separately as a univariate distribution, to be below $5\sigma$ in each band.

3.3 Problems with the photographic images images

Through the initial analysis of the two YSO subsets we were able to recognise a major source of contamination in our variable star candidates list. We found that in many objects the non-detection of the source in the SSS catalogues was probably due to problems with the plate images rather than true variability of the source. These issues were most apparent in the areas with high extinction, and were more evident in the $B$ and $R$ SSS images. In order to define whether the variability was driven by this issue, we first queried Gaia DR2 and the SSA to select sources within a box size of $3\arcmin\times 3\arcmin$ around the YSO being analysed. Then we crossmatched both catalogues to determine the number of objects detected both in SSS and Gaia, $\mathrm{N}_{\rm sss}$. In individual filters, YSOs found in problematic areas showed a low value for the ratio $\mathrm{N}_{\rm sss}/\mathrm{N}_{\rm Gaia}$, with $\mathrm{N}_{\rm Gaia}$ the total number of Gaia sources in the $3\arcmin\times 3\arcmin$ box. Therefore, if the individual bands driving the variability of the YSO showed $\mathrm{N}_{\rm sss}/\mathrm{N}_{Gaia}<0.45$ then this object was classified as unlikely to be a variable star.

3.4 Visual inspection of variable star candidates

We next used the results from our subsamples to inform a search for high-amplitude variable stars using all the 15404 YSOs in our sample. The method yielded 4815 objects that fulfilled the condition established by Equation 2. This number reduced to 1576 objects when we imposed the conditions of Section 3.2. From the latter, 501 objects were flagged as being in the list due to the problems with photographic plates explained in Section 3.3. Thus, we were left with 1075 variable YSO candidates.

To confirm or discard its variability, each of the 1075 variable candidates was analysed individually. The analysis included the inspection of the $3\arcmin\times 3\arcmin$ cutout images from the photographic plates provided by the SSA, as well as comparing with the more recent images (when available) provided by the Panoramic Survey Telescope and Rapid Response System (Pan-STARRS, Chambers et al., 2016) and the SkyMapper Southern Sky Survey (Wolf et al., 2018). In addition we also searched for additional photometry from publicly available catalogues through the VizieR catalog access tool (Ochsenbein et al., 2000). The list of the individual surveys found in Vizier that were used in our analysis is given in Appendix C. This final step revealed further sources of contamination that were not considered before the cleaning of the final sample. Most of these related to problems with bright objects, crowding or incorrect magnitudes in SSS catalogues.

After inspection of the individual candidates, we found 139 true high-amplitude variable stars. In each case, the additional photometry from other surveys allowed us to classify the objects according to their light curves. We note that this is based only on a handful of epochs, therefore in many cases our classification is uncertain and is thus followed by a ? sign. This classification aims to select those objects where variability resembles that of long timescale outbursts. An example of the photometry obtained for each high amplitude variable is presented in Table 1. The full light curve data are presented as supplementary material.

The list of 139 high-amplitude variable stars selected with our method is shown in Table 2. Column 1 shows the designation given by us. Since the objects are ordered by amplitude, the number gives an indication of their level of variability. Columns 2 and 3 mark the right ascension and declination of the object respectively. The SIMBAD identification and type are shown in columns 4 and 5, whilst column 6 gives the observed ($R$) amplitude taking into account only the SSS and Gaia epochs (after transforming the Gaia $G$, SSS $B$ and $I$ into $R$-band magnitudes). We note that these are approximate values. In column 7 we present a classification for the behaviour of the light curve of the YSO (Section 3.4). Column 8 provides a final classification of the light curve, which is obtained either by the analysis done in this work (see Section 5) or by information found in the literature (e.g. V26 is a known UXor), whilst column 9 shows the references to such classification. In column 10 we mark whether the object is a known variable star, whilst column 11 shows whether the SSS vs Gaia variability relates to the known variability of the source of if our detection relates to previously unknown variability of the source. Finally column 12 shows the source that was used to classify the YSOs as being in the class II stage (see Section 2.2.2).

Our sample contains 66 new variable stars, which represents 48$\%$ of our sample. From the known variable stars (72 objects) we find that in 17 objects the SSS vs Gaia variability has not been reported previously. We note that upon further inspection objects V5, V12, V24, V32, V45 and V84 are found to be more likely non-YSOs. Therefore, although these are shown in Table 2, we did not include these objects in the following analysis of variability in the class II YSOs.

4 The distribution of class II amplitudes

The top plot of Figure 4 shows the histogram of the amplitudes of our high-amplitude variable stars as determined from the SSS and Gaia epochs. We note that in the histogram we did not include the objects that were found to be likely non-YSOs. In addition we also did not include the three long-lasting outbursting YSOs that are more likely class I YSOs (see Section 5.4).

In the histogram we can see that we do not detect objects with amplitudes below 1 mag. In addition the sample appears incomplete below 2 magnitudes.

To obtain the distribution of amplitudes for low amplitudes we used the dataset of $i$-band photometry of Cep OB3b presented in Sergison et al. (in prep). We took the objects classified as Class II on the basis of their Spitzer photometry (Allen et al., 2012), which are in Sergison et al’s primary sample. Importantly, variability has played no part in the selection of these stars. We took one long observation from each of 2004, 2005, 2007 and 2013 as the closest analogue we could obtain to the sampling in this paper. For each star which had “good” photometry (as defined by Sergison et al) for all four datapoints we then calculated a full amplitude.

To compare both samples we needed to estimate the fraction of stars per unit amplitude (d$f$/d$A$). However, given the observed amplitudes of our YSO variables, we considered the possibility that some objects in our class II YSO sample were too faint to have been detected in Gaia, even if they showed large amplitude outbursts. We found that 3712 YSOs had no Gaia DR2 counterparts, but this was found to be consistent with the fact that they were also faint or not detected in SSS, so they were flagged as non-variable stars.

There is no detailed selection function for Gaia, but the star counts turn over at about $G\simeq 20.5$ mag (Gaia Collaboration et al., 2018). Hence, objects with $G\leq 22.5$ mag showing variability with $\Delta R\geq$2 mag should have been detected in Gaia DR2. To estimate an approximate $G$ magnitude for the 3712 YSOs that are not in Gaia DR2, we obtained sloan-$r$, $J$ and $K$ photometry of the YSOs by crossmatching with catalogues from Pan-STARRS, SkyMapper, the SDSS photometric catalogue (Alam et al., 2015), DENIS (Epchtein et al., 1994), 2MASS Cutri et al. (2003), UKIDSS GPS (Lucas et al., 2008) and the VISTA Variable in the Via Lactea Survey (Minniti et al., 2017). Objects that were detected with $r<24$ mag would have been detected in Gaia if they had shown high-amplitude variability. Objects without $r$-band detections, do have $J$, and $K$ photometry, so we use the information from YSOs that were detected in Gaia DR2 to determine a relation between the $G-J$ and $J-K$ colours. The latter fit was then used to estimate an approximate value of $G$ for the objects that are not originally detected in Gaia DR2.

This method showed that 1318 YSOs had approximate $G$-band magnitudes that were fainter than 22.5 mag, so we believe these would have not been detected even if they had shown high-amplitude variability. Thus for the remainder of our analysis of variability, our total sample of class II YSOs reduces to 14086 objects.

The comparison of d$f$/d$A$ between our sample and that of Sergison et al. (see Figure 5) confirms that our sample is not complete below amplitudes of 1 magnitude. Both samples agree remarkably well at 2 magnitudes. The figure also shows that at $\Delta m>3.5$ mag the number of objects does not fall as expected from the rate of decline observed at lower amplitudes. Instead, the distribution reaches a plateaux, which could point to YSOs at these amplitudes being part of a different population of variable stars.

5 Selecting high amplitude accretion-driven variability

We detected a large number of high-amplitude variable YSOs. The next step consisted of determining what is the physical mechanism driving these large changes.

5.1 Physical mechanism

Variability is one of the defining characteristics of pre-main-sequence stars (e.g. Joy, 1945). This occurs on a broad range of timescales and amplitudes due to various physical mechanisms. Periodic modulation of the stellar flux resulting from the rotation of cool spots in the stellar photosphere is mostly observed in weakly accreting YSOs (the so-called weak-lined T Tauri stars or WTTS) and its characterised by low level variability in the order of 0.1 mag (Herbst et al., 1994). The irregular (although sometimes periodic) variability from short-lived hot spots at the stellar surface due to accretion can reach can 2-3 mags in extreme cases (see e.g. Grankin et al., 2007). Luminosity dips caused by circumstellar dust obscuration can last days to months and with amplitudes that are effectively limitless as they depend of the optical depth of the dust obscuring the central star (see e.g Carpenter et al., 2001). Finally, changes in the accretion rate from the disc onto the central star can also cause variability, with timescales of days to hundreds of years, in YSOs (see Section 1).

The amplitude of the variability of stars in our sample gives us an insight into the physical mechanism that drives the variability in these YSOs. We have previously discussed that at amplitudes below 2 magnitudes the sample is incomplete. However, this is not a problem for the purposes of this study as we are interested in the most extreme cases of accretion-related variability which we expect to have amplitudes greater than this level. We would not expect to detect YSOs where variability is caused by cool spots in the stellar photosphere, as amplitudes due to this effect are not expected to reach more than a few tenths of a magnitude. In addition our sample is comprised of class II YSOs, where this effect is not expected to dominate. The lower number of YSOs at amplitudes below 2 magnitudes and the non-detection at below 1 mag, allows us to conclude that we are not observing variability due to cool spots.

Figures 4 and 5 also show that the number of variable YSOs drops dramatically at amplitudes larger than 3.5 mag, but there appears to be a secondary peak at around 4.5 mag, although the position of the peak is highly dependent on the choice of bin size. However the cumulative distribution shown in Figure 4 does reveal a bump at around 4 mag. Thus the variable YSOs with amplitudes larger than 3.5 magnitudes could be a different population from the objects described above. Considering the photometric characteristics of known long-lasting YSO outbursts, it is at these amplitudes that we would expect to find the majority of objects that can be classified as showing long-term, accretion-driven outbursts.

5.2 Selecting long timescale variability

The selection of objects where variability is driven by extreme, long-term changes in the accretion rate was difficult. In many cases using just the data from SSS and Gaia (5 epochs of observations) could give a false impression that the object is an eruptive YSO. Therefore, we made use of all the available photometry (see Section 3.4) to discard or confirm a classification as long-lasting YSO outbursts. In many cases the classification was aided by a more exhaustive literature search.

In this analysis we found that the majority of YSOs with amplitudes between 1 and 3 mag (101 out of 108 objects) are characterised by light curves with fading events or irregular variability (i.e. going from bright to faint states repetitively across their light curves). The variability in this group is likely to be explained by hot spots due to accretion, representing extreme cases of this type of variability (see Grankin et al., 2007) or from dust obscuration, similar to the fading events observed in AA Tau (Bouvier et al., 2013) or in RW Aur (e.g. Bozhinova et al., 2016). The remaining 7 YSOs (V36, V51, V53, V56, V68, V85 and V94) were classified as a potentially having long lasting outbursts based on their light curves.

From the 28 objects with amplitudes of $\Delta R\geq 3$ mag (variable stars V1 to V4, V6 to V11, V13 to V23, and V25 to V31), 7 were classified as showing long-lasting outbursts. Amongst the 28 objects we also find objects with apparent long-lasting fading events, again likely related to extinction, and objects with irregular variability. However, given the larger amplitudes, variability due to hot-spots is less likely to explain this kind of variability. Instead these are more likely explained by short-term accretion events that appear as irregular given the long baseline of the light curves.

The variable star V20 (2MASS J20521294$+$4420534) illustrates the issues of classifying YSOs based only a handful of epochs. Figure 6 shows the light curve of the source using the SSS and Gaia epochs, in addition we also compare the $R$ and $I$ images from SSS. Based only on this information, this object appears to have gone into outburst after the second epoch of $R$-band observations, becoming $\simeq$3.2 mag brighter during the I SSS observations and remaining bright until the Gaia observations. However, the addition of literature data and inspection of images provided by the Palomar Transient Factory (PTF Law et al., 2009) reveals that the faint $R$ magnitudes are probably part of the “normal” variability of the source, which seems to be characterised by irregular, short-term variability (see the bottom images of Figure 6).

A second example is V3 (V1219 Cyg), which could be easily mistaken as showing a 4.7 mag outburst between the first and second epochs of $R$ SSS observations (see Figure 7), in addition the object drives an H2 outflow (Makin & Froebrich, 2018). However, further analysis reveals that this source has been classified in the past as showing repetitive high-amplitude fading events rather than a long timescale outburst (see Stecklum & Linz, 2013).

In both examples the light curves, given only the SSS and Gaia data, could easily be classified as long-term outbursts. We were able to discard such classification only through the inclusion of additional data from the literature.

Therefore in trying to find objects with long-lasting outbursts, we discovered a large population of high-amplitude variable YSOs that cannot be classified as showing long-term outbursts (see Section 5.3). Changes in the accretion rate could still be the underlying cause of such variability, but objects that show repetitive and/or short-term (lasting less than $\simeq$10 yr) outbursts, such as V18 or V8 (the known short-term eruptive YSO ASASSN13-db, see Sicilia-Aguilar et al., 2017), are not included in the subsequent analysis. However, we emphasise that these are an interesting class of objects in their own right.

In addition, many objects show apparent long-lasting fading events. This type of variability is explained by obscuration of the central star by structures in the disc located at large distances and above the midplane of the disc (see e.g. Bouvier et al., 2013). This is interesting in the context of discs and planet formation. The migration of giant planets due to interactions with the disc (also known as type II migration, Ida & Lin, 2004) can open up gaps in the disc which could lead to scaleheight bumps (Baruteau & Papaloizou, 2013) or to a misalignment between the inner and outer disc (Lubow & Martin, 2016), that could lead to fading events. We note, however, that mechanisms such as dusty outflows could also lead to these type of events (Davies et al., 2018).

Again, these objects are interesting findings, but further analysis of them is beyond the scope of this paper.

5.3 Outbursts

To determine an outburst rate from our analysis, we needed to select only objects where a classification as a long-lasting YSO outburst was very likely. In this sense YSOs where changes in the accretion rate might be the physical mechanism driving the variability, but which show repetitive and/or short-term outbursts are not included (see above). Additionally we did not include objects where lack of data did not allow a firm classification. For example, five out of seven objects with $\Delta R\leq 3$ mag and that were classified as potentially having long-lasting outbursts, were not included in our final sample as there was not enough data that would allow to confirm such classification.

Collating literature and database information in this way allowed us to classify 9 objects as being part of the eruptive YSO class showing long-lasting outbursts (the 7 objects with $\Delta R\geq 3$ mag classified as having long-lasting outbursts in Section 5.2 and two objects with $\Delta R<3$ mag, V36 and V51). These include 6 known eruptive YSOs, V1 (LkHA 190, most commonly known as V1057 Cyg, e.g. Hartmann & Kenyon, 1987), V2 (V2493 Cyg, most commonly known as HBC722, Semkov et al., 2010), V6 (V350 Cep, Ibryamov et al., 2014a), V10 (V733 Cep, Reipurth et al., 2007), V31 (UCAC2 30330609, most commonly known as V960 Mon, Hillenbrand, 2014) and V36 (V582 Aur, Ábrahám et al., 2018). There are 3 previously unknown eruptive YSOs, V4 (2MASS J02335340$+$6156501), V9 (2MASS J08410676$-$4052174) and V51 (WRAY 15-488). A more extensive explanation for the classification of the previously unknown objects is presented as individual notes in Appendix A.

5.4 Class I contamination

Due to the effects of geometry, reddening or disc inclination there might not be a direct correlation between the YSO class and its evolutionary stage (see e.g. Robitaille et al., 2006). The classification of objects in Section 2.2.2 was done using near- to mid-infrared data, which lead to a classification as class II YSOs, however, longer wavelengths could reveal the presence of envelopes (see also Section 7.3 for further discussions on the contamination from class I YSOs).

Therefore, before we determined the outburst rate from our selected YSOs, we took into consideration the possible contamination from YSOs at earlier evolutionary stages in our sample (class I or flat-spectrum sources). The concern arises from the fact that if the outbursts are $\simeq$10 times more common during the class I stage (see Section 6.1.1), then even a low contamination from class I YSOs would result in a high fraction of our outbursts sources being from the embedded phase.

To study this effect we performed a deeper analysis of the likely evolutionary stage of the YSOs in our long-term outburst sample. From the 9 YSOs in this sample, 8 were selected from the Marton et al. (2016) catalogue and given a class II YSO classification using the Koenig & Leisawitz (2014) criteria, whilst V4 is selected from the near-infrared excess as determined from its 2MASS colours. We determined the 2–24 $\mu$m spectral index, $\alpha$, for the observed spectral energy distribution of the burst sample, using 2MASS and WISE magnitudes. Table 3 lists the values determined for the sample. We see that $\alpha\leq-0.3$ for four objects in the sample is consistent with the definition for class II YSOs of Greene et al. (1994) . For five objects, the spectral index indicates an earlier evolutionary stage.

For objects where $\alpha$ is consistent with a class II stage, we find further evidence to support the idea that the envelope has dissipated and that these are pure disc-bearing objects. The interferometric 13CO and C18O survey of Fehér et al. (2017) determines that V10 (V733 Cep) and V2 (HBC722) have evolved beyond the class I and flat-spectrum stages. Ábrahám et al. (2018) determines that the pre-outburst spectrum of V36 (V582 Aur) is that of a classical T Tauri star. There is no literature information in the case of V4, however the spectral index of $\alpha=-0.4$ is an upper limit and makes this object a very likely class II YSO.

In the case of objects with $\alpha>-0.3$ we find that V1 (V1057 Cyg) shows evidence for a large and dense envelope (Fehér et al., 2017), whilst the SED of V31 (V960 Mon) suggests the presence of an envelope in Kóspál et al. (2015). V6 (V350 Cep) is classified as a protostar candidate in Muzerolle et al. (2004) (based also on the 70$\mu$m magnitude of the YSO) and Miranda et al. (1994) establish that the object is embedded in nebulosity. Thus, these 3 YSOs were not included in our calculation of the outburst rate (see Section 6).

The value of spectral index for V9 and V51 is not within the class II limits of the Greene et al. (1994) classification. However, the value is not far from those of class II YSOs, and we note that we did not consider the effects of reddening in the calculation of $\alpha$, which would likely place the objects in the class II definition. In addition is hard to determine whether the slope of the SED has changed due to the outburst. In the case of V9 we find that from observations prior to outburst, the object is an H$\alpha$ emission line star from Pettersson & Reipurth (1994). The authors provide a measure of the strength of the H$\alpha$ line on a scale from 1 (faint) to 5 (very strong), where V9 was listed with I(H$\alpha)=3$. We can estimate the equivalent width (EW) this corresponds to using objects from Pettersson & Reipurth (1994) with measured EWs which also have I(H$\alpha$) classes (see their table 3). YSOs that show the same value for the intensity of the emission line, have EWs that go between 10 and 100Å  with an average of 38Å. Thus, is very likely that the progenitor of V9 is in fact a disc-bearing (class II) YSO. Therefore the object was retained in the outburst sample for the remainder of our analysis. V51 was also originally classified as an emission line star prior to outburst (Wray, 1966). In addition the object was classified as a classical T Tauri star (from it SED during outburst, see Appendix A.3) by Gregorio-Hetem & Hetem (2002). Therefore we included both objects in our sample of class II outbursts. Thus we have a final sample of 6 large-amplitude, long-timescale accretion driven class II YSO outbursts.

6 Outburst Rate

To determine the timescale of high-amplitude, long-lasting outbursts events during the class II stage, we first obtained the baseline of observation for each YSO in our sample. Therefore, we queried the SSA database for the epoch of observation for the corresponding B, R1, R2 and I photographic plates. The baseline of observations for each star was given by the difference between the Gaia DR2 epoch (set at 2015.5) and that of the oldest photographic plate observation. The total time covered was calculated as the simple sum of all of the baselines, which yields 782959 years or a mean baseline of 55.6 years.

We can place formal limits on the “waiting time" ($\tau$) between outbursts, provided that our model is that all YSOs have quiescent periods of order centuries or longer, punctuated by outbursts where the star rises by more than four magnitudes on a timescale of around a decade or less, and then remains bright for decades. If that is the case, then the probability of observing $k$ rises, given an event rate for the whole sample of $R$ yr*-1*, and a total observation time $t$ is given by the Poisson distribution

[TABLE]

This is formally derived by considering the binomial distribution and approximating for low event rates. We are interested not in the event rate for the whole sample, but in the waiting time for a single star $\tau$=$1/(RN)$. Substituting this into Equation 3, we then used Bayes’ theorem to derive the probability density function $p(\tau|k)$ of a particular waiting time $\tau$, given $k$ as

[TABLE]

We then evaluated $P(k)$ in the normal way as

[TABLE]

where we assumed that the prior is independent of $\tau$ (an assumption we shall return to later), and hence solved the integral using the $\Gamma$-function. This yielded the result that

[TABLE]

Differentiating this to find the turning point gave the unsurprising result that the most probable value of $\tau$ is given by

[TABLE]

We then calculated the confidence limits as the values of $\tau$ which enclosed the most likely 68 percent of the probability given by Equation 6. We achieved this by integrating Equation 6 over all values of $\tau$ where $p$ exceeded a certain value. We then decreased that value of $p$ until the integral reached 0.68, at which point the extreme values of $\tau$ are our 68 precent confidence limits. Our observed values of $N$=14 077, $t$=55.6 years and $k$=6 gave $\tau$=$131^{+91}_{-48}$ kyr (68 percent confidence).

To test the effect of priors, we repeated the calculation using $p(r)\propto 1/\tau$, which gave

[TABLE]

and the most probable value of $\tau$ as

[TABLE]

These yielded $\tau$=$112^{+68}_{-38}$ kyr (68 percent confidence). The prior is equivalent to assuming that equal decades of $\tau$ have equal probabilities of containing $\tau$. This is normally thought a more reasonable prior for a quantity which must be greater than zero, and so this is our preferred result, although the prior does not have a great impact.

6.1 Comparison with literature data

The number we obtain in the above analysis represents the first observational determination that long-lasting outbursts do in fact occur during the class II stage and the frequency of such events. To our knowledge, there is only one other observational study, from Scholz et al. (2013), that provides an outburst recurrence timescale for YSOs. However, our analysis differs from Scholz et al. (2013) as the authors do not divide their YSO sample into different classes.

In this subsection we compare our results with that of Scholz et al. (2013) as well as comparing it with other observational and theoretical evidence that give some estimate of the frequency of YSO outbursts.

6.1.1 Scholz

et al. (2013)

With the aim of determining the outburst recurrence timescale in YSOs Scholz et al. (2013) searched for young eruptive variables by comparing Spitzer with WISE photometry, thus providing a 5 year baseline. The authors performed their analysis for about 8000 YSOs divided into two samples that were analysed separately. Sample A contains 4000 YSOs selected from c2d and clusters surveys (Gutermuth et al., 2009), with additional objects found in the NGC2264, Taurus and North American/Pelican Nebulae regions. Sample B comprises YSOs found in the Robitaille et al. (2008) catalogue of intrinsically red objects, in which Scholz et al. find $\simeq 4000$ YSOs.

In their analysis Scholz et al. (2013) did not divide objects into different YSO classes, but we require such a division to compare their data with our results. According to the authors, sample A contains about 2800 class II YSOs. In this sample they detected one very likely outburst candidate, V2492 Cyg, which is a class I YSO (Kóspál et al., 2011). Given the lack of class II outbursts, we cannot use the same arguments from Section 6 as the integral of Equation 8 over $\tau>0$ is not finite for k=0. All we can say is that the outburst rate is probably longer than the monitored YSO years from sample A, i.e. $5\times 2800=14$ kyr. From sample A, we can only estimate the outburst rate during the class I stage. To obtain this we have used the same $1/\tau$ prior we preferred in Section 6, but also note that without this prior we cannot obtain a confidence limit as the integral of Equation 6 over $\tau>0$ is not finite for $k=1$. This yields a timescale of $3.0^{+12.7}_{-1.9}$ kyr for the class I stage.

In sample B Scholz et al. (2013) find two outburst candidates, 2MASS J16443712-4604017 (hereafter SFW13 YSO1) and 2MASS J15111357-5902366 (hereafter SFW13 YSO2). Using the multi-epoch $K_{\rm s}$ photometry from the VISTA Variables in the Via Lactea survey and its extension (VVV and VVVx respectively, see e.g. Saito et al., 2012; Bica et al., 2018) we determined that SFW13 YSO2 showed high-amplitude recurrent variability over the 2010-2018 period, whilst SFW13 YSO1 has remained at a bright state since 2010 until at least the latest epoch of VVVX (May 2018, see Figure 8), making the YSO a true eruptive variable. Using the Koenig & Leisawitz (2014) classification criteria, we determined that SFW13 YSO1 is a class II YSO from its WISE colours. However the object shows a rising SED (Figure 8) with an spectral index, $\alpha=0.12$ that is more consistent with a class I YSO. Thus, we assume the latter classification for the object.

An analysis of near- to mid-infrared colours of candidate YSOs from red objects in Robitaille et al. (2008) shows that $\simeq 40$% can be classified as class II YSOs. Assuming the same fraction for sample B in Scholz et al. (2013) implies that around 1600 objects are class II YSOs. Since the only outburst in sample B, SFW13 YSO1, is more likely a class I YSO, we can only estimate that the recurrence timescale in the class I stage, using the $1/\tau$ prior, of $6.0^{+25.3}_{-3.9}$kyr. Similarly to the results of sample A, we can only say that the outburst recurrence timescale during the class II stage is longer than $\simeq$8 kyr.

The lack of class II outbursts in the samples of Scholz et al. (2013) agrees with the results from our work, i.e. given the outburst rate determined by us we would not expect to observe an outburst for a sample of a few thousand class II YSOs over a 5 yr baseline.

The analysis of the Scholz et al. (2013) yields another important result, as samples A and B show that during the class I stage, outbursts due to variable accretion appear to be $\approx$10 times more frequent than during the class II phase.

6.1.2 Theoretical estimate

Bae et al. (2014) find that gravitational instability(GI)-induced spiral density waves can heat the inner disc and trigger magneto-rotational instabilities (MRI) that will lead to accretion outbursts (or GI+MRI mechanism). They find that this process is more efficient at earlier stages when the envelope dominates the disc, with the number of outbursts decreasing as the star approaches the post-infall phase. Based on figure 3 of Bae et al. (2014), Hillenbrand & Findeisen (2015) estimate an outburst rate of $8\times 10^{-5}$ year*-1* star*-1* during the class I stage, and $3\times 10^{-6}$ year*-1* star*-1* during the class II stage, or a recurrence timescale of 12.5 and 333 kyr respectively. The value for the class I stage is consistent with our class I estimates. The theoretical estimate for the class II stage agrees with our results in the fact that is much longer than that of the class I stage, and that outburst occur over scales longer than 100 kyr. However, it is still longer than the recurrence timescale from our study, and is formally rejected at a 95.3% confidence level.

6.1.3 Other observational evidence

The recurrence timescale from our study is larger than the time between ejection events of $\simeq 1000$ years determined from the observed gaps between H2 knots along jets from YSOs (Ioannidis & Froebrich, 2012; Froebrich & Makin, 2016; Makin & Froebrich, 2018) and the 100 year timescale observed in MNors by (Contreras Peña et al., 2017a). However this is not surprising as H2 jets likely trace the accretion related events during the class I phase (Ioannidis & Froebrich, 2012) and that the young eruptive variables in Contreras Peña et al. (2017a) are biased towards class I YSOs.

Fischer et al. (2019) studied the large amplitude mid-infrared variability for a sample of 319 protostars (YSOs with class designations 0, I or flat) in Orion, by comparing Spitzer versus WISE photometry. Based on the recovery of two outbursts over the 6.5 yr baseline Fischer et al. (2019) determine an outburst rate of 1000 yr with a 90% confidence interval of 690 to 40300 yr. This agrees with our results that suggest that outbursts are more frequent towards YSOs at earlier evolutionary stages the class II stage.

6.1.4 Gaia17bpi and the rate from the Gaia alerts

At the time of writing, Hillenbrand et al. (2018), based on photometry from Gaia, reported the discovery of a new YSO outburst that has the characteristics of the most powerful events. The outburst occurred around 2014 and is still ongoing, so it may well be a long-lasting event. The YSO prior to outburst was faint at $r=22$ mag and it is difficult to classify the source from its SED, but it could correspond to a class II YSO (see figure 2 of Hillenbrand et al., 2018).

The discovery of a potential class II YSO outburst over the 3 year baseline of Gaia does not contradict our results. Given an outburst rate of 112 kyr, to observe an outburst over a 3 year baseline would imply that we are observing a sample of approximately 40000 YSOs. This is conceivable if we take into account that there is a large population of YSOs that are too faint for current and past surveys. The number of YSOs increases as we go towards fainter magnitudes (see e.g. Contreras Peña et al., 2017a) and Gaia17bpi itself was in fact unknown prior to its outburst. Hillenbrand et al. (2018) followed up the YSO as a Gaia alert (Hodgkin et al., 2013) was triggered on a source that was found in proximity of a star forming region.

The outburst of Gaia17bpi is remarkable as the analysis of Hillenbrand et al. (2018) shows that the outburst started first at mid-infrared wavelengths, with the optical outburst occurring at least a year later. This is consistent with instabilities the lead to an “outside-in” type outburst, which is predicted by the outburst mechanisms that could explain outbursts during the class II stage (see Section 8).

7 Caveats

7.1 Contamination by other classes of variable stars

The process of planet formation is more likely to be affected by long-lasting rather than short-term accretion events. In this sense we needed to be careful not to include objects that go through sudden episodes of large changes in the accretion rate, but with a short duration. For example, V8 is the known short-term eruptive YSO ASASSN13-db (Sicilia-Aguilar et al., 2017). When inspecting its long-baseline light curve the object did not resemble a long-lasting outburst, but instead appeared to be irregular, i.e. going from bright to faint states repetitively across time. As a result we are confident that similar objects that suffer short-term accretion events were not included in our final sample of long-lasting YSO outbursts.

The classification of high-amplitude variable stars as showing long-lasting outbursts, based only on the photometric properties of the sample, also has the risk of selecting objects with light curves that can mimic an accretion-related outburst, especially if the light curve sampling is sparse (as discussed in Section 5.2). For example is not hard to imagine that sparse sampling of objects brightening from high-amplitude, long-lasting eclipses by inhomogeneities at long distances in the accretion disc (RW Aur, AA Tau) or a precessing disc in a binary system (e.g. KH 15D, Aronow et al., 2018) could resemble a YSO outburst.

Contemporaneous multi-colour monitoring and/or spectroscopic follow-up are usually used to determine that the dramatic change in the brightness of the system is due to large changes in the accretion rate. The sample of 6 long-lasting YSO outbursts that were used to determine the outburst rate in Section 6 contains 3 known YSO outbursts that have been confirmed as such via spectroscopic follow up (V2, V10, V36). Hence in this work we have identified three previously unknown outbursting YSOs that are included in our calculation. Of these, V51 has literature data confirming that this object is an eruptive YSO (see Appendix A.3). Recent spectroscopic monitoring shows that V4 has the spectroscopic characteristics of objects at the high accretion rates of long-lasting YSO outbursts (Contreras Peña et al. 2019, in prep.). V9 is the only object lacking such an spectroscopic confirmation.

However, as discussed in Section 5.2 we tried to avoid including objects where identification as long-lasting YSO outburst was uncertain due to bad sampling and/or lack of additional data from the literature. As a result we are likely excluding objects that could mimic strong YSO outbursts.

7.2 Non-YSO contamination

The sources of contamination in the classification of YSOs arises from red extragalactic sources, such as obscured active galactic nuclei (AGNs), and Galactic sources, mainly reddened giant stars or dust-enshrouded asymptotic giant branch (AGB) stars (see e.g Robitaille et al., 2008; Povich et al., 2013). However, we do not expect the contamination in our sample to be high, as most of our objects arise from samples where contaminants had already been dealt with. For example, Marton et al. (2016) estimate a less than 1$\%$ contamination in their sample of YSOs, whilst this number is not larger than 2$\%$ in the samples of Megeath et al. (2012) and Gutermuth et al. (2009).

We could expect a higher percentage of contamination from the sample of objects that are classified using ALLWISE photometry. Koenig & Leisawitz (2014) establish that a contamination of $7\%$ is expected in objects classified as class II YSOs from ALLWISE photometry. If we assume that objects classified from ALLWISE or 2MASS photometry alone (6419 objects) suffer this high contamination, and that objects selected from other methods suffer a 1% contamination, then we expect $\simeq 492$ or 3.5% of objects are contaminants in our sample.

7.3 Mis-classification

The likely evolutionary stage for the majority of our sample was determined using the Koenig & Leisawitz (2014) criteria. Inspection of their figure 5 shows that the area in the colour-colour plane that defines the class II YSOs could be contaminated by objects at earlier evolutionary stages, i.e. flat-spectrum and class I sources. Given the much higher outburst rate of lass I sources compared with class II (Section 6.1.1), even a small percentage of misclassification would lead to the detection of class I YSOs in our outburst sample, an effect that we observed in our results (see Section 5.4).

To estimate this contamination, we obtained ALLWISE photometry for a sample of YSOs from Megeath et al. (2012) that also have measured values of the spectral index, $\alpha$. Then we applied the Koenig & Leisawitz (2014) criteria to determine the likely evolutionary stage of these YSOs. Finally, we compared the latter with the stage determined from $\alpha$. We found that 7% of objects classified as class II YSOs from the Koenig & Leisawitz (2014) criteria are classified as being at earlier evolutionary stages from the value of the spectral index. As we will show in the following, this estimate is consistent with the number of class I YSO outburst that contaminate our sample in Section 5.4

7.4 Effects of contamination and mis-classification on the outburst rate

We can now examine the effect of all of the uncertainties discussed in Sections 7.1 to 7.3 on our outburst rate. If we first only consider the effect of the lack of spectroscopic confirmation as an accretion outburst for V9 (Section 7.1), then our outburst sample reduces from 6 to 5 objects. This implies an increase of the outburst rate, from 112 to 130 kyr.

Considering possible contamination of non-YSOs, our sample reduces to 13589 objects (see Section 7.2). If we further consider a 7% contamination from class I YSOs (951 objects, see Section 7.3) then we are left with 12638 class II YSOs. If we use a sample of 6 YSO outbursts, the outburst rate in class II YSOs reduces from 112 to 100 kyr. If we omittV9, the outburst rate increases again to 117 kyr.

Thus considering all of the possible contamination and exclusion effects, the outburst rate changes to between 100 and 130 kyr. This is a change that is within the 68% confidence intervals we estimated in Section 6.

We can also use the contamination of class I YSOs to determine another approximate estimate of the outburst rate in the class I stage. In Section 5.4 we found that 3 YSO outbursts are more likely to be class I YSOs, if we consider a total sample of 1043 class I YSOs implies an outburst rate of $13^{+15}_{-6}$ kyr (68 percent confidence, assuming the $1/\tau$ prior). This number agrees remarkably well with that estimated from the Scholz et al. (2013) YSO sample.

7.5 Selection criteria

The value for the recurrence timescale derived from our sample of class II YSO outbursts is subject to some caveats from the initial selection. As we are interested in a statistical characterisation of variability, we should not use variability as a selection criterion. However, is clear that some of the objects in our initial catalogue were initially identified because of their variability, primarily the FU* and or* classes. However given the need of IR data to identify class II objects, we believe that they would have also have been identified as YSOs in these surveys, regardless of their variability. Hence we left them in the sample.

8 Outburst mechanism

The results from our analysis show for the first time that the suggestion that the frequency of the outbursts is lower in the class II stage than during the class I stage is correct. This has implications for the physical mechanism likely to be driving the outbursts during class II phase.

As noted above GI+MRI instabilities become less efficient in disc dominated systems (Bae et al., 2014). Due to GI, and provided that the cooling time is shorter than the dynamical timescale, disc pressure is not able to support the collapsing gas against its own gravity, leading to disc fragmentation (see e.g., Vorobyov & Basu, 2005; Machida et al., 2011). The fragments are later transported to the inner disc via gravitational interactions with the spiral arms (Machida et al., 2011), and are likely to trigger mass accretion bursts (Vorobyov & Basu, 2010). This process is more efficient during the embedded phase. Mass accretion outbursts could occur during the class II phase, but it is not clear whether these fragments survive beyond the embedded stage (see Audard et al., 2014). Zhu et al. (2010) find that for the lower infall rates of class II objects, layered MRI turbulence might be able to accumulate mass and trigger the MRI outbursts. However, these bursts are likely of lower amplitude and of shorter duration (Audard et al., 2014).

Lodato & Clarke (2004) propose that the interaction of a massive planet with the disc can trigger thermal instabilities in the outer disc. In this scenario, the migration of the planet opens up a gap in the disc, and the inner disc is emptied out. Due to the tidal effect induced by the planet, material will pile up at larger radii. When the density reaches a critical value, the thermal instability is triggered. A final possibility is that the gravitational force of a companion star may perturb the disc and thus enhance accretion (Bonnell & Bastien, 1992). This idea is supported by the fact that a number of known eruptive variables are members of binary systems (see e.g., Fedele et al., 2007; Reipurth et al., 2010; Caratti o Garatti et al., 2015). In some of these, both stars in the system are found to have characteristics of strong accretors (e.g. systems AR6A+B, RNO1B/C, Aspin & Reipurth, 2003).

Therefore, the planet-induced thermal instability and binary models are most likely to explain the existence of large accretion events during the class II stage. However, we note that the fact that more than one mechanism can explain the outbursts has implications for our interpretation of the outburst rate. As noted by Scholz et al. (2013), if outbursts are triggered by two or more different mechanisms not all of which might apply to every YSO (such as outbursts due to encounters in multiple systems) would imply that the frequency of outbursts is highly variable amongst young stars, and it would not be valid to extrapolate the outburst rate for a single YSO from our estimated class II recurrence timescale.

9 Summary

In this work we compared the all-sky digitised photographic plate surveys provided by SuperCOSMOS with the latest data release from Gaia (DR2) for a large sample of YSOs. The mean baseline of 55 years implies that we monitored $\simeq 800,000$ YSO years.

We retrieved 139 high-amplitude variables with $\Delta R\geq 1$ mag. We believe we are complete for amplitudes larger than 2 magnitudes. In Section 5 we determined that most of the variable stars are found with amplitudes between 1 and 3 mag, where variable extinction or extreme changes due to hot-spot variability are the likely mechanisms driving the variability in these YSOs. For amplitudes larger than 3 mag we also find objects with apparent long-lasting fading events, likely related to extinction, and objects with irregular variability. At these larger amplitudes, the latter are more likely explained by short-term accretion events that appear as irregular given the long baseline of the light curves.

YSOs with this extreme variability, but not related to variable accretion, are interesting in their own right. For example, long lasting fading events could be explained by extinction due to inhomogeneities at large distances in the disc, which can relate to planet-formation processes.

Objects where variability could be due to changes in the accretion rate, but with outburst duration that is less than 10 years, were not included in the calculation of the recurrence timescale of outbursts in the class II stage. This is because short duration outbursts are less likely to have an impact on the formation and evolution of protoplanets Hubbard (see e.g. 2017b).

We classified 9 YSOs as showing long-term accretion driven outbursts, with 3 of them representing new discoveries (Section 5.3). We discussed the possible contamination by class I YSOs and we determine that 3 objects in this group are more likely at a younger evolutionary stage. Therefore, 6 YSOs are left in our final sample of class II YSO outbursts (Section 5.4).

In Section 6 we determined, for the first time, that long-term accretion driven outbursts occur during the class II stage of young stellar evolution. The frequency of these events was found at $\tau$=$112^{+68}_{-38}$ kyr (68 percent confidence). Our results imply that planet formation models must take into account the likely effects that these episodic accretion events can have on the formation and evolution of planetary systems.

Through the analysis of the YSO sample of Scholz et al. (2013) (see Section 6.1.1) as well as possible contamination of class I YSOs in our sample (Section 7.4), we have estimated the recurrence timescale of accretion driven outbursts at this earlier evolutionary stage. We found that episodic accretion events are $\simeq$10 times more frequent than during the class II stage, in agreement with theoretical expectations.

Acknowledgements

This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). This research has made use of the NASA/ IPAC Infrared Science Archive, which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration. Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. The contributions of C.C.P. and T.N. were funded by a Leverhulme Trust Research Project Grant and of S.M. through a Science and Technology Facilities Council (STFC) studentship. We are grateful to N. Hambly for his help in understanding the SSS database.

Appendix A Individual notes

In the following section we will describe the three YSOs that are previously unknown as eruptive variables. The individual subsections include a description of the photometry and any information that we were able to gather from the literature and that allowed us to include these objects as YSOs displaying long-term outbursts.

A.1 V4

V4 (2MASS J02335340$+$6156501) was originally classified as a YSO by Panwar et al. (2014) based on its near-infrared excess from 2MASS data. There is no additional literature information for this source. The near-infrared $JHK_{\rm s}$ colours from 2MASS indicate that this is likely a class II YSO, hence its inclusion in our sample.

The classification as a high-amplitude variable for this source results from its Gaia DR2 magnitude being brighter than all of the SSS bands. The object is not visible in the images of the 2 epochs of $R$-band from SSS. Although the SSS catalogues yield a detection with $R=19.3$ this in fact corresponds to a star located 2″south west from V4. The variable star becomes visible in the $B$ and $I$ SSS observations and continues to get brighter as confirmed from IPHAS, Pan-STARRS and PTF observations (see Figure 9).

A.2 V9

V9 (2MASS J08410676-4052174) is an H$\alpha$ emission line star (Pettersson & Reipurth, 1994), selected from the Marton et al. (2016) catalogue and classified based on its WISE colours. This object was selected beacuse its Gaia magnitude appeared brighter than expected compared with all of the 4 SSS bands. Photometry from Vizier confirms this is as a long-term outburst as the object appears to be faint in all SSS epochs with an eruption likely occurring between the second R epoch from SSS and DENIS observations. The object now appears to be slowly declining, as shown by APASS, VPHAS+ observations and the light curve from the Bochum Galactic disc survey (Hackstein et al., 2015). Figure 10 shows the light curve of V9 as well as the comparison between SSS and VPHAS$+$ images of the source.

A.3 V51

V51 (WRAY 15$-$488) is an H$\alpha$ emission line star (Wray, 1966), which we selected from the Marton et al. (2016) catalogue and classified based on its WISE colours. The object is part of the 139 high-amplitude variable star catalogue as the $R$-band epochs of SSS and Gaia DR2 are brighter than the earlier B and I SSS observations. Photometry collected from Vizier confirms its nature as a long-term outburst as the object becomes $\simeq$2.5 magnitudes brighter between 1980 (I-band SSS) and 1984 (first epoch in $R$ from SSS) and has remained at approximately the same level ever since, as confirmed by data from the All Sky Automated Survey (ASAS, Pojmanski, 2002) and SkyMapper (see Figure 11). Evidence that the object was at a quiescent state when it was classified as an H$\alpha$ emission line star by Wray (1966) arises from the catalogue of Wackerling (1970) which contains the data from Wray (1966) and where V51 is listed with V$=13.1$, approximately 3 magnitudes fainter than the current brightness of the object.

This object has spectroscopic observations that confirm its nature as an eruptive YSO. Carmona et al. (2010), based on high resolution FEROS optical spectra observed in 2004, find that the object shows an H$\alpha$ P Cygni profile, with an spectrum consistent with that of an F-type giant star. The latter classification is characteristic of eruptive YSOs during outburst (Connelley & Reipurth, 2018). Previous spectra, also observed during outburst, from Suárez et al. (2006) and Gregorio-Hetem et al. (1992) are also consistent with an F-type spectral type. V51 also shows broad Na I D $5890$Å absorption (Suárez et al., 2006), another characteristic of YSO outbursts, as well as Li I absorption (Gregorio-Hetem et al., 1992), evidence of the youth of the system. Remarkably during these spectra (which are observed at earlier epochs of the outburst than that of Carmona et al., 2010) the H$\alpha$ line is seen purely in emission.

Appendix B Light curves

Appendix C Vizier Tables

This appendix contains the list of surveys used to add photometry to the light curves of our high-amplitude variable YSOs.