The Outburst of the Young Star Gaia19bey

The Gaia $G$-band photometric data were obtained from the Gaia alert data base and form the basis for the light curve in Figure 2. ³³3http://gsaweb.ast.cam.ac.uk/alerts From the start of the measurements on 2014 November 30 to 2016 August 4, Gaia did not return significant detections of our object, and from the faintest reported measurements we treat these non-detections as upper limits at $G$ = 20.5.

Many of the bright data points of the optical light curve of Gaia19bey are archival data from the ATLAS project described by Tonry et al. (2018). ATLAS usually takes more than one image of any given region of the sky in each suitable night to follow fast-moving asteroids. For the light curve points in Figure 2, we have median-combined the individual measurements for each night when the photometric zero points and sky brightness were stable. From the small number of epochs when closely coinciding Gaia and ATLAS or Pan-STARRS data points were available, the magnitude offset in each filter was determined and the other light curve points in that filter were shifted accordingly, so that the full light curve is effectively in Gaia-$G$ magnitudes. We did not have the data to also match second order (color) effects. The ATLAS data points have fairly large errors in the fainter phases of the light curve and we therefore include only the brighter (effectively $G<$ 18.0) ATLAS “Orange” data points in Figure 2. The current outburst is the largest observed in the past decade, but smaller variations are seen in the older Pan-STARRS data.

The Panoramic Survey Telescope and Rapid Response System (Kaiser et al., 2010, Pan-STARRS) is a wide-field optical imaging sky survey project. The PS1 3$\pi$ survey (Chambers et al., 2016) patrolled the entire sky north of Declination.$=-30^{\circ}$ at multiple epochs between June 2009 and April 2014 and in 5 filters ($grizy_{\rm P1}$) (Tonry et al., 2012). The 5$\sigma$ single-epoch depths are down to 22.0, 21.8, 21.5, 20.9, and 19.7 mag, respectively and photometric calibration uncertainty of the survey is $\sim 0.01$ mag. We extracted the photometric data from the PS1 third full reduction data set through the PS1 internal database system. Additionally, we ingested some post-PS1 3$\pi$ survey data covering our target, hence the time baseline extends up to May 2017. The additional data were also calibrated by a reference catalog generated from the previous reductions data (Magnier et al., 2016). For the light curve in Figure 2, we have shifted the magnitudes in any of the PS1 filters to the Gaia $G$ photometry, when closely coincident epochs were available. This shifting is just the first order of adjusting the different filter data and does not include color correction because coincident color information was not available. These adjusted Pan-STARRS data give an approximate light curve covering the years from 2009 to 2017 effectively in Gaia $G$ magnitudes, when neither Gaia nor ATLAS data are available.

To complement the ATLAS data, we have downloaded ⁴⁴4https://irsa.ipac.caltech.edu archival $r$-band photometry of Gaia19bey from the Zwicky Transient Facility (ZTF) (Bellm et al., 2019) archive (Masci et al., 2019). These data reach about one magnitude fainter than the ATLAS data and cover the outburst maximum very well. Similar to the procedure for ATLAS and Pan-STARRS, we shifted the ZTF magnitudes by a fixed amount to match the Gaia $G$ photometry.

We have included archival data from the WISE (Wright et al., 2010) and NEOWISE (Mainzer et al., 2014) missions downloaded from the NASA/IPAC Infrared Science Archive ⁴ in the light curve (Middle panel in Figure 2). Most of the WISE data in the W1 and W2 bands on Gaia19bey are saturated, but the fainter W1 data points prior to the present outburst are unsaturated and confirm the brightness minimum that is otherwise only indicated by just two Pan-STARRS measurements. The W1 band data points above saturation still contain some useful information in that they confirm the basic shape of the light curve during the outburst.

