Spectral Evolution and Radial Dust Transport in the Prototype Young Eruptive System EX Lup^*

We obtained N-band spectra of EX Lup using VISIR (Lagage et al. 2004) on the European Southern Observatory's (ESO) Very Large Telescope (097.C-0639, PI: P. Ábrahám). The observations were performed on 2016 August 20 with an exposure time of 2300 s, and a slit size of 1'' × 3174, which provided a spectral resolution of R ∼ 350. We observed HIP 78820 as our standard star, but to correct for airmass differences we also used HD 178345, that was measured during the same night close to the zenith, from the ESO archive. For the basic data reduction and the extraction of the spectra we run the ESO VISIR spectroscopic pipeline. In order to correct for telluric features, we divided the target spectrum by a spectrum derived via interpolation between the two standard star measurements, one obtained at higher and one at lower airmass than EX Lup. Flux calibration was carried out by multiplying by the model spectra of the standard stars.

EX Lup was observed in the post-outburst phase with the mid-infrared interferometer MIDI (Leinert et al. 2003), as part of the program 091.C-0668 (PI: S. Antoniucci). In total, 17 interferometric observations were carried out on 2013 May 27/28, using the U2–U3 and U3–U4 baselines, with baseline lengths ranging from 31 to 62 m. The MIDI data are publicly available, and have been already published in an online interferometric atlas of MIDI observations of low- and intermediate-mass young stars (Varga et al. 2018). The reduced data consist of 7.5–13 μm total and correlated spectra, spectrally resolved visibilities, and differential phases. The visibilities imply that the object was barely resolved at most baselines. Both the total and the correlated spectra show the 10 μm silicate emission feature.

Since all MIDI spectra from 2013 exhibited very similar spectral shapes, we averaged them to increase the signal-to-noise ratio (S/N). Special care was needed to average correlated spectra, as data taken at different baselines and position angles are generally not comparable. Before the averaging process, we calculated Gaussian sizes from each measurement at 10.7 μm in the same manner as in Varga et al. (2018). By combining the resulting sizes with the baseline position angles we could estimate the disk inclination (${62}_{\ -{20}^{^\circ }}^{\circ +{7}^{^\circ }}$) and position angle (${81}_{\ -{39}^{^\circ }}^{\circ +{4}^{^\circ }}$). These numbers are roughly consistent with i = 32°–38° and PA = 65°–78°, derived from ALMA maps for the outer disk (Hales et al. 2018). Assuming elliptical symmetry, we determined an effective baseline length for each observation as if we would observe a face-on circular disk. Since in such a configuration the results do not depend on the position angle, we averaged all correlated spectra within two bins of effective baseline of 21–31 and 43–51 m. These average spectra represent disk regions inside 4.1–2.8 au and 2.0–1.7 au radii, respectively (adopting the Gaia DR2 distance of 157 pc; Gaia Collaboration et al. 2018). Additionally, we averaged all total spectra, in order to have a spectrum representative of the whole circumstellar disk.

EX Lup was observed with the IRS (Houck et al. 2004) on board the SST (Werner et al. 2004) at six epochs between 2004 and 2009 (Table 1). Apart from the fourth epoch on 2008 May 2, the observations were performed using the short-low module (5.2–14.5 μm) of the low-resolution (R = 60–120) spectrograph. For the first and the two last epochs also, low-resolution spectra using the long-low module (14–35 μm) were obtained. In addition to the low-resolution observations, for epoch one to four, spectra were obtained using the short-high (9.9–19.5 μm) and long-high (18.7–37.2 μm) modules of the high-resolution (R = 600) spectrograph. All observations, with the exception of epoch four, used a Pointing Calibration and Reference Sensor (PCRS) peak-up, to accurately center EX Lup at the correct field-of-view position of the IRS instrument. Observations from four epochs were already published in Juhász et al. (2012), but here we reprocessed them with the latest calibration files.

