© ESO 2020

1. Introduction

The young stellar object (YSO) phase represents a very important stage in the life of a star and influences its subsequent evolution. YSOs can be divided into four classes (Class 0, I, II, and III), where Class I to Class III are defined by their spectral index (α) measured from the near- to mid-infrared (NIR to MIR) portion of the spectrum (Lada 1987). Class 0 stars are normally observed only at millimetre (mm) and radio wavelengths and represent the earliest phase when over 50% of the mass is still contained in an envelope surrounding the protostellar core. Although Class I YSOs are deeply embedded, they are nevertheless observable in the NIR. Since they are still strongly accreting and generating powerful outflows, it is possible to study both accretion and ejection processes at this relatively early stage through IR spectroscopy and imaging using state-of-the-art ground-based telescopes.

Young stars have been known to exhibit episodic variability in their accretion and ejection over the course of their evolution (see e.g. Audard et al. 2014, and references therein). It is important to note that in this case, an increase in accretion is usually associated with an increase in luminosity as more material is accreted, producing strong shocks onto the stellar photosphere and additional radiation. Two evident forms that this variability can take are FU Orionis-type outburts (FUors, named after the prototype FU Ori) and EXor outbursts (named after EX Lupi), which were first discovered in the context of optical observations (Herbig 1966, 1977, 1989) and since then have had their definition broadened to include more embedded types of young stars (e.g. Connelley & Greene 2010).

FUors are YSOs that exhibit accretion bursts of several orders (3–4) of magnitude, reaching Ṁ_acc ∼ 10⁻⁴ M_⊙ yr⁻¹ for a relatively short timescale (∼10² years) and might accrete up to ∼30 − 40% of their final mass during these bursts over the course of their formation (Fischer et al. 2019). It is believed that these kinds of dramatic bursts preferentially occur during the early stages of star formation when mass is still falling onto the disc from an envelope (e.g. Vorobyov & Basu 2015), even though the phenomenon was first discovered in pre-main sequence (PMS) stars.

EXor bursts are phenomena that occur over shorter timescales (∼1–2 years) and are less violent (Ṁ_acc increases of 1–2 orders of magnitude typically up to 10⁻⁷ − 10⁻⁶ M_⊙ yr⁻¹) than their FUor counterparts (Audard et al. 2014). The brightness of these objects can increase by a few magnitudes over mere months, according to photometric observations (e.g. Audard et al. 2010). Their frequency is also higher than FUors, with bursts occurring potentially only a few years apart (Herbig 2008). Similar to FUors in the quiescent state, most EXors are optically observable classical T Tauri stars. However, there is evidence that earlier stage protostars also exhibit episodic bursts (see, e.g. Audard et al. 2014, and references therein), and it has been found that, in Class I YSOs, this eruptive variability is at least an order of magnitude more common than in Class II YSOs (Contreras Peña et al. 2017). Certainly the increased use of IR observations has helped to shed light on this phenomenon at earlier phases in stellar evolution.

Here, we investigate the variability of a single object (IRS 54) over 9 years using NIR spectroscopic and photometric data. IRS 54 (YLW52) is located in the Ophiuchus star-forming region at a distance of ∼137 pc (Sullivan et al. 2019). It is a Class I very low-mass star (VLMS) (M_⋆ ∼ 0.1 − 0.2 M_⊙) of estimated spectral type M (Garcia Lopez et al. 2013, hereafter GL13) with a bolometric luminosity of L_bol = 0.78 L_⊙ (van Kempen et al. 2009). Observations of this YSO have revealed an accretion disc and a H₂ molecular jet (Khanzadyan et al. 2004, GL13) typical of protostars at an early evolutionary phase (see, e.g. Lee 2020, and references therein). IRS 54 in fact is one of the lowest luminosity sources where an H₂ outflow has been spatially resolved (GL13). Moreover, it is an ideal candidate for studying variability due to its edge-on disc geometry that allows us to view the red- and blueshifted components of its outflow, and because multi-epoch spectra and imaging are available spanning almost a decade.

2. Observations and data reduction

The Class I protostar IRS 54 was observed over four epochs (2005, 2010, 2013, and 2014) in the NIR, as reported in Table 1. Epochs 2005 and 2013 were obtained with the Very Large Telescope (VLT) at the European Southern Observatory (ESO) Paranal Observatory in Chile using the Infrared Spectrometer and Array Camera (ISAAC, Moorwood et al. 1998). ISAAC employed medium spectral resolution (R ∼ 10 000, see Table 1), a slit width of 0.3″, and a slit length of 120″, with a pixel scale of 146 milli-arcseconds (mas). The K band data in 2005 cover a larger wavelength range than in 2010, because two contiguous spectral segments were acquired in 2005. The seeing values for each individual night are included in Table 1. To correct for the atmospheric response, telluric standard stars were also observed (see Col. 6 in Table 1).

The 2014 data were acquired over two separate nights with the VLT using the Spectrograph for Integral Field Observations in the Near Infrared (SINFONI; Eisenhauer et al. 2003), an integral field unit (IFU). SINFONI observations in the H and K bands had a pixel scale of 100 mas, with a corresponding field of view of 3″  ×  3″, and a spectral resolution of ∼4000. The seeing measurements for each night are included in Table 1. As with the ISAAC data, to correct for atmospheric effects, telluric standard B-type stars were observed (see Table 1).

The data reduction was completed with the GASGANO¹ data file organiser to run the standard SINFONI pipeline recipes. These were used to apply dark and bad pixel masks, flat field correction, optical depth correction, and a wavelength calibration using either OH lines (in the case of the H band) or arc lamps (in the case of the K band, where there were not enough strong OH lines) to the data cubes. A spectrum was then extracted from the telluric cube using IRAF² (Image Reduction and Analysis Facility). The region to extract (about the central region on the cube) was determined using CASA³ (Common Astronomy Software Applications) Viewer. Hydrogen-recombination lines were manually removed from the spectrum of the telluric standard star before this spectrum was used to correct for telluric absorption. The H band reduction of the 22 May 2014 SINFONI data required a further manual sky subtraction because OH line residuals were present in the datacube. This was done by selecting a region of sky in the field of view with little to no emission from the source and subtracting this from the science cube. The resulting spectra were extracted on source.

