<figure><img src="image_1.png" /></figure>

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

P. Ábrahám¹, L. Chen¹, Á. Kóspál¹,² †, J. Bouwman², A. Carmona³, M. Haas⁴, A. Sicilia-Aguilar⁵, C. Sobrino Figaredo⁴, R. van Boekel², and J. Varga¹,⁶

¹Konkoly Observatory, Research Centre for Astronomy and Earth Sciences, Hungarian Academy of Sciences, Konkoly-Thege Miklós út 15-17, 1121 Budapest, Hungary; abraham@konkoly.hu

²Max Planck Institute for Astronomy, Königstuhl 17, D-69117 Heidelberg, Germany

³Université de Toulouse, UPS-OMP, IRAP, Toulouse F-31400, France

⁴Astronomisches Institut, Ruhr-Universität Bochum, Universitätsstrasse 150, D-44801, Bochum, Germany

⁵SUPA, School of Science and Engineering, University of Dundee, Nethergate, Dundee, DD1 4HN, UK

⁶Leiden Observatory, Leiden University, P.O. Box 9513, NL2300, RA Leiden, The Netherlands

Received 2019 August 11; revised 2019 October 25; accepted 2019 October 27; published 2019 December 17

## Abstract

EX Lup is the prototype of a class of pre-main-sequence eruptive stars defined by their repetitive outbursts lasting several months. In 2008 January–September EX Lup underwent its historically largest outburst, brightening by about 4 mag in visual light. In previous studies we discovered ongoing silicate crystal formation in the inner disk during the outburst, but also noticed that the measured crystallinity fraction started decreasing after the source returned to the quiescent phase. Here we present new observations of the 10 µm silicate feature, obtained with the MIDI and VISIR instruments at Paranal Observatory. The observations demonstrate that within five years practically all crystalline forsterite disappeared from the surface of the inner disk. We reconstruct this process by presenting a series of parametric axisymmetric radiative transfer models of an expanding dust cloud that transports the crystals from the terrestrial zone to outer disk regions where comets are supposed to form. It is possible the early Sun also experienced similar flare-ups, and the forming planetesimals might have incorporated crystalline silicate material produced by such outbursts. Finally, we discuss how far the location of the dust cloud could be constrained by future *James Webb Space Telescope* observations.

*Unified Astronomy Thesaurus concepts: Pre-main sequence stars (1290)*

## 1. Introduction

EX Lup is the prototype of EXors, a class of young pre-main-sequence eruptive stars defined by their repetitive outbursts lasting several months (Herbig 1977; Herbig et al. 2001; Herbig 2007). The flare-ups represent periods of temporarily increased accretion from the circumstellar disk onto the star, possibly in a process similar to the more energetic outbursts of the FU Orionis-type objects (Hartmann & Kenyon 1996). In 2008, EX Lup underwent its historically largest outburst, brightening by about 4 mag in visual light (Figure 1). Our group obtained a 5–37 µm spectrum with the InfraRed Spectrograph (IRS) of the *Spitzer Space Telescope* (SST) on 2008 April 21, shortly after the peak of the outburst. Comparing our spectrum with a similar pre-outburst observation from 2005, we discovered that the originally amorphous silicate grains were transformed to crystalline particles in the inner disk due to the outburst heat (Ábrahám et al. 2009). This was the first direct observation of ongoing silicate crystallization in a celestial object. The presence of crystalline silicates was confirmed by our VLTI/MIDI interferometric observations obtained in 2008 June and July (Juhász et al. 2012).

The outburst of EX Lup ended in 2008 September. In order to document any subsequent changes in the silicate features, we obtained post-outburst *Spitzer* spectra on 2008 October 10 and 2009 April 6 (Figure 2). We expected that the narrow

crystalline features in the 10 µm emission (the sharp peak at 10 µm and the shoulder at 11.3 µm) become relatively stronger after the outburst. It is because in the high state both a central crystallized zone, and an outer amorphous region contribute to the mid-infrared silicate emission. In the low luminosity phase, however, the area that is warm enough to emit at 10 µm shrinks, and may partly or completely overlap with the zone where surface crystallization had occurred during the outburst ($r < 0.7$ au Ábrahám et al. 2009). But our observations clearly showed that the degree of crystallinity in the 10 µm feature decreased after the outburst, while strong peaks related to cold forsterite grains developed in the 30 µm range (Juhász et al. 2012).