Figure 3 shows the spectral energy distribution (SED) of Gaia19bey based on ground-based and space-based survey data obtained from VizieR. The infrared and far-infrared data are from the compilation of catalog data by Abrahamyan et al. (2015, and references therein) and include data from WISE (Wright et al., 2010), MSX (Egan & Price, 1996), AKARI (Murakami et al., 2007), and IRAS (Neugebauer et al., 1984). The 2MASS data are from the point source catalog (Skrutskie et al., 2006). The photographic data points (DSS) are from the guide-star catalog of Lasker et al. (2008). In the optical and near-infrared we have compiled all available individual photometric values from Pan-STARRS PS1, the 2MASS and UKIDSS surveys, and NEOWISE. All these data were taken without any knowledge of the light variations of Gaia19bey and do not represent the full range of brightness variations that Gaia19bey may have experienced in the past, but they give an impression of the typical variations of the short end of the SED with its variability. Since Gaia has covered this object from invisibility, i.e. G $>$ 20.5, early in its mission to the current maximum, those data points represent a good lower limit of the total variability this object has experienced in the past. Also, we note that Gaia19bey must have been in a previous, but unrecognized outburst when the 2MASS photometry was obtained on 1998 November 3 since the 2MASS catalog photometry is brighter than the maximum of infrared photometry measured during the present outburst. For wavelength longer than 4.6 $\mu$m, we do not have multi-epoch data from individual missions, but some of the observed scatter may still be due to variability.

Deep infrared images were obtained in the $J$, $H$, $K$ bands and the narrow-band S(1) line filter with the Wide-Field Camera (WFCAM) described by Casali et al. (2007) on UKIRT on 2020 June 20. The exposure time of individual frames was 1.0 s in J, H, and K with 5 coadds per position and a total of 72 dithered frames taken, resulting in a total integration time of 360s. We have verified that the detector did not saturate in the K band on this relatively bright object. In the S(1) filter, the individual frame exposure time was 20s, and the total integration time in the 72-point dither pattern was 1440 s. The continuum subtracted image in Figure 4 is the S(1) line image with 1.2 times the K-band in the 1 s individual exposures image subtracted for optimal subtraction of field star. The stars in the group of highly reddened stars near Gaia19bey had higher $K$-band flux than the average field stars and are therefore oversubtracted in Figure 4. It is noteworthy that the UKIDSS image of Gaia19bey had an exposure time of 5 sec and was in the non-linear regime of the detector, despite the fainter magnitude at the time. ⁴⁴4http://casu.ast.cam.ac.uk/surveys-projects/wfcam/technical/linearity To account for this possible systematic effect, we estimate a larger error for this data point, as indicated in Table 2.

We obtained an optical spectrum of Gaia19bey with the “SuperNova Integral Field Spectrograph” (SNIFS) described in Lantz et al. (2004) at the UH 88” telescope on 2019 December 9 (MJD 58826). SNIFS is a dual-channel instrument, but we only show the red part of the spectrum in Figure 5 because the blue channel spectrum was too faint to show any features. The spectral resolution is approximately R $\approx$ 1300. The spectrum has not been corrected for telluric absorption, and the absorption feature of O₂ at $\approx$ 760 nm (Rudolf et al., 2016) is clearly visible, but does not affect the interpretation of this spectrum. The line identifications are based on the study of Herbig AeBe stars by Hamann & Persson (1992).

We obtained infrared spectroscopy covering parts of the wavelength range 1.0 - 2.5 $\mu$m during two observing nights at the Keck 1 telescope, and using director’s discretionary time on the NASA Infrared Telescope Facility (IRTF). All infrared spectra are shown in Figure 6.

The light curve (Figure 2) shows a broad maximum from MJD 57500, when Gaia began recording significant detections, to MJD 59000, when the light curve appears to have faded back to within half a magnitude of the pre-outburst levels below $G$ = 20.5. The light curve showed a brief spike up to the maximum of $G$ = 17.09 on 2019 June 1 (MJD 58605) that had a total duration of less than one year and led to the Gaia alert at a magnitude of $G$ = 17.49. While the one-year brightness spike is of the same duration as typical EXor outbursts, the overall outburst duration above $G$ = 20.5 (the Gaia detection limit) of $\approx$ 4 yr is longer than those found in classical EXors, but much shorter than that of FUors. The longer duration is in line with many of the higher luminosity eruptive events discovered and studied at infrared wavelengths by Lorenzetti et al. (2012), who called them “Newest EXors” and that Contreras Peña et al. (2017b) used to tentatively define the new MNOR class.