ISAAC spectroscopic data were reduced in the standard way using IRAF. Wavelength calibration relied on the OH atmospheric lines in each frame. Spatial distortion and curvature caused by the long slit were corrected using the calibration file STARTRACE. An average wavelength accuracy of 2 Å was achieved. As in the case of the SINFONI data, hydrogen-recombination lines were removed from the telluric standard spectra before telluric corrections were applied.

The photometric data obtained with ISAAC (epochs 2005 and 2013) were reduced with IRAF. Flat-fielding of the raw data, sky subtraction, and cosmic ray corrections were all performed. Approximately five nearby stars in the field of view were used to flux calibrate the final science images using their known 2MASS catalogue values. However, the J band flux calibration was completed using the J band acquisition image as photometric science images were not available.

When calculating the line velocities, the spectra were corrected to the parent cloud velocity of ∼ 3.5 km s⁻¹ (Wouterloot et al. 2005; André et al. 2007).

The 2010 VLT/SINFONI archival data were taken from GL13.

3. Results

3.1. Morphology

Detected features in IRS 54 include a disc, jet, and illuminated outflow cavity walls (see Fig. B.2 of GL13). The geometry of the system is such that the disc is seen roughly edge-on (GL13). This geometry poses challenges from an observational perspective, specifically in viewing the inner disc where most of the accretion activity takes place. Nevertheless, the edge-on disc configuration of IRS 54 also provides good conditions in which to trace its bipolar jet back to the source (see below).

Images were generated from the SINFONI data for the H₂ and [Fe II] (1.644 μm) emission lines in the H and K bands. Figure 1 shows both the blue- and redshifted continuum-subtracted images of the emission from the [Fe II] line at 1.644 μm observed in 2014, where four spectral channels were averaged from 1.6430 μm (−171 km s⁻¹) to 1.6438 μm (−25 km s⁻¹, see Fig. 1a) and from 1.6440 μm (11 km s⁻¹) to 1.6443 μm (66 km s⁻¹, see Fig. 1b). This emission traces the jet of the YSO and is extended with respect to the source position, which is indicated with the black and white contours in Figs. 1a and 1b, respectively. It is spatially asymmetric about the central source, with much stronger blueshifted than redshifted emission. In summary, the IRS 54 jet predominately emits [Fe II] in the blueshifted lobe.

Fig. 1. Panel a: blueshifted component of the [Fe II] emission in the H band in IRS 54 from the SINFONI 2014 data. Four spectral channels were averaged from 1.6430 μm (−171 km s−1) to 1.6438 μm (−25 km s−1). Panel b: redshifted component of the [Fe II] emission in the H band in IRS 54 from the SINFONI 2014 data. Four spectral channels were averaged from 1.6440 μm (11 km s−1) to 1.6443 μm (66 km s−1). For reference, the centre black (a) and white (b) contours represent the continuum of the source taken at levels of 0.1, 0.3, 0.5, 0.7 and 0.9.

In contrast, most of the H₂ emission comes from the redshifted lobe: it traces not only the bright redshifted jet but also what appear to be cavity walls that straddle the source. Figure 2 shows the red- and blueshifted continuum-subtracted images of IRS 54 at the H₂ 1-0 S(1) emission line in the K band. In Fig. 2a four spectral channels were averaged from 2.12087 μm (−125 km s⁻¹) to 2.12160 μm (−22 km s⁻¹), and in Fig. 2b four spectral channels were averaged from 2.12185 μm (30 km s⁻¹) to 2.12258 μm (134 km s⁻¹). The blue contours represent the location of the central source and its continuum emission. The behaviour of this H₂ emission traces a different spatial component of the jet than that of the [Fe II] emission. The redshifted component of the jet is primarily radiating H₂ at 2.122 μm. The molecular jet was already observed to be asymmetric by GL13, with a redshifted molecular jet component and also possibly a blueshifted atomic jet component. Here, we observe this asymmetry as well and also observe the atomic component. Our observations therefore adhere to the morphology sketch presented by GL13 (their Fig. B.2).

Fig. 2. Panel a: blueshifted component of the H2 emission in the K band in IRS 54 from the SINFONI 2014 data. Four spectral channels were averaged from 2.12087 μm (−125 km s−1) to 2.12160 μm (−22 km s−1). Panel b: redshifted component of the H2 emission in the K band in IRS 54 from the SINFONI 2014 data. Four spectral channels were averaged from 2.12185 μm (30 km s−1) to 2.12258 μm (134 km s−1). For reference, the centre blue contours represent the position of the continuum of the source taken at levels of 0.2, 0.4, 0.6 and 0.8 of the continuum flux.

3.2. Spectroscopy and IFU on source

Figure 3 shows the on-source flux calibrated spectra (J, H, and K bands) along with the main lines detected from our observations (epochs 2005, 2013, and 2014) and from 2010 archival data (GL13). It is apparent from Fig. 3 that not only the flux intensities of both line and continuum have changed from one epoch to another, but also the shape of each spectral energy distribution (SED). Between 2005 (blue), 2010 (red), and 2013 (orange) the flux increased and the shape of the SED went from being approximately flat (especially in the K band) to having a much steeper slope. The SED of 2014 (green) receded to a flux below that of 2010 (red), becoming less steep in the K band. Because 2010 data were only available in the K band, it is impossible to say definitively whether this was the case in the J and H bands as well. Emission lines at different epochs have been identified and labelled in Fig. 3; Table 2 provides a list of the main lines detected along with their full width at half maximum (FWHM), fluxes, radial velocities and full width at zero intensity (FWZI). These quantities are analysed further in the coming sections to derive visual extinction (A_V) and mass accretion rates at different epochs.