In order to disentangle the effects potentially responsible for the decay of crystalline peaks around 10 µm, Juhász et al. (2012) performed radiative transfer and turbulent mixing calculations. Their results excluded that vertical mixing replaced the freshly formed crystals by amorphous particles from the interior of the disk. As an alternative mechanism, they speculated about outward transportation of the silicate crystals by a radial wind driven by the outburst. The appearance of forsterite peaks in the 15–30 µm range in the late phase of and after the outburst would then be interpreted as the arrival of crystalline particles at larger distances from the star where they cool down and emit at longer wavelengths, supporting the outward radial transport scenario.

The quick drop of crystallinity after an outburst could have been suspected because of an earlier observational result, too. In 1955–56, EX Lup had already produced a major eruption, with a peak visual brightness of 8.5 mag, very close to the

<footnote>* Based on observations collected at the European Organisation for Astronomical Research in the Southern Hemisphere under ESO programmes 091.C-0668 and 097.C-0639.</footnote>
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Light curves of EX Lup. V-band observations before 2010 and visual brightness estimates are from the AAVSO database (http://www.aavso.org). V-band data points from 2016 and later are from the ASAS-SN survey (Shappee et al. 2014; Kochanek et al. 2017). The *i*-band data were observed with the RoBoTT telescope at the Universitätssternwarte Bochum (Section 2.4). The 3.4 and 4.6 µm photometry was taken by the WISE satellite (Wright et al. 2010) and published in the AllWISE Multiepoch Photometry Table and in the NEOWISE-R Single Exposure (L1b) Source Table. Additional 3.6 and 4.5 µm data were taken with the Spitzer Space Telescope and were published in Kóspál et al. (2014). For better visibility, the infrared light curves are shifted by the values listed in the upper right corner. Vertical dashed lines mark the epochs when mid-infrared spectra were observed (see also Table 1).</figcaption></figure>

maximum of the outburst in 2008 (Herbig 1977). The comparable peak luminosities of the two outbursts imply that similar amounts of crystalline silicates must have formed. However, by 2004–2005 (the dates of two pre-outburst Spitzer spectra) the crystalline spectral features vanished, setting an upper limit of 50 yr for their disappearance. Aiming to outline the fate of the crystalline particles after the 2008 outburst and to understand the physical processes in action, here we present new *N*-band observations of EX Lup and model the disappearance of crystalline particles by their transportation to more distant disk regions.

**2. Observations and Data Reduction**

**2.1. Mid-infrared Spectroscopy with VLT/VISIR**

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" × 31"/74, 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.

**2.2. Archival VLTI/MIDI Observations**

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°+7°-20°) and position angle (81°+4°-39°). 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.

**2.3. Archival Spitzer Spectroscopy**

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
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Left and center; continuum subtracted spectra of EX Lup at different epochs. Black curves are observations, red curves show our models. Right: density distributions (unit: g cm⁻³) of the EX Lup system in the modeled epochs.</figcaption></figure>

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
<table><caption>**Table 1**
Log of Mid-infrared Spectra of EX Lup Analyzed in This Paper</caption><thead><tr><th>Date</th><th>Telescope</th><th>Instrument/<br>Spectral Module</th><th>AOR/<br>Prop. ID</th></tr></thead><tbody><tr><td>2004 Aug 30</td><td>Spitzer</td><td>ch0, ch1, ch2, ch3</td><td>5645056</td></tr><tr><td>2005 Mar 18</td><td>Spitzer</td><td>ch0, ch1, ch3</td><td>11570688</td></tr><tr><td>2008 Apr 21</td><td>Spitzer</td><td>ch0, ch1, ch3</td><td>27039232</td></tr><tr><td>2008 May 2</td><td>Spitzer</td><td>ch1, ch3</td><td>27063808</td></tr><tr><td>2008 Oct 9</td><td>Spitzer</td><td>ch0, ch2</td><td>28476672</td></tr><tr><td>2009 Apr 6</td><td>Spitzer</td><td>ch0, ch2</td><td>28476416</td></tr><tr><td>2013 Apr 27/28</td><td>VLTI</td><td>MIDI</td><td>091.C-0668</td></tr><tr><td>2016 Aug 20</td><td>VLT</td><td>VISIR</td><td>097.C-0639</td></tr></tbody></table>