Continguous coverage of the light curve started around MJD 55000 when Pan-STARRS 1 (PS1) began the PS1 3$\pi$ sky survey. Until the Gaia mission started reporting significant detections of Gaia19bey, these and a few unsaturated $W1$ data points from the WISE and NEOWISE missions are the only available observations of this object. The earliest Pan-STARRS measurements, combined from all the filters, between MJD 55000 and 55300 ($\approx$2009) showed a brightness near $G$ = 23. Around MJD 55500, the brightness increased by 2 magnitudes to $\approx$21 mag, and this rise was partly confirmed by WISE. The light curve between MJD 56500 and the start of the Gaia data at MJD 57500 is only poorly covered since PS1 had completed the 3$\pi$ sky survey and was now concentrating on asteroid searches and only covered the position of Gaia19bey on two occasions. This flux minimum was, however, observed by the NEOWISE-R mission with four additional data points. These data points were, in fact, the only NEOWISE $W1$ data points that were not saturated and therefore are a reliable confirmation of the minimum prior to the outburst when the brightness dropped again to $\approx$23 mag. In the brighter phase at approximately 21 magnitude, the object showed additional variations of about one magnitude and timescale of $\approx$ 500 days. The outburst amplitude depends on how the quiescent brightness is defined. Prior to the present outburst, our Pan-STARRS photometry suggests a minimum of $G$ $\approx$ 22.5 around MJD 57000, even though that time interval is poorly covered with only two Pan-STARRS observations. Relative to this minimum, the outburst amplitude is 5.5 mag in the $G$ band.

We have no spectroscopic data from that time period and from optical photometry alone cannot distinguish variations due to accretion variations from variations causes by variable extinction in the light path. Given that the morphology of the forbidden line emission (Figure 7) strongly suggests a disk seen nearly edge-on, variability due to extinction variations is very likely to occur and may be the dominant reason for those brightness variations. Between the start of the PS1 observations and the first Gaia detections, the data indicate initially a faint phase, a brighter phase of 1000 days, and again fainter magnitudes. This is, of course, insufficient to determine that those changes are periodic. However, if we assume some semi-periodicity, it is likely that the latest, much brighter outburst coincided with a period of lower extinction, even though it was not caused solely by lower extinction. In this working assumption of a combination of extinction variations and a stronger variation of accretion, we assume that the baseline for the accretion outburst is the brighter level of the extinction variations in the light curve at $\approx$ 21 mag. Calculating the outburst relative to this mean gives an outburst amplitude of $\approx$ 4 mag, well in the range of other eruptive events in YSOs.

Some evidence for a prior major outburst comes from the 2MASS catalog (Skrutskie et al., 2006) that shows $J$, $H$, and $K_{s}$ magnitudes (Table 1) brighter than those seen in the present outburst, and from the fact that the UKIDSS survey recorded much fainter $J$, $H$, and $K$ values that we assume represent the quiescent state of Gaia19bey. The UKIDSS $K$-band data point is brighter than the saturation limit reported by Lucas et al. (2008) and is therefore not reliable.

While the early parts of the Gaia light curve prior to MJD 58200 appear to be smooth, starting shortly before MJD 58300, through the maximum, and in the declining phase we see several brief spikes in brightness with typical duration of $\approx$ 10 to 50 days and $G$ amplitudes of order 0.5 mag, the earlier of which are very well confirmed by the high-cadence ZTF observations. The durations and amplitudes of these brief spikes are similar to those reported by Semkov (2006); Semkov & Peneva (2012) in the declining phases of the recent outbursts of V1647 Ori. Since these minor spikes in brightness occur in the phase where our spectra indicate on-going accretion activity, they are most likely associated with fluctuations in the magnetospheric accretion process.

The SED is a typical “flat spectrum” distribution (Greene et al., 1994). The flat part of the SED extends from the WISE 2 band (4.6 $\mu$m) out to 160 $\mu$m from AKARI, characterizing Gaia19bey as a YSO in transition from a deeply embedded infrared object (Class I) to a moderately embedded Class II. At optical and near-infrared wavelengths, the SED rises steeply, indicating substantial dust obscuration in the line of sight. In Section 4.5, we will determine the extinction to the immediate vicinity of the star to be A_V = 12 mag. Figure 3 also shows that the observed outburst amplitude is strongly wavelength dependent, being much larger at optical wavelengths than in the infrared.