In the case of a low-resolution mode, the data reduction process started from the droopres intermediate data product. For the high-resolution data we used the rsc data product as a starting point. All data products were processed through the Spitzer Science Center pipeline version S18.18.0. For the spectral extraction and flux calibration we used the data reduction packages developed for the c2d and feps legacy programs (Lahuis & Boogert 2003; Bouwman et al. 2008; Carpenter et al. 2008). The background emission in the low-resolution data has been subtracted using associated pairs of imaged spectra from the two nod positions, also eliminating stray light contamination and anomalous dark currents. Pixels flagged by the data pipeline as being "bad" were replaced with a value interpolated in the dispersion direction from an elongated eight pixel perimeter surrounding the flagged pixel. The spectra were extracted using a six pixel and five pixel fixed-width aperture in the spatial dimension for the short-low and the long-low modules, respectively. The low-level fringing at wavelengths >20 μm was removed using the irsfringe package (Lahuis & Boogert 2003). The spectra were calibrated with a position dependent spectral response function derived from IRS spectra and MARCS stellar models for a suite of calibrators provided by the Spitzer Science Center. To remove any effect of pointing offsets in the short-low module data, we matched orders based on the point-spread function (PSF) of the IRS instrument, correcting for possible flux losses (Swain et al. 2008). We estimate the accuracy of our flux estimates to be at the level of 1%.

For the high-resolution data we applied an optimal source profile extraction method which fits an analytical PSF derived from sky-corrected calibrator data and an extended emission component, derived from the cross-dispersion profiles of the flat-field images, to the cross-dispersed source profile. It is not possible to correct for the sky contribution in the high-resolution spectra, subtracting the two nod positions as with the low-resolution observations, due to the small slit length. The fourth epoch observations had a dedicated background observation which we used to subtract the background emission. For the other high-resolution observations, we used the background estimate from the source profile fitting extraction method to remove the background emission. For correcting "bad" pixels we used the irsclean package. We further removed low-level (∼1%) fringing using the irsfringe package. The flux calibration for the high-resolution spectrograph has been done in a similar way as for the low-resolution observations. For the relative spectral response function, we also used MARCS stellar models and calibrator stars provided through the SSC. The spectra of the calibration stars were extracted in an identical way to our science observations using both extraction methods. As with the low-resolution observations, we also corrected for possible flux losses due to pointing offsets. We estimate the absolute flux calibration uncertainty for the high-resolution spectra to be ∼3%, slightly higher than that of the low-resolution observations. Given that the flux calibration of the short-low module of the IRS spectrograph is the most precise, for those epochs where these were available, the flux of the high-resolution data was scaled to the flux of the low-resolution spectra. All scaling factors are consistent with the expected uncertainties.

EX Lup was monitored inbetween 2009 February 12 and September 28 with the 15 cm RoBoTT telescope at the Universitätssternwarte Bochum near Cerro Armazones.⁷ The observations were carried out in the i-band (λ_eff = 752 nm, zero-magnitude flux f₀ = 3631.0 Jy). The good weather conditions allowed us almost daily data sampling; in total there were 134 observing nights. The light curve of EX Lup is constructed by using 20 calibration stars on the same exposures. The absolute photometric calibration was achieved by comparing EX Lup to Landolt standard fields which were observed during the same period of time with the same telescope. More details on the telescope, observations, and reduction steps can be found in Haas et al. (2012).

In order to check whether the weakening of the crystalline silicate features could be related to changes in the irradiation of the disk, i.e., in the luminosity of the central star, we constructed optical and mid-infrared light curves of EX Lup between 2007 and 2018 (Figure 1). The data imply that no other eruption of amplitude similar to the one in 2008 has occurred in this period. Localized brightness fluctuations, however, appeared several times.

