THE 2008 OUTBURST OF EX Lup—SILICATE CRYSTALS IN MOTION

The optical brightness variations during the course of the outburst were monitored by amateur and professional astronomers. We collected visual brightness estimations from Albert Jones, the AAVSO⁷ database, and V-band CCD observations from the All Sky Automated Survey (ASAS)-3 database (Pojmanski 2002). For the visual magnitude estimation we used an uncertainty of 0.3 mag, which was the average difference between EX Lup and the comparison stars. From the ASAS-3 database we used only the highest quality measurements (quality flag "A"), for which the uncertainty of the measured V-band magnitude was less than 0.1 mag. The optical light curve is presented in Figure 1.

Zoom In Zoom Out Reset image size

We observed EX Lup with the Gamma-Ray Burst Optical and Near-Infrared Detector (GROND), which is a seven-channel simultaneous optical–infrared imager at the ESO 2.2 m telescope, on 2008 April 20. Optical images, covering a field of view of 54 × 54, were obtained in the g', r', i', and z' bands. The total observing sequenced consisted of four dithered exposures with 1 s integration time. At near-infrared wavelengths (JHK_S) images were obtained in 4 (JH) or 24 (K_S) dither positions, with a total of 40 s exposure time. The field of view was 10' × 10'. The optical images were reduced with standard IRAF tasks, and the final magnitudes were determined with aperture photometry and adopting the nominal zero point of each optical band from Greiner et al. (2008). Near-IR magnitudes were determined by averaging the photometric zero points derived for several stars in the field using the Two Micron All Sky Survey (2MASS) catalog. Because the near-infrared images were defocused, a large aperture was adopted in order to include all fluxes. The estimated photometric accuracy is 2.5%, 2.5%, and 5% in the J, H, and K_S bands, respectively, while in the optical the uncertainty is significantly larger (15%) due to the lack of dedicated calibration measurements.

EX Lup was observed with the Wide Field Imager (WFI) on the 2.2 m telescope in La Silla on 2008 April 20, with about one hour difference with respect to the GROND observations. The field of view of WFI spans 33' × 34' with small holes between the eight CCD detectors. EX Lup was observed using standard Johnson UBV filters. A total of three exposures were taken with each filter, with individual exposure times of 4 s for U and 0.5 s for BV. The data were reduced using standard IRAF procedures.

The photometry was done as standard aperture photometry with a 15 pixel aperture. The calibration and aperture corrections were obtained from observations of the Landolt fields SA107 and SA110 (Landolt 1992) observed at six different airmasses, resulting in zero-point errors <0.01 mag.

Optical photometries corrected for emission line contamination are also listed in Table 1. Correction was done by convolving optical spectra (A. Sicilia-Aguilar et al. 2012, in preparation) with the WFI and the GROND filter curves and calculating the ratio of the total observed emission over the continuum. The emission line corrected photometries were used during the SED fitting, since our model did not include optical emission lines.

EX Lup was observed on 2008 April 20 using the Son of ISAAC (SOFI) at the ESO 3.5 m New Technology Telescope (NTT). We obtained near-infrared images using the narrowband filters NB1.215, NB1.71, and NB2.195. The observations were performed using the template img_obs_AutoJitter, taking five offsets, and repeating the exposure several times at each jitter position. The high brightness of the target prevents to use common broadband filters. Individual frames were taken using the minimum possible exposure time of 1.2 s. The total exposure time per filter was 30 s in all cases. The data were reduced using the SOFI pipeline provided by ESO together with the Gasgano file manager. Photometric calibration was obtained by comparison with field-of-view stars included in the 2MASS catalog. The zero points were derived using between six and eight stars with accurate 2MASS magnitudes; the standard deviation was about 0.03 mag in all cases.

In addition, SOFI was used to obtain a spectrum using the high-resolution grism covering the K band with a slit width of 1''. A total exposure time of 560 s was dedicated to the EX Lup spectrum, split into eight individual frames of 70 s each. A close solar-type star (HIP78466) was observed to correct for telluric absorption. The data were also reduced using the SOFI pipeline.

The SOFI and the GROND near-infrared photometries agree with each other within the uncertainties in the H and K bands. In the J band the difference is about 3.8σ, which indicates a real difference, though with low significance. Such a difference can be caused by short-timescale (several hours) low-amplitude intrinsic photometric variability.

The total spectral line contamination in EXors is usually weaker at near-infrared wavelengths compared to the optical (see also, e.g., Lorenzetti et al. 2007). To verify this in the case of EX Lup we used the SOFI K-band spectrum for the NB2.195 filter. Since for the remaining two filters we did not have simultaneous spectroscopic data, we used the spectra from Kóspál et al. (2011) for the NB1.71 and NB1.215 filters, respectively. Although the data in Kóspál et al. (2011) were taken three months later than our photometry, the optical brightness of the source was about the same at both epochs and we assumed that the strength of the line emission was also similar on these dates. Since the observed line contamination was less than 1% in all three filters, which is significantly less than the photometric accuracy, we did not correct the photometries for line contamination.

EX Lup was observed with the Infrared Spectrograph (IRS; Houck et al. 2004) on board the Spitzer Space Telescope at four epochs: 2005 March 18 (PID: 3716), 2008 April 21 (PID: 477), 2008 October 10 (PID: 524), and 2009 April 7 (PID: 524). The data of the first two epochs were published in Ábrahám et al. (2009) and we refer to that paper for the details of the data reduction procedure of those spectra.

At the last two epochs EX Lup was measured using only the low-resolution modules (R = 60–120). Our spectra are based on the droopres products processed through the S15.3.0 version of the Spitzer data pipeline. First, the background 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 from an 8 pixel perimeter surrounding the flagged pixel. The spectra were extracted using a 6.0 pixel and 5.0 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 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, we matched orders based on the point-spread function of the IRS instrument, correcting for possible flux losses.

