© ESO, 2013

1. Introduction

Circumstellar disks are a natural outcome of the star formation process and are known to dissipate over time (~10⁶–10⁷  yr; Haisch et al. 2001) by means of stellar winds or by photoevaporation caused by either the radiation of the central star or by an external source of radiation (e.g., Hollenbach et al. 2000; Clarke et al. 2001; Alexander & Armitage 2007), by accretion onto the central star (Hartmann et al. 1998), by grain growth and fragmentation (e.g., Dullemond & Dominik 2005; Dominik et al. 2007; Natta et al. 2007; Birnstiel et al. 2011; Sauter & Wolf 2011; Garaud et al. 2013; Ubach et al. 2012), and by the formation of planets (e.g., Pollack et al. 1996; Boss 2002; Armitage 2007, 2010; Fortier et al. 2013). The way in which the disk material dissipates has important implications for the possibility of planet formation.

During the last decades several theories to explain the formation of planets in circumstellar disks have been proposed, for example, the core accretion-gas capture scenario (Pollack et al. 1996; Papaloizou & Terquem 2006; Lissauer & Stevenson 2007). Sub-micron-sized dust particles, which are strongly coupled to the motion of the gas in the disk via gas drag, agglomerate through low-velocity collisions to eventually form planetesimals of several kilometers in size (e.g., Beckwith et al. 2000; Dominik et al. 2007; Natta et al. 2007). During this coagulation process, larger particles decouple from the turbulent gas motion, settle toward the disk midplane and radially drift toward the inner parts of the disk as a result of the impact of stellar gravity and gas drag (e.g., Weidenschilling 1977; Barrière-Fouchet et al. 2005; Fromang & Papaloizou 2006; D’Alessio et al. 2006). The growth beyond planetesimals occurs via direct collisions, with an increasing role for gravitational focusing as masses become larger, leading to the gravitational agglomeration of these bodies to rocky planets (Safronov & Zvjagina 1969). A phase of runaway growth occurs, followed by an oligarchic growth phase (Kokubo & Ida 1998; Thommes et al. 2003).

There are, however, many unresolved problems in the picture of this formation process. For example, it is uncertain how dust grains overcome the radial drift barrier (Weidenschilling 1977). This occurs as it is expected that meter-sized bodies migrate toward the star in a timescale that is shorter than the timescale for further growth which effectively stops further grain growth. Another example is the fragmentation barrier, where the boulders are destroyed at typical collision speeds, halting the dust growth at centimeter to meter sizes (Benz 2000; Blum & Wurm 2008; Brauer et al. 2008). Moreover, the dust coagulation process, the first stage of planet formation, must also overcome the charging barrier (Okuzumi 2009) and bouncing barrier (Zsom et al. 2010) which affect the later stages, too.

Besides the hypothesis that planetesimals may form from pairwise collisional growth of smaller bodies alone, there is the hypothesis that, when the dust particles have reached cm-scale as a consequence of this collisional growth, the planetesimals may form from the gravitational fragmentation of a dense particle sub-disk near the equatorial plane resulting from the dust settling. Protoplanetary disks may support regions within which turbulence acts to locally enhance the ratio of solids to gas. As a result, patches of the disk may become dense enough to trigger a gravitational instability and collapse into planetesimals (Safronov & Zvjagina 1969; Goldreich & Ward 1973; Armitage 2007; Chiang & Youdin 2010). This suggestion entirely bypasses the size scales that are most vulnerable to radial drift and would allow the radial drift and fragmentation barrier to be overcome.

With the current observational capacities it is not possible to directly observe large dust boulders (>1 m) and planetesimals, but signs and implications of their formation can be detected. The stratified structure resulting from grain growth and settling has a significant impact on observable quantities of the disk (e.g., Dullemond & Dominik 2004). The spectral energy distribution (SED) provides constraints on the dust mass and grain-size distribution. Scattered light images at near-infrared (NIR) and mid-infrared (MIR) wavelengths probe the surface layers of the optically thick disks and reveal properties of small dust grains through wavelength-dependent opacity (e.g., Watson et al. 2007). Continuum observations at submillimeter and millimeter (mm) wavelengths are mainly sensitive to the bulk of the disk mass closer to the midplane and play an important role in constraining disk structure properties, for example, the spatial dust distribution of larger grains (Andrews & Williams 2005, 2007). The analysis of scattered light or re-emission images or the SED only provides a limited view of a disk leading to ambiguities in the derived disk and dust properties (e.g., Chiang et al. 2001). To precisely study the physical processes like dust evolution and to strongly reduce model degeneracies, it is essential to combine resolved images at different wavelengths and a SED with a good wavelength coverage in a multi-wavelength observational and modeling approach.

In this paper we study the circumstellar disk of the prominent Butterfly Star, also known as IRAS 04302+2247, using this approach. This Class I young stellar object (YSO; e.g., Adams et al. 1987; Lada 1987) is located in the Taurus-Auriga molecular cloud complex at a distance of ~ 140 pc (Kenyon et al. 1994). It is surrounded by an edge-on seen flared circumstellar disk, so the central star is not visible directly (inclination i = 90° ± 3°; Padgett et al. 1999, 2001; Wolf et al. 2003). As shown by interferometric observations at submm and mm wavelengths, the disk exhibits a radius of about 300 AU (Wolf et al. 2003, 2008). Several attempts have been undertaken to model the structure and physical conditions in the circumstellar environment of this object (e.g., Lucas & Roche 1997; Wolf et al. 2003; Stark et al. 2006). In particular, Wolf et al. (2003) were the first to base their modeling on scattered light images and resolved mm maps and, therefore, were the first to consider wavelength regimes tracing different physical processes in different regions of the circumstellar environment. They concluded that the grains in the outer parts of the circumstellar environment are comparable to those from the interstellar medium (ISM), while dust grains have grown via coagulation by several orders of magnitude in the much denser parts of the disk. However, the observational data they presented did not allow them to constrain the spatial dependence of the dust grain properties in the disk interior.

We present new observations and a coherent multi-wavelength model for the circumstellar disk of the Butterfly Star that accounts for spatially resolved data sets ranging from the NIR to mm wavelengths as well as for the well-sampled SED of this object. These observations provide different, complementary detailed perspectives on the disk structure and dust properties. Here, coherent means that the model is capable of reproducing all considered observations obtained at different wavelengths at the same time. This new modeling became necessary as new observations are not properly reproduced by the model of Wolf et al. (2003, 2008). Furthermore, these new observations with much higher angular resolution promise to allow one to put stronger constraints on the disk structure and dust grain properties.