Fig. 3. Spectral Energy Distribution (SED) of J, H, and K bands over the four epochs (2005, 2010, 2013, and 2014). We note that the J band has two overlapping spectra from 2005 (blue) and 2013 (orange).

The spectra from IRS 54 display multiple hydrogen-recombination lines, the brightest of which are the Brγ and Paβ emission lines in the K and J band, respectively. These lines are primarily accretion signatures (e.g. Muzerolle et al. 1998). Forbidden iron ([Fe II] 1.257 μm and 1.644 μm) and molecular hydrogen (H₂ 2.122 μm) emission lines, which trace the jet, are also visible. The [Fe II] emission is also useful in understanding physical properties of the surroundings of the YSO, such as the amount of foreground extinction in the observations. Also present are the R(0–15) and P(1–9) (between 2.31 and 2.37 μm) rotational lines (J) of the v = 2 − 0 CO band head, which trace relatively mild temperatures (a few hundred Kelvin). These lines are seen in absorption in the 2010 and 2014 data (see Fig. 4). Many of these line measurements are below 3σ, however the signal-to-noise ratio is higher in the 2014 spectrum, making these absorption features easier to identify in this epoch. The 2010 epoch also contains these CO lines, less clearly, but just barely visible in absorption. Notably, a star of spectral type M would have CO photospheric absorption lines, including the high rotational lines (i.e. those that actually pile up, producing the band heads), which typically trace gas at a few thousand Kelvin. The high-J lines are not seen here. Moreover, such a young source as IRS 54 should present very high veiling and thus photospheric lines should not be detected. Therefore these features in absorption most likely originate from the outer disc or from the envelope, seen against a much hotter inner disc gas. These features are observed in emission during outbursts of EXors (Audard et al. 2014).

Fig. 4. Low J lines from the v = 2 − 0 CO band head (2.31–2.37 μm) in the 2010 (red) and 2014 (green) data, seen in absorption. The vertical lines indicate the different lines.

3.3. Variability in visual extinction

[Fe II] transitions can be used to determine visual extinction (A_V), however, uncertainties are still prevalent when estimating the radiative transition probabilities used in the calculation (Giannini et al. 2015). The ratio of two bright NIR lines, 1.644/1.257 μm, is useful because they originate from the same upper level and are optically thin. The line ratio 1.644/1.320 μm can similarly be used to calculate A_V, however the signal-to-noise ratio of the 1.320 μm line in our observations is below 3σ. Because these transitions originate from the same upper level, their theoretical intensity ratio depends not on the physical conditions in the emission region, but on the frequencies and transition probabilities. The observed ratio is:

$$\begin{array}{c}\hfill \frac{I_{{\lambda }_1}}{I_{{\lambda }_2}}=\frac{A_{\mathit{ij}}{\lambda }_2}{A_{\mathit{ik}}{\lambda }_1}{10}^{-(A_{{\lambda }_1}-A_{{\lambda }_2})/2.5},\end{array}$$(1)

where I_λ1 and I_λ2 are the observed intensities of the two [Fe II] lines; A_ij and A_ik are Einstein coefficients of 4.65 s⁻¹ ([Fe II] at 1.644 μm) and 4.83 s⁻¹ ([Fe II] at 1.257 μm) representing the transition rates for each line (taken from Nussbaumer & Storey 1988); and λ₁ and λ₂ are the wavelengths of the lines. The following equation, along with an extinction law (namely, Rieke & Lebofsky 1985), allows for the calculation of A_V:

$$\begin{array}{c}\hfill A_{\lambda }=\alpha {\lambda }^{\beta }{A}_V;\alpha =0.42,\beta =-1.75,\end{array}$$(2)

where A_λ is the extinction at a specific wavelength λ.

To study if and how the visual extinction changes with time, line ratios were measured at different epochs. Visual extinction (A_V) was calculated for the ISAAC data (2005 and 2013 epochs) using Eq. (1) for the [Fe II] 1.644/1.257 μm line ratio from the spectra extracted on-source, using the flux values reported in Table 2. To increase the signal-to-noise ratio and thus reduce the uncertainties, we combine the fluxes of the different velocity components. As the SINFONI data (2010 and 2014 epochs) did not contain J band observations, a measurement of the integrated flux of the [Fe II] 1.257 μm line was not possible. The measured A_V value changes from 15 ± 1 mag in 2005 to 24 ± 1 mag in 2013 (as reported in Col. 2 of Table 3). These values are similar to those found in other studies towards Class I protostars (Davis et al. 2011). It is important to note that these values represent lower limits for the extinction on source, because the [Fe II] lines originate in the jet.

3.4. Accretion-Tracing Lines

The Brγ and Paβ lines are understood to primarily trace the accretion activity of the young star, rather than outflow activity further out (e.g. Muzerolle et al. 1998; Calvet et al. 2004; Alcalá et al. 2017). An empirical correlation exists between the accretion luminosity (L_acc) and the luminosity of these hydrogen-recombination lines (L_line). As expected of lines tracing the same processes, both of these lines follow similar trends of increasing flux in the accretion burst between 2005 and 2013. In 2014, a decrease is seen in the Brγ line (Fig. 5). The values of the fluxes of these lines can be found in Table 2.

Fig. 5. Line profiles of the (panel a) Brγ 2.166 μm and (panel b) Paβ 1.282 μm hydrogen-recombination lines. Both lines are continuum-subtracted and have been corrected to the cloud velocity and normalised to the 2013 peak flux densities. In (a), the ISAAC spectral data (epochs 2005 and 2013) have been smoothed to match the lower spectral resolution of the SINFONI data (epoch 2014). In (b), this was not necessary as only ISAAC data contained J band observations.