We obtained photometry at 70 μm and a low-resolution spectrum of EX Lup, using the Multiband Imaging Photometer for Spitzer (MIPS; Rieke et al. 2004). The 70 μm imaging was performed in photometry mode. The data reduction was started with the Basic Calibrated Data (BCD) products generated by the pipeline (version S17.2) developed at the Spitzer Science Center. We performed column spatial filtering and time median filtering on the BCD images following Gordon et al. (2007). The filtered BCD data were co-added and corrected for array distortions with the SSC MOPEX (MOsaicking and Point source Extraction; Makovoz & Marleau 2005) software. Outlier pixels were rejected using a 3σ clipping threshold. The output mosaic was resampled to 4'' pixel⁻¹.

We used a modified version of the IDLPHOT routines to detect the source on the final image and to obtain aperture photometry with an aperture radius of 16'' and with sky annulus between 18'' and 39''. The aperture correction, appropriate for a 60 K blackbody, was taken from Gordon et al. (2007). The final uncertainties were computed by quadratically adding the internal error and the absolute calibration uncertainty for which we adopted 7% (see the MIPS Data Handbook).

The low-resolution far-IR (55–95 μm; R ∼ 15–25) spectrum of EX Lup was obtained in the SED mode of MIPS. Six observing cycles were performed, each consisting of 10 s long on- and off-source exposures. The on-source and the background positions were separated by 1''. We began our data reduction with the BCD images (pipeline version S17.2) and MOPEX was used to perform further processing steps (combination of data, background removal, application of the dispersion solution). The spectrum was extracted in a 5 pixel wide aperture. The aperture correction factors were taken from Lu et al. (2008).

EX Lup was observed with the Mid-Infrared Interferometric Instrument (MIDI; Leinert et al. 2003), the N-band instrument on the Very Large Telescope Interferometer (VLTI; Glindemann et al. 2000). Spectrally resolved interferometric visibilities were obtained between 8 and 13 μm with a spectral resolution of R ≈ 35. Observations were performed at two epochs, during the nights starting 2008 June 21 (baseline U2U3, projected length 37.4 m, position angle 576 E of N) and 2008 July 16 (baseline U1U4, projected length 121.1 m, position angle 734 E of N). Data were obtained in "high sense" mode, i.e., each interferometric observation in which the beams coming from both telescopes are combined coherently is followed by photometric observations in which first only the light coming from telescope A, then from telescope B, is measured, and the intensities of the individual beams are determined.

Calibration stars were observed to determine the system response and atmospheric transparency. During the night of 2008 June 21 the calibrator HD 149447 (K6III, θ = 4.68 ± 0.05 mas; Cohen et al. 1999) was observed immediately following the EX Lup observation. During the night of 2008 July 16 the same calibrator was observed immediately before and after the EX Lup observations. We used the EWS package to reduce the data, and verified that reductions using the MIDI Interactive Analysis yield consistent results.

Low-resolution (R ≈ 250) N-band spectroscopy of EX Lup was performed with the VLT Imager and Spectrometer for mid-Infrared (VISIR) at the Very Large Telescope (VLT). We observed the source during two epochs, first around 2008 July 25 (for which data from the nights of the July 24, 25, and 27 were combined), and second in the night of 2008 August 28. Standard chopping and nodding techniques were applied to correct for the high celestial and instrumental background inherent to mid-IR observations.

During the July observations, HD 149447 was observed as telluric and flux calibrator immediately following each observation of EX Lup, at essentially identical airmass. During the observations in August we used HD 196171. Both calibrators are taken from the catalog of Cohen et al. (1999) that provides spectral types and angular diameters for each star, which we used to compute Kurucz model spectra to calibrate the absolute flux levels.

EX Lup was observed at 870 μm using the 295-element Large Apex Bolometer Camera (LABOCA; Siringo et al. 2009) on the 12 m Atacama Pathfinder Experiment (APEX) radio telescope (Güsten et al. 2006). The observation was performed on 2008 April 21 in continuous integration mode using a spiral pattern. This observing mode provides a fully sampled map of the total field of view of LABOCA. The total on-source observing time was 2.2 hr. The data processing was performed by the Bolometer array Analysis Software (BoA) package following the data reduction steps (correction for the atmospheric opacity, flux calibration, removal of correlated noise using an iterative method, and removal of spikes) outlined in the BoA User Manual (version 3.1). For the flux calibration we used the sources IRAS16293, B13134, and G34.3. The co-added final map had pixels with sizes of 6''. Aperture photometry was used to extract flux from the map using the position of EX Lup.

The brightest ever recorded outburst of EX Lup was announced in 2008 February (Jones 2008). The light curve, composed from visual estimations and V-band observations, is presented in Figure 1. Unfortunately the beginning of the outburst was not observed; one can only extrapolate for the starting date. The last optical measurement (ASAS) in the quiescent phase was on 2007 October 10 (∼110 days before the peak brightness), when brightness was V = 12.7 mag. We estimated the total length of the outburst to be about 300 days, which is an average value within the class.

The brightness of EX Lup at the peak of the outburst was about V = 8 mag. Taking the quiescent brightness to be V ≈ 13 mag (Sipos et al. 2009) one gets a V-band amplitude of 5 mag for the outburst. With this value, the 2008 outburst of EX Lup belongs to the strongest ones observed in EXors. This is an unusually high value compared also to the historical record of EX Lup. The only comparable outburst was observed in the 1950s.