The observations and data reduction are introduced in Sect. 2. In Sect. 3, we characterize the data set that is the basis for our modeling. A detailed description of our model can be found in Sect. 4 and the entire modeling process is described in Sect. 5. The results of this study are presented in Sects. 6 and 7 and their implications are discussed in Sect. 8. Finally, Sect. 9 contains a summary and concluding remarks.

2. Observations and data reduction

The data set used in this study includes both spatially resolved images at mm, submm, and NIR wavelengths, as well as a well-sampled SED. In this section these observations and the data reduction procedures are discussed.

2.1. Millimeter observation

The Butterfly Star was observed with the IRAM Plateau de Bure interferometer (PdBI, Lucas 1991) in its two most extended configurations (A and B). The dual polarization receivers were tuned near 220 GHz (λ = 1364 μm). The wideband WIDEX correlator provided a total bandwidth of 4 GHz in each of the two polarizations. The covered frequency range included lines of the ¹³CO and C¹⁸O isotopologues, as well as H₂CO and SO. Spectral line results will be reported in a separate paper, but the continuum data discussed here was generated by avoiding line contamination.

Observations were performed on February 3 and 4, 2010 (Aq configuration, two baselines of 730–760 m, excellent rms (baseline-based) phase noise <35°), January 22 (Aq configuration), and February 10, 2011 (B configuration, up to 450 m baselines, rms phase noise <45°). Although the repeated configurations produced similar results, data from January 22, 2011 was not included in the final analysis because of its higher phase noise. With very good observing conditions (<1 mm water vapor, going up to 1.5 mm only for the last 2 h of February 10, 2011), single sideband system temperatures ranged between 100 to 140 K. Phase noise yields an effective seeing of approximately 0.15″.

Phase calibration was performed using the nearby quasars 0400+258 and 0507+179, and flux calibration is referred to MWC 349, with an assumed flux density of 1.7 Jy at this frequency, with an overall calibration accuracy of order 10%.

With natural weighting, the angular resolution of the data is 0.59″ × 0.35″ at a position angle (PA) of 31.1°. The shortest baseline is large, 79 m including projection effects, so that structures larger than about 4″ could be substantially resolved out. However, given the source elongation and orientation and the beam shape, deconvolution is efficient in recovering most structures.

The total flux density is about 101 ± 10 mJy. The expected thermal noise in the map is 0.18 mJy/beam, but because of phase noise, the dynamic range is limited and the effective noise level in the deconvolved map is ~0.4 mJy/beam.

2.2. Submillimeter observation

Observations of the Butterfly Star with the Submillimeter Array (SMA, Ho et al. 2004) were carried out on January 9, 2006 using the upper and lower sidebands at 330 GHz and 340 GHz, respectively (λ = 894 μm). The angular resolution of the data is 0.67″ × 0.53″ at PA =  −76.5°. The SMA observations and the data reduction were described in detail by Wolf et al. (2008).

2.3. Near-infrared observations

The NIR observations of the Butterfly Star were obtained with the Near Infrared Camera and Multi-Object Spectrometer (NICMOS, Thompson et al. 1998) on the Hubble Space Telescope (HST) on August 19, 1997. The NICMOS NIR Camera, NIC2, and the filters F110W (λ_c = 1.13 μm), F160W (λ_c = 1.60 μm), F187W (λ_c = 1.87 μm), and F205W (λ_c = 2.07 μm) were used. The data reduction and flux calibration were described in detail by Padgett et al. (1999). Figure 3 shows a composite image of the Butterfly Star at NIR wavelengths.

2.4. Spectroscopic data

The Butterfly Star was observed with the low-resolution modules of the Spitzer Space Telescopes (SST) InfraRed Spectrograph (IRS) on October 1, 2007 and March 1, 2004 for short-low (SL) and long-low (LL), respectively. The LL observation was obtained as part of the IRS guaranteed time observations (PI: Houck), while the deeper SL observations were obtained as part of the open time program 30765 (PI: Stapelfeldt). The total exposure time was 587 s for the SL modules and 12.6 s for the LL modules. The two-dimensional droop-corrected frames obtained using pipeline version S16.1.0 were co-added for each nod position, and the two nods were pairwise differenced to remove the background. One-dimensional spectra were extracted using apertures of 3 pixels for the SL modules, while 4 and 5 pixel apertures were used for the LL2 and LL1 modules, respectively. The one-dimensional spectra were flux-calibrated using the spectral response function included in the SST pipeline.

3. Basis for modeling

The basis of our modeling forms the (sub)mm and NIR maps, as well as the SED, including the SST/IRS spectrum. This data set is briefly characterized in this section.

3.1. Images

The interferometric continuum maps at 1.3 mm and 894 μm are shown in Figs. 1 and 2. The emission seen in these maps is thermal re-emission of the dust. While the source is unresolved vertically to the disk in both maps, it is well-resolved in radial direction. The circumstellar disk seen in both maps has a north/south extension of ~600 AU. While the brightness structure of the 1.3 mm map shows a clear maximum at the center, a local minimum is visible at this position in the 894 μm map. The nature of this local minimum is discussed in Sect. 8.1.

The emission seen in the NIR images is scattered light from the central star (see Fig. 3 and Fig. 1 in Wolf et al. 2003). Their appearance is dominated by a totally opaque linear feature that bisects the scattered light structure. This feature is interpreted as a circumstellar disk seen edge-on. Because of the large optical depth of this disk, no point source is detected at any of the observed NIR wavelengths. The dependence of the height of the disk on the wavelength is clearly visible. The bipolar scattered light nebula, which extends ~900 AU in the direction north/south, also shows a complex morphology with approximately equal brightness between the eastern and western lobes. As this appearance reminds one of a butterfly, the name Butterfly Star has been established (Lucas & Roche 1997). Comparing the (sub)mm maps with the NIR images it can be seen that the position and orientation of the elongated (sub)mm structure perfectly fits that of the dark lane in the NIR.

We rotate the observed images by the position angle of the image y-axis and the position angle of the Butterfly Star (PA = 175°, as simulated images are symmetric we use PA =  −5°) in order to align the major axis of the dust lane with the vertical axis, as it is for the simulated images.

3.2. Spectral energy distribution

The photometric measurements used in this modeling are presented in Table 1. In addition, the SST/IRS spectrum in the range from 12 μm to 38 μm is considered. In our disk modeling we do not include the scattered-light part of the SED below 12 μm as it is not an objective of this study to reproduce the highly-structured envelope seen in the NIR images (see Sect. 5.2).