Figure 5a shows the continuum-subtracted line profiles of the Brγ (2.166 μm) line emission present in the data from the epochs 2005, 2010, 2013, and 2014 normalised to the peak intensity of the 2013 epoch. Due to the different spectral resolution of ISAAC and SINFONI, the ISAAC spectra were re-sampled to match the lower resolution of the SINFONI data for comparison purposes. Subsequently, this line was found to have an average FWZI of ∼765 km s⁻¹. The intensity of the Brγ line flux increased by about a factor of five from 2005 (blue) to 2013 (orange), and then in 2014 (green) dropped to an intensity of only ∼40% of the peak of 2013 (orange). As a tracer of accretion processes, this decrease in Brγ emission suggests that accretion decreased dramatically between the 2013 and 2014 epochs, a much sharper change than the peak increase from 2005 to 2013. The integrated flux and radial velocity of the Brγ emission across all four epochs can be found in Table 2.

In Fig. 5b, the Paβ (1.282 μm) line profiles from two epochs are shown normalised to the peak intensity of the 2013 epoch in the J band. It is apparent that the change in flux from 2005 to 2013 was much less pronounced at this wavelength range than in the K band (see Brγ in Fig. 5a); the intensity in 2005 was ∼70% that of the peak intensity in 2013. In comparison the Brγ line in 2005 only reached ∼25% of the peak intensity of that in 2013. However, it is worth noting that the line profiles shown in Fig. 5 are not dereddened. Once dereddened (see Figs. A.1 and A.2), the differences in flux between the 2005 and 2013 epochs are readily visible.

3.5. Outflow-tracing lines

Jets provide an environment where shocks can break up dust grains releasing, for example, refractory elements like Fe into the gas phase (see, e.g. Nisini et al. 2002; Nisini 2008). Two [Fe II] lines in particular (1.644 μm and 1.257 μm) are understood to be tracing the outflow activity in IRS 54 (Connelley & Greene 2014). They are consistent in their trends of flux increase and eventual decrease, as for the Brγ and Paβ lines, as seen in Fig. 6. This trend suggests that the YSO ejection activity close to the source follows the same path as that of accretion.

Fig. 6. Panels a and b: line profiles of [Fe II] at 1.257 μm and 1.644 μm. Panel c: line profile of molecular hydrogen (H2) 2.122 μm. All lines are continuum-subtracted and have been corrected to the cloud velocity. In (b) and (c), the ISAAC spectral data (epochs 2005 and 2013) have been smoothed to match the lower spectral resolution of the SINFONI data (epoch 2014); in (a), this was not necessary as only ISAAC data contained J band observations. We note the different velocity components present. The redshifted HVC in (b) cannot be separated as it is blended with the LVC. In (b) the Br12 (1.641 μm) line is also shown to highlight how the 2014 epoch (green line) changes differently in the outflow-tracing [Fe II] emission from the accretion-tracing Br12 line.

These forbidden emission lines can exhibit multiple and often complex velocity components (Davis et al. 2001). The high-velocity component (HVC) is generally associated with the jet at higher velocities on larger scales, while the low-velocity component (LVC) originates from a more compact and dense region at the base of the jet (see, e.g. Garcia Lopez et al. 2009). Both [Fe II] lines clearly show a HVC and a LVC, which are both blueshifted. A redshifted component of the HVC is also visible in the 1.257 μm line at ∼+100 km s⁻¹. In the 1.644 μm line, the emission at ∼+100 km s⁻¹ could potentially be the redshifted component of its HVC, but it is strongly blended with the LVC. Our difficulty in separating the different velocity components is not unexpected given the (almost edge-on) geometry of the disc.

In addition to atomic emission, molecular hydrogen (H₂) emission is also detected. This traces dense molecular gas of relatively low excitation (n_H2 ≥ 10⁵ cm⁻³, T ∼ 2000 K; Caratti o Garatti et al. 2006). The brightest transition detected is the 1–0 S(1) line at 2.122 μm (other H₂ lines are also present in the data, but their signal-to-noise ratio values are much lower). Figure 6c shows the line profile at different epochs. Notably, between 2013 and 2014 the intensity dropped to ∼20% of the peak, below even that of the earliest (2005) epoch. This trend can be seen in the [Fe II] 1.644 μm line as well, but not in the Bracket or Paschen hydrogen-recombination lines (whether this trend is also seen in the [Fe II] 1.257 μm line cannot be confirmed due to the lack of J band data in 2014). This is consistent with our expectations that the H₂ and [Fe II] emission are tracing different processes in IRS 54 than the Brackett and Paschen lines. The H₂ and [Fe II] lines trace the outflow, while the Brackett and Paschen lines trace accretion activity closer to the star and inner disc⁴⁵.

3.6. Photometry

In order to put our observations in perspective and investigate the variability of IRS 54 across a broader time frame and at different wavelengths, we combine our SINFONI and ISAAC photometry with archival and literature photometric data for IRS 54 obtained by 2MASS (2 Micron All Sky Survey), the SQIID (Simultaneous Quad-Color Infrared Imaging Device) at Kitt Peak National Observatory, DENIS (Deep Near Infrared Survey), NSFCAM (NASA), the Anglo-Australian Telescope (AAT), and the Wide-field Infrared Survey Explorer (WISE). These data can be seen in Fig. 7 and Table 4, where their respective sources are cited. The WISE data (MIR) were used to obtain an idea of how the object has varied at these longer wavelengths where extinction is less influential. J, H, and K bands all show a similar trend in luminosity variability, while the MIR observations (bands W1 and W2 at 3.4 and 4.6 um, respectively) show the brightness peaking in 2010. Their steep rise suggests that the maximum has happened around the epochs of 2010 or 2011 rather than in 2013. However, due to the gaps in photometry between 2010 and 2013, it is impossible to define when the maximum of each light curve takes places. As monitoring data were available from NEOWISE between 2014 and 2019, it can be seen from the photometry in 2014 that W1 and W2 magnitudes dropped with respect to 2010, and after that there is a smooth secondary maximum followed by an erratic dimming of the source.

Fig. 7. Photometry of IRS 54. The present archival data are described in Table 4. Our data from SINFONI and ISAAC are also included.