The shape of the light curve during the 2008 outburst of EX Lup showed some peculiar features. The observations began about three weeks before the peak of the outburst. The brightness of the source rose steeply until it achieved a peak brightness of V = 8 mag. The star remained this bright only for a few days, then it faded to about 9.5 mag. Then the average brightness decreased by another 2 mag in the next 200 days. One of the most interesting properties of the light curve was the quasi-periodic brightness variation seen on top of the plateau-like slow fading. The amplitude of these fluctuations increased with time, while their period remained roughly constant, about 35 days.

The SEDs of EX Lup on 2008 April 21 and in the quiescent phase are compared in Figure 2. The outburst-to-quiescent flux ratio was the highest in the U-band and it decreased gradually with wavelength. Although EX Lup was already fainter on 2008 April 21 than at the peak of the outburst, at optical wavelengths it was still more than an order of magnitude brighter in the outburst compared to the quiescent phase. The shape of the SED in the optical domain is consistent with blackbody emission at ∼6500 K, which is far higher than the effective temperature of the star in quiescent phase (3800 K; Gras-Velázquez & Ray 2005). Although the absolute flux levels in the GROND i and z filters are somewhat higher than what would be expected for the aforementioned blackbody emission, the i − z color fits well to the modeled emission. The discrepancy in the absolute flux level can be explained by the relatively high calibration uncertainty in the GROND optical bands. In the near-infrared domain the flux increased only by a factor of about 2.5 in the outburst compared to the quiescent phase. Moving to longer wavelengths the ratio of the outburst-to-quiescent fluxes increased again to about an order of magnitude at 5.5 μm, which is the shortest wavelength of the Spitzer IRS spectra. Longward of 5.5 μm the ratio of outburst-to-quiescent fluxes decreased with wavelength.

Zoom In Zoom Out Reset image size

The brightening was certainly caused by the increased accretion rate onto the star during the outburst, indicated by the increasing luminosity in the accretion indicator lines. In Figure 3, we compare the spectral region around the Br γ line in the outburst and in the quiescent phase from Sipos et al. (2009). It can be seen that the luminosity of the Br γ was far higher in the outburst than in the quiescent phase. Using the relationship of Muzerolle et al. (1998) we derived an accretion rate of 2 × 10⁻⁷ M_☉ yr⁻¹ in the outburst.⁸ This is about three orders of magnitude higher than the value in the quiescent phase found by Sipos et al. (2009). We want to note that different lines were used to derive the accretion rate in the outburst (Br γ) and in the quiescent (Pa β), since none of the accretion indicator lines were measured both in the quiescent and in the outburst. The usage of different accretion indicators might introduce some additional uncertainty when a comparison of the accretion in the outburst and in the quiescent is done. This uncertainty, however, does not affect our conclusion that the 2008 brightening of EX Lup was certainly caused by the increased accretion rate onto the star during the outburst, indicated by the increasing luminosity in the accretion indicator lines.

Zoom In Zoom Out Reset image size

In this section we will describe the mid-infrared spectra measured toward EX Lup in our observing campaign. In order to avoid confusion, we will use the following naming convention. The ensemble of features between 8 and 13 μm will be called the "10 μm feature complex." The "10 μm feature" term will be used for the forsterite band peaking at 9.9 μm. The Spitzer IRS spectra of EX Lup are presented in Figure 4. In order to show also the weaker dust features, a spline continuum was fitted to and then subtracted from each spectrum. The continuum subtracted spectra are shown in Figure 5.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The pre-outburst spectrum from 2005 March was already analyzed by Sipos et al. (2009). They showed that the dust in the disk atmosphere of EX Lup is dominated by small amorphous grains with a mass-weighted average grain size of 0.5 μm. They also showed that the abundance of crystalline silicates in the disk atmosphere is negligible.

The Spitzer IRS spectrum from the outburst, taken on 2008 April 21, showed remarkable changes. First of all, the absolute flux level increased by about an order of magnitude at the shortest wavelength of the spectra (5.5 μm) and by about a factor of five at the longest wavelength (∼35 μm). The amount of brightening decreased with wavelength shortward of about 20 μm, but it was completely wavelength independent longward of 20 μm. As reported by Ábrahám et al. (2009) the spectrum showed a small but narrow and sharp feature at 10 μm on top of the 10 μm feature complex (see Figure 5). A shoulder at 11.2 μm and another feature at 16 μm also appeared in the outburst spectrum which were not present in the pre-outburst spectrum. Ábrahám et al. (2009) identified the source of these new bands as crystalline forsterite. It was also reported that no crystalline feature was seen in the outburst spectrum longward of 20 μm.

The third Spitzer IRS spectrum, taken on October 10 after the end of the outburst, showed further changes. As can be seen in Figure 4, the absolute flux level decreased compared to the outburst spectrum. The flux level was still higher than in the quiescent spectrum by about a factor of 1.7, but this factor was wavelength independent in the whole Spitzer IRS wavelength domain. The shape of the spectral features between 2008 April and 2008 October also changed. In Figure 5 it can be seen that the 10 μm feature complex became broader and the 11.2 μm shoulder became stronger, relative to the underlying broad amorphous feature, in October compared to the April spectrum. The most remarkable changes occurred, however, at longer wavelengths. Apart from the narrow 16 μm band of forsterite, new bands appeared in the 2008 October spectrum at approximately 23 μm, 28 μm, and 33 μm. These wavelengths correspond to the peak wavelengths of the strongest forsterite bands in this wavelength interval (see, e.g., Koike et al. 2003 for a list of forsterite band positions). Thus, we identified the source of these bands to be crystalline forsterite.