Figure 4 shows the complete SED of the Butterfly Star consisting of the data shown in Table 1, the SST/IRS spectrum ranging from 5 μm to 38 μm, as well as HST/NICMOS, SST/IRAC, and WISE data points. The SED of our target is typical for a Class I YSO. It peaks close to 100 μm and shows a 10 μm absorption feature. This silicate feature and the steep slope of the IRS spectrum beyond 10 μm require a high inclination angle which is consistent with the spatially resolved NIR images. In Fig. 5 the fit to the (sub)mm slope of the SED of the Butterfly Star can be seen which yields the millimeter spectral index α_mm =  −dlog (F_λ)/dlog (λ) = 2.48 ± 0.07. This value is significantly smaller than that found for small ISM-sized dust grains (~3.7), indicating grain growth in the disk (see, e.g., Natta et al. 2007).

4. Model description

4.1. The disk

Our model for the Butterfly Star consists of a parameterized disk. We describe this disk with a radially and vertically dependent density distribution based on the work of Shakura & Sunyaev (1973) which can be written as (1)Here, z is the cylindrical coordinate with z = 0 corresponding to the disk midplane and r_cyl is the radial distance from this z-axis (see, e.g., Wolf et al. 2003; Stapelfeldt et al. 1998, see Fig. 6 for illustration). The parameter ρ₀ is determined by the total dust mass m_dust of the disk. The quantity r₀ is a reference radius with r₀ = 100 AU, whereas h is the vertical scale height and a function of r_cyl(2)The parameters α and β describe the radial density profile and the disk flaring, respectively. The disk extends from an inner cylindrical radius r_in to an outer one r_out. Together with the quantity h₀, these five parameters are used to adjust the disk structure in order to fit the data. This ansatz has already been successfully applied to other edge-on circumstellar disks, such as HH 30 (Madlener et al. 2012), CB 26 (Sauter et al. 2009), IM Lupi (Pinte et al. 2008), HV Tauri C (Stapelfeldt et al. 2003), and HK Tauri (Stapelfeldt et al. 1998). Integrating Eq. (1) along the z-axis yields the surface density (3)

4.2. Dust

In our model we focus on radiative transfer through the dust that dominates the transport of radiation and thermal structure of the disk (Chiang & Goldreich 1997). The density structure described in Eq. (1) is given by the gas in the disk. Because the dust in the disk is coupled to the gas, its density distribution is also described by this equation. The dust grain properties in our model can be divided into three groups: the shape of the dust grains, their chemical composition, and their size distribution.

4.2.1. Grain shape

We assume the dust grains to be homogeneous, spherical, non-aligned, and non-oriented particles (Mie theory), although they are expected to feature a much more complex and fractal structure. As shown by Voshchinnikov (2002), shape, chemical composition, and size of dust grains cannot be determined separately, but only in combination. Hence, we restrict our model to the simplest but least ambiguous shape possible.

4.2.2. Chemical composition

To model the chemical composition of the dust grains we use the homogeneous mixture of smoothed astronomical silicate and graphite with a mean density of ρ_grain = 2.5 g cm^-3 and optical properties described by Weingartner & Draine (2001). Relative abundances of 62.5% astronomical silicate and 37.5% graphite are used (Draine & Lee 1984; Weingartner & Draine 2001). This grain model has been employed successfully in the modeling of HH 30 (Madlener et al. 2012) and CB 26 (Sauter et al. 2009), for example. Given the index of refraction, the optical attributes of a homogeneous sphere can be calculated using Mie theory. We extrapolate the complex refractive indices to a wavelength of 2.7 mm to cover all available observational data. This is applicable since both the real and the imaginary part of the refractive index show asymptotic behavior in this wavelength regime. Because graphite is a highly anisotropic material, it is necessary to treat the two possible alignments of the electric field to the crystal axis independently using the so-called approximation for graphite spheres. That means, if Q_ext is the extinction coefficient, then (4)explains the dependence of Q_ext on the components of the dielectric tensor (ϵ_∥ and ϵ_⊥) of the electric field parallel and perpendicular to the crystallographic axis. Draine & Malhotra (1993) showed that this approximation is sufficiently accurate for extinction curve modeling.

4.2.3. Grain sizes

The sizes of the dust grains in our model are distributed according to a power law of the form (5)Here, a is the dust grain radius and n(a) the number of dust grains with a specific radius. For a_min = 5 nm and a_max = 250 nm this grain-size distribution becomes the commonly-known MRN distribution of the ISM by Mathis et al. (1977). To properly model the different grain sizes in a dust grain mixture, an arbitrary number of separate dust grain sizes within a given interval [a_min : a_max] has to be considered. Wolf (2003a) showed that the observables resulting from radiative transfer (RT) simulations considering each grain species separately are close to those based on weighted mean dust grain parameters of the dust grain mixture. Thus, we use weighted mean values for the optical properties of the dust grain mixture. Moreover, since the separate processes (e.g., grain growth, dust settling, grain-grain interaction, mixing processes) and their mutual influence during the evolution of the circumstellar environment are still rather poorly understood and also to reduce the number of free parameters necessary to describe the dust grain mixture, we assume no spatial dependence of its properties and assume the power law distribution (Eq. (5)) to be valid throughout the circumstellar disk.

4.3. Heating sources

Stellar heating is the primary heating source for the circumstellar disk. The radiation of the central star heats the dust which then, in turn, re-emits at longer wavelengths. The star, which is represented by a single black body, is characterized in our model by its effective temperature T_⋆ and bolometric luminosity L_⋆. Owing to obscuration by the circumstellar disk, the illuminating source of the system cannot be observed directly, leading to weak constraints on these two parameters from observation.

Viscous heating of the disk is neglected and we ignore accretion and turbulent processes in the disk.

5. Means of modeling

The goal of the modeling process is to find a coherent model for the circumstellar disk of the Butterfly Star explaining all the described observational data. Using continuum RT simulations we compute observables from the model that are then compared to the observed data to find the best-fit model. In this section we describe the means of modeling in detail.

5.1. Radiative transfer

For our continuum RT simulations we use the program MC3D (Wolf et al. 1999; Wolf 2003b). It is based on the Monte-Carlo method and solves the continuum radiative transfer problem self-consistently. It makes use of the temperature correction technique described by Bjorkman & Wood (2001), the absorption concept introduced by Lucy (1999), and the enforced scattering scheme proposed by Cashwell & Everett (1959). The optical properties of the dust grains (scattering, extinction and absorption cross sections, scattering phase function) and their interaction with the radiation field is calculated using Mie theory. Multiple and anisotropic scattering is considered.

In order to derive a spatially resolved dust temperature distribution, the model space has to be subdivided into volume elements inside which a constant temperature is assumed. Both the symmetry of the density distribution and the density gradient distribution have to be taken into account. We use a spherical model space, centered on the illuminating star and an equidistant subdivision of the model in θ-direction, while a logarithmic radial subdivision is applied to resolve the temperature gradient at the very dense inner region of the disk.