Multi-epoch archival data in the J, H, and K bands show a large variability of up to two magnitudes during the decade preceding our 2005 observations. The observed variability in IRS 54 seems therefore to be episodic rather than being a single event we happened to witness with our observations. 2MASS NIR archival data indicate that in 1999 the object was ∼1 mag brighter than the peak seen in this study. However, it is worth noting the differing spatial resolution and likely contamination from the surrounding nebulosity in the region in the 2MASS data. Nevertheless, as there are also DENIS data (in J and K bands) from the 1999 epoch that agree with the 2MASS data points, we can trust that IRS 54 brightened around 1999 and then a quiescent state followed until the more recent surge we observed between 2010 and 2013.

3.7. Accretion variability

The mass accretion rate (Ṁ_acc) can be derived from measuring the release of accretion energy (L_acc) in the form of UV continuum and line emission (Gullbring et al. 1998). The relation between L_acc and Ṁ_acc is expressed in Eq. (3) (Gullbring et al. 1998; Hartmann et al. 1998):

$$\begin{array}{c}\hfill L_{\mathrm{acc}}=\frac{G{M}_{\star }{\dot{M}}_{\mathrm{acc}}}{R_{\star }}(1-\frac{R_{\star }}{R_{\mathrm{in}}}),\end{array}$$(3)

where R_⋆ is the radius of the YSO and R_in is the inner (truncation) radius. We adopted the values M_⋆ ∼ 0.15 M_⊙ and R_⋆ ∼ 2 R_⊙ (GL13). The truncation radius is assumed to be R_in = 5 R_⋆ (Gullbring et al. 1998).

To calculate L_acc, we used the empirical relation between line luminosity and accretion luminosity given by Alcalá et al. (2017) for the hydrogen-recombination lines Paβ (1.282 μm) and Brγ (2.166 μm). The observed line flux is reddened by the dust in the YSO system. The degree to which this reddening occurs can be calculated using the following equation:

$$\begin{array}{c}\hfill F_{\lambda }^0=F_{\lambda }^{\mathrm{obs}}{10}^{A_{\lambda }/2.5},\end{array}$$(4)

where $F_{\lambda }^0$ is the line flux emitted from the source corrected for extinction, F_obs is the measured line flux (affected by extinction), and A_λ is the extinction measured at a specific wavelength. It is especially important to take this reddening into account when studying these early stages of star formation, as protostars are deeply embedded in their parent cloud and the value of A_λ is not only non-negligible but can be quite large. With the extinction-corrected flux and assuming d = 137 ± 5 pc (Sullivan et al. 2019), we derive the line luminosities and use the following relation to derive L_acc:

$$\begin{array}{c}\hfill log\left(\frac{L_{\mathrm{acc}}}{L_{\odot }}\right)=alog\left(\frac{L_{\mathrm{line}}}{L_{\odot }}\right)+b,\end{array}$$(5)

where a and b are the parameters derived in Alcalá et al. (2017) and depend on the particular line being measured: the Brγ line has corresponding values of a = 1.19 ± 0.10 and b = 4.02 ± 0.51 and the Paβ line has corresponding values of a = 1.06 ± 0.07 and b = 2.76 ± 0.34.

The derived Ṁ_acc values can be found in Table 3. We derive an average from these two lines to obtain a value of Ṁ_acc = (1.7 ± 0.5) × 10⁻⁸ M_⊙ yr⁻¹ in 2005 and (2.6 ± 0.5) × 10⁻⁷ M_⊙ yr⁻¹ in 2013. This shows that the mass accretion rate increases by one order of magnitude between 2005 and 2013. Notably, the A_V increases in tandem with the Ṁ_acc, although an increase in extinction is normally associated with a decrease in flux. However, in the case of IRS 54 the accretion burst is so strong that the flux increase is apparent even with the extinction dampening it.

To determine A_V, we used the [Fe II] line ratio, which originates from the jet rather than directly from the accretion region, where one would expect the HI lines to originate. As the jet is further from the surface of the star, the extinction would only decrease from the central source. However, as we extracted spectra on source, the extinction measured is on a part of the jet that is very close (within 1″–1.2″ i.e. within ∼150 au of the source) to the central object, and is therefore a reasonable estimate. Even so, our values of A_V should be taken as lower limits. Extinction dampens the apparent accretion luminosity, and therefore the apparent Ṁ_acc. As a lower limit for A_V has been calculated, the actual Ṁ_acc may be higher than reported. This is explored further in Sect. 4.

For the 2010 archival data, GL13 estimate extinction on source using a different method from ours, exploiting the measured versus expected ratio of H₂ lines. This leads to an estimate of A_V of ∼30 mag. Using this value, GL13 estimate Ṁ_acc ∼ 3 × 10⁻⁷ M_⊙ yr⁻¹ from the relation between L_line and L_acc derived by Calvet et al. (2004). In order to compare Ṁ_acc in 2010 with what was found in the other epochs (2005 and 2013), we recalculated it using the same Alcalá et al. (2017) relation between L_line and L_acc used throughout this study. The result is Ṁ_acc = 4.3 × 10⁻⁷ M_⊙ yr⁻¹, which is slightly higher than that found in GL13. This result would then strengthen the idea of the burst maximum being closer to the 2010-2011 epochs rather than the 2013 one.

4. Discussion

4.1. Accretion and extinction variability

From our observations, it is clear that IRS 54 underwent significant changes in luminosity, accretion, and extinction during the period from 2005 to 2014. We interpret this sharp increase in flux as an accretion burst that peaked between 2010 and 2013.