The last Spitzer IRS spectrum was taken on 2009 April 7. The absolute flux level decreased further compared to the first post-outburst spectrum from about five months earlier. Longward of 20 μm the flux decrease was independent of the wavelength. The flux variation showed, however, clear wavelength dependence shortward of approximately 20 μm, where the brightness decreased the most. The flux at 5.5 μm decreased to about half of the pre-outburst value, while longward of 20 μm the flux in the last IRS spectrum was approximately 30% above the pre-outburst value. The shape of the spectral features longward of 13 μm remained unchanged within the uncertainties (see Figure 5), while the 10 μm silicate feature complex clearly changed compared to the previous epoch. It can be seen in Figure 5 that the 10 μm forsterite feature on top of the broad amorphous feature was still present, but the shoulder at 11.2 μm became much weaker in 2009 April than it was in 2008 October. The general shape of the 10 μm feature seems to be somewhere between the outburst and pre-outburst spectra.

Due to the difference in the noise levels between the outburst and the post-outburst spectra, we cannot exclude the possibility from the Spitzer IRS spectra alone that weak crystalline features (of the same strength as in the post-outburst spectra) longward of 20 μm were already present in the outburst spectrum. In the post-outburst spectra the peak of the 24 and 28 μm features is 9 and 6 times the corresponding flux uncertainty above the underlying continuum, respectively. The absolute noise level of the outburst spectrum is, however, worse than the post-outburst spectra by a factor of about 3.4. If we measured the post-outburst spectra with the same noise level as we had in the outburst spectrum, the detection of the 24 and 28 μm features would be about 2.6 and 1.8 times the corresponding flux uncertainty, respectively, i.e., they could be hidden by the noise. We will further discuss this question in Section 4.2.1

In the Spitzer IRS spectra, the disk of EX Lup remains unresolved allowing us to investigate the integrated spectrum of the disk as a whole only. To obtain information about the spatial distribution of the newly formed crystals we used MIDI on the VLTI. In our first measurement in the U2U3 configuration (2008 June 22) with a projected baseline length of 37.4 m EX Lup remained spatially unresolved. The second measurement in the U1U4 configuration on 2008 July 16 (projected baseline length 121.2 m), however, resolved the source. In the following we will use the naming convention of van Boekel et al. (2004) and we will refer to the correlated spectrum as "inner disk spectrum." We also subtracted the correlated spectrum from the total spectrum and the resulting uncorrelated emission will be called the "outer disk spectrum." In Figure 6 (left) we present the inner disk spectrum of EX Lup. The spectrum peaks at 11.3 μm which corresponds to the wavelength of the strongest forsterite peak in the 8–13 μm region. The inner disk spectrum of EX Lup resembles that of HD 142527 which is known to originate in regions with ≳ 90% crystallinity (van Boekel et al. 2004), suggesting also very high crystallinity in the EX Lup inner disk spectrum. The outer disk spectrum of EX Lup is shown in Figure 6 (right). It is broader compared to the pre-outburst total spectrum and peaks at about 10.2 μm. The feature is smooth and lacks any substructure, which would be the signature of crystalline silicates. Interestingly, the outer disk spectrum of EX Lup is very similar to the 10 μm feature complex of FU Ori (Quanz et al. 2007). Quanz et al. (2007) have shown that the 10 μm feature of FU Ori is known to emerge from a region where significant grain growth has taken place, indicating also the presence of large (few μm) grains in the outer disk atmosphere of EX Lup.

Zoom In Zoom Out Reset image size

Two spatially unresolved spectra were taken with the VISIR instrument on VLT on 2008 July 25 and on 2008 August 28 (see Figure 7). Although the VISIR spectra have better signal-to-noise ratio than those taken by the MIDI instrument, the quality of the spectra is not as good as that of Spitzer IRS, allowing only moderate analysis of the spectral features. In 2008 July, at the epoch of the first VISIR measurement the source was at about the same brightness level as in 2008 April, when the outburst Spitzer IRS spectra was taken. The absolute flux level of the two spectra agrees well with each other, but the peak-to-continuum ratio is lower in the VISIR spectrum than in the Spitzer IRS (1.6 versus 1.4). The 10 μm feature complex seems to be broader in the VISIR spectrum too. The broader and shallower features are usually interpreted as indication for the presence of larger, micron-sized grains. Similar to the MIDI data the VISIR spectrum points toward the presence of large grains too. The 2008 August VISIR spectrum was taken at about the same optical brightness level as the first Spitzer IRS post-spectrum. Due to the lower signal-to-noise ratio of the VISIR spectrum no significant difference can be found in the shape or the strength of the Spitzer and VISIR spectra above the noise level.

Zoom In Zoom Out Reset image size

Ábrahám et al. (2009) reported the in situ formation of forsterite crystals via annealing in the disk atmosphere of EX Lup during its 2008 outburst. The crystalline features we found in the 20–35 μm wavelength interval are also associated with crystalline forsterite, supporting the conclusion of Ábrahám et al. (2009). It is interesting that no features of other crystalline dust species, apart from forsterite, can be seen in the mid-infrared spectra. For instance, enstatite (the magnesium-end member of the pyroxene solid solution series), which is frequently observed in protoplanetary disks (see, e.g., van Boekel et al. 2005; Watson et al. 2009; Juhász et al. 2010), is not present in the spectra of EX Lup. Laboratory annealing experiments of amorphous silicates show that forsterite will form independently of the starting stoichiometry if the grains are small (<1 μm) and porous (Fabian et al. 2000). For larger (few μm) and compact particles, annealing of amorphous grains with pyroxene stoichiometry results in the formation of enstatite. Using the activation energies from Fabian et al. (2000) the annealing time of enstatite in the 1000–1300 K temperature interval is about an order of magnitude longer than that of forsterite, but it decreases rapidly with increasing temperature. At 1200 K the calculated annealing time of enstatite is on the order of a minute. It is unlikely that dust grains in the outburst are heated up to 1200 K then cooled down below 1000 K on a one-minute timescale to prevent the formation of enstatite. The only possibility for the formation of forsterite without enstatite production is that the parent amorphous grains were sub-micron-sized porous grains. Indeed, the analysis of the Spitzer IRS spectrum of EX Lup from the quiescent phase show the presence of amorphous grains with a mass-averaged grain size of 0.5 μm (Sipos et al. 2009).