To determine the luminsity of Gaia19bey during the outburst, we integrated over the available F$\nu$ values and selected the flux at the outburst maximum for the Gaia, 2MASS, and WISE wavelengths, while using the available archival data for longer FIR wavelengths, even though these measurement were not taken during an outburst. Assuming the minimal possible distance of 1.4 kpc, this method gives a lower limit to the luminosity during outburst of $\approx$182 $L_{\odot}$, placing Gaia19bey in the luminosity range of typical Herbig Ae/Be stars (Hillenbrand et al., 1992). Considering the uncertainty in the distance, an order of magnitude higher luminosity is within the range of possibilities.

Herbig (2008) has found that eruptive young stars of the EXor type occupy a characteristic locus in the $J-H/H-K$ color diagram (their Fig. 10), with redder colors than what main-sequence stars obscured by interstellar extinction would have. Lorenzetti et al. (2012) has shown a similar color-color diagram including all known EXors at the time. While classical EXors, including the prototypical EX Lupi, are only moderately reddened, many of the recently discovered deeply embedded stars with EXor-like spectra are much redder and occupy the upper right area in Figure 8. All these deeply embedded EXors for which $J$, $H$, and $K$ photometry in different phases of their outburst is available change their position in the color-color diagram in the sense of being bluer when bright. The variation path in the color-color diagram is flatter than the interstellar extinction vector. If the variations were to be interpreted as solely the effect of extinction variations, this would indicate a less wavelength dependent extinction law, consistent with larger grains than are prevalent in the interstellar medium. However, we will show in Sections 4.4 and 4.5 that the outburst of Gaia19bey is not solely caused by extinction variations. The newly discovered Gaia19bey shares the color-color diagram locus and its variation with other deeply embedded outburst stars studied by Lorenzetti et al. (2012) and is closely similar to ESO-H$\alpha$ 99, studied by Hodapp et al. (2019) While these “new EXors” (Lorenzetti et al., 2012) share the basic characteristic of having an emission line spectrum with classical EXors, their outbursts last longer, they are of higher luminosity, and their locus in the color-color diagram is markedly redder.

In the optical light curve in Figure 2 we have indicated the times when the spectra were taken with the same color coding as used in Figures 5 and 6 for the spectra. The first infrared spectrum (shown in red in Figure 5), taken on 2019 June 29 (MJD 58664), was taken just one month after the maximum brightness of Gaia19bey during its present outburst at a $G$ magnitude of 18.0. The second spectrum (shown in green in Figure 6) on 2019 October 23 (MJD 58779) shows an intermediate stage of the brightness decline ($G$ = 18.5) after the maximum when the brightness was rising again to a minor temporary spike. The emission lines were fainter at that time than in the first spectrum both in absolute terms and relative to the continuum. The change in line strength (Table 3) was more pronounced for the Na I and Mg I doublet lines than for the hydrogen lines. The last spectrum shown in blue in Figure 6, taken with SPEX on 2020 March 6 (MJD 58915) when the Gaia magnitude of 20.0 mag was about 3 magnitudes below the maximum at optical wavelengths, has changed fundamentally, showing an almost pure continuum spectrum with emission lines being an order of magnitude fainter than on 2019 October 23. The most pronounced emission lines in this spectrum are the forbidden shock-excited lines presumably emerging in the outflow at some distance from the star, and apparently not directly affected by the variability.

We have measured the equivalent width of several emission lines (Table 3) and plot them against the continuum flux at 2190 nm in Figure 9. This particular wavelength was chosen to be free of detectable emission lines and therefore representative of the continuum component. In this plot, extinction variation due to dust absorption would not change the equivalent width of emission lines shining through that extinction, leading to horizontal lines in this plot. In contrast to this, the equivalent widths strongly decline with declining continuum brightness of the object, indicating that the strong photometric outburst from MJD 57500 to 59000 was caused by a substantial increase in magnetospheric accretion, leading to an EXor-type emission line spectrum.