Note. EX Lup underwent a large outburst in 2008 January–August, thus the first two spectra represent a pre-outburst state, the ones in 2008 April–May were taken close to the peak of the outburst, while the remaining spectra are post-outburst observations.

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.

## 2.4. Optical Photometric Monitoring

EX Lup was monitored in between 2009 February 12 and September 28 with the 15 cm RoBoTT telescope at the Universitätssternwarte Bochum near Cerro Armazones.<fnref n="7" /> The observations were carried out in the *i*-band ($\lambda_{eff}$ = 752 nm, zero-magnitude flux $f_0$ = 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).

# 3. Results

## 3.1. Brightness Evolution of EX Lup

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 $\Delta i \sim (0.67 \pm 0.06) \times \Delta 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).

<footnote n="7">⁷ https://www.astro.ruhr-uni-bochum.de/Astrophysik/all_info_en.html</footnote>
Apart from these two events, EX Lup exhibited no obvious brightness variations in the last decade.

## 3.2. Mid-infrared Spectroscopy

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.

## 4. Discussion

Analyzing the long wavelength parts of the *Spitzer* spectra of EX Lup, Juhász et al. (2012) reported the appearance of crystalline silicate features at 23, 28, and 33 µm after the end of the outburst (see the 2008 October 10 and 2009 April 6 panels in Figure 2). This was a strong argument against the post-outburst destruction of the crystals, either by amorphization due to high-energy radiation or by being accreted into the star. The result was also not compatible with the vertical mixing scenario, because if the fresh crystals were rapidly mixed down below the disk surface then the appearance of any long-wavelength forsterite peak would also not be expected. Thus, Juhász et al. (2012) proposed a qualitative scenario in which the crystallized grains are transported outward, and after cooling down they emit radiation primarily at wavelengths longward of 10 µm. In the following we will adopt and perform a quantitative study of this hypothesis, by fitting the crystalline features, and determining the mass and actual radial distance of the crystals for the different epochs.

<table><caption>**Table 2**
Parameters of the Best-fit Quiescent Model</caption><thead><tr><th></th><th></th><th></th></tr></thead><tbody><tr><td colspan="3"><strong>System Parameters</strong></td></tr><tr><td>Distance (d)</td><td>157</td><td>pc</td></tr><tr><td>Inclination (i)</td><td>35</td><td>°</td></tr><tr><td>Visual extinction (A<sub>V</sub>)</td><td>0</td><td></td></tr><tr><td colspan="3"><strong>Stellar Parameters</strong></td></tr><tr><td>Temperature (T<sub>*</sub>)</td><td>3800</td><td>K</td></tr><tr><td>Luminosity (L<sub>*</sub>)</td><td>0.54</td><td>L<sub>☉</sub></td></tr><tr><td colspan="3"><strong>Disk Component</strong></td></tr><tr><td>Inner radius (R<sub>in</sub><sup>d</sup>)</td><td>0.3</td><td>au</td></tr><tr><td>Outer radius (R<sub>out</sub><sup>d</sup>)</td><td>150</td><td>au</td></tr><tr><td>Outer scale height (h<sub>out</sub><sup>d</sup>)</td><td>0.12</td><td></td></tr><tr><td>Power-law index for scale height (q)</td><td>0.09</td><td></td></tr><tr><td>Power-law index for surface density (p)</td><td>-1</td><td></td></tr><tr><td>Mass (M<sub>dust</sub>)</td><td>2.5 × 10<sup>-4</sup></td><td>M<sub>☉</sub></td></tr><tr><td colspan="3"><strong>Inner rim curvature:</strong></td></tr><tr><td>Scale height of the rim (h<sub>rim</sub>)</td><td>0.05</td><td></td></tr><tr><td>Ending point of the rim (R<sub>rim</sub>)</td><td>0.6</td><td>au</td></tr><tr><td colspan="3"><strong>Dust composition</strong></td></tr><tr><td>Diameter (a)</td><td>0.1</td><td>µm</td></tr><tr><td colspan="3"><strong>Dust species:</strong></td></tr><tr><td>Amorphous carbon</td><td>5</td><td>%</td></tr><tr><td>Amorphous olivine</td><td>95</td><td>%</td></tr><tr><td colspan="3"><strong>Cloud Component</strong></td></tr><tr><td>Inner radius (R<sub>in</sub><sup>c</sup>)</td><td>1.2 (2008 Oct)<br>1.5 (2009 Apr)</td><td>au<br>au</td></tr><tr><td>Radial width (R<sub>out</sub><sup>c</sup> - R<sub>in</sub><sup>c</sup>)</td><td>0.4</td><td>au</td></tr><tr><td>Center height<sup>a</sup> (h<sub>center</sub><sup>c</sup>)</td><td>0.6</td><td></td></tr><tr><td>Vertical size<sup>a</sup> (Δ h<sub>cloud</sub>)</td><td>1.2</td><td></td></tr><tr><td>Mass (M<sub>cloud</sub>)</td><td>9.5 × 10<sup>-11</sup></td><td>M<sub>☉</sub></td></tr><tr><td colspan="3"><strong>Dust composition</strong></td></tr><tr><td>Diameter (a)</td><td>0.1</td><td>µm</td></tr><tr><td colspan="3"><strong>Dust species:</strong></td></tr><tr><td>Amorphous carbon</td><td>5</td><td>%</td></tr><tr><td>forsterite</td><td>95</td><td>%</td></tr></tbody></table>