The RT is simulated at 107 wavelengths; 100 of them are logarithmically distributed in the wavelength range [0.05 μm, 2000 μm] and the remaining 7 wavelengths are distributed in the range [894 μm, 3000 μm]. The quantities we derive with MC3D from the model are

• 1.

  (sub)mm maps at 894 μm and1.3 mm;

• 2.

  NIR images at 1.13 μm, 1.60 μm, 1.87 μm, and 2.07 μm; and

• 3.

  the SED.

5.2. Model fitting properties

The simulated (sub)mm maps are convolved with the point spread function (PSF) that is described by an elliptical Gaussian function based on the beam major and minor axis as well as the beam position angle from the corresponding observation (see Sect. 2). To account for the rotation of the observed maps, the beam position angle of the PSF is adjusted by the same angle. For the simulated NIR images we use the corresponding PSF obtained with the PSF modeling tool Tiny Tim v.7.4 (Krist et al. 2011; Krist & Hook 2004) for the convolution. Additionally, the convolved simulated maps have the same pixel scale as the observed data. This way we ensure that the simulations are comparable to the observations.

There are primary features that we want to reproduce with our model. These also determine the criteria for the best-fit model. For the SED we aim to reproduce the thermal re-emission spectrum over almost three orders of magnitude from the MIR down to the mm wavelength regime. For resolved images the issue is not as simple as for the SED. Our model is rotationally symmetric and thus does not account for any related asymmetry that is seen in the observations. We consider the maps of the disk to be two different groups, the (sub)mm maps and the NIR maps.

• Maps in the (sub)mm: although the two (sub)mm maps are simplystructured they provide two decisive features that constrain ourmodel and that we want to reproduce, the peak flux density and thespatial brightness distribution (mainly the radial brightnessdistribution with constraints for the vertical extent). The uv data is, ingeneral, used to analyze the (sub)mm data because this avoids thenon-linear deconvolution step. Here, as the source is quite strongand barely resolved in one direction, the deconvolution isstraightforward and we use image plane comparison.

• Maps in the NIR: the four images in the NIR show more structures and details than the (sub)mm maps. In addition to the disk, which appears as a dark lane, a very complex and highly-structured envelope surrounds the disk. It is not an objective of this study to reproduce this envelope with its wavelength-dependent morphology, so we restrict ourselves to reproducing the width of the dust lane and, consequently, the wavelength dependence of the width of the dust lane. This is the main information provided by the NIR images because it constrains the vertical opacity structure of the disk.

5.3. Quality of the fit

For each comparison between model and observation on the above points, we determine a goodness of fit value that characterizes how well the model fit matches the observational data and so the best-fit model is represented by the smallest goodness of fit value of the investigated parameter space.

SED: the goodness of fit of the SED is given by (6)Here, N is the number of observed wavelengths, ω_i are the observed flux densities and μ_i the synthetic flux densities at these wavelengths. The observational uncertainties are considered by σ_i. As the re-emission part of the SED of the Butterfly Star, which we want to fit, is better sampled in the wavelength range from 12 μm to 38 μm thanks to the SST/IRS spectrum than at longer wavelengths up to 2.7 mm, we include weighting factors g_i in the calculation of the . They are choosen in such a way to ensure that the better and the worse sampled parts of the SED are equally considered in the calculation, i.e., that the flux measurements ranging from 12 μm to 38 μm have in total the same weighting as the data points in the range from 70 μm to 2730 μm altogether. A natural weighting fit would rely excessively on a narrow range of wavelengths and could possibly miss gross features which have a much lower S/N. By using an adjusted weighting, we ensure that we do not sacrifice the overall agreement for tiny details in the SED.

Maps: for every map the goodness of fit is calculated by (7)where μ is the modeled data, ω the observed data, and σ the observational uncertainty. Therefore, for the (sub)mm maps every pixel j whose flux density in the observed data is larger than the background noise in a square box with a side length according to the diameter of the object and centered on the object is taken into account, which lets us rewrite Eq. (7) as (8)In this case, N is the number of pixels taken into account in this box and σ represents the background noise of the observed map. Instead of a map-based goodness of fit there is also the possibility of choosing an approach where only the profile along the disk axis is taken into account. However, this method would not allow for any direct constraints on the disk scale height. Therefore, we decided to choose the approach described above.

To get a total goodness of fit that characterizes a model fit we combine all (9)Based on the we get from Eq. (9), we give our modeling errors, i.e., constraints on the model parameters, as the range where we can alter the parameter values without changing by more than 10%. Allowing for a larger variation of than 10% gives generally worse results.

5.4. Parameter space study

The modeling process is divided into two steps. We simultaneously fit the (sub)mm data for the first time and figure out if it is possible to find a model fit that reproduces these data using the described model (Sect. 4) with a single dust grain-size distribution. Based on the results of this first modeling step (see Sect. 6.2) we modify our model setup and also include the SED and the NIR images in the fitting process.

On the basis of the model described in the previous sections, the parameter space we deal with is defined by several free parameters. These adjustable model parameters are used to reproduce the key features of our observations as described in Sect. 5.2. The volume of the parameter space consists of six and nine free parameters in the first and second modeling step, respectively. For our study, the selection and range of these parameters is based on modeling of other similar objects and previous modeling efforts (see, e.g., Sauter et al. 2009; Wolf et al. 2003). To constrain the model parameters as strongly as possible and to find the parameter set of the best-fit model we use a grid search based method.

6. First modeling step

We first consider a disk model with no spatial variations of the dust grain properties.

6.1. Model parameters

The six adjustable model parameters we use in the first step are the following:

• The exponents α and β that describe the radial density profile of thedisk and the disk flaring, respectively (see Eqs. (1)and (2)). In our modeling, both parameters are treated asindependent quantities. We vary α in the range of [1.0, 4.0] and β in the range of [1.0, 2.0].

• The scaling factor h₀ of the disk’s scale height (see Eq. (2)). Because the value of the reference radius r₀ is 100 AU, h₀ equals the scale height at a radial distance from the star of r_cyl = 100 AU (h₀ = h(100  AU)). We consider h₀ in the range of [1.0, 25.0] AU.

• The dust mass m_dust of the disk. The dust mass range we probe is [1.0 × 10^-5, 1.0 × 10^-2] M_⊙.

• The inner disk radius r_in. Inner holes in protoplanetary disks are a common prediction of theories of disk evolution and are observed in various circumstellar disks (e.g., Alexander & Armitage 2009; Andrews et al. 2011; Gräfe et al. 2011; Espaillat et al. 2012). The inner radius is varied in the range of [0.1, 50.0] AU. The lower limit is approximately the dust sublimation radius.