The optical spectrum taken with SNIFS at the UH 88” telescope (Figure 5) was taken on 2019 December 9 (MJD 58826) at around a $G$ magnitude of 19, between the last two infrared spectra both in time and in object brightness. It clearly shows the forbidden line of [O I] at 630.0 nm and weakly indicated the 636.3 nm  line The H$\alpha$ emission is strong and the blended [S II] 671.6/673.1 nm are weakly indicated. In combination these lines are characteristic of Herbig-Haro shocks. In addition, the Ca II triplet is very strong in emission as is often found in T Tauri stars as discussed, e.g., by Muzerolle et al. (1998). The ratio of these 3 lines is expected to be 1:9:5 proportional to their (gf) values in the optically thin case, as discussed by Azevedo et al. (2006) for the case of classical T Tauri stars. In Gaia19bey, the line ratio is much closer to unity, indicating that the conditions in the line emitting region are close to optically thick. This is similar to the pattern of line strength found by Hamann & Persson (1992) and Hamann & Persson (1992) for T Tauri and Herbig AeBe stars in that the Ca II line at 849.8 nm has the highest line flux. The only absorption line detected in the optical spectrum is the Li I line at 670.8 nm (the insert in Figure 5). The lithium line is an indicator of young age in stars since in convective atmospheres, it is being transported sufficiently deep into the stellar interior to be destroyed by nuclear reactions, as described in the review by Pinsonneault (1997) and the references therein.

We have used two infrared emission lines of hydrogen to measure the line luminosity and deduce an estimate for the extinction and the accretion luminosity. For Herbig AeBe stars, less obscured analogs to Gaia19bey in terms of luminosity, Fairlamb et al. (2017) have recently measured both UV continuum emission, a direct indicator of accretion onto the stellar surface, and the line luminosities of many optical and near-infrared emission lines, and have given empirical relationship between these line luminosities and the total accretion luminosity. Gaia19bey is in the luminosity range of Herbig AeBe stars, and we used these L_acc vs. L_line relationships to derive L_acc for two strong hydrogen emission lines (Br$\gamma$ and Br 11-4) in our Keck OSIRIS spectra from 2019 June 29 and 2019 October 23. As a side note, Fairlamb et al. (2017) do not give the calibration for Br 10-4, so that we could not use this strong line. Used naively, their calibration procedure ignores the effect of extinction. We can get an estimate of the extinction, and correct the accretion luminosity for extinction by assuming that the two emission lines, when properly corrected for extinction, should give the same accretion luminosity. We used the continuum extinction law measured by Boogert et al. (2011) in lines of sight through dense molecular cores as the best available estimate of the extinction through the edge-on protostellar disk of Gaia19bey. For the 2019 June 29 data when the emission lines were strong, we solved for the extinction, expressed as A_K when the accretion luminosity determined by each line becomes the same and obtained an A_K = 1.59 $\pm$0.6 and $log(L_{acc})$ = 1.39 $\pm$0.35. For the 2019 October 23 data, when the emission lines were substantially weaker and the determination of the extinction much more difficult, we assumed the same extinction as on 2019 June 29: A_K = 1.59, and took the average of the accretion luminosities measured in Br$\gamma$ and Br 11-4 as the best estimate listed in Table 4.

The uncertainties in these measurements are large, partly due to the faintness of the lines, the less-than-perfect subtraction of the hydrogen absorption lines in the spectrum of the telluric absorption standard star, and due to uncertainty over the extinction law. Despite all this, we conclude that shortly after the maximum of the light curve, Gaia19bey has a similar accretion luminosity and mass accretion rate than the higher end of the distribution of mass accretion rates in Herbig AeBe stars, whose accretion has been studied by Mendigut´ıa et al. (2011). The accretion rate is much lower than that of FU Orionis outbursts where dM/dt = $1.\times 10^{-4}$ M_⊙yr^-1 is possible (Audard et al., 2014). For an approximate conversion of the A_K extinction into the commonly used A_V, we use A_V = 7.5 * A_K, corresponding to R_V = 5 (Cardelli, 1985). Our extinction measurement of A_K = 1.59 corresponds to A_V = 12 $\pm$ 5.