All crystalline silicate features, so far identified in protoplanetary disks, were associated with iron-poor, magnesium-rich crystals (see, e.g., van Boekel et al. 2005; Bouwman et al. 2008; Sargent et al. 2009; or Henning & Meeus 2011 for a review). An important consequence of the absence of iron is the dramatic decrease of the optical and near-infrared opacity of the grains. Iron-poor silicates are, therefore, cooler than ferromagnesian silicates at the same distance from the central star. Comparing the positions of the observed forsterite features in the spectrum of EX Lup with those of synthetic and natural olivines from Koike et al. (2003), we estimated the iron content of the lattice to be less than 10% with respect to magnesium. Such iron-poor crystals can be formed from iron-poor amorphous silicates, but also from ferromagnesian silicates if the annealing occurred under reducing conditions. Laboratory experiments showed that in the latter case the reduced iron will remain in thermal contact with the resulting crystals in the form of metallic inclusions (Davoisne et al. 2006). Such iron inclusions do not affect the positions of the forsterite bands while the grain as a whole will have higher optical to near-infrared opacity than that of isolated, pure forsterite crystals (see, e.g., Ossenkopf et al. 1992).

Besides iron, amorphous carbon is the other important opacity source in protoplanetary disks. Carbon, as well as iron, has a featureless opacity curve and high optical to near-infrared opacity (Jaeger et al. 1998). Sipos et al. (2009) used amorphous carbon (20% in terms of the total dust mass) in their radiative transfer (RT) modeling of the quiescent-phase SED of EX Lup. Amorphous carbon, however, is not stable above 1000 K as it oxidizes quickly and forms CO (Gail 2001). Since annealing is expected to occur above 1000 K, we do not expect amorphous carbon to be in thermal contact with silicate crystals.

Based on the above constraints, we propose two scenarios to be tested for the formation of forsterite crystals in the disk of EX Lup. In scenario A, the parent grains were porous ferromagnesian silicates in thermal contact with carbon. As the outburst set in and the temperature rose above 1000 K, the silicate grains crystallized. The formation of CO from amorphous carbon lowered the partial pressure of oxygen providing reducing conditions under which iron-poor crystals can form from ferromagnesian silicates. To explain the observed variations of the strength of forsterite features, an additional mechanism, e.g., radial/vertical mixing, is required.

In scenario B, the parent amorphous silicate grains were iron-poor sub-micron-sized porous grains in thermal contact with carbon. As the temperature rose above 1000 K at the beginning of the outburst, the silicate grains were annealed to crystalline forsterite, while the amorphous carbon oxidized and formed CO. Due to the low opacity of the crystals (in the absence of carbon) their temperature decreased. This resulted in the weakening of the forsterite features in the 10 μm region while in the same time the features in the 20–35 μm region became stronger, which might explain the observed spectral variations.

We modeled the SED, the mid-infrared spectra, and the interferometric visibilities simultaneously using the two-dimensional RT code RADMC (Dullemond & Dominik 2004a) and a (1+1)-dimensional turbulent mixing code (Dullemond & Dominik 2004b). The goal of our modeling was to find answers to the following three questions. Are the freshly produced crystals isolated iron-free particles or are they in thermal contact with iron? Which physical process was responsible for the observed variation of the mid-infrared solid-state features? What kind of constraints can we put on the mechanism of the outburst? We calculated, therefore, seven series of models to simulate the two crystal formation scenarios (series A-, and B-), with (-AD-) and without (-NAD-) an optically thick accretion disk and with different strength of the vertical turbulent mixing (-α0.01, -α0.1, -α1.0). Details of the modeling procedure are described in the Appendix. We now briefly describe the modeling strategy in the different model series.

• 1.  
  Series A-AD-α1.0, A-AD-α0.1, A-AD-α0.01. First, we started with a disk containing only amorphous silicates and amorphous carbon. The parameters of the starting disk model were taken from Sipos et al. (2009) and are summarized in Table 2. The inner radius of the disk was taken to be 0.3 AU (instead of 0.2 AU as in Sipos et al. 2009) to get a better agreement between the model and the quiescent-phase Spitzer IRS spectrum in the 5–10 μm wavelength interval. After the accretion rate has been set, we ran the RT code to obtain the temperature structure at the beginning of the outburst. Amorphous dust grains were replaced by forsterite grains with iron inclusions wherever the temperature in the disk exceeded 1000 K. Then the turbulent mixing code was run to simulate vertical mixing over a 10 day period. Turbulent mixing was treated as a diffusion process using the α prescription (Shakura & Sunyaev 1973; see also the Appendix) for the diffusion coefficient. The strength of the turbulence was controlled via the value of α, for which we used 0.01, 0.1, and 1.0. The last step was to update the density structure according to the results of the mixing code. The whole modeling cycle was then repeated. Accretion heating was simulated with an optically thick accretion disk as well as a hot spot on the stellar surface.
• 2.  
  Series A-NAD-α1.0, A-NAD-α0.1, A-NAD-α0.01. Same as the A-AD- models, but no accretion disk was present in the models. The total accretion luminosity was assumed to be radiated away from the hot spot.
• 3.  
  Series B-AD. We started the modeling in the same way as in the A-AD- models. Amorphous dust grains were replaced with isolated iron-free forsterite crystals wherever the temperature in the disk exceeded 1000 K. Amorphous carbon was also removed from the dust phase above 1000 K. In this model series we did not use the vertical mixing code, we only updated the accretion rate and calculated models on each 10th day.