The Ṁ_acc between 2005 and 2013 (Table 3) increases by a factor of ∼20, while A_V increases by nine magnitudes (flux increases by ∼4000 in the V band and by ∼2.4 in the K band). While both of these are significant changes, it is clear that the accretion appears to have the larger effect on the SED of IRS 54 between this time period because we see an overall increase in luminosity, especially in the K band where the continuum flux increases by a factor of approximately six. The photometry also reflects this increase in flux during the burst (Fig. 7), supporting the idea of an accretion burst. The MIR data show that there is a large increase of over a magnitude in 2010 followed by a decrease and with subsequent small fluctuations from 2014 to 2019. The MIR shows a sharper increase in flux within 2010 than seen in the NIR which is likely due to extinction affecting these latter wavelengths to a lesser extent. To measure the effect of visual extinction on a given wavelength, the ratio A_V/A_λ can be used. In the case of the MIR, A_V/A_L is ∼17, and A_V/A_M is ∼43 (Rieke & Lebofsky 1985). Namely, A_V = 15 mag and A_V = 24 mag (the values found in this study) translate to 0.88 mag and 1.41 mag in the L band, and to only 0.35 mag and 0.56 mag in the M band, respectively.

As the extinction was seen to increase between 2005 and 2013, it is important to note that the spectra shown in Fig. 3 are reddened due to extinction. A_V was thus calculated in the data where appropriate emission lines were available for the analysis (see Sect. 3.3). Between 2005 and 2013, we find a large increase in A_V (from ∼15 mag to ∼24 mag). As line and continuum fluxes are also seen to increase in this time frame, this change cannot be accounted for solely by variable extinction. All increases of line fluxes that we measure, especially in the K band, are in fact more pronounced in the de-reddened spectra (see Figs. A.1 and A.2). The de-reddened spectra indicate that the 2010 epoch was close to the peak of the burst, as its flux is similar to that of 2013 in the K band, when both are corrected for extinction. However, it is important to remember that the A_V in 2010 was estimated in GL13 using H₂ lines rather than the [Fe II] lines which were used for 2005 and 2013, and therefore we must be cautious when comparing these values.

Both extinction and accretion can affect the shape of the SED of a protostar, but in different ways. An increase in accretion would result in a more pronounced increase in luminosity at shorter wavelengths. For example, in the NIR, it would increase the flux in the J band more than in the K band, effectively flattening the SED. Alternatively, an increase in extinction would result in a steeper slope of the SED in this wavelength range, decreasing the observed flux more at the shorter wavelengths (J band) than at the longer wavelengths (K band). Our results indicate that the change in SED shape observed is not due solely to one or the other of these phenomena, but to a combination of the two. We see an increase in steepness in the slope with increased extinction, but also an increase of flux in the H and especially the K bands. This combination of both accretion and extinction is representative of the complex processes at work and how they are related to one another in IRS 54. The combined effect of these two processes is quite unique, because it tends to flatten the light curves more at shorter wavelengths. This can be seen qualitatively in Fig. 7, where the J band light curve is much flatter with respect to the K, W1, or W2 bands.

A possible explanation for this tandem increase of both extinction and accretion is that the increase in accretion and ejection lifts a large amount of dust from the disc, which crosses the line of sight and therefore produces more extinction. The edge-on geometry of the system supports this interpretation as any dust lifted as a result of an accretion or ejection burst would easily intersect the line of sight between the observer and the source. Therefore, this system demonstrates an accretion or ejection burst activity that also increases the visual extinction along the line of sight. A similar increase of extinction was observed for RW Aur, whose photometric and polarimetric variability were explained by the presence of dust in the disc wind (Dodin et al. 2019; Koutoulaki et al. 2019).

The mass accretion rates derived in this study are likely lower limits, as the method used to determine A_V utilised outflow-tracing lines, which originate further from the source where the extinction is lower. As can be seen in Table 3, the values of Ṁ_acc derived from the two different lines are slightly different. Incorporating a higher A_V value into the calculation of Ṁ_acc would provide a higher mass accretion rate. The accuracy of our A_V estimates is thus investigated further by plotting Ṁ_acc derived from the Paβ and Brγ lines and the difference in accretion luminosities measured using the Paβ and Brγ lines (ΔL_acc) as a function of A_V for the 2005 (Fig. 8 top panel) and 2013 (Fig. 8 bottom panel) epochs. The Ṁ_acc measured from each line should theoretically agree if the A_V was measured accurately. As can be seen in Fig. 8, the agreement between Paβ and Brγ in 2005 is strongest in the low extinction limit, and it diverges for A_V values larger than ∼20 mag. In the 2013 case, the difference between Ṁ_acc calculated using either line is largest at A_V values ∼25 mag, which is where our measurement of A_V ∼ 24 mag resides. However, this strengthens our assumption that this is a lower limit for extinction, as values approaching 30 mag provide better agreement. The red star in Fig. 8 represents the value for A_V that minimises the disparity between the Ṁ_acc values derived from Brγ and Paβ measurements, which is at A_V∼ 18.5 mag in 2005 and ∼29.5 mag in 2013. This is most likely close to the actual extinction toward the source in 2013. These A_V values correspond to Ṁ_acc ∼ 3.1 × 10⁻⁸ M_⊙ yr⁻¹ in 2005 and Ṁ_acc ∼ 7.3 × 10⁻⁷ M_⊙ yr⁻¹ in 2013, respectively.

Fig. 8. Ṁacc and ΔLacc vs. AV for the 2005 (top panel) and 2013 (bottom panel) epochs. The blue line represents the measurement of Ṁacc using the Brγ line luminosity and the orange line represents the measurement of Ṁacc using the Paβ line luminosity. The red stars represent the value of AV where the Brγ and Paβ de-reddened fluxes best agree in calculating the Ṁacc.

4.2. Jet variability and asymmetry

Examining the H₂ line emission maps in Fig. B.1, the 2014 epoch exhibits lower luminosity in line emission than that of the 2010 epoch. We sought to investigate whether the luminosity in the region of the jet was also changing. In principle, the jet would not be expected to change instantaneously with a change in accretion. Although accretion and ejection are linked, the observed jet emission originates from further out in the system and consequently takes longer to reflect an increase in accretion. Therefore, between 2010 and 2014, where accretion was observed to change, the jet would not necessarily change accordingly. To determine this, spectra from multiple regions of the same size across the SINFONI image were extracted in 2010 and 2014 and compared. Unexpectedly, we found the flux of the regions extracted along the jet to vary between the two epochs. A possible explanation is that these regions are contaminated by scattered light, which does reflect the changes in accretion luminosity. The cavity walls surrounding the central source are also illuminated by this scattered light and their flux also increases during the burst.