A remarkable brightness variability pattern can be seen in early 2009, a few months after the end of the large outburst, in both the V and i bands. It began with a deep dimming in late 2008, followed by a rapid brightening of ∼0.5 mag in the V band in early 2009, and then by a gradual fading until the fall of 2009. The brightness increase corresponds in time to a moderate rise of the accretion rate (a factor of few to 10 with respect to the minimum quiescence value) detected spectroscopically by Sicilia-Aguilar et al. (2015). Due to the lack of multiwavelength measurements, it is unclear whether a drop in extinction along the line of sight had also contributed to this brightness increase. The gradual fading during 2009 is documented in two optical bands (Figure 1), and the relationship between these two light curves could be fitted with a first order polynomial as Δi ∼ (0.67 ± 0.06) × ΔV. The slope term agrees, within the formal uncertainties, with the variability amplitude ratio predicted by the interstellar extinction law for these two wavelengths (0.65; Cardelli et al. 1989). Thus, the available observations cannot exclude that both accretion and extinction processes have been involved in the variability processes in 2009, but the data do not allow us to determine their respective contributions. Variable extinction could have been caused by some rearrangement of the inner disk structure as a consequence of the outburst.

There are indications for a brightening also in early 2011, but it is less documented, having only a few infrared data points on the rising part and only visual estimates close to the peak. Simultaneous spectroscopy also reveals a slight increase in the accretion rate during this date (Sicilia-Aguilar et al. 2015), which nevertheless stays more than one order of magnitude below the 2008 outburst accretion rate value. Similar short-term variations in the accretion rate of EX Lup are found throughout its observed history (e.g., Lehmann et al. 1995). Apart from these two events, EX Lup exhibited no obvious brightness variations in the last decade.

The list of mid-infrared spectroscopic observations available for EX Lup is presented in Table 1 (we did not include the VLTI/MIDI observations from 2008 June and July, since the present study focuses on the post-outburst phase and the outburst period is well represented by the higher S/N Spitzer spectra). Our new VISIR spectrum from 2016, together with the average total spectrum of the MIDI observations from 2013, are plotted in Figure 2. For comparison, we also overplotted all Spitzer spectra. Following Juhász et al. (2012), we subtracted a spline fit continuum from each Spitzer spectrum. The VISIR and MIDI observations included only a short continuum on both sides of the 10 μm silicate feature, thus we subtracted a simple linear trend.

In Figure 2 both the VISIR and MIDI spectra exhibit very similar spectral shapes. The peaks and shoulders that were prominent in the Spitzer spectra due to the presence of crystalline silicates are not visible in these two recent ground-based spectra. Their spectral profiles are reminiscent of the pre-outburst Spitzer observations in 2004 and 2005. The spectrum from 2005 March 18 was modeled by Sipos et al. (2009), who concluded that the emitting dust consisted of small amorphous silicates. Thus, here we may conclude that by 2013, and also later in 2016, any signature of the crystalline forsterite grains formed in 2008 had disappeared from the 10 μm feature of EX Lup.

Due to the different baseline combinations, the MIDI interferometric observations offer a possibility to compare the silicate spectra in three different radial zones in the post-outburst period. We found that all the total, and the two correlated spectra (representative of disk regions at r < 2 au and r < 4 au, respectively; Section 2.2), exhibit 10 μm silicate features characteristic of pristine, amorphous interstellar grains. While in 2008 the MIDI correlated spectra still outlined crystalline spectral features (Juhász et al. 2012), forsterite grains disappeared from the whole inner disk by 2013.

We performed a quasi-static radiative transfer modeling of EX Lup in its pre-outburst, outburst, and post-outburst states, using the RADMC3D code (Dullemond et al. 2012). For the pre-outburst (quiescent) phase, we used a disk model based on Sipos et al. (2009), with only minor modifications due to the updated knowledge of disk inclination from ALMA (Hales et al. 2018). Following Sipos et al. (2009), we introduced a curved inner rim, approximated by a vertical step in the radial density distribution. More details about the disk model are presented in the Appendix and in Table 2. For the post-outburst phase, we supplemented the disk with a tenuous axisymmetric cloud above the disk surface. The disk component was fixed to the best-fit quiescent model (Figure 2).