The adaptive optics OSIRIS data cubes were used to extract continuum subtracted images of the emission in the [Fe II] line at 1644 nm, and the H₂ 1–0 S(1) line at 2122 nm in Figure 7. For the S(1) line, we have two images in the 50 mas and 100 mas per spaxel spatial scales. The [Fe II] emission is extended over $\approx$ 1 ″and shows bifurcating band of low emission across the position of the unresolved continuum source. This morphology strongly suggests that we look at a region of extended [Fe II] emission in an edge-on disk system with the disk plane oriented at P.A. $\approx$ 125°. Similar morphologies were found to be common in young outflow sources when imaged at high spatial resolution with HST by Padgett et al. (1999). Morphological features related to outflow activity are also seen in most other deeply embedded EXors and FUors. The S(1) line shows very different features. At the 50 mas scale, the S(1) emission is only slightly extended to the NW, in the same direction as the dark lane in [Fe II]. The wider view in the 100 mas scale, the right panel in Figure 7, shows what looks like the inner walls of an outflow cavity. The western side of this outflow cavity is not included in our field of view, but in combination with the absorption lane seen in the [Fe II] image (left panel in Figure 7) it is consistent with the axis of the outflow cavity being oriented at P.A. 215°.

We have recorded both the [Fe II] lines at 1256 and 1644 nm, from which, in principle, a value for the extinction can be calculated. However, the lines are very faint in our spectra, and for the OSIRIS spectra, the close proximity of the [Fe II] 1644 nm with the H Br 12 emission line, and consequently the excess noise associated with the correction of the telluric standard absorption line, is problematic. We have therefore decided to forego this exercise. In any case, as Figure 7 shows, the [Fe II] 1644 nm emission arises in an extended region around the star and a dark lane bisecting this extended emission indicates absorption in an edge-on disk, as discussed in Section 4.5.

The presence of an edge-on disk is also consistent with the flat-spectrum characteristics of its SED, which indicates the presence of a substantial amount of relatively cold dust around the object. In young stars, such cold dust is almost inevitably in the form of a disk.

We have imaged the wider field around Gaia19bey with the UKIRT WFCAM in search of more shock-excited outflow features and the continuum-subtracted image (Figure 4) does indeed show a H₂ jet with several knots and shock fronts. While this jet lies in the general direction of the opening of the outflow cavity seen in the adaptive optics image (Figure 7), it does not align well with Gaia19bey, and, more importantly, the morphology of the shock fronts suggest that this jet is moving in the general SW to NE direction, i.e., in the general direction towards Gaia19bey. We cannot determine the source of this jet with certainty, but a plausible candidate source is the WISE object J204008.40+464708.0, 6.0 arcmin SW of the jet, that is detectable out to the WISE band 4, and is associated with faint S(1) line emission closer to it. This object is not included in the SIMBAD data base and no records exist in the literature. The important point is that Gaia19bey is not associated with a large scale H₂ S(1) jet, outside of the features seen in the adaptive optics images very close to the star. As was mentioned in Section 2, there is also no H$\alpha$ emission that could be morphologically associated with Gaia19bey.

Both the optical spectrum (Figure 5) and the first two of the infrared spectra in Figure 6 show the characteristic emission lines seen in very young EXors. The hydrogen emission lines are generally assumed to originate in the accretion funnels in the magnetosphere of the accreting star, where material flows onto the star along magnetic field lines. The CO bandhead emission is also produced in the hot innermost parts of the accretion disk by irradiation from the accretion hotspots near the surface of the star (Calvet et al., 1991). In the infrared, the emission lines are on a strong continuum indicating substantial emission from hot dust probably irradiated and heated by the accretion hotspots.

While not two eruptive YSOs have identical characteristics, it is useful to put Gaia19bey in the context of similar objects. Gaia19bey is a fairly high luminosity, deeply embedded object in a distant star-forming region. An interesting comparison case is V346 Norma, for which Kóspál et al. (2020) have recently published a detailed spectrum and have discussed changes of the spectrum over the past several decades. V346 Nor has undergone a long-lasting outburst from 1980 to 2010, a relatively brief minimum, and for most of the past decade, a rise back close to the previous outburst level. While the timescales of its variations are an order of magnitude longer than for Gaia19bey, the dip in brightness before the rise to the present maximum appears qualitatively similar.