For the amorphous grains we applied the dust model of Sipos et al. (2009), which was used to model the quiescent-phase SED. During the modeling we used two grain models for forsterite crystals. For isolated iron-free forsterite crystals we used mass absorption coefficients calculated from the optical constants using Distribution of Hollow Spheres (DHS; Min et al. 2003) scattering theory. We used the optical constants of Sogawa et al. (2006) longward of 2 μm and supplemented them with that of Scott & Duley (1996) shortward of 2 μm. For composite grains, where forsterite and iron are in thermal contact with each other, we calculated the optical constants of the aggregate from those of the constituents first, using effective medium theory (APMR mixing rule; Min et al. 2008). For the constituents of the aggregate we used the DHS polarizability. Finally, Mie theory was applied to calculate the mass absorption coefficient of the aggregate. The same optical constants were used for forsterite and for the isolated crystals, while for iron we used the optical constants of Ordal et al. (1988) between 0.66 μm and 285 μm and extended it with the optical constants of Pollack et al. (1994) shortward of 0.66 μm and longward of 285 μm.

In Figure 8, we compared the observed SED of EX Lup in the outburst to model calculations. Although our models reproduced the general shape of the observed SED, there were quantitative differences. The largest differences between our models and the observations occurred in the 1–8 μm wavelength interval (A-AD models) and in longward of 60 μm. The former is related to the presence or absence of an optically thick accretion disk inward of 0.3 AU. We will discuss this problem in Section 4.2.2.

Zoom In Zoom Out Reset image size

The most likely explanation why our models overestimated the observed fluxes longward of 60 μm is that the disk of EX Lup was not in an equilibrium state. RADMC, as well as other RT codes used to model protoplanetary disks, calculates the temperature structure of the disk in radiative equilibrium. The emission longward of 60 μm is dominated by the disk interior which is optically thick to the heating radiation of the star and the disk atmosphere above it. The high optical thickness prevents the disk interior to react to the variable irradiation luminosity quickly. Since the optical depth in the far-infrared to millimeter wavelengths is lower than at shorter wavelengths the observer can see deeper into the disk observing those regions that are still not accommodated to the changed irradiation luminosity. The optically thin, upper layers of the disk can, however, react to the variable irradiation luminosity almost immediately.

4.2.1. Radial and Vertical Mixing of Solids

The first question we tried to answer with our modeling was which of the two crystal formation scenarios, outlined in Section 4.1, was at work in the disk of EX Lup. In Figure 9, we compared the calculated mid-infrared spectra of models A-AD⁹ and B-AD. As can be seen in Figures 9(a)–(c), the positions of forsterite bands in the observed and modeled spectra matched very well in models A-AD. The strengths of the forsterite peaks are dependent, however, on the strength of the vertical mixing. The positions of the forsterite peaks in the B-AD models were shifted systematically toward shorter wavelengths compared to the observations (see Figure 9(d)). The shift of the peak position was certainly related to temperature/RT effects since in the mass absorption coefficients there was no difference in the peak positions of isolated forsterite crystals and forsterite–iron aggregates. Moreover, spectra of B-AD models showed also stronger forsterite features, especially in the post-outburst epochs, compared to the observed spectra. Thus, our modeling suggests that the freshly produced forsterite crystals had to be in thermal contact with iron, which excludes scenario B as a crystal formation mechanism.

Zoom In Zoom Out Reset image size

The next question we addressed was whether vertical mixing can explain the observed variation in the shape of the mid-infrared solid-state emission features. Without considering vertical mixing, the apparent crystallinity in spatially unresolved spectra depends on the ratio between the size of the crystallized region and the total (crystalline + amorphous) silicate emission zone seen at 10 μm. The former depends on the irradiation luminosity at the peak of the outburst and does not change with time. The total silicate emission zone, however, scales linearly with the actual irradiation luminosity (see, e.g., Kessler-Silacci et al. 2007). If the irradiation decreases, the total silicate emission would shrink accordingly. Since the size of the crystalline region remains unchanged, the apparent crystallinity increases. So far we assumed that in the crystalline region the value of crystallinity is 1.0, i.e., every dust grain is in crystalline form and this does not change with time. The effect of vertical mixing is to decrease the crystallinity with time in the crystalline zone by stirring up amorphous grains from deeper layers of the disk and, at the same time, mixing down crystals. In this way vertical mixing might compensate the effect of the shrinking total silicate emission zone.

In Figures 9(a)–(c), we compared the observed spectra in the 7.5–13.5 μm region with those of models A-AD with different values of α. In our weakest mixing case (α = 0.01, Figure 9(a)), our models predicted stronger forsterite features than observed. In contrast, models with the strongest possible mixing (α = 1.0, Figure 9(a)), predicted too weak crystalline features compared to our observations. The best agreement between our models and the observed spectral features in the 7.5–13.5 μm was achieved using α = 0.1. The strength of the 10 μm feature complex is well reproduced for the outburst (2008 April) and the second post-outburst spectra (2009 April), but our model overestimates the feature strength compared to the spectrum taken in 2008 October. The disagreement between the model and the observed spectra is about 35%. Since the emission at the 10 μm is optically thin, it reacts to the changing irradiation immediately and the infrared radiation is proportional to the irradiation flux. This means that the 35% difference between our models and the observed flux level translates to a difference of ∼0.3 mag difference in the optical. Since we had no optical photometry at the time when the Spitzer IRS spectrum was taken, only a few days after, the difference could be caused by the extrapolation to the optical brightness used in our modeling. In Figures 10(a)–(c), we compared our models to the observed spectra in the 13.5–36 μm wavelength interval. Although our models with α = 0.1 matched the observed spectra in the 7.5–13.5 μm region well, the observed forsterite features longward of 20 μm are far stronger than in any of our models independently of the value of α (see Figure 10).

Zoom In Zoom Out Reset image size