**Notes.** Those in italics were kept fixed during the modeling.

<sup>a</sup> We assume the cloud to extend vertically from $z/r = h_{center}^{c} - \Delta h_{cloud}/2$ to $z/r = h_{center}^{c} + \Delta h_{cloud}/2$.

### 4.1. Modeling

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_{\text{out}} - r_{\text{in}} = 0.4$ au, where $r_{\text{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_{\text{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 \times 10^{23}$ g, equivalent to $3.2 \times 10^{-5} M_{\oplus}$. The obtained radius $r_{\text{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_{\text{in}}$. We found that it is about $r_{\text{in}} \ge 3$ 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_{\text{in}} = 0.3$ 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.

## 4.2. Transportation of the Crystals

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_{\text{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

<figure><img src="image_4.png" /><figcaption>**Figure 3.** Predicted radial locations of the silicate crystals in the EX Lup system. Black rectangle marks the time interval of the outburst (x-axis) and the 0.3–0.7 au zone where crystallization occurred during the outburst (y-axis). Filled black dots are the *Spitzer*, MIDI, and VISIR observations. For the arbitrarily selected epoch of a future *James Webb Space Telescope (JWST)* observation three possible radii of the expanding dust crystalline cloud are marked. Asterisks denote the approximate location of the snowline.</figcaption></figure>

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.

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\alpha$ and the Balmer lines, but also in the CaII H and K lines and in the FeII 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\alpha$ and CaII 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 CaII 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 \text{ km s}^{-1}$ (Banzatti et al. 2019). In outburst, however, velocities were measured at $-12 \text{ km s}^{-1}$ (Banzatti et al. 2019) and in the range of $-15$ to $-18 \text{ km s}^{-1}$ (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 \times 10^{-11} M_{\odot}$ or $1.9 \times 10^{23} 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 \approx 1.2 \text{ au}^2$. We assume that crystallization happened in the hot disk surface where the vertical optical depth $\tau \le \tau_0 \approx 1$ for the stellar radiation. Then the vertical column density of this surface layer is $\tau_0/\kappa_{\text{abs}}$, where $\kappa_{\text{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 = \tau_0/\kappa_{\text{abs}} \times A \approx \tau_0 \times 1.1 \times 10^{23} g$. This is consistent with the total crystal mass used in our post-outburst models of ~$1.9 \times 10^{23} 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 \times 10^{23} g$ is equal to ~$10^4$ 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.

### 4.3. Detection of the Crystalline Dust Cloud with JWST

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_{\text{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 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.

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
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Left and center: continuum subtracted model spectra of EX Lup at different radii as seen by JWST. Right: density distributions of the EX Lup system at the modeled radii.</figcaption></figure>

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.

**5. Summary and Conclusions**

The large outburst of EX Lup in 2008 crystallized the surface of its inner disk. The strength of the crystalline signatures in the 10 µm silicate emission peak, however, started weakening soon after the eruption, and our new observations demonstrate that by 2013 no crystalline features were visible in the 10 µm spectrum. Juhász et al. (2012) claimed that the >20 µm parts of the *Spitzer* spectra taken in 2008–09 exhibit emission features from cold forsterite grains. The appearance of these long-wavelength indicates that the crystalline particles formed in the outburst were not destroyed by high-energy radiation, not accreted into the star, or not mixed down to the disk midplane, but they could have been transported to the outer disk.

Here we perform a quantitative analysis of this scenario, by fitting at each epoch the spectral shape of the 10 µm emission feature by a parametric axisymmetric radiative transfer model of an expanding dust cloud of mainly crystalline silicates. Our modeling suggests that in this scenario the expanding cloud must be at a stlocentric radius larger than 3 au by 2013, roughly corresponding to the snowline radius in the system. A detailed dynamical modeling of the trajectory of the expanding cloud and the paths of the individual crystalline particles is out of the scope of the present paper. Our results, however, suggest that the freshly created crystals could have reached the distance of the water snowline, and if a fraction of them fell back onto the disk they could have been incorporated into forming comets. The crystalline mass of the expanding cloud of

1.9 × 10<sup>23</sup> g (~10,000 Hale–Bopp-like comets) provides an upper limit for the amount of crystals transported to the outer disk in a single outburst.

Our scenario of an expanding cloud could only be validated by new >20 µm observations as the crystals have cooled down by now and radiate at longer wavelengths. We present a simulation of a JWST observation, and demonstrate that with its extreme sensitivity JWST/MIRI will be able to detect the forsterite crystals out to 5 au radius. The successful detection of the expanding motion will support the conclusion that episodic outbursts of young stars may contribute to the build-up of the crystalline dust component ubiquitously seen in comets.

This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 Research and Innovation programme under grant agreement No. 716155 (SACCRED). This publication makes use of data products from the Near-Earth Object Wide-field Infrared Survey Explorer (NEOWISE), which is a project of the Jet Propulsion Laboratory/California Institute of Technology. AllWISE makes use of data from WISE, which is a joint project of the University of California, Los Angeles, and the Jet Propulsion Laboratory/California Institute of Technology, and NEOWISE, which is a project of the Jet Propulsion Laboratory/California Institute of Technology. WISE and NEOWISE are funded by the National Aeronautics and Space Administration.

**Facilities:** VLT:Melipal (VISIR), VLTI (MIDI), *Spitzer* (IRS), OCA:RoBoTT.

**Software:** ESO VISIR spectroscopic pipeline, *Spitzer* Science Center pipeline (vS18.18.0), *Spitzer* irsfringe package (Lahuis & Boogert 2003), RADMC3D (Dullemond et al. 2012), JWST Exposure Time Calculator V1.3 (Pontoppidan et al. 2016).

**Appendix**

**Pre-outburst Disk Model**

Our pre-outburst disk model is based on Sipos et al. (2009), with only minor modification due to the updated knowledge of disk inclination (Hales et al. 2018).
For the dust disk, we assumed a power-law radial distribution for its surface density $\Sigma$,

$$ \Sigma(r) = \Sigma_{\text{in}} \left(\frac{r}{R_{\text{in}}}\right)^p, \quad R_{\text{in}} < r < R_{\text{out}}, \qquad (1) $$

where $\Sigma_{\text{in}}$ is the surface density at the inner radius $R_{\text{in}}$, $p$ is the power-law exponent, and $R_{\text{out}}$ is the outer radius of the dust disk. We assumed the vertical dust density distribution to be a Gaussian function, and the resulting density structure of the disk is

$$ \rho(r, z) = \Sigma(r) \frac{1}{\sqrt{2\pi}H} \exp\left(-\frac{z^2}{2H^2}\right), \qquad (2) $$

where $\rho(r, z)$ denotes the dust density as a function of $r$ and the height $z$ above the midplane, and $H$ is the scale height. We assumed that the dependence of $H$ on $r$ is also a power law, except for the region close to the inner disk rim,

$$ h(r) = \frac{H(r)}{r} = \begin{cases} h_{\text{rim}}, & R_{\text{in}} < r < R_{\text{rim}} \\ h_{\text{out}} \left(\frac{r}{R_{\text{out}}}\right)^q, & R_{\text{rim}} < r < R_{\text{out}}, \end{cases} \qquad (3) $$

where $h(r)$ is the dimensionless scale height, $h_{\text{in}}$ is the dimensionless scale height at the inner radius, and $q$ is the power-law index.

## ORCID iDs

Á. Kóspál https://orcid.org/0000-0001-7157-6275

## References

* Ábrahám, P., Juhász, A., Dullemond, C. P., et al. 2009, Natur, 459, 224

* Banzatti, A., Pascucci, I., Edwards, S., et al. 2019, ApJ, 870, 76

* Bouwman, J., Henning, T., Hillenbrand, L. A., et al. 2008, ApJ, 683, 479

* Cardelli, J. A., Clayton, G. C., & Mathis, J. S. 1989, ApJ, 345, 245

* Carpenter, J. M., Bouwman, J., Silverstone, M. D., et al. 2008, ApJS, 179, 423

* Dullemond, C., Juhasz, A., Pohl, A., et al. 2012, RADMC-3D: A multi-purpose radiative transfer tool, Astrophysics Source Code Library, ascl:1202.015

* Fang, M., Pascucci, I., Edwards, S., et al. 2018, ApJ, 868, 28

* Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2018, A&A, 616, A1

* Haas, M., Hackstein, M., Ramolla, M., et al. 2012, AN, 333, 706

* Hales, A. S., Pérez, S., Saito, M., et al. 2018, ApJ, 859, 111

* Hanner, M. S., Lynch, D. K., & Russell, R. W. 1994, ApJ, 425, 274

* Hartmann, L., & Kenyon, S. J. 1996, ARA&A, 34, 207

* Herbig, G. H. 1977, ApJ, 217, 693

* Herbig, G. H. 2007, AJ, 133, 2679

* Herbig, G. H., Aspin, C., Gilmore, A. C., Imhoff, C. L., & Jones, A. F. 2001, PASP, 113, 1547

* Houck, J. R., Roellig, T. L., van Cleve, J., et al. 2004, ApJS, 154, 18

* Jager, C., Mutschke, H., & Henning, T. 1998, A&A, 332, 291

* Juhász, A., Dullemond, C. P., van Boekel, R., et al. 2012, ApJ, 744, 118

* Kochanek, C. S., Shappee, B. J., Stanek, K. Z., et al. 2017, PASP, 129, 104502

* Koike, C., Chihara, H., Tsuchiyama, A., et al. 2003, A&A, 399, 1101

* Kóspál, Á, Mohler-Fischer, M., Sicilia-Aguilar, A., et al. 2014, A&A, 561, A61

* Lagage, P. O., Pel, J. W., Authier, M., et al. 2004, Msngr, 117, 12

* Lahuis, F., & Boogert, A. 2003, in SFChem 2002: Chemistry as a Diagnostic of Star Formation, ed. C. L. Curry & M. Fich (Ottawa: NRC Press), 335

* Lehmann, T., Reipurth, B., & Brandner, W. 1995, A&A, 300, L9

* Leinert, C., Graser, U., Richichi, A., et al. 2003, Msngr, 112, 13

* Liffman, K., Cuello, N., & Paterson, D. A. 2016, MNRAS, 462, 1137

* Min, M., Hovenier, J. W., & de Koter, A. 2005, A&A, 432, 909

* Pontoppidan, K. M., Pickering, T. E., Laidler, V. G., et al. 2016, Proc. SPIE, 9910, 991016

* Quanz, S. P., Henning, T., Bouwman, J., et al. 2007, ApJ, 668, 359

* Rieke, G. H., Wright, G. S., Böker, T., et al. 2015, PASP, 127, 584

* Ruíz-Rodríguez, D., Cieza, L. A., Williams, J. P., et al. 2017, MNRAS, 468, 3266

* Shappee, B. J., Prieto, J. L., Grupe, D., et al. 2014, ApJ, 788, 48

* Sicilia-Aguilar, A., Fang, M., Roccatagliata, V., et al. 2015, A&A, 580, A82

* Sicilia-Aguilar, A., Kóspál, Á., Setiawan, J., et al. 2012, A&A, 544, A93

* Sipos, N., Ábrahám, P., Acosta-Pulido, J., et al. 2009, A&A, 507, 881

* Swain, M. R., Bouwman, J., Akeson, R. L., Lawler, S., & Beichman, C. A. 2008, ApJ, 674, 482

* Varga, J., Ábrahám, P., Chen, L., et al. 2018, A&A, 617, A83

* Wells, M., Pel, J.-W., Glasse, A., et al. 2015, PASP, 127, 646

* Werner, M. W., Roellig, T. L., Low, F. J., et al. 2004, ApJS, 154, 1

* Wright, E. L., Eisenhardt, P. R. M., Mainzer, A. K., et al. 2010, AJ, 140, 1868