A spectrum of V346 Norma obtained in 1983 by Reipurth (1985) showed only an absorption spectrum. Less than one year later, in 1984 Graham and Frogel 1985 obtained a spectrum showing weak Li absorption, but also a strong H$\alpha$ emission line with a P Cygni profile. Later in the outburst of V346 Nor, in 1993, an infrared spectrum by Reipurth et al. (1997) detected emission of Br$\gamma$, NaI (2206 nm) FeI (2240nm) and CO bandhead emission, similar to what we have observed for Gaia19bey during its outburst. However, the latest infrared spectrum from 2015 (Kóspál et al., 2020) do not detect any of these lines, and the CO bandhead may be absent or weakly in absorption. Their spectrum is dominated by a multitude of [FeII] and H₂ emission lines originating in shocks, possibly indicating the emergence of a new HH object or jet. Otherwise, the spectrum is essentially a continuum, similar to what we observe in Gaia19bey in the declining phase of its present outburst.

Another object with characteristics similar to Gaia19bey is the eruptive variable PV Cep. The infrared spectrum obtained by Caratti o Garatti et al. (2013) shows many emission lines, both permitted lines, in particular of atomic hydrogen, similar to those seen in Gaia19bey and probably originating in magnetospheric accretion. In addition, many forbidden lines of [FeII] and H₂ were detected and indicate outflow activity, similarly, but with stronger forbidden lines than in Gaia19bey. However, Kun et al. (2011) have demonstrated that the brightness variations in PV Cep cannot solely be explained by accretion variations, but that variations in the extinction also play an important role.

The eruptive young object V1647 Ori, sometimes referred to as McNeil’s Nebula (reference) and being the prototype of the “MNOR” class of deeply embedded high luminosity EXor that Contreras Peña et al. (2017b) have defined, shared many characteristics with Gaia19bey:

1) have SEDs of class I or are flat-spectrum sources, 2) show outburst durations $>$1.5 yr but shorter than those usually associated with FUors, 3) show spectroscopic characteristics of eruptive variables, i.e. CO emission or absorption, and 4) usually have 2122 nm H₂ emission. Also, their sample contains several objects of high bolometric luminosity in the range of hundreds to thousands of $L_{\odot}$, the same range as Gaia19bey, considering the substantial uncertainty of its distance. 5) Similar to at least one object in their sample, Gaia19bey changed from an emission spectrum to a continuum during the declining phase of the outburst.

Aspin et al. (2009) derived a mass accretion rate of $4.6\times 10^{-6}M_{\odot}yr^{-1}$ during the outburst of V1647 Ori from the Br$\gamma$ line, and Aspin et al. (2008) obtained $1.0\times 10^{-6}M_{\odot}yr^{-1}$ during quiescence. They also point out the substantial role that extinction variations in the scattering path of the reflection nebula must play in explaining the light curve of this object. While we do not detect a reflection nebula near Gaia19bey in continuum, the prevalence of extinction variations in other YSOs strongly suggests that such variations play a role in explaining brightness variations in phases of the light curve when magnetospheric accretion, evident by emission lines, has already subsided.

The spectral changes observed in Gaia19bey are similar to those observed by Guo et al. (2020) in at least two YSOs of the MNOR type in their sample. Their object VVVv374 changed from an emission line $K$-band spectrum to absorption of Br$\gamma$ over the course of roughly one year, associated with a $\approx$ 0.5 mag increase in brightness during a longer brightening phase. Their VVVv662 went from Br$\gamma$ emission to a pure continuum spectrum after a 1 magnitude drop in $K$-band magnitude.

Changes in the spectrum of Class I protostars were studied by Connelley & Greene (2014) in a sample of non-eruptive objects. Even without major outburst events, changes in emission line ratios were routinely found. However, only one object in their sample, IRAS 03301+3111 showed a strong increase in the near-infrared continuum and a transition from CO in absorption to emission over the course of 3 years.