• The maximum grain radius a_max. Grain growth is a major issue for protoplanetary disks as it is the first step toward the formation of planets. For this upper radius of the grain-size distribution we consider nine different values: [0.25, 1.0, 10, 30, 50, 75, 100, 300, 1000] μm.

For the minimum grain radius a_min we use 5 nm for both modeling steps. The outer radius of the disk r_out is fixed at 300 AU for both steps based on the (sub)mm observations and the modeling results from Wolf et al. (2003). On the basis of the mass of the star of ~1.7   M_⊙ resulting from observations of the gas kinematics (Dutrey et al., in prep.) and pre-main-sequence evolutionary tracks (e.g., Hillenbrand & White 2004), we found that L_⋆ = 5   L_⊙ and T_⋆ = 4500  K characterize the central star best. These values are compatible with the radiative transfer simulations and are used in the whole modeling. According to the precisely determined inclination i of the circumstellar disk of the Butterfly Star of 90° ± 3° (Sect. 1), we fix the inclination at 90° in the first modeling step. Because of its location in the Taurus-Auriga molecular cloud complex, a distance d of the object of 140 pc is adopted throughout the modeling process. An illustration of the density structure of a model fit from step one is seen in Fig. 6. Table 2 shows an overview of the fixed parameters used in the modeling process.

6.2. Results

The results presented in this subsection are based on the model and applied assumptions outlined in Sects. 4 and 5.

Using the described modeling setup of the first step, which is mainly characterized by a single dust grain-size distribution, it is not possible to find a coherent model that properly reproduces the (sub)mm data. Figures 7 and 8 show the best-fit results for a grain-size distribution with small ISM-sized dust grains (a_max = 0.25 μm) and larger dust grains (a_max = 100 μm), that yield the motivation for the second modeling step. Dust grains with sizes comparable to those found in the ISM (a_max = 0.25 μm) are neither capable of reproducing the brightness structure of the (sub)mm maps nor do they allow us to fit the shallow mm slope of the SED (Fig. 7, see also Fig. 5 and Sect. 3.2). The radial brightness distributions of the (sub)mm data in Fig. 7 show an insufficient flux density at almost all regions in the disk. Only in the outermost disk regions the small dust grains yield a flux density that is in agreement with the observation. From the SED in Fig. 7 it can be seen that the MIR and far-infrared (FIR) part is reproduced using ISM-sized dust grains, but that it is not possible to account for the (sub)mm part of the SED. These results strongly suggest the need for larger dust grains in the disk, particularly in the region close to the central star.

We find that the brightness structure of the (sub)mm maps is reproduced best using a dust grain-size distribution with a_max = 100 μm (Fig. 8). Furthermore, the large grains fit the (sub)mm part of the SED better than the small ones while particles significantly smaller or larger than a_max = 100 μm are inapplicable. In contrast, it is not possible to account for the IR part of the SED using the large dust grains. Moreover, they yield a flux density that is too high, especially in the outer disk regions as can be seen in the (sub)mm radial brightness profiles (Fig. 8). This indicates that small dust grains, especially in the outer parts of the disk, are still needed in the model.

Nevertheless, with no combination of the used variable parameters it is possible to account for the proper shape of the radial brightness profiles. We conclude, that to properly fit all the observational data and to find a coherent model that explains all these data, we get from the first modeling step that we need large dust grains with a_max ~ 100 μm as well as small ISM-sized dust grains at the same time in the disk model. Additionally, from the disk scale height h(100  AU) that amounts to 6 AU and 12 AU using the small and the large dust grains, respectively, and from the described brightness profiles (Figs. 7 and 8), we conclude that we need the large grains to be located around the disk midplane with a small vertical extent and concentrated toward the central star.

In addition, no need for a hole in the inner disk is found from the modeling. As a consequence we fixed the inner disk radius in the second modeling step at 0.1 AU which is approximately the dust sublimation radius.

7. Second modeling step

The results of the first modeling step show that it is not possible to properly fit the (sub)mm data at the same time and that at least two dust grain-size distributions are necessary, with larger dust grains concentrated around the disk midplane.

7.1. Model parameters

Owing to the modification of our model setup based on the results of the first modeling step (Sect. 6.2), we use a set of nine adjustable parameters in the second step instead of six. In contrast to the first step, we consider two different dust grain-size distributions in the disk model. Their location in the disk can be seen in Fig. 9. They differ in the maximum grain size with a_max,ld > a_max,sd (ld: large dust, sd: small dust). The parameter a_max,sd is set to 0.25 μm which corresponds to the dust grain-size distribution found in the ISM (Mathis et al. 1977) and that we expect in the outer regions and upper layers of the circumstellar disk of the Butterfly Star (Wolf et al. 2003). Because it is not necessary to vary the inner disk radius r_in we fix it at 0.1 AU. The minimum grain radius a_min for both grain-size distributions, the outer radius of the disk r_out, the distance d of the object, L_⋆, and T_⋆ are the same as in the first modeling step (see also Table 2). The parameters α, β, h₀, and m_dust are used again as free parameters. In addition, we use the following adjustable parameters in the second step of our modeling:

• The two quantities ζ_ld and r_out,ld that describe the vertical and radial extentof the region where the grain-size distribution with the larger dustgrains is located (see Fig. 9). The parameter ζ_ld isradially dependent in the same way as in Eq. (2). Weemphasize that ζ_ld is not a scale height in the sense of h₀.In our model the circumstellar disk is described by one densitydistribution (see Eq. (1)). Therefore, the verticalextent described by ζ_ld has a sharp upper limit. The range we considerfor ζ_ld (at 100 AU) is [1.0, 25.0] AU andfor r_out,ld it is [50, 300] AU.

  Fig. 10 Radial brightness profile at 1.3 mm based on the observation and our best-fit model. For details see Sect. 7.2.

• A scaling factor sm that adjusts the distribution of the dust mass of the disk within the two different dust regions. This parameter is affected by the other free disk parameters, i.e., if the total number of dust grains is changed in either of the two different dust regions by any of the adjustable parameters, the ratio of the dust mass of these two regions is modified although this scaling factor is not varied.

• The inclination i of the circumstellar disk. Despite the small uncertainty of the inclination measurement, here we allow the inclination to vary in the range of 90° ± 3°.

• The maximum grain radius of the dust grain-size distribution containing the larger dust grains a_max,ld. We consider three different values: [50, 100, 150] μm.

Using this model setup in the second modeling step we are accounting for grain evolution and its dependence on the radial and vertical position in the circumstellar disk.

7.2. Results