Checking the light curves of Figure 1, it is a reasonable assumption that there were no strong flare-ups since 2008 that could have produced new crystals (the size of the region heated above 1000 K has always been smaller than the inner radius of the disk). To comply with the information that crystallization was an episodic process in a limited period, we adopted an expanding hollow geometry for the crystal-rich component, with the following assumptions. We postulate that all crystal grains were created during the outburst, within the radial range between 0.3 and 0.7 au (Ábrahám et al. 2009). This cloud of crystals started expanding during the outburst, pushed away by a stellar or disk wind (for a discussion on the possible origin of this wind, see Section 4.2) and we assume that the crystalline dust cloud continues expanding after the end of the outburst. We prescribe the cloud structure to have a constant width of ${r}_{\mathrm{out}}-{r}_{\mathrm{in}}=0.4\,\mathrm{au}$, where r_in is increasing with time. We assume that the mass and dust composition of the expanding cloud is constant during the post-outburst phase. The best-fitting value of the cloud mass will be determined by our modeling. Therefore, its optical depth, and in turn the fraction of stellar light that the cloud can absorb and re-emit in the infrared, is decreasing approximately as $\propto {r}_{\mathrm{in}}^{-2}$. While the adopted geometry is an ad-hoc assumption, the cloud is mainly optically thin, thus its integrated emission does not depend on the actual spatial distribution of the dust particles. At each epoch of observation we will fit the shape of the 10 μm emission feature and deduce the actual radius of the dust cloud, from which we can constrain the expansion velocity between the epochs. While we will derive some kinematical information, our study is not a dynamical modeling, rather a quasi-static sequence of equilibrium models.

Based on these assumptions, first we fitted the cloud parameters to match the Spitzer spectra on 2008 October 10 and 2009 April 6. In the best-fit models, the expanding shell consists of 95% forsterite and 5% amorphous carbon. The optical properties of amorphous carbon were calculated with a distributed hollow sphere model (Min et al. 2005), assuming grain size of 0.1 μm, using complex refractive index from Jager et al. (1998). For the optical properties of forsterite we used the results from a laboratory test (Koike et al. 2003), following Juhász et al. (2012) and Ábrahám et al. (2009). The crystalline mass in the cloud is 1.9 × 10²³ g, equivalent to 3.2 × 10⁻⁵ M_⊕. The obtained radius r_in is 1.2 au for 2008 October 10, and 1.5 au for 2009 April 6.

In order to reproduce the observations in 2013 and 2016, we also run models with increasingly larger r_in. We found that it is about ${r}_{\mathrm{in}}\geqslant 3\,\mathrm{au}$ where the crystalline features in the 10 μm peak fade below the detection limit due to the low temperature of the grains. We also fitted the Spitzer observations during the outburst (2008 April 21 and May 2), although with some different assumptions. Since the crystal formation was still ongoing at these two dates, we allowed the total forsterite mass to be a free parameter, but fixed the inner radius of the crystalline cloud to ${r}_{\mathrm{in}}=0.3\,\mathrm{au}$. In addition, during the outburst we included an accretion heating term in our model. The best-fitting models are plotted in Figure 2, together with the actual density distribution. In 2013 and 2016 we could not determine the exact location of the dust cloud. Thus, in the figure we plotted models with radii of 3 au for both epochs, although these are lower limits only.

Our modeling suggested that during the ∼6 months that separated the two post-outburst Spitzer measurements (2008 October 10 and 2009 April 6), the expanding cloud moved from ${r}_{\mathrm{in}}=1.2$ to 1.5 au. This corresponds to an expansion velocity of ∼3 km s⁻¹. If we assume that after the end of the outburst the crystalline dust cloud was not further accelerated, a higher threshold for the actual location of the cloud is set by a constant velocity expansion of 3 km s⁻¹ as marked with the straight line in Figure 3. For a lower threshold, we take into account that after 2013 the cloud must have been at least 3 au from the star (this result, derived from the spectral fitting, is also supported by our interferometric observations in Section 3.2). The likely position of the inner radius of the expanding forsterite cloud is delimited by the two curves (shadowed area). Note that by quadratically adding the Keplerian orbital velocity of the cloud at 1.2 au to the 3 km s⁻¹ radial expansion velocity the result is ≈17.5 km s⁻¹, which is significantly lower than the escape velocity of the system at this radius, ≈24.3 km s⁻¹. It implies that the expanding cloud will not leave the EX Lup system, but the trajectories of the particles will cross the circumstellar disk at some distances.