Two explanations could be used to solve the discrepancy between the observed and modeled crystalline features longward of 20 μm. A solution could be that, in reality, the size of the crystallized region was higher than in our models. This can be due to a higher peak luminosity of EX Lup in the outburst or due to a lower crystallization temperature. In this case, the outburst spectrum could have contained some weak forsterite emission features longward of 20 μm, which were hidden by the noise in the measured spectra. The uncorrelated MIDI spectrum, however, excludes this possibility. We tested this scenario with our two-dimensional RT model. If the crystallized region is larger than about 1.5–1.8 AU, the calculated uncorrelated spectra of our models show strong forsterite bands, which are not observed in the MIDI spectra. Moreover, models which reproduced the strength of the forsterite features longward of 20 μm in the spatially unresolved post-outburst spectra overestimated the strength of the crystalline features in the 10 μm region. Using our modeling framework with vertical mixing, we were not able to reproduce the strength of the crystalline features both shortward and longward of 20 μm in the post-outburst spectra, and the spatially resolved mid-infrared interferometry simultaneously.

The other explanation would be that forsterite crystals were not present in the disk atmosphere outward of 1.5–1.8 AU in 2008 July, but they were transported to a few AU by 2008 October. The transport of crystals over a 2 AU distance within three months implies an average radial velocity of about 38 km s⁻¹. Radial mixing is far too slow for such a transport velocity. Instead, it is likely that in this scenario the crystals are transported by, e.g., a wind-driven transport mechanism. Indeed, optical and near-infrared spectroscopic observations suggest the presence of a wind during the outburst of EX Lup (Aspin et al. 2010; Goto et al. 2011; Kóspál et al. 2011; A. Sicilia-Aguilar et al. 2012, in preparation). The quantitative modeling of this scenario is, however, beyond the scope of this paper.

4.2.2. Outburst Mechanism

EXor outbursts are usually considered to be a "down-scaled" version of FUor outbursts, suggesting the same triggering mechanism of the outbursts in the two classes (Hartmann & Kenyon 1996). The most widely accepted model of FUor outbursts is the self-regulated thermal instability (TI; Bell & Lin 1994). Although being the most widely accepted model, even TI cannot explain all observational properties of FUors (see, e.g., Zhu et al. 2007). In the case of TI, the released viscous energy in the disk midplane is trapped in the disk above a certain threshold surface density value. The radiative trapping is caused by the opacity of H⁻ ions in the partially ionized disk regions. The effective viscosity as well as the mass transfer rate in the disk within the ionized radius increases. The ionization front propagates outward in the disk up to a certain radius then retreats again. The gas inward of the ionization front is accreted onto the central star.

In the following we will investigate the TI model in the context of our observations. TI requires an optically thick accretion disk (due to the presence of H⁻ ions) which dominates the SED of FU-Ori-type objects at optical and near infrared wavelengths (see, e.g., Zhu et al. 2007). The comparison of our A-AD and A-NAD models to the observed SED of EX Lup showed that models without an optically thick accretion disk fit the SED better in the 1–8 μm wavelength range (Figure 8). In order to quantify the goodness of the SED fit we calculated and compared the χ² of the fit for all model series. The number of tuned parameters were only two, the value of α (except for B-AD models) and whether or not an optically thick accretion disk was present within the inner hole of the dust disk. The number of fitted data points were 1435. The χ² values for the B-AD model was 3409. The A-AD-α0.01, A-AD-α0.1, and A-AD-α1.0 models resulted in a χ² of 1423, 1372, and1290, respectively. Finally, the A-NAD-α0.01, A-NAD-α0.1, and A-NAD-α1.0 models provided a χ² of 533, 470, and 492, respectively. It can be easily seen that the goodness of fit is very sensitive to the presence (or rather to the lack) of an optically thick accretion disk, while only marginal dependence is visible on the value of α. Independent of the value of the viscosity parameter, the best-fit model family is the A-NAD, while the absolute best-fit model is A-NAD-α0.1. Mid-infrared visibilities of models without an optically thick accretion disk gave also better agreement with the observations (see Figure 11). This suggests that there was either no optically thick accretion disk present during the outburst of EX Lup or the emission at optical wavelengths came from a narrow ring of hot gas in the inner edge of the gas disk.

Zoom In Zoom Out Reset image size

One of the main differences between FUor and EXor outbursts is the timescale over which the outburst occurs. The TI model alone, without any additional perturbation, is not able to reproduce the rise time in the case of FU Ori and of V1515 Cygni (Clarke et al. 2005). In these cases the brightening happened over about a year, which is an order of magnitude longer than in the case of EX Lup (few weeks). The relatively short duration time of the outburst of EX Lup puts constraint also on the outer radius, from where gas can be accreted onto the star. The mass transfer in the disk occurs by viscous diffusion, the timescale of which can be estimated using

$$\begin{equation} t_{\rm visc}\approx \frac{r^2}{D}\approx \frac{r^2}{H_p^2}\frac{1}{\alpha \Omega }, \end{equation} \tag{ 1 }$$

where D is the diffusion coefficient (D = αc_sH_p), r is the radius, H_p is the pressure scale height of the gas, and Ω is the Kepler frequency. We used Equation (1) to estimate the outer radius from which material can be accreted within 10 months. We assumed that r/H_p = 10 and for α we used the highest theoretically possible value, α = 1.0. The resulting radius is ≈0.12 AU, which is well within the fitted inner radius of the dust disk (0.2–0.35 AU; Sipos et al. 2009). This implies that the reservoir of material for the 2008 outburst should have been present in the inner hole of the dust disk. This calculation implies that the hole in the disk, observed in the quiescent-phase SED, should not be empty, but that it is filled up with gas, which is optically thin in the continuum. On the other hand, if thermal instability takes place in the disk, the disk should be optically thick due to the presence of H⁻ ions, which should be visible in the SED.