The results presented in this subsection are based on the model and applied assumptions outlined in Sects. 4 and 5.

With the model setup of the second modeling step, which is mainly characterized by two different dust grain-size distributions, it is possible to find a model fit that explains all the observational data taken into account and therefore a coherent model for the circumstellar disk of the Butterfly Star. The values of the parameters of our best-fit model and the constraints on the model parameters resulting from our multi-wavelength modeling can be found in Table 3. Only in the small parameter space described by the lower and upper limits of the listed parameters is it possible to find model fits that satisfactorily reproduce all observational data, i.e., does not change by more than 10% compared to the best-fit value. Figures 10 and 11 show the radial brightness profiles of the (sub)mm observations overlayed with the corresponding counterparts based on the best-fit model. The optical depth τ at 894 μm amounts to 1.75 and at 1.3  mm to 0.64. This means that the disk is optically thick at submm wavelength and optically thin at mm wavelength. The effective specific dust opacity is determined to be 2.58   cm²   g^-1 and 0.96   cm²   g^-1 at 894 μm and 1.3  mm, respectively (based on the dust model outlined in Sect. 4.2). The observed and modeled re-emission SED is represented in Fig. 12. In Table 4 the width of the dust lane at the point of minimum separation between the two lobes in both the observed and modeled scattered light maps is listed (see also Fig. 3). The width of the dust lane stays almost constant across the whole radius of the disk such that small deviations are covered by the listed uncertainties. Overall, our best-fit model is in very good agreement with all considered observations, especially taking into account the range of wavelengths and the variety of observations analyzed in the fitting procedure.

In Fig. 10 it can be seen that our model almost perfectly fits the 1.3 mm observation. The overall accuracy, i.e., how good the modeled profile matches the observed one in terms of the observed background noise, is 1.39   σ_mm ± 0.98   σ_mm. Although the central brightness minimum in the submm observation (Fig. 11) is not reproduced by our model we do see a clear flattening of the profile which is definitely caused by an optical depth effect (see also Sect. 8.1). With an accuracy of 0.96   σ_submm ± 0.77   σ_submm, our model also provides a satisfying fit to the submm observation. Furthermore, the modeled SED is in very good agreement with observations from the MIR over the FIR to the mm regime, meaning that almost all modeled flux densities are within the uncertainties of the observed flux densities (Fig. 12). In particular, it reproduces the (sub)mm flux densities and hence the mm spectral index. The width of the dust lane in the modeled NIR images is in agreement with the observed width taking into account the measurement uncertainties (Table 4). That means that the four modeled NIR images nicely reproduce the observed width of the dust lane and therefore also its wavelength dependence which was the goal for this wavelength range.

8. Discussion

8.1. Brightness minimum in the SMA observation

The central brightness minimum seen in the observed 894 μm map was first described in Wolf et al. (2008) where two explanations for this minimum are discussed, an inner disk hole and an optical depth effect. When comparing the beam size of the two (sub)mm maps (see Sect. 2) it can be seen that the beam size of the 1.3  mm observation is significantly smaller than that of the 894 μm map. This means that any disk hole that would be capable of producing such a local brightness minimum in the submm observation would also show up in the 1.3 mm observation. However, at 1.3  mm no sign of a central brightness minimum or even a flattening of the profile was found. Therefore, this rules out the possibility that the local minimum seen in the 894 μm map is caused by the lack of emitting dust, i.e., an inner disk hole. Moreover, when comparing both maps and taking the beam sizes into account, we can conclude that the disk of the Butterfly Star is optically thick even at submm wavelength. Thus, the central brightness minimum in the submm observation is produced by an optical depth effect (for details see Wolf et al. 2008). This is also supported by the optical depth at both wavelengths (see Sect. 7.2) resulting from our best-fit model. We emphasize that the observed brightness minimum has a low significance of ~1.4σ and ~1.8σ with respect to the maxima to the left and right of the center, respectively. As our best-fit model shows a clearly flattened profile in the submm that accounts for the optical depth effect and is in very good agreement with the observation (Sect. 7.2), it is a satisfying result.

8.2. Dust evolution

We found a coherent model capable of explaining all considered observations. While it was not possible to achieve this goal with a model using a single dust grain-size distribution we succeeded by using two different distributions. Both distributions differ in the size of the maximum particle radius a_max and in the spatial distribution in the disk. We found that for a_max values of 0.25 μm and 100 μm for the smaller and larger dust grain-size distribution, respectively, fit the data best. The grain size we determine for the larger dust in the disk is almost three orders of magnitude larger than upper grain sizes given for the ISM in the literature. This strongly indicates that grain growth is taking place in the circumstellar disk of the Butterfly Star. However, this dust grain size is still one order of magnitude smaller than the maximum grain size derived from modeling other young circumstellar disks, such as for IM Lupi (Pinte et al. 2008). There, a maximum size of a few mm has been found. In the framework of our model we exclude the strong influence of such large grains on the (sub)mm observations of the disk of the Butterfly Star. Our model clearly shows the need for small ISM-sized as well as larger (a_max ~ 100 μm) dust grains in this disk. In particular, the smaller grains are located in the upper layers and outer region of the disk and the larger ones in the inner disk region. The findings on the grain sizes are only valid in the context of the assumed dust properties (see Sect. 4.2). Assuming a different distribution of the sizes of the dust grains may lead to different results. But what in general can be concluded is that the larger dust grains have a different spatial distribution than the smaller ones.

Comparing the disk scale height of our model and the corresponding vertical extent of the region where the large dust grains are located at 100 AU yields that the large dust grains are settled toward the disk midplane. Furthermore, comparing the radial extent of the large dust region with the outer disk radius it can be seen that the large dust grains are not distributed over the whole radial extend of the circumstellar disk, but are concentrated toward the star.

We conclude that besides grain growth in the disk of the Butterfly Star, we find evidence for vertical settling and radial segregation of the dust in this disk. These three mechanisms are the key issues in the process of planetary core formation, which start the formation process of planetesimals. Our coherent model not only gives strong constraints on the dust properties in this disk, but we can also provide quantitative constraints on the vertical and radial distribution of the phase with large dust grains for the first time.

8.3. Disk mass

Our best-fit model shows that we need a dust mass of M_dust = 9    ×    10^-4   M_⊙ to account for the observed flux densities. This value is well constrained (see Table 3) and stronger variations of the dust mass would have a significant influence on the re-emission SED, for example. The dust mass of the region where the large grains are distributed is 2.2 × 10^-4   M_⊙ and that of the small grains is 6.8 × 10^-4   M_⊙. Assuming a typical dust-to-gas ratio of (10)we get a total disk mass of ~ 9 × 10^-2   M_⊙ and determine a maximum absorption opacity of 1.97   cm²   g^-1 at 1.3  mm. This is in good agreement with the results of Wolf et al. (2003) and very similar to other YSOs such as IM Lupi (Pinte et al. 2008), DL Tau, GO Tau, and HL Tau (Andrews & Williams 2007).