It is clear that the jet of IRS 54 has both an atomic and a molecular component, as seen in Figs. 1 and 2 where the H₂ line emission is predominantly redshifted and the [Fe II] line emission is predominantly blueshifted. This implies different excitation conditions and possibly different velocities of the jet material ejected in the two lobes. This asymmetry could be due in part to the inhomogeneity of the interstellar medium (ISM) in the region, or is potentially an effect of a misalignment in the magnetic fields of the protostar, but investigating this issue is beyond the scope of this study.

It is also interesting to note that during the 2014 epoch dimming, the flux from the jet-tracing lines (H₂ and [Fe II]) dropped to a flux lower than that of 2005. In contrast, the flux from accretion-tracing lines (Paβ, Br12 and Brγ) during 2014 dropped only to a flux between that of 2005 and 2010, which reflects how the continuum changed during this period across all bands (J, H, and K). This difference in how the flux changed in 2014 is most obviously seen in Fig. 6b where Br12 line is adjacent to the [Fe II] 1.644 μm emission line. This implies that the processes of accretion and ejection are not behaving synchronously following the burst, as we would expect considering they originate from different regions.

4.3. A new EXor object?

The increase of mass accretion rate seen in IRS 54 is consistent with that observed in EXor objects. In addition, EXor burst timescales and frequencies are also in line with our observations of IRS 54. As our observations are at NIR wavelengths, the observed increase in luminosity of IRS 54 is lower than what would be observed in typical, much less embedded EXors, which are usually observed at optical wavelengths and are therefore expected to be larger. If one considers this and the fact that IRS 54 is a very low-mass star, its luminosity increase is likely on par with typical Exor-type bursts. The brightness of IRS 54 increases by more than one magnitude in the K band (and varies by more than two magnitudes over 20 years in the J, H, and K bands) and Ṁ_acc increases by at least an order of magnitude.

We also note a similar behaviour of EX Lup and IRS 54 light curves in the MIR in terms of strength, duration and shape (see Fig. 1 of Ábrahám et al. 2019, and our Fig. 7). As for IRS 54, EX Lup W1 and W2 light curves (during the 2011 burst) have a similar shape, showing a steep rise at the beginning and a smoother flickering decline lasting a few years, as indeed seen in IRS 54 (see top panels of Fig. 7). In both objects the MIR brightness increases by a couple of magnitudes.

Both the duration (a few years) and intensity of the burst hint at an EXor-type event, although some of the typical disc spectral features are missing, such as the CO band head lines in emission and Na lines. This might be due to the system geometry of IRS 54, which prevents us from seeing the inner disc, as these signatures originate from that region. On the other hand, the hydrogen lines detected in IRS 54 that trace accretion are likely seen in scattered light. Indeed, the low-J lines from the v = 2 − 0 CO band head are seen in absorption, suggesting a cold gas (with temperature of a few hundred Kelvin), likely originating from the outer disc atmosphere, absorbing emission coming from a hotter (inner) region emitting CO. This might possibly hint at CO band heads in emission, typical of EXor bursts, being present in the inner gaseous disc of IRS 54. While it is not possible at this time to definitively decipher whether or not IRS 54 is an EXor-type object, the issue is worthy of further investigation, as this source would be the first very low-mass protostar where an EXor-type burst is observed.

5. Conclusions

In the course of this study, data obtained with ISAAC and SINFONI over four epochs (2005, 2010, 2013, and 2014) were reduced and analysed to assess variability in IRS 54. This was done over the wavelength range 1.24 μm to 2.45 μm (J, H, and K bands). Significant changes in flux were found between the epochs, reflecting an increase between 2005 and 2013 and a drop in 2014. The lightcurves in the W1 and W2 bands show a similar behaviour, showing a steep rise in 2010 and possibly a secondary maximum after 2014. This increase in luminosity is accompanied by a burst in the mass accretion rate, Ṁ_acc, which increases from ∼1.7 × 10⁻⁸ M_⊙ yr⁻¹ in 2005 to ∼2.6 × 10⁻⁷ M_⊙ yr⁻¹ in 2013. This burst is consistent with the photometric data that we have been gathering during these same epochs. Going back to archival data from 1999, the photometry illustrates that this protostar went through a previous change in luminosity between 1999 and 2005 of approximately two magnitudes in the J, H, and K bands. These two large changes in flux suggest these bursts may be episodic.

Specific emission lines were analysed for this variability and to calculate Ṁ_acc and A_V, which demonstrated an increase in tandem with one another. Maps of the jet-tracing emission (H₂ and [Fe II]) were also generated to examine how the emission varied spatially about the central star of the YSO. The [Fe II] was found to be emitting predominantly from the blueshifted component of the jet, and the H₂ was found to be illuminating the cavity walls with scattered light, while primarily originating from the redshifted component of the jet. This asymmetry is notable and may help understand the inner mechanism at work in the YSO. Two of the possible causes for this asymmetry in emission from the jet are that the interstellar medium in the region may be inhomogeneous or that there may be misaligned magnetic fields in the protostar.

Examining the SED of each epoch and its variability, we deduce that the changes it exhibits reflect a combination of both an increase in accretion and extinction. A possible explanation for this tandem increase in these two parameters is that the increased accretion and ejection activity during the burst lifts up material into the line of sight and obscures the YSO.

The timescales of the burst seen in IRS 54 are reminiscent an EXor-type object, and its increase in luminosity is exceptionally large over a short period of time, especially considering it is a VLMS. Further investigation of IRS 54 as a potential EXor-type object is warranted, as if found to be true, IRS 54 would be the lowest mass Class I source observed to have this type of violent bursts in accretion and ejection.