The mass of gas accreted onto the star during this outburst was about 10⁻⁷ M_☉, and the calculation above implies that all of this gas should be located within 0.12 AU from the central star. If we assume that such a high-amplitude outburst, like that in the 1950s or in 2008, consumes all the available gas within 0.12 AU, a mass transfer rate of 2 × 10⁻⁹ M_☉ yr⁻¹ or higher is required to maintain the observed amplitude and frequency (one in every 50 years) of these large eruptions. The estimated mass transfer rate is a lower limit, as we neglected all the smaller amplitude outburst, which occurred since the 1950s, but it is more than an order of magnitude higher than what Sipos et al. (2009) estimated. The required mass transfer rate in the disk for TI to operate was found to be 5 × 10⁻⁷ M_☉ yr⁻¹ or higher (Bell & Lin 1994), which is far higher than the estimated accretion rate in the quiescent phase of EX Lup. If the outbursting region is fed with a lower rate, thermal instability does not occur, the accretion onto the star is stable and steady. The equilibrium curves calculated by Bell & Lin (1994) also suggests that between 5 × 10⁻⁷ M_☉ yr⁻¹ and about 10⁻⁵ M_☉ yr⁻¹ there is no stable solution for the accretion. If thermal instability sets in, the disk quickly flips to the outburst phase and the accretion rate increases to above 10⁻⁵ M_☉ yr⁻¹.

The mass transfer rate in the disk is, however, not necessarily the same as the accretion rate onto the star. The accretion rate onto the star can be lower than the mass transfer rate in the disk due to, e.g., interaction between the stellar magnetic field and the accreting gas (see, e.g., Romanova et al. 2002, 2004; D'Angelo & Spruit 2010), which causes the gas to pile up in the disk. The gas lines, used as accretion indicator, are formed in the funnel flow between the inner radius of the gas disk and the surface of the star (e.g., Muzerolle et al. 1998), and therefore they only measure the accretion rate onto the star. Thus, the mass transfer rate in the disk can be significant even if the accretion rate derived from, e.g., the Br γ line is zero. In FUor outbursts, all gas piled up over 10⁴ years (estimated time between two subsequent FUor outbursts) is accreted onto the star since the accretion onto the star is as high as 10⁻⁴ M_☉ yr⁻¹. In contrast, the accretion rate of EX Lup in the outburst was three orders of magnitude lower. If the mass transfer rate in the disk of EX Lup was 5 × 10⁻⁷ M_☉ yr⁻¹ (required for the TI to operate), EXor outbursts cannot accrete all piled up mass from the disk onto the star. In this case gas accumulates continuously in the inner regions of the disk and at some points the system should go into a FUor-type outburst, meaning that EXors are in fact FUors in the "quiescent" phase. Since EX Lup has a disk mass of 0.02 M_☉ and no evidence for a massive envelope has been found so far, the disk of EX Lup would be entirely consumed by only very few FUor eruptions (assuming $\dot{M}$ = 10⁻⁴ M_☉ yr⁻¹ and an outburst duration of 50 yr), which makes this scenario unlikely. Our study suggests that the 2008 outburst of EX Lup is caused by increased accretion onto the central star as indicated by the increase of the Br γ line luminosity in the outburst compared to the quiescent phase. On the basis of the observational data, we conclude, however, that it is unlikely that the self-regulated thermal instability is the responsible mechanism triggering the outburst.

The recent work of D'Angelo & Spruit (2010) could provide an alternative explanation for the outburst mechanism of EX Lup. If the magnetic truncation radius of a strongly magnetized star is close to the co-rotation radius the accretion onto the star will not be continuous anymore, but becomes episodic. After the end of an outburst, in the quiescent phase, the inner radius of the disk is located just outside of the co-rotation radius. The accretion onto the star is low and material piles up at the inner edge of the disk. The inner radius of the disk moves inward and when it crosses the co-rotation radius a burst of accretion occurs that empties the piled-up reservoir of material and moves the inner edge of the disk outside of the co-rotation radius again. The interesting property of this model is the shape of the accretion rate versus time curve (see e.g., Figure 4 middle panel in D'Angelo & Spruit 2010). The outburst starts with a steep rise of the accretion rate which is followed by a slow plateau-like fading with quasi-periodic oscillations on top of the slow decline. This is very similar to what is observed in the optical light curve of EX Lup, which is thought to be closely related to the accretion rate (i.e., hot spot on the stellar surface). Not only the shape of the light curve, but also the timescales predicted by these simulations are in qualitative agreement with the observations. For the parameters that are the closest to that of EX Lup (see Figure 5 in D'Angelo & Spruit 2010) the simulations of D'Angelo & Spruit (2010) predict that the duration of the outbursts is less than the viscous timescale at the co-rotation radius, while the time between two subsequent outbursts is about 10 times the viscous timescale. To compare the observations to the simulation of D'Angelo & Spruit (2010) we estimated the viscous timescale at the co-rotation radius for EX Lup. Sipos et al. (2009) measured the v · sin  (i) of EX Lup to be 4.4 ± 2 km s⁻¹. Using a stellar mass and radius of 0.6 M_☉ and 1.6 R_☉ (Gras-Velázquez & Ray 2005), respectively, and assuming 20° inclination and α = 0.1 the resulting viscous timescale is about three years. According to the models of D'Angelo & Spruit (2010) the duty cycle of the outbursts would be about 30 years while a single outburst would last less than three years. Although these numbers are close to the observed values, the agreement could only be considered qualitative, due to the large uncertainty in the value of α and the unknown inclination angle. Despite these uncertainties, the remarkable similarity of the observed light curve of EX Lup and the accretion rate curve predicted by D'Angelo & Spruit (2010) suggests that a further quantitative comparison between this model and observations of EXors should be done.