Zoom In Zoom Out Reset image size

The best candidate for driving the motion of the silicate crystals from the inner disk outwards is a wind, either from stellar or disk origin. Optical spectroscopic observations of EX Lupi in both quiescence (Sicilia-Aguilar et al. 2015) and outburst (Sicilia-Aguilar et al. 2012) confirm the presence of wind features. The wind signatures are ubiquitous during the outburst maximum, appearing not only in Hα and the Balmer lines, but also in the Ca ii H and K lines and in the Fe ii emission lines. At least two wind components with different densities and temperatures are observed, peaking at ∼−200 and ∼−100 km s⁻¹ in outburst (Sicilia-Aguilar et al. 2012), and their strength decays quickly as the accretion rate decreased. The presence of wind in quiescence is less ubiquitous, and manifests itself mostly in the Hα and Ca ii H and K lines, which show rapidly variable absorption features at velocities ∼−100 to −20 km s⁻¹ that are correlated with the accretion rate, as expected for wind-related features (Sicilia-Aguilar et al. 2015). The slowest component of the wind is also detected in the Ca ii H and K lines. These wind velocities would be nevertheless too high compared to those required for the models, although lacking information on the direction of the wind and considering that EX Lup is observed at a low inclination angle, the component of the wind parallel to the disk surface could be significantly smaller than the radial velocity component.

Very low-velocity disk winds have been identified in EX Lup via forbidden line emission by Fang et al. (2018) and Banzatti et al. (2019). These works revealed significant differences in the observed wind velocities in quiescence and in outburst. In quiescence, a low-velocity component was detected at −1.5 km s⁻¹ (Banzatti et al. 2019). In outburst, however, velocities were measured at −12 km s⁻¹ (Banzatti et al. 2019) and in the range of −15 to −18 km s⁻¹ (Fang et al. 2018). Very slow winds and slow outflows have also been detected in FUors (Ruíz-Rodríguez et al. 2017).

The main question is whether these highly variable wind components would be able to induce motion of the silicate grains formed within the disk. While determining the origin and physics of possible wind mechanisms in the EX Lup system is beyond the scope of this paper, there are claims in the literature that winds can have a significant role in radial transportation of grains at distances of 1 au. A possible example is the recently proposed scenario for the radial transportation of grains (in particular the calcium-aluminium-rich inclusions, CAIs) in the early protosolar disk by Liffman et al. (2016), which is based on the interaction of the stellar magnetic field with the inner disk rim.

According to Figure 3, by now the crystals are at 3–7 au from the central star. From the temperature profile of the EX Lup disk (Ábrahám et al. 2009, Figure 4 in Supplementary material) we can conclude that the water snowline, corresponding to ∼160 K, is at about 5 au at the disk surface, and much closer in the midplane. Thus, we suggest that the crystalline particles, produced in the large outburst in 2008, already approached, and in the near future potentially cross the water snowline. If part of them would fall back onto the disk, they could mix with the icy grain population there, and become incorporated into the forming planetesimals, comets. This scenario would help to solve the long-standing mystery of high crystalline fraction observed in pristine comets in the solar system (e.g., Hanner et al. 1994).

In order to judge the potential impact of an outburst, like that of EX Lup in 2008, on the mineralogical evolution of the disk, we can estimate the mass of forsterite crystalline particles in our model and compare it with typical comet masses. In our model, the total amount of forsterite in the expanding cloud is 9.5 × 10⁻¹¹ M_⊙ or 1.9 × 10²³ g. As an alternative check of this number, we can compute the expected amount of freshly formed crystals on the disk surface during the outburst. Assuming that crystallization took place in a region between 0.3 and 0.7 au radii, the surface area of crystal production is A ≈ 1.2 au². We assume that crystallization happened in the hot disk surface where the vertical optical depth τ ≤ τ₀ ≈ 1 for the stellar radiation. Then the vertical column density of this surface layer is τ₀/κ_abs, where κ_abs is the absorption coefficient. Assuming that annealing is efficient in this hot surface and all amorphous silicate grains are turned into crystals, the total crystallized silicate mass in this surface is therefore M = τ₀/κ_abs × A ≈ τ₀ × 1.1 × 10²³ g. This is consistent with the total crystal mass used in our post-outburst models of ∼1.9 × 10²³ g, giving an independent support to the mass estimate in our model, and to the assumption of an optically thin dust cloud. The crystalline mass of 1.9 × 10²³ g is equal to ∼10⁴ times the mass of comet Hale–Bopp, thus the outward transportation of processed crystalline material might contribute significantly to the building blocks of cometary material.