The derived disk mass depends on the adopted assumptions about the chemistry and shape of the dust grains (see Sect. 4.2). Fractal grain shapes and the spatial orientation of every dust grain are not taken into account as this would only introduce further free parameters without providing a qualitative improvement of our understanding of the disk properties and distribution of large grains. The complex shapes of the dust grains would have implications on the light scattering and absorption behavior of the grains (e.g., Wright 1987; Lumme & Rahola 1994). Besides that, it also can be thought of dust grains with almost the same absorption cross section as spherical dust grains but with much less mass. Such fluffy particles with a porosity up to 90% are discussed by Voshchinnikov et al. (2007).

Furthermore, we point out that the disk mass is proportional to the grain density ρ_grain and the number of dust grains. Within our study we can constrain the disk structure, the grain size, and their number, but not the density of one grain. For our investigation of the circumstellar disk, we used ρ_grain = 2.5   g   cm^-3, but a more porous, fractal structure of the dust grains may result in a smaller value. In turn, this will alter our estimate for the disk mass by the same factor. However, our general results on the different spatial distributions of the smaller and larger dust grains are not affected by these assumptions.

8.4. Disk structure

Our multi-wavelength modeling allows us to quantitatively constrain most of the geometrical parameters of the circumstellar disk of the Butterfly Star. A flared geometry with a scale height of ~10 AU at a reference radius of 100 AU is required. The determined flaring index β of ~1.14 and scale height are comparable to the values for the YSO IM Lupi by Pinte et al. (2008). Sauter et al. (2009) also found for the circumstellar disk in the Bok globule CB 26 a scale height of ~10 AU and with β ~ 1.4 a slightly larger flaring index. For the circumstellar disk of HH 30, Madlener et al. (2012) found a flaring index comparable to what we find for the disk of the Butterfly Star, but a larger scale height (h ~ 15 AU). Our best-fit model has a midplane temperature of 20 K at 100 AU. Figure 13 shows the calculated midplane temperature as a function of the radial distance from the center. The determined value for T₁₀₀ is very similar to other YSOs such as IM Lupi (Pinte et al. 2008), CB 26 (Sauter et al. 2009), DG Tau, DG Tau-b, HL Tau (Guilloteau et al. 2011), and HH 30 (Madlener et al. 2012).

The hydrostatic scale height H is derived by (11)where k_B is Boltzman’s constant, G is the gravitational constant, μ = 2.33   g   mol^-1 (Ruden & Pollack 1991) is the mean molecular weight, and m_p is the proton mass. This yields a value of H = 6.9 AU at a radius of 100 AU. To compare this to the vertical extent of the region where the large dust grains are located (ζ_ld) an equivalent scale height at a radius of 100 AU h_ld can be determined by (12)For our best-fit model this results in a scale height for the large dust of h_ld = 3.7 AU and for the constraints on h₀ and ζ_ld in a range of 2.5 AU≤ h_ld ≤ 4.3 AU. This is substantially smaller than ζ_ld and substantially smaller than the hydrostatic scale height H which strongly indicates that the large dust grains in the circumstellar disk of the Butterfly Star are vertically settled.

8.4.1. Surface density

In Fig. 14 the derived surface density of our best-fit model for the circumstellar disk of the Butterfly Star as well as the upper and lower limit can be seen. The surface density exponent (see Eq. (3)) is found to be in the range of [0.04, 0.48]. These values are lower than those predicted by most theoretical models of disks (p ~ 1; Bell et al. 1997). Furthermore, studies of large samples of circumstellar disks in the Taurus-Auriga and Ophiuchus-Scorpius star formation regions by, e.g., Kitamura et al. (2002), Andrews & Williams (2007), and Guilloteau et al. (2011) show that values of p less than ~1 are common, but they are in general larger than 0.5. Therefore, our result for p is slightly smaller than that found for most other objects. We determine the surface density to be in the range of [3.3, 15.2]   g   cm^-2 and [2.9, 3.6]   g   cm^-2 at 5 AU and 100 AU, respectively, which is consistent with the previously mentioned surveys.

8.4.2. Disk stability

To estimate if the circumstellar disk of the Butterfly Star is self-gravitating, we determined the ratio μ_sg of the Kepler time τ_K and the local free-fall time τ_ffHere, Ω_K is the Keplerian angular velocity, G is the gravitational constant, and M_⋆ is the stellar mass. If μ_sg ≪ 1 the disk is dominated by the gravitational potential of the central star and therefore Keplerian, i.e., it is non-self-gravitating. In contrast, if μ_sg ≫ 1 the disk is dominated by the gravitational potential of the local mass distribution and so becomes near-Keplerian self-gravitating to fully self-gravitating with increasing μ_sg. Figure 15 shows the radial dependent ratio μ_sg and it can be seen that in our model the disk of the Butterfly Star is non-self-gravitating at all radii. This also implies that gravitational instabilities in the disk are very unlikely.

To check this, we made use of the Toomre criterion (Toomre 1964) that describes when a disk becomes unstable through gravitational collapse. The Toomre Q parameter is defined as (16)where c_s is the local sound speed. In terms of Q, a disk is unstable to its own self-gravity if Q < 1, and stable if Q > 1. As Fig. 16 clearly shows, in our model Q > 1 throughout the entire disk, which means that the circumstellar disk of the Butterfly Star is gravitationally stable at all radii as is expected by the determination of μ_sg.

8.5. Comparison with theoretical models