Acknowledgments

The authors wish to thank the anonymous referee for their comments. C. S., A. C. G., P. Mc G., and T. P. R. would like to acknowledge funding from the European Research Council under Advanced Grant No. 743029, Ejection, Accretion Structures in YSOs (EASY). R. G. L. acknowledges support from Science Foundation Ireland under Grant No. 18/SIRG/5597. R. F. acknowledges support from the Chalmers Initiative on Cosmic Origins (CICO) postdoctoral fellowship.

References

Appendix A: De-reddened fluxes of the accretion and outflow-tracing lines

Fig. A.1. De-reddened flux of the accretion-tracing (Paβ and Brγ) lines from IRS 54. These plots were generated using the extinction values found in Table 3 for the 2005 and 2013 epochs. The 2010 epoch uses the extinction value found in GL13.

Fig. A.2. De-reddened flux of the outflow-tracing ([Fe II] and H2) lines from IRS 54. These plots were generated using the extinction values found in Table 3 for the 2005 and 2013 epochs. The 2010 epoch uses the extinction value found in GL13.

Appendix B: H₂ emission line maps from 2010 and 2014

Fig. B.1. Continuum-subtracted SINFONI data from the 2010 (left panel) and 2014 (right panel) epochs taken by summing five spectral channels centred on the H2 emission line at 2.122 μm in the K band. The blue contours represent the spatial position of the continuum of the source taken at levels of 0.2, 0.4, 0.6, and 0.8 of the continuum flux.

All Tables

All Figures

Fig. 1. Panel a: blueshifted component of the [Fe II] emission in the H band in IRS 54 from the SINFONI 2014 data. Four spectral channels were averaged from 1.6430 μm (−171 km s−1) to 1.6438 μm (−25 km s−1). Panel b: redshifted component of the [Fe II] emission in the H band in IRS 54 from the SINFONI 2014 data. Four spectral channels were averaged from 1.6440 μm (11 km s−1) to 1.6443 μm (66 km s−1). For reference, the centre black (a) and white (b) contours represent the continuum of the source taken at levels of 0.1, 0.3, 0.5, 0.7 and 0.9.

In the text

Fig. 2. Panel a: blueshifted component of the H2 emission in the K band in IRS 54 from the SINFONI 2014 data. Four spectral channels were averaged from 2.12087 μm (−125 km s−1) to 2.12160 μm (−22 km s−1). Panel b: redshifted component of the H2 emission in the K band in IRS 54 from the SINFONI 2014 data. Four spectral channels were averaged from 2.12185 μm (30 km s−1) to 2.12258 μm (134 km s−1). For reference, the centre blue contours represent the position of the continuum of the source taken at levels of 0.2, 0.4, 0.6 and 0.8 of the continuum flux.

In the text

	Fig. 3. Spectral Energy Distribution (SED) of J, H, and K bands over the four epochs (2005, 2010, 2013, and 2014). We note that the J band has two overlapping spectra from 2005 (blue) and 2013 (orange).
In the text

	Fig. 4. Low J lines from the v = 2 − 0 CO band head (2.31–2.37 μm) in the 2010 (red) and 2014 (green) data, seen in absorption. The vertical lines indicate the different lines.
In the text

Fig. 5. Line profiles of the (panel a) Brγ 2.166 μm and (panel b) Paβ 1.282 μm hydrogen-recombination lines. Both lines are continuum-subtracted and have been corrected to the cloud velocity and normalised to the 2013 peak flux densities. In (a), the ISAAC spectral data (epochs 2005 and 2013) have been smoothed to match the lower spectral resolution of the SINFONI data (epoch 2014). In (b), this was not necessary as only ISAAC data contained J band observations.

In the text

Fig. 6. Panels a and b: line profiles of [Fe II] at 1.257 μm and 1.644 μm. Panel c: line profile of molecular hydrogen (H2) 2.122 μm. All lines are continuum-subtracted and have been corrected to the cloud velocity. In (b) and (c), the ISAAC spectral data (epochs 2005 and 2013) have been smoothed to match the lower spectral resolution of the SINFONI data (epoch 2014); in (a), this was not necessary as only ISAAC data contained J band observations. We note the different velocity components present. The redshifted HVC in (b) cannot be separated as it is blended with the LVC. In (b) the Br12 (1.641 μm) line is also shown to highlight how the 2014 epoch (green line) changes differently in the outflow-tracing [Fe II] emission from the accretion-tracing Br12 line.

In the text

	Fig. 7. Photometry of IRS 54. The present archival data are described in Table 4. Our data from SINFONI and ISAAC are also included.
In the text

	Fig. 8. Ṁacc and ΔLacc vs. AV for the 2005 (top panel) and 2013 (bottom panel) epochs. The blue line represents the measurement of Ṁacc using the Brγ line luminosity and the orange line represents the measurement of Ṁacc using the Paβ line luminosity. The red stars represent the value of AV where the Brγ and Paβ de-reddened fluxes best agree in calculating the Ṁacc.
In the text

	Fig. A.1. De-reddened flux of the accretion-tracing (Paβ and Brγ) lines from IRS 54. These plots were generated using the extinction values found in Table 3 for the 2005 and 2013 epochs. The 2010 epoch uses the extinction value found in GL13.
In the text

	Fig. A.2. De-reddened flux of the outflow-tracing ([Fe II] and H2) lines from IRS 54. These plots were generated using the extinction values found in Table 3 for the 2005 and 2013 epochs. The 2010 epoch uses the extinction value found in GL13.
In the text

	Fig. B.1. Continuum-subtracted SINFONI data from the 2010 (left panel) and 2014 (right panel) epochs taken by summing five spectral channels centred on the H2 emission line at 2.122 μm in the K band. The blue contours represent the spatial position of the continuum of the source taken at levels of 0.2, 0.4, 0.6, and 0.8 of the continuum flux.
In the text