It was a long-standing puzzle that while the extra outburst heat provides optimal conditions for the transformation of pristine amorphous silicate grains into crystalline ones, no sign of crystallization in the outbursts of young eruptive stars was ever seen (see e.g., the infrared spectroscopic study of a sample of FUors by Quanz et al. 2007). Our infrared spectroscopic observations of EX Lup's outburst in 2008 revealed crystalline silicates and their formation in a young eruptive star. From that, one would naively expect that a high observed level of crystallinity can be used as a proxy for the erupting nature, and identify dormant eruptive stars in their quiescent phase. However, the results presented in this paper demonstrate that—at least in an individual case—crystallinity disappears on a timescale of a few years. Thus, instead of the high level of crystallinity, perhaps it is a rapid temporal variation of the degree of crystallinity that could signal a recent overlooked eruption. Detecting it would require a multiepoch high S/N spectral monitoring in the mid-infrared of a larger sample of pre-main-sequence stars, for which the Mid-infrared Instrument (MIRI, Rieke et al. 2015) on board James Webb Space Telescope (JWST) will be the best instrument in the coming years.

Due to its unprecedented sensitivity, JWST will also be able to observe EX Lup, and possibly localize the actual position of the crystalline dust cloud. In Figure 3 we marked the epoch of a possible JWST observation set arbitrarily to 2022 July 20. The location of the crystalline cloud is unknown at this epoch, thus we calculated the expected spectral shapes for three different radii within the expected range (shaded area in Figure 3). The closest possible distance of r_in = 3 au is the most favorable case for detection, the intermediate distance at 5 au corresponds to the snowline, and the largest distance at r = 9.68 au is consistent with the constant expansion velocity in Figure 3.

We calculated density distributions and integrated spectra for all three radii following the same method as described in Section 4.1. Then we used the JWST Exposure Time Calculator V1.3 (Pontoppidan et al. 2016), to simulate observations with the MIRI Medium Resolution Spectrograph (MRS; Rieke et al. 2015; Wells et al. 2015), assuming an effective integration time of 1000 s for each of the three dichroic settings, to cover the full wavelength range between 5 and 28 μm. The spectral resolution of the resulting synthetic MIRI spectra was reduced by binning the original full resolution of MRS (R = 2000–3000) to be comparable with the the existing EX Lup observations (R = 127). This step increased the S/N of the resulting spectra by a factor of ∼20. As a next step, we subtracted a continuum from each of the three spectra in the same way as we did for the Spitzer observations. The resulting spectra and the corresponding density distributions are plotted in Figure 4.

Zoom In Zoom Out Reset image size

All three model spectra at the different radii exhibit amorphous 10 μm peaks, which is not surprising as the same was found already at earlier epochs when the dust cloud was still closer to the star (Figure 2). However, weak forsterite peaks at longer wavelengths, emitted by cold crystalline grains, may be detectable for JWST/MIRI. From Figure 4 we can conclude that if the expanding shell is situated at 3 au, then we have a clear detection of the 16, 19, and 24 μm forsterite bands. The intermediate distance at R = 5 au seems to be the detection limit for the forsterite bands exhibiting some marginal features, while at 9.68 au none of the bands are visible any more. The results imply that the long-wavelength forsterite features are detectable with JWST by up to three times larger distances than the radius of the most distant Spitzer observation. If the shell is between 3 and 5 au in 2022, then JWST will determine its actual position, outlining the further path of the expansion motion. A non-detection will probably signal a relatively high expansion velocity and an actual location of >5 au. Both results will provide constraints on the effectiveness of the radial mixing in the EX Lup disk.