We compared our results qualitatively with theoretical models for the evolution of dust in circumstellar disks by Garaud et al. (2004), Garaud & Lin (2004), Barrière-Fouchet et al. (2005), and Garaud (2007), for example. The Butterfly Star is characterized as a Class I YSO and is surrounded by a circumstellar disk that is optically thick even at submm wavelengths and by a substantial envelope. This implies that this object is rather young and potential evolution of dust is in its early phase. The main points of current models on the first phase of planet formation are summarized below. The gas in the protoplanetary disk is partially pressure supported, meaning that both the centrifugal force and the gas pressure counteract the gravity. Because of an outward decreasing gas pressure gradient, the gas orbits at sub-Keplerian velocity around the star. Small dust particles are well-coupled to the gas. As they do not experience the same radial pressure gradient as the gas they orbit with the Keplerian velocity and therefore feel a net inward force causing them to drift inward. As a consequence of processes such as Brownian motion and turbulence in the disk, the dust grains agglomerate as a result of collisions and form aggregates of ever increasing size and mass. In the context of our model we found evidence for the early stage of this grain growth phase in the circumstellar disk of the Butterfly Star with grain radii of up to ~100 μm. This is significantly larger than what is commonly found in the ISM (~0.25 μm). With the increase of the importance of gravity resulting from the continuing grain growth, the particles decouple from the pure gas motion and settle toward the disk midplane. Because they have a faster velocity than the gas does, the larger particles see a headwind that saps their angular momentum and causes them to spiral in toward the star. This radial drift occurs on a time scale that is much shorter than the disk lifetime. From our multi-wavelength modeling of the disk of the Butterfly Star we also found quantitative evidence for the settling of larger dust grains toward the disk midplane while smaller particles are retained in the upper layers of the disk. We found the larger dust grains to be concentrated toward the star which is in agreement with the predicted inward spiraling. Moreover, theoretical models predict that the gas drag efficiency varies according to the grain size, where intermediate-sized particles (100 μm to 10 cm) experience the strongest perturbation to their movement. The largest particles that we have found based on our model exhibit a radius of ~100 μm. Therefore, our findings for radial segregation of the dust is not unexpected. Theoretical work by Garaud et al. (2013) suggests the coexistence of two particle populations, one with larger particles and one with smaller particles, at different radial positions in the protoplanetary disk in contrast to a single continuous population from small to large sizes. Our need for two different dust grain-size distributions to properly model the properties of the disk of the Butterfly Star agrees well with this prediction.

We find that our results are consistent with current models for the evolution of dust in circumstellar disks. Our findings of grain growth, vertical settling, and radial segregation in the disk of the Butterfly Star are common predictions of theoretical models.

8.6. Observability with ALMA

Figures 17 and 18 show simulated ALMA observations at 1.3 mm and 894 μm, respectively. They are based on the modeled (sub)mm maps of our best-fit and are created using the CASA simulator (Jaeger 2008). Both figures represent the use of the full ALMA array with an observing time of 2 h, in which Fig. 17 reflects configuration 15 and Fig. 18 configuration 13. The difference between the simulations shown on the left and the right in both figures is that for the right one only one single dust grain-size distribution with a_max = 100 μm is used whereas the simulation on the left exactly represents our best-fit model. The dust re-emission in the region that is encased by the 20% and 40% contour level in Fig. 17, left and Fig. 18, left, respectively, is strongly dominated by the larger dust grains. Because of the much higher angular resolution compared to the SMA and the PdBI, the radial extent of this region is cleary visible. The difference in the radial brightness structure between the presence of only one grain-size distribution or, as we have found in our study, two different grain-size distributions in the circumstellar disk of the Butterfly Star can be clearly seen in both figures. Thus, observations with the full ALMA array will allow us to distinguish between the presence of one or two different dust grain-size distributions and let us check our current model of the disk of the Butterfly Star. Most likely, these observations will allow us to give further constraints on the properties of this prominent circumstellar disk.

9. Summary and conclusion

We have compiled a high-quality data set for the circumstellar disk of the Butterfly Star, spanning a wide range of wavelengths. We obtained images in the NIR and in the (sub)mm regime, as well as photometric data and a spectrum. We have constructed a detailed model that allows the interpretation of all of these observations with one single set of parameters. A systematic analysis of the parameter space allows us to establish strong constraints on all the parameters of the model. The conclusions we obtain are the following:

• 1.

  The disk structure is well constrained. The disk extends froman inner radius of 0.1 AU up to 300 AU. The scaleheight of the disk is ~ 10 AU at 100 AU and varies witha flaring index of ~ 1.14. The surface density exponent is found to bein the range of [0.04, 0.48], which is slightly smaller than therange found in general for other circumstellar disks (Kitamuraet al. 2002; Andrews& Williams 2007; Guilloteauet al. 2011). The midplanetemperature is determined to be 20 K at100 AU.

• 2.

  The dust mass is found to be ~ 9 × 10^-4   M_⊙ under the assumption of spherical grains and ρ_grain = 2.5   g   cm^-3. The dust mass of the region where the large grains are distributed is ~ 2.2 × 10^-4   M_⊙ and that of the small grains is ~ 6.8 × 10^-4   M_⊙. With a typical dust-to-gas ratio of 1/100, the total disk mass amounts to ~ 9 × 10^-2   M_⊙, compared to a mass of the central star of ~ 1.7   M_⊙. Based on our constraints on the model parameters, the disk of the Butterfly Star is found to be non-self-gravitating and is, according to the Toomre criterion, gravitationally stable at all radii.

• 3.

  The disk of the Butterfly Star has already entered the first phases of planetary formation. We found quantitative evidence for grain growth, vertical settling, as well as radial segregation in this disk. It is not possible to find a coherent multi-wavelength model of the disk using a single grain-size distribution. To do so, the use of at least two different distributions is necessary. In the outer parts and upper layers of the disk, small ISM-sized dust grains can be found. In the inner disk regions, larger dust grains with a maximum particle radius of up to ~ 100 μm are located. We constrain this region where the larger grains are distributed to a radial extent of [170, 210] AU and a vertical extent of [4, 7] AU at 100 AU. The equivalent scale height of this region amounts to a range of [2.5, 4.3] AU and for our best-fit model to 3.7 AU which is substantially smaller than the vertical extent and the hydrostatic scale height (H = 6.9 AU), thus providing strong support for vertical settling of the large dust grains in the disk. Our findings of dust evolution in this disk are consistent with current theoretical models.

• 4.

  The millimeter spectral index indicates that the dust grains in the disk midplane have already grown to sizes larger than those found in the ISM, supporting our results.

• 5.

  Although each individual observation provides valuable information on the disk, it is necessary to combine a large set of independent observations (SED and spatially resolved images) from different wavelength regimes in a multi-wavelength study. This approach allowed us to derive qualitatively new conclusions which were not obvious on the basis of individual data sets alone and to strongly reduce model degeneracies.

The large number of observations in different wavelength domains, as well as our coherent multi-wavelength modeling, make the Butterfly Star one of the best studied protoplanetary disks so far. With the upcoming completion of the next-generation interferometer ALMA, observations using this facility will allow us to test our findings of dust evolution, will let us further constrain the radial and vertical structure of this fascinating disk, and will refine our understanding of the evolution of protoplanetary disks.

Acknowledgments

This work is supported by the DFG through the research group 759 “The Formation of Planets: The Critical First Growth Phase” (WO 857/4-2) and project WO 857/12-1. We are grateful to the anonymous referee for providing useful suggestions that greatly improved this paper.

References

All Tables

All Figures