© ESO, 2011

1. Introduction

Herbig Ae/Be stars are intermediate-mass pre-main-sequence stars first described as a group by Herbig (1960). They are characterized by the presence of a circumstellar gas and dust disk (e.g. Waters & Waelkens 1998), which is the remnant of the star formation process. There is by now a wealth of observational evidence suggesting that these disks are the site of planet formation. For instance, imaging at optical, infrared and millimeter wavelengths shows the presence of disk gaps and/or inner holes (e.g. Grady et al. 2005; Verhoeff et al. 2011), pointing to clearing of substantial parts of the disk. The composition of the dust in the disk strongly differs from that of dust in the interstellar medium. Large, millimeter sized dust grains have grown in the disk and settled to the mid-plane, depleting the disk of small, sub-micron sized particles. In addition, crystalline silicates are detected in the mid-infrared spectra of many disks, pointing to substantial grain processing (e.g. Bouwman et al. 2001; van Boekel et al. 2005). The recent discovery by direct imaging of exo-planets orbiting intermediate-mass stars (e.g. Marois et al. 2008; Kalas et al. 2008; Lagrange et al. 2010) provides strong support for the interpretation that disks surrounding Herbig Ae/Be stars are in the process of planet formation.

However, many basic questions related to planet formation need still to be clarified. For example: what is the main channel for planet formation as a function of stellar mass? How does planet formation affect the composition of the gas and dust and its spatial distribution? Which observed properties of disks can be used as a signpost of planet formation, and can we derive an empirical evolutionary sequence towards mature planetary systems? In order to answer these questions, multi-wavelength observations of the spatial structure and of the gas and dust composition of proto-planetary disks are needed. In this study, we focus on the spatial distribution of crystalline silicates in the Herbig Be star HD 100546. This star has been observed with Herschel/PACS as part of the DIGIT¹ open time key programme that aims to address these questions (Sturm et al. 2010).

Crystalline silicates are not found in interstellar space (Kemper et al. 2004) but are abundant in a significant fraction of proto-planetary disks (e.g., Olofsson et al. 2009; Juhász et al. 2010). They must therefore have been formed in situ, and trace the thermal and chemical history of the grains in the disk. They can form by thermal annealing of amorphous silicates (Fabian et al. 2000), or by direct condensation from the gas phase. Both processes require high temperatures, respectively above the glass temperature (~1000 K) or near the dust evaporation temperature (~1500 K). These temperatures are typically found close to the star.

Crystals are also observed much farther away from the star than expected based on their temperature (e.g., Olofsson et al. 2010) and are abundant in solar system comets, requiring efficient radial mixing or local production in shocks (Harker & Desch 2002), collisions (Morlok et al. 2010), parent-body processing (Huss et al. 2001) or stellar outbursts (Ábrahám et al. 2009). Additionally, crystals can also be destroyed by stellar winds (Glauser et al. 2009). However, the correlation of crystallinity with other stellar parameters such as age and luminosity in a large sample remains unclear (e.g. Watson 2009; Oliveira et al. 2010; Juhász et al. 2010). Detailed studies of individual disks are therefore needed in which spatially resolved information and spectroscopy are combined in order to establish the spatial distribution, abundance and chemical composition of crystalline silicates in relation to the overall disk geometry.

In this paper we focus on HD 100546, a Herbig Be star of spectral type B9.5Vne with one of the strongest observed crystalline dust features (Malfait et al. 1998). Apart from the unusual mineralogy, the spectral energy distribution (SED) of HD 100546 is characterized by a large mid-infrared excess and a modest near-infrared excess. This SED is explained by a disk gap, which creates a frontally illuminated wall at the far edge of the gap with a temperature of  ~200 K (Bouwman et al. 2003). The existence of the gap has been confirmed by spatially revolved spectroscopy of the gas as well, placing the wall at 13 AU (e.g. Grady et al. 2005). The small near-infrared excess is caused by a tenuous inner disk (Benisty et al. 2010).

The mid-infrared spectrum of HD 100546 has been observed with ISO (Malfait et al. 1998), TIMMI2 (van Boekel et al. 2005) and Spitzer (Juhász et al. 2010). The ISO observations revealed very pronounced crystalline silicate features, similar to those of comet Hale-Bopp (Crovisier et al. 1997). Comparison to laboratory spectra indicated that the features between 11 and 30 μm can be described by a dust population with a single temperature of 210 K (Malfait et al. 1998), coinciding with the disk wall. The longer wavelength features (33 and 69 μm) require an additional cold component (~50–70 K) to explain the feature strength.

The ISO spectrum was analysed in more depth by Bouwman et al. (2003), who study the forsterite abundance and spatial distribution with a radiative transfer model that assumes the dust is optically thin. They find a crystalline mass fraction that increases from 2% in the inner disk to 19% in the outer disk – which is inconsistent with radial mixing models that predict a decreasing fraction with radius.

The TIMMI2 spectrum was analysed using a similar optically thin model, but without taking into account the geometry, by van Boekel et al. (2005), who find a forsterite mass fraction of 5.0% in the 10 μm range. The forsterite features have been observed more recently with Spitzer by Juhász et al. (2010). They find a forsterite mass fraction of 5.6% from the 10 μm spectral range and 5.1% from the 20 μm spectral range.

In addition to the strength, the shape and central wavelength of the features also contain information on the temperature of the emitting dust, as these are temperature-dependent (Koike et al. 2006). This is especially true for the 69 μm feature where the feature broadens and shifts to longer wavelengths at higher temperatures (Bowey et al. 2002; Suto et al. 2006, Fig. C.1).

The 69 μm feature of HD 100546 has been observed in great detail with Herschel/PACS as part of the DIGIT open time key programme, and analysed by Sturm et al. (2010) using an optically thin approach. The observed central wavelength of the feature is 69.2 μm, whereas cold (50 K), iron-free forsterite² has a peak at 69.0 μm. Sturm et al. (2010) propose two different scenarios for the required wavelength shift of 0.2 μm.

In the first scenario, warm iron-free forsterite creates the observed shape. Using a weighted sum of the opacities, they find that pure forsterite with a temperature of 200–150 K has a peak at 69.2 μm and matches the feature shape. This scenario is consistent with the inferred temperature of the disk wall, but strongly overpredicts the features at shorter wavelengths. It is therefore inconsistent with the ISO and Spitzer data, unless optical depth effects play a significant role.

In the second scenario, Sturm et al. (2010) proposed an alternative explanation: colder forsterite, with a temperature around 50–70 K creates the 69 μm feature, and an admixture of a few percent iron shifts the peak longwards to 69.2 μm. Such a model does not overpredict the short wavelength features, consistent with the earliest model from Malfait et al. (1998). Optical constants for forsterite with a few percent iron have not (yet) been measured in the lab, and a compositional fit could not be made. The inferred wavelength shift follows from interpolation of measurements with 0% and 10% iron content (see Sturm et al. 2010).

The two scenarios lead not only to a different chemical composition of the forsterite, but also to a different location of emission within the disk. Iron-free forsterite originates from the disk wall, whereas iron-containing forsterite comes from colder material farther out – possibly near the disk midplane. Because forsterite has a large number of features over a broad wavelength range, regions with different optical depths are probed at different wavelengths, providing extra diagnostic power. To shed more light on the chemical composition and spatial distribution of forsterite in the disk, we need to explore how feature strengths and shapes over the entire wavelength range are influenced by these optical depth effects.

To do this, we need to use a 2D radiative transfer code, described in Sect. 3, to construct a physical and mineralogical model of the dust surrounding HD 100546. We construct a structural hydrostatic model in Sect. 4, which is consistent with a new VISIR mid-infrared image (Appendix A). In Sect. 5 we use this structural model to study the spatial distribution of crystalline forsterite, by fitting the forsterite feature strengths over a broad wavelength range and in particular the shape of the 69 μm feature.

We will discuss our results in the context of previous work in Sect. 6, as well as address the question of the origin of the forsterite and correlation with the disk gap. We will summarize our conclusions in Sect. 7.

2. Observations

2.1. SED and spectroscopy

To construct the observed SED of HD100546 (see Fig. 2), we use photometric data from van den Ancker et al. (1997) combined with the ISO spectrum from Malfait et al. (1998). We do not use the Spitzer (Juhász et al. 2010) and PACS (Sturm et al. 2010) spectra for constructing the SED because the ISO spectra encompasses the wavelength range of both, and the current flux calibration of the PACS spectrum is not yet better than that of ISO. There is no significant offset between the ISO and Spitzer spectra that could influence our analysis. We do make use of the increased sensitivity of PACS and Spitzer for calculating the integrated strength of the forsterite features.

3. Radiative transfer and dust model

3.1. Radiative transfer model

The dust radiative transfer code used in this paper is MCMax (Min et al. 2009), a 2D Monte Carlo code. It is based on the immediate re-emission procedure from Bjorkman & Wood (2001), combined with the method of continuous absorption by Lucy (1999). The code has been benchmarked against other radiative transfer codes for modelling protoplanetary disks (Pinte et al. 2009) and has been successfully applied for modelling spatially and spectrally resolved observations (e.g. Verhoeff et al. 2010). The radial grid around the inner radius and disk wall is refined to sample the optical depth logarithmically. The SED, images and flux contributions (see Sect. 5.3) are calculated by integrating the formal solution to the equation of radiative transfer by ray-tracing.

In addition to radiative transfer, the code explicitly solves for the hydrostatic vertical structure of the disk, with the implicit assumption that the gas temperature is set by the dust temperature³. Therefore it only needs the radial distribution of dust to construct the structural model. For the radiative transfer that sets the dust temperature – and therefore its vertical structure – also the stellar and dust properties are required.

The radial distribution of the dust in the disk is defined by the following parameters: the inner and outer radius (R_in and R_out) as well as the gap location R_gap,in and R_gap,out; the total dust mass (M_dust) and its distribution across the disk, the surface density profile (SDP); and finally an additional dust depletion factor for the inner disk f_inner. The stellar spectrum is described by a Kurucz model with effective temperature T_eff, luminosity L _∗ ⁴ and mass, and is set at a distance d. These parameters are summarized in Table 2.

Although the stellar parameters have recently been updated by Tatulli et al. (2011), our fit parameters provide an equally good fit to the stellar photosphere. Adopting the lower luminosity of 26 L_⊙ would have only a small impact on the model, as it requires a higher continuum opacity in the optical to achieve the same dust temperatures and fit the SED.

3.2. Dust model

We compute the optical properties of our grains using a distribution of hollow spheres (DHS, Min et al. 2005a), which has been shown to represent mid-infrared feature shapes very well (Juhász et al. 2010).

We will restrict ourselves to a model with three dust components that dominate the opacity from the optical to the far-infrared: amorphous silicates, crystalline forsterite and carbon. Amorphous silicates produce the 10 and 20 μm emission features which are very prominent in HD 100546. We will use the composition and optical constants of interstellar silicates towards the galactic center as derived by Min et al. (2007), which are summarized in Table 1. Crystalline forsterite dominates the mid infrared opacity only at wavelengths where it has strong features, most notably at 11.3, 23.5, 33.5 and 69 μm. We will use the optical constants from Suto et al. (2006), as these have been measured at different temperatures. At optical and near-infrared wavelengths, the opacity is dominated by continuum sources without spectral features in the mid and far-infrared. For this component we use carbon, with optical constants from Preibisch et al. (1993). Note that probably not all of the continuum-opacity sources are amorphous carbon, and that also different materials with similar optical properties could contribute, for example metallic iron.

Dust grains in protoplanetary disks are most likely in the form of mixed aggregates. Lacking a good effective medium theory for DHS, we calculate the opacities for every dust species from the optical constants as if they were separate particles, and then force them to be in thermal contact. Although this approach does not have as good a physical basis as an effective medium theory, it does allow an accurate representation of feature shapes. In addition, opacities do not have to be recalculated when the composition changes, making it easier to conduct large parameter studies and make models with spatially varying dust compositions, including temperature dependent opacities. The implications of the dust grains being in thermal contact will be discussed in Sect. 6.6.

Other dust properties – grain size range and composition – have to be constrained from the observations and are discussed in Sect. 4.2.

4. Disk structure

Before we can model the spatial distribution of forsterite, we need to construct a structural model for the disk that fits the available observations. The disk of HD 100546 has been well-studied, and the disk geometry is already well-constrained. In the next section we will describe these observational constraints. In Sect. 4.3 we fit the SED to fill in most of the remaining parameters of our disk model. The only parameter important for this study that can not be well constrained from the SED is the surface density profile (SDP). We will assume the same SDP as previous modeling attempts, which is consistent with VISIR imaging (see Appendix A.3). Figure 1 shows a Q band image of the final disk model, illustrating its geometry with an outer disk (Sect. 4.1.1), gap (Sect. 4.1.2) and inner disk (Sect. 4.1.3).

4.1. Observational constraints

4.1.1. Outer disk

The (outer) disk of HD 100546 was first resolved in scattered light by Augereau et al. (2001), and has since then been imaged at different wavelengths (e.g. Liu et al. 2003; Ardila et al. 2007). Its position angle and inclination have been measured using different techniques, indicating an inclination between 40 and 50 degrees from face-on, and a position angle east of north between 130 and 145 degrees (e.g. Pantin et al. 2000; Ardila et al. 2007). We adopt the inclination of 42 degrees and a position angle of 145 degrees.

The outer radius of the dust disk is difficult to determine. Scattered light is observed up to at least ten arcseconds (1000 AU at a distance of 103 pc, Ardila et al. 2007), yet it is not clear whether this belongs to the disk or to a remnant envelope, and the dust disk could be significantly smaller (e.g. Ardila et al. 2007). Modelling of rotational CO lines are consistent with a gas-rich disk with an outer radius of 400 AU (Panić et al. 2010). We therefore adopt an outer radius of 1000 AU, but with an exponential cutoff that sets in at R_exp = 350 AU (Hughes et al. 2008).

The total dust mass of the outer disk is 7.2  × 10^-4 M_⊙, as has been determined by Henning et al. (1998) from 1.3 mm measurements, which probe mostly big grains in the disk midplane. The mass distribution in the outer disk is not well constrained. Various authors have adopted an SDP of Σ(r) = r^−p with p = 1 (e.g. Panić et al. 2010; Benisty et al. 2010), a value commonly used for protoplanetary disks. We will also adopt this value, which is consistent with our VISIR image in Appendix A.3

4.1.2. Gap size

The SED of HD 100546 (Fig. 2) is dominated by a relatively small near-infrared excess between 2 and 8 μm, followed by a steep rise towards the mid-infrared at 10 μm. This SED has been attributed to a disk gap around 10 AU – possibly cleared out by a proto-planet – which creates a frontally illuminated wall at the far edge of the gap (Bouwman et al. 2003). This wall has a temperature of 200 K, and is responsible for the huge mid-infrared excess compared to other Herbig Ae/Be stars. The gap manifests itself in mid-infrared imaging, where the source is much larger than other Herbig Ae/Be stars, and its size is independent of wavelength (Liu et al. 2003).

The gap is also seen in the distribution of the circumstellar gas (e.g. Acke & van den Ancker 2006; Brittain et al. 2009; van der Plas et al. 2009), and tight constraints have been put on the size of the disk gap using spatially resolved spectroscopy by Grady et al. (2005). We adopt their value for the gap outer radius of R_gap,out = 13 AU, which is the location of the disk wall.

4.1.3. Inner disk

The inner disk of HD 100546 causes a near-infrared excess that is relatively small compared to other Herbig Ae/Be stars (Bouwman et al. 2003), indicating that the dust is depleted (Dominik et al. 2003) or settled. Benisty et al. (2010) have resolved the disk using near-infrared interferometry, constraining the inner radius to R_in = 0.25 AU. The authors also find that the inner disk is very tenuous, and has a low dust mass of only 3  × 10^-10 M_⊙. The AMBER observations cannot resolve the outer radius, but they put it at R_gap,in = 4 AU. We adopt these values for the inner and outer radius.

4.2. Mineralogy

Before determining the remaining model parameters using SED fitting we need to specify our dust composition. Although the peculiar forsterite mineralogy of HD 100546 is the main focus of this paper, forsterite itself plays only a minor role in determining the disk structure and fitting the SED. The main reason for this is that it has a very low opacity in the optical compared to other dust species such as amorphous silicates and carbon. The latter species therefore determine the heating/cooling balance of the dust, and set the disk’s temperature and vertical density structure. Leaving out the forsterite does not affect the shape of the modeled Q band radial profile, and changes to the overall SED are small (Fig. 2). The shape of the 10 micron feature is better fitted with forsterite included, but this does not affect the overall SED shape (see Fig. 2). We will therefore constrain the disk geometry in this section using only the other two components, amorphous silicates and carbon. By doing so, the derived geometry is independent of the forsterite mineralogy, and we can study the abundance and spatial distribution of forsterite in Sect. 5 without refitting the SED.

The dust composition is determined by the mass fraction of carbon, and follows from our SED fit in Sect. 4.3. Both Bouwman et al. (2003) and Juhász et al. (2010) show that the small grains responsible for the mid-infrared dust emission features in HD 100546 cannot be larger than a few μm, though larger grains could be present at higher optical depths not probed by Spitzer. The small grains are most important for setting the disk density and temperature structure because of their high opacities (Meijer et al. 2008), and we adopt a grain size distribution between 0.1 and 1.5 μm for all grains. We do not fit the SED in the (sub)mm because it probes a separate population of large, millimeter-sized grains (Benisty et al. 2010), which do not influence the strength of mid-infrared features. A summary of the mineralogy can be found in Table 1.

4.3. SED fitting

With the disks radii (R_in, R_gap,in, R_gap,out, R_out) and opacities constrained, the next step is to constrain the dust mass (M_disk), inner disk depletion (f_inner) and carbon fraction by SED fitting. The SED is relatively insensitive to the SDP, but some limits can be placed on it using spatially resolved VISIR imaging (see Appendix A). For now, we will adopt p = 1. The total dust mass for HD 100546 has been measured at millimeter wavelengths (Henning et al. 1998). However this mass probes mostly the big, millimeter sized grains in the midplane, while we are interested in the dust mass in small micron-sized grains responsible for mid-infrared and far-infrared features. We therefore have to refit the dust mass to the optically thin part of the SED at  ~60–200 μm. We find a dust mass of 10^-4  M_⊙, consistent with the dust mass in small grains from Dominik et al. (2003) and Benisty et al. (2010).

We proceed by fitting the inner disk (0.25–4 AU), which has been done before with a similar 2D radiative transfer code by Benisty et al. (2010). They use a parametrized vertical structure with a fixed scale height and flaring exponent. In our model, the vertical structure is computed under the assumption of vertical hydrostatic equilibrium. Therefore we can not freely adjust these parameters, so instead we fit the inner disk SED by varying its dust mass. We do this by assuming the mass distribution in the inner disk follows the same power law as the outer disk (p = 1) – as would be the case for a disk without a gap – and reduce the surface density⁵ to lower the optical depth and fit the near-infrared flux in the same way as Benisty et al. (2010). We fit a dust mass of 6  × 10^-9 M_⊙, consistent with the value found by Benisty et al. (2010). This corresponds to a surface density normalization of Σ_1   AU = 2.3  × 10^-3 g/cm² . This is a reduction of a factor f_inner = 200 with respect to an extrapolation inward of the SDP of the outer disk, which has Σ_1   AU = 0.45 g/cm².

The modelled inner disk is not completely optically thin in the radial direction at optical wavelengths, as can be inferred from the shadow seen on the disk wall in Fig. 1. Even when we take into account that the star is not a point source (Espaillat et al. 2010), the shadow reaches up to a height above the midplane of around  ~1 AU at the disk wall. The shadow only disappears when the inner disk has an increased scale height such that it effectively becomes a halo (Mulders et al. 2010; Verhoeff et al. 2011), as has been suggested for HD 100546 by Vinković et al. (2006).

Even though the wall at 13 AU is partly shadowed, the illuminated part still intercepts a significant fraction of the stellar light. It is frontally illuminated by the star, and has a temperature of 200 K. Such a structure fits the mid-infrared part of the SED without any increase in scale height beyond hydrostatic equilibrium – which is often needed to account for extreme near-infrared excesses in Herbig Ae/Be stars (Acke et al. 2009). The steep rise in the SED between 8 and 10 μm is caused by this wall, but it also contains information on the dust composition. At 8 μm the opacity is dominated by carbon, at 10 μm by amorphous silicates. From the 10/8 μm flux ratio, we estimate a carbon/silicate mass ratio of 5%. Increasing the carbon fraction “fills” the gap in the SED at 8 μm: with 10% carbon the flux at this wavelength is already significantly overpredicted. Because carbon combusts (see, e.g., Lee et al. 2010) at the same temperatures needed to anneal silicates (Gail 2001), the relatively low carbon abundance could be linked to the high crystallinity.

4.4. Final model

Model parameters of the final disk model used for the analysis of the forsterite features can be found in Table 2. Note that this model is in hydrostatic equilibrium, and therefore does not have a predefined scale height. The calculated scale height increases with height above the midplane due to the temperature gradient, but can be fitted with a power law in radius (H = H_100   AU × (r/100   AU)^β) at both the surface and midplane for comparison. The outer disk is well described by a flaring exponent of β = 1.4 and a scale height that increases from H_100   AU = 13 AU in the midplane to H_100   AU = 21 AU at the disk surface and in the disk wall. The inner disk is fitted with β = 1.3 and H_100   AU = 13 AU in the midplane to H_100   AU = 25 AU in the disk rim. The midplane temperature follows T =  ~150 × (r/AU)^-0.25 for the outer disk, and is a factor of 2 higher in the inner disk.

5. Disk mineralogical model: Forsterite

In this section we will study the abundance and spatial distribution of forsterite in the disk of HD 100546, using the disk structural model described in the previous section. The aim is to distinguish between the two scenarios for the forsterite composition described by Sturm et al. (2010) (warm iron-free forsterite or cold iron-containing forsterite), which have also been discussed in the introduction.

5.1. Optical depth effects

We explore in detail the effects of optical depth and the gapped disk geometry of HD 100546 on the observed feature strengths of forsterite. Throughout the paper, we will assume that the forsterite is in thermal contact with the other dust species (see Sect. 4.2). Figure 3 shows the temperature structure in the outer disk from the model derived in Sect. 4.4. Figure 3 also shows the vertical τ = 1 surfaces – as seen from the observer to the emission point – of the continuum at 24, 33 and 69 μm. These surfaces mark the boundaries below which the dust can no longer contribute substantially to the feature at that particular wavelength.

At 69 μm, the vertical τ = 1 surface lies below the disk midplane (i.e., the disk becomes optically thin), allowing most of the disk to contribute to the 69 μm feature. At shorter wavelengths, these surfaces lie above the disk midplane out to radii of 50–100 AU. For the 150 and 200 K components, a significant part lies below the τ = 1 surface, especially when taking into account that the density increases towards the midplane. Most of the warm forsterite that can contribute to the 69 μm feature is therefore hidden from view at shorter wavelengths.

To illustrate the influence of these optical depth effects on the measured feature strengths of forsterite, we compare the continuum subtracted spectra (see Appendix B) of our model (described in Sect. 5.2.4) with the optically thin model from Sturm et al. (2010) in Fig. 4. Both models are fitted to the feature shape and strength at 69 μm, and contain only 150 K and 200 K forsterite. The optically thin approach clearly overpredicts the strength of the short wavelength features in the 10 to 40 μm range, while the 2D radiative transfer does not overpredict these because the bulk of the warm forsterite is hidden at such short wavelengths.

The gapped geometry of HD 100546 plays a key role here: because part of the wall at 13 AU is directly illuminated by the central star, it heats regions close to the midplane to temperatures of 150 to 200 K. In a disk without a gap, these regions would be much colder and these high temperature regions would only be found in the midplane at much smaller radii (~4 AU). There, the surface density would be high enough also for the τ = 1 surface at 69 μm to be above the midplane. Without a gap, the optical depth effects would not be so pronounced as the hot forsterite would be hidden from view at 69 μm as well.

This first analysis shows that while optically thin models with iron-free forsterite overpredict feature strengths at shorter wavelengths, a 2D radiative transfer model that takes optical depth effects into account does not. Additionally, none of the longer wavelength features (33, 69 μm) require an additional cold (<150 K) component to fit the feature strengths. This indicates that the shape of the 69 μm feature can be fitted with only warm forsterite, suggesting an iron-free composition (the first scenario from Sturm et al. 2010). We will therefore first pursue a disk model that explains all feature strengths and shapes with iron-free forsterite. Afterwards, we will discuss models that include a fraction of iron. We describe how we implement the temperature-dependent opacities of iron-free forsterite into our radiative transfer code in Appendix C.

5.2. Forsterite spatial distribution

To determine the spatial distribution of the forsterite in the outer disk, we use the disk structural model from Sect. 4.4. We replace a small fraction of the amorphous silicates by forsterite, and compare continuum subtracted model spectra with the data. All discussed models can be found in Fig. 5. These models are fitted to the integrated flux at 69 μm, which is  ~9 Jy μm.

We will show that models with constant abundance throughout the outer disk severely underpredict short wavelength features due to optical depth effects (model G0). We therefore run an extensive parameter study where the radial distribution of forsterite peaks towards the disk wall. We consider two parametrizations: an abundance gradient (model series G), predicted by radial mixing models (Gail 2001), and a (narrow) ring of constant abundance (model series R), which makes it easier to control the minimum temperature of the forsterite. In addition, we also ran a series of models where the forsterite is only located in the disk wall and at the disk surface (model RW), and not in the midplane. Below we discuss the different models.

5.2.1. Constant abundance throughout the disk

In contrast with the result from Bouwman et al. (2003), models with a constant forsterite abundance throughout the outer disk consistently fail to reproduce all observed feature strengths. This failure is due to optical depth effects, which suppress short wavelength features (Sect. 5.1). A very low abundance of 2% fits the integrated 69 μm flux (Fig. 5, model G0), but underpredicts shorter wavelength features, most notably shortwards of 33 μm. An abundance increase of at least a factor of 10 is necessary in the disk wall to fit the flux at 11 μm. To achieve this without overpredicting the 69 μm feature by the same factor, we need to concentrate the forsterite towards the disk wall.

5.2.2. Abundance gradient

First, we concentrate the forsterite towards the disk wall by introducing an abundance gradient that gradually decreases outward. Radial mixing models predict that if radial mixing transports crystalline material from the inner to the outer disk, it follows a smooth distribution (Gail 2001), which is well described by a power law (Stevenson 1990). We parametrize the forsterite abundance as f(r) = f_13   AU   (r [AU] /13   AU)^−ξ, where free fitting parameters are the abundance in the disk wall (f_13   AU) and the steepness of the gradient (ξ). The grid scanned for these parameters is summarized in Table 3. We start by assuming the forsterite is well-mixed in the vertical direction, i.e. there is no vertical abundance gradient.

With a gradient of ξ = 1, we can increase the forsterite abundance in the wall to 10% without overpredicting the 69 μm flux. This is still not high enough to explain the short wavelength features, which originate mainly in the disk wall (G1). Steeper gradients of ξ = 2 or 3 concentrate the forsterite closer to the disk wall, and can reach higher abundances for the same 69 μm flux. They reach abundances of respectively 30% (G2) and 40% (G3), which is sufficient to fit the feature strengths at shorter wavelengths. With even steeper gradients, the abundance drops off so rapidly even in the wall that short wavelength features are again underpredicted (not plotted).

Although models G2 and G3 can explain all feature strengths, the shape of the 69 μm feature does not fit the observations well. The peak is located at 69.0 μm, indicating that cold, 100 K forsterite dominates the peak location. Although an increasing contribution at 69.2 μm is seen with steeper gradients, it never dominates the feature shape. These solutions therefore require some iron to shift the feature towards longer wavelengths. However, even then the feature shapes do not match very well (Fig. 5).

The reason why the cold forsterite is so prominent at 69 μm is that even with a steep gradient a few percent forsterite remains far out in the disk. Such a low abundance provides a significant contribution at 69 μm, while it does not add a lot of flux to the short wavelength features (see also model G0). We will therefore focus on a series of models with a forsterite abundance that is non-zero only near the disk wall.

5.2.3. Abundance ring

To prevent a low percentage of cold forsterite in the outer disk from dominating the shape at 69 μm, we choose a second parametrization where the forsterite is located in a narrow ring starting at the disk wall. The abundance of forsterite is constant within the ring, between 13 and r_ring AU, and is zero outside this ring. Both the abundance in the ring and its outer radius r_ring are free fitting parameters. The explored parameter space is summarized in Table 3, a selection of models is displayed in Fig. 5.

Because the midplane is close to optically thin at 69 μm, avoiding a strong contribution from cold (50–100 K) forsterite requires a narrow ring of only a few AU. A ring with outer radius of 14 AU (R14) requires a high abundance of 60% to fit the integrated feature strength at 69 μm. At such high abundances, features become optically thick and flatten off, producing a 69 μm feature that is much too broad. Additionally, the short wavelength features flatten off, and become too weak compared to the observations.

Lower abundances require broader rings to produce enough emission at all wavelengths. Rings with abundances of 40 and 30% require outer radii of 20 (R20) and 25 AU (R25) respectively. At these radii, midplane temperatures are down to 100 K (Fig. 5), and the 69 μm feature is dominated by this component with a peak at 69 μm. In contrast to the gradient models G1 and G2, these rings do not contain a 50 K component, since these temperatures occur farther out in the disk. Therefore the feature shape matches much better than models with a gradient, although still some iron is required to shift the feature towards longer wavelength. At even lower abundances (R60, 10% and 60 AU), the feature strengths at shorter wavelengths cannot be reproduced.

5.2.4. Only warm forsterite

In order to reproduce the 69 μm feature with iron-free forsterite, we need to exclude contributions from 100 K forsterite. However, ring models that do so (R14) do not produce enough flux at all wavelengths because the emitting region is too small. Larger rings include warm forsterite from the disk surface, but also colder forsterite from the midplane (R20 and R25). We therefore run the same series of ring models, but now with forsterite only in regions that are 150 and 200 K, disabling the contribution from colder forsterite. This model is easily compared to the optically thin fit of Sturm et al. (2010) (Fig. 4), which also has only these components. We implement this by replacing the opacities for 50, 100 and 295 K forsterite with that of amorphous silicates. A ring with an outer radius of 30 AU and abundance of 60% provides a good fit to the short wavelength features (RW30). The shape of the 69 μm feature matches the observed one well.

Looking at Fig. 3, it becomes clear that such a model with only warm forsterite has an abundance that varies with height above the midplane. Outside of the disk wall (>15 AU), a vertical temperature gradient is present from a warm surface to a cold midplane. The forsterite in this region must therefore be located only at the disk surface, with an abundance that decreases to zero near the disks midplane. The forsterite is present in the disk wall and near the disk surface farther out, but not in the midplane. This scenario agrees well with the geometry proposed by Bouwman et al. (2003), where forsterite is produced in the disk gap and transported into the disk wall and outer disk surface.

5.2.5. Best fit to 69 μm feature

The parameter study from the previous section yields a number of good fits to the integrated feature strengths, but only 2 models fit the shape of the 69 μm feature well: model R20, with 40% forsterite between 13 and 20 AU and model RW30, with 60% forsterite between 13 and 30 AU, but only in the “warm” regions with a temperature of 150 and 200 K. Figure 6 shows the modeled feature shape, as well as the contribution of the individual temperature components.

Model R20 peaks at 69.10 μm, and the model has been shifted by 0.1 μm to obtain a good fit. We interpret this shift as an increase in the iron fraction from 0 to 0.3%, following the interpolation for low iron content by Sturm et al. (2010), which gives a peak shift of 3.13%/μm. The peak shows dominant contributions from both 200, 150 and 100 K, and is therefore slightly broader than the optically thin fit from Sturm et al. (2010), but well within the error bars. The higher temperature components dominate the feature by integrated emission, but the 100 K component is narrower and is equally important for the peak position.

Model RW30 peaks at 69.15 μm, and provides a good by-eye fit without a wavelength shift. Fit quality increases with an additional shift of 0.05 μm, corresponding to an iron fraction of less than 0.2%. The relative contribution of the 200 and 150 K components are 57 and 42% in contrast with 62 and 27% for the optically thin model (Sturm, priv. comm.). This also explains the small wavelength shift required, as the 150 K component has its peak at 69.10 μm.

The feature of model RW30 is also a little bit too broad, especially at the blue side of the feature. The reason for this broad feature is the high forsterite abundance of 60%, at which the feature itself starts to become optically thick. A lower abundance would require a larger region with a temperature of 200 K, which is difficult to achieve in the current model.

5.3. Location of forsterite emission

To visualize where the different forsterite features originate in the disk, we calculate the contribution to the features at every location in the disk. Although the temperature and optical depth also give some information about where the features originate (Fig. 3), it does not take into account that the density decreases with disk height and radius, as well as the increase in opacity in the wavelength of the feature and inclination effects.

To calculate these contributions, the ray tracer stores the individual contributions to the total spectrum in every grid cell. We then fit a 3rd order polynomial to the continuum, subtract it, and calculate the integrated feature strength. The contributions of the individual temperature components to the 69 μm feature are calculated in a similar fashion.

Because we view the disk at a 42° angle, the upper half of the disk (z/r > 0) is in our direct line of sight, while the bottom half (z/r < 0) is mostly obscured from view. This results in an asymmetric profile around the midplane (z/r = 0) in Figs. 7 and 8, while the forsterite abundance itself is symmetric around the midplane in all our models.

5.3.1. Short wavelength features

The contributions to the integrated flux of the features at 11.3, 23.5, 33.5 and 69 μm are plotted in Fig. 7 for best-fit models R20 and RW30. All features have a strong contribution from the disk wall, approximately half of the total integrated feature strength at all wavelengths. The remaining emission originates at larger radii at or near the disk surface. Features at longer wavelength can be produced by colder material, and have stronger contributions from larger radii. Due to the lower continuum opacity at longer wavelengths, features also originate from deeper within the disk, especially at larger radii where surface densities become lower.

The shortest wavelength feature at 11.3 μm originates in the warmest regions of the outer disk, mostly in the disk wall. Its surface component extends up to only a few AU outwards of the disk wall, and is located at the disk optical surface for both models. At 24 and 33 μm, the emission from the disk surface extends farther out, and its contribution to the total features increases up to more than 50%. As the optical depth of the continuum decreases at these wavelengths, features start originating from deeper in the disk. Both models look similar at 24 μm but start to deviate from each other at 33 μm, as here the cold 100 K forsterite component is not present, truncating most of the emission after 17–18 AU.

At 69 μm, the optical depth of the continuum is so low that features originate from regions all the way down to the midplane, even in regions just behind the disk wall (14 AU). The total contribution from outside the disk wall is more than 50% (see also Fig. 9), and peaks towards the midplane with only a very small contribution from the disk optical surface. The model with only warm forsterite doesn’t extend as far out because the cold forsterite component in the midplane is not present.

The detailed analysis of the feature location confirms the simple analysis from the beginning of this section: towards longer wavelengths, the optical depth decreases and features originate from regions deeper in the disk which are hidden from view at shorter wavelengths. Compared to Fig. 3, regions of emission don’t reach down as far as the vertical τ = 1 surface. The main reason is that disk’s inclination is (by definition) not taken into account in calculating the vertical τ = 1 surface. When viewing the disk at an inclination angle i, the τ = 1 surface of an inclined disk is raised by approximately a factor  ~1/cos(i) with respect to a face-on disk.

5.3.2. The 69 μm feature

The location of emission of the different temperature components of the 69 μm feature are displayed in Fig. 8 for model R20 and RW30. These are calculated similarly to the shorter wavelength features – continuum subtracting the spectrum in every grid cell – but for each of the five temperature components separately. For model R20, the cumulative radial contribution to the 69 μm feature is plotted in Fig. 9.

For model R20, the 295 K component originates from the disk wall (between 13 and 14 AU), except from near the midplane. This is because the temperature in the midplane is slightly lower than at the surface, an effect of the shadow of the inner disk (see Fig. 1). The 200 K component originates only for a small part in the disk wall. Most of the emission comes from the upper half of the disk between 13 and 15 AU, and peaks below the disk surface. At 150 K, a fraction of the emission comes from the shadowed part of the disk wall, and the bulk from near the midplane below the surface out to 17 AU. The 100 K component does not originate in the disk wall, but from the midplane outwards of 14 AU. The coldest, 50 K component originates from even larger radii(>18 AU) in the midplane.

For Model RW30, the 200 K and 150 K components originate in more or less the same regions as for model R20, while the 50 K, 100 K and 295 K components are not present.

6. Discussion

6.1. Forsterite abundance in comparison with previous work

The silicate mineralogy of HD 100546 has been addressed in a number of studies (e.g. Bouwman et al. 2003; Juhász et al. 2010), with the underlying assumption that emission features arise in an optically thin surface layer in the disk. Assuming a homogeneous composition for the emitting region, the mass fraction of forsterite in small dust grains is found to be around 5–6% (van Boekel et al. 2005; Juhász et al. 2010). Because this number can be derived from the mid-infrared spectrum directly by measuring crystalline and amorphous feature strengths, we will refer to this as the apparent crystallinity. It is clear that this number is not neccessarily representative for the mass ratio in the entire disk, but an extrapolation is of course tempting. Taking a separate dust composition for the inner and outer disk, Bouwman et al. (2003) find a forsterite abundance in the depleted inner disk of 2%, compared to 19% in the outer disk. Since the latter dominates the total mass, the total mass fraction is also 19%, a factor of 4 higher than the apparent crystallinity.

Our analysis shows that the crystalline forsterite is concentrated towards the disk wall, with locally very high abundances (40–60%). Although uncertainties in the grain properties and model assumptions can lead to small errors in the derived abundance, the abundances we find are significantly higher than those quoted above. These high abundances are confined to a very narrow region (from 13 AU up to 20–30 AU). Outside of this region – where the bulk of the mass is located – the forsterite abundance is consistent with zero. Therefore, the overall forsterite mass fraction – the mass in crystalline forsterite grains compared to the total disk mass in small dust grains – is around 0.5–0.8% (see Table 4). This is a factor of 10 lower than quoted above.

The reason why such a small amount of forsterite shows up so prominently in the SED is a typical case of window-dressing: the wall reprocesses approximately a third of the total stellar energy captured by the disk, but represents only a few percent of the total disk mass. The wall therefore acts as a display case, showing a strong apparent crystallinity from only a small amount of forsterite with a high local abundance. This may also explain the lack of correlations found in larger samples between apparent crystallinity and other disk parameters.

6.2. Forsterite abundance in comparison with solar system comets

It has been noted before that the mid-infrared spectrum of comet Hale-Bopp shows a striking resemblance to that of HD 100546 (e.g. Malfait et al. 1998). The derived crystallinities are in the range of 25% to 50% (Brucato et al. 1999; Wooden et al. 1999; Harker et al. 2002). It should be noted that the degree of crystallinity depends on the assumed dust model, see e.g., Hanner & Zolensky (2010). With the inclusion of larger grains (>μm), crystallinities tend to be lower, around 7% for Hale-Bopp (e.g. Min et al. 2005b; Bouwman et al. 2003). Because we measure the crystallinity of small dust grains in HD 100546, we can compare the local abundance in the disk wall of 40–60% with the crystallinity of small grains in comets.

Also other longer period comets with high crystallinities have been observed, for example COMET C/2001 Q4 has a crystallinity of  ~30% (Wooden et al. 2004). The surface of shorter period comets might be subject to amorphization by solar irradiation, but ejecta from the inside comet 9P/Tempel released during the Deep Impact experiment show that these comets can also have a high crystallinity of 20–40% (Harker et al. 2007).

Such high crystallinities at these large radii are not expected from radial mixing models ( 2004), and not observed in typical protoplanetary disks. Because the proposed formation region of comets is beyond the snowline, it coincides with the location of increased crystallinity in the disk wall of HD 100546, making it an ideal place for comet formation.

6.3. Crystallinity-gap correlation

Both the high abundance and spatial distribution indicate a correlation between the disk geometry and the strong crystalline silicate features in the disk of HD 100546. If the disk gap were filled with amorphous material – so that an illuminated wall is no longer present – the spectral features of forsterite would be weaker by roughly a factor of 3. This might indicate that the strength of crystalline features – the apparent crystallinity – does not reflect the overall crystal abundance, but rather its spatial distribution and the presence of a disk gap.

Another star which shows a similar geometry and apparent crystallinity is RECX 5 (Bouwman et al. 2010). This opens up the possibility that a disk gap-crystallinity correlation also exists for other protoplanetary disks. However, not all disks with gaps show crystalline dust features (Brown et al. 2007). A reason for this might be that the crystalline features only show up if the gap is located in the right temperature range. A good census of crystallinity and gaps in the right radius/temperature range is needed to establish the presence of such a correlation. The strength of crystalline emission features can be well determined from mid-infrared spectroscopy (e.g. Juhász et al. 2010). But the characterization of disk gaps from the SED is less trivial, and requires spatially resolved observations (Varnière et al. 2006).

Aside from the apparent crystallinity of HD 100546 – which is one of the highest among Herbig Ae stars (e.g. Juhász et al. 2010) – also the local abundance of forsterite at tens of AU is unusually high (40–60% at 13–20 AU). Other Herbig Ae/Be stars show lower forsterite abundances in the outer (2–20 AU) disk, around 10–15% derived with an optically thin approach (van Boekel et al. 2004), with the exception of HD142527, which also has a prominent disk gap (Verhoeff et al. 2011). If all Herbig Ae/Be stars would have the same crystallinity as HD 100546 in the 13–20 AU range – and an equal or higher abundance farther in – they would produce much stronger spectral features than observed. Only if the abundance would increase outwards would it be possible to produce the observed crystallinities in the 0–5% range with the same high local abundance at 13 AU. However, this is in sharp contrast with observations (van Boekel et al. 2004) and radial mixing models (Gail 2001, 2004), that show a crystallinity that is constant or decreases outwards.

Whether the spatial distribution of forsterite around HD 100546 is truly unique among Herbig Ae/Be stars requires an analysis of disks with lower crystallinity as well. The Herschel open time key program DIGIT will observe the 69 μm feature in a number of Herbig Ae/Be stars, allowing for a characterization of the spatial distribution of forsterite in these disks in a similar fashion. If there is a mechanism that produces a high crystallinity at 13–20 AU independent of a disk gap, it will be identified in this sample.

6.4. Forsterite origin

How can such a high abundance of forsterite originate in the disk wall, and is it related to the formation or presence of the disk gap? There are two chemical pathways of producing forsterite: by condensation from the gas phase or by annealing from amorphous olivines. Both require high temperatures above the dust evaporation temperature (~1500 K) or the glass temperature (~1000 K) respectively. These high dust temperatures exist only in the inner regions of protoplanetary disks, but not at 13 AU. We can think of two main scenarios, radial transport or in-situ production.

In the first scenario, forsterite is produced in the inner disk, and transported outwards to 13 AU. The current dust mass in the inner disk is much lower than the total amount of forsterite in the disk wall between 13 and 20 AU. It therefore seems unlikely that radial mixing from the current inner disk can produce such a high crystallinity in the disk wall. As discussed before, it is unlikely that radial mixing has produced a locally high abundance before the opening of the disk gap, as it would have shown up in the mid-infrared spectra. This means that radial transport must have taken place during the formation of the disk gap and depletion of the inner disk. Radiation pressure on the crystalline grains can transport them across the gap and blow them into the disk wall, depleting the inner disk. Also a different mechanism, photophoresis, can transport grains outwards efficiently when the inner disk depletes (Petit et al. 2006). The depleted dust mass from the inner disk (200 times the current inner disk mass, Sect. 4.3) is large enough to explain the crystalline silicates in the disk wall, but it requires extremely efficient radial transport and almost a 100% crystallinity in the inner 4 AU before the opening of the gap.

In the second scenario, forsterite is produced locally in the disk wall. Because the equilibrium temperature of the dust in the wall (200 K) is far below the dust evaporation or glass temperature, another mechanism must then be responsible for heating the amorphous dust grains so they can crystallize:

• Flash heating of grains in a stellar outburst (Ábrahámet al. 2009). However, HD 100546shows currently no signs of such outbursts. Such crystals areobserved to disappear quickly after formation, which is not thecase for HD 100546, where crystals persist overmore than a decade of observations.

• Shocks in the disk midplane can locally heat the grains above the glass temperature (Harker & Desch 2002). From scattered light images, HD 100546 is known to possess two spiral arms in the outer disk (Grady et al. 2001), which might be related to shocks in the disk midplane that could create forsterite (Harker & Desch 2002). If a companion is responsible for clearing the disk gap, it likely induces a similar spiral wave at the location of the disk wall (e.g. Nelson & Papaloizou 2003), producing forsterite locally.

• Amorphous silicates can be crystallized in parent bodies large enough to retain the heat from decay of radioactive elements or accretion (e.g. Huss et al. 2001). Micron-size forsterite grains can be created in a collisional cascade of such objects, that can take place in the disk gap or wall (see Bouwman et al. 2003 for a discussion on this topic.) The amount of forsterite found in the disk wall of HD 100546 requires the destruction of an Earth-like planet or several hundred Pluto-sized objects.

All plausible scenarios for forsterite formation (shock waves, parent body processing and radial transport during gap formation) are connected to the presence or formation of the disk gap. Therefore we conclude that the anomalously high forsterite abundance in the disk wall of HD 100546 is – one way or the other – connected to the presence of the gap.

6.5. Forsterite in the inner disk

The mass fraction of forsterite in the inner disk is hard to constrain from our models. The contribution to mid-infrared features in the model is low, and shows a significant contribution only at 11 μm. Regions near the inner rim that are above a temperature of 1000 K should be fully crystalline, but even if the silicates in the entire inner disk are fully crystalline, the contribution is only a few Jy μm at 11 μm. The contribution to longer wavelengths is negligible because the inner disk is too warm, and optical depth effects do not play a role because it is also (vertically) optically thin (Benisty et al. 2010). Therefore we cannot put any limits on the crystallinity of the inner disk with the current model and observations.

6.6. Thermal contact

We assume that the different dust species in our model (amorphous silicates, forsterite, carbon) are in thermal contact. However, Bouwman et al. (2003) found with an optically thin approximation that the bulk of the forsterite is at a temperature of  ~70 K while the other dust species are much warmer, leading them to the conclusion that forsterite cannot be in thermal contact with the other dust species. Because the 2D radiative transfer does not need this cold component due to optical depth effects (see Sect. 5.1), the forsterite has the same temperature as the other dust species in the disk wall, which means they can be in thermal contact.

If the forsterite in the outer disk is forced to be out of thermal contact with the other dust species, it does not reach temperatures high enough to explain the shortest wavelength features of forsterite. The reason for this is the very low opacity of forsterite at optical wavelengths compared to the other dust species, especially carbon. Therefore the forsterite is not heated directly by the central star, but only by mid-infrared radiation from the disk surface. It will therefore have a temperature similar to dust species in the midplane (50–100 K), which is too cold to explain any short-wavelength features. Even if the forsterite is moved in to the inner disk, where it would normally be hotter, it falls short of explaining the strength of the short wavelength features (see also Sect. 6.5).

We can therefore exclude that the dust species are not in thermal contact. Bouwman et al. (2003) argue that the dust in HD 100546 is more primitive than in comet Hale-Bopp because it is out of thermal contact, and that we see the creation of such comets rather than their destruction. Since we exclude that the dust species are not in thermal contact, this implies that the dust in HD 100546 is more similar to comet Hale-Bopp than thought before, and not necessarily more primitive. This opens up again the possibility that we are witnessing the destruction of forsterite-rich comets like Hale-Bopp in HD 100546, rather than their formation.

6.7. Iron content

Finally, we use the shift in wavelength required to fit the shape of the 69 μm feature to estimate the iron content of the silicates. For large changes in the iron fraction, the central wavelength shifts linearly with iron content (Koike et al. 2003). The wavelength shifts between 0 and 10% have not been measured, but relies on the interpolation by Sturm et al. (2010), which assumes the shift in this regime is also linear.

We also assume that the shape of the 69 μm feature does not change with a small change in iron content. Although this has not been measured in the lab for forsterite with low iron fractions, small changes in iron fraction at higher iron content have been measured (Koike et al. 2003). For example, between 9.3 and 8.6%, the shape does not change significantly. The feature broadens a little bit between 0 and 8% iron content, but this change is much smaller than the broadening due to a temperature change (see Fig. C.1).

Therefore we expect the shape of the feature not to change significantly for the small changes in iron content (<0.3%) found in this paper. Laboratory measurements of forsterite with low iron content are necessary to confirm this.

7. Conclusion

We have studied the spatial distribution of crystalline silicate dust (forsterite) in the gapped disk of Herbig Be star HD 100546. We have used a 2D radiative transfer code to model the dust geometry, and study the effects of optical depth on the forsterite feature strengths in the mid and far-infrared, as well as the shape of the temperature-dependent 69 μm feature observed with Herschel/PACS. Our conclusions are:

• 1.

  2D radiative transfer effects need to be taken into account when interpreting spectral features over a broad wavelength range. In the case of HD 100546, the bulk of the warm (200–150 K) forsterite that contributes to the 69 μm feature is hidden from view at shorter wavelengths (10–40 μm). An additional cold forsterite component (~50–70 K) is therefore not necessary to explain the strength of the 69 μm feature.

• 2

  The forsterite of HD 100546 is located in the disk wall at 13 AU, with extremely high local abundance. Based on mid-infrared features strengths and the shape of the 69 μm feature, we arrive at two best-fit models:

  • a)

    a model with 40% forsterite between 13 and20 AU.

  • b)

    a model with 60% forsterite between 13 and 30 AU but restricted to regions where the temperatures are around 150 and 200 K. This excludes the cold midplane outside of the disk wall.

• 3.

  All forsterite features show a strong contribution from the disk wall, but also probe a region farther out. The 69 μm feature from this region originates from below the disk surface. The shorter wavelength features from this region originate at the disk surface.

• 4.

  The mass fraction of forsterite compared to the total disk mass is low – around 0.5–0.8% – a factor of 10 lower than the apparent crystallinity derived from and optically thin analysis of relative feature strengths in the mid-infrared. The wall at 13 AU acts as a display case for showing a high local forsterite abundance.

• 5.

  The crystalline forsterite in the disk of HD 100546 is in thermal contact with the other dust species.

• 6.

  We cannot rule out forsterite that is completely iron free. Our best fits indicate an iron content of  <0.3% of the crystalline olivines.

• 7.

  Due to a degeneracy in the modelling between iron content and dust temperature, it is not straightforward to use the 69 μm feature as an independent temperature indicator.

• 8

  A disk model in hydrostatic equilibrium (with T_gas = T_dust) can self-consistently explain the SED and mid-infrared imaging. Our best fit has a surface density power law  ∝ r^-1 and an inner disk depletion factor of 200. Additionally, continuum opacity sources – such as carbon – must have a low abundance to fit the low flux at 8 μm.

• 9.

  The apparent crystallinity of dust as measured from the ratio of crystalline and amorphous components of the mid-infrared features is not necessarily a good measure of the real crystallinity, the mass fraction of small dust grains in the disk that is crystalline. Both disk geometry and spatial distribution of the crystals can significantly distort the picture, and a one-to-one correlation between crystal feature strength and disk state is not to be expected. However, in our case, the detailed analysis shows a highly localized region with high abundance, which may very well inform us about special processing in this region, possibly linked to planet formation.

Online material

Appendix A: VISIR imaging

During the night of March 26 2005, Q-band imaging was performed using the VLT Imager and Spectrometer for the mid-IR (VISIR, Lagage et al. 2004). For photometric calibration and PSF determination, the standard stars HD 50310 and HD 150798 were observed before and after the science observation respectively. Standard chopping and nodding was employed to remove the atmospheric background emission. The imaging was performed with the Q2 filter (λ_c = 18.72 μm) in the small field mode (pixel field of view = 0.075″). The observing conditions were fair, airmass  ~1.4, and the optical seeing  ~0.8″. The achieved sensitivity was 80 before and  ~100 mJy/10σ1h after the science observation. The full width at half maximum (FWHM) of the point spread function (PSF) was  ~0.5″. The data reduction was performed with a dedicated pipeline, which corrects for various instrumental signatures (see Pantin et al. 2008, 2009; Pantin 2010). The reduced image is presented in Fig. A.1

A.2. Image analysis

The VISIR image (Fig. A.1) shows a resolved disk that is almost completely spherical, consistent with previous imaging at this wavelength (Liu et al. 2003). Deviations from spherical symmetry are down to the 10% level, which means that inclination must be small. Due to the relatively large size of the PSF at this wavelengths, it is difficult to get a direct constraint on the inclination out of the image. Comparison with a set of inclined, PSF-convolved model images show that it is consistent with an inclination less than  ~40 degrees and a position angle along the south-east to north-west.

Due to the increased sensitivity with respect to previous observations, the image shows resolved emission up to  ~1.4″, corresponding to a physical radius of 145 AU at 103 pc. To characterize this emission, we construct a surface brightness profile. Because the disk is not seen face-on, we construct the surface brightness profile from the image by averaging the surface brightness over an elliptical annulus – rather than a circular one. The shape of the annuli correspond to a disk inclination of 42° and a position angle of 145°. The semi-major axis of this ellipse corresponds to the distance from the star at which the radiation is emitted (radius). Annuli are 0.075″ wide, the size of one pixel on the VISIR chip. They are thus much smaller than the width of the VISIR PSF.

Comparing the measured radial intensity profiles to that of the VISIR PSF⁶ (Fig. A.2), we see that the disk is clearly resolved. The profile itself is not easily characterized by simple power-laws as is often the case for scattered light images, due to the different emission mechanism (thermal versus scattered light) as well the relatively large size of the PSF.

Instead, we recognize two components. The profile up to 0.5″ is well described by a Gaussian with a FWHM of 0.73″. After quadratic subtraction of the PSF, we derive a spatial scale of 0.33″ for this central component, corresponding to a physical scale of FWHM = 34 AU at a distance of 103 pc. This number agrees well with the value of 34 ± 2 found by Liu et al. (2003) with direct imaging at the same wavelength. This component therefore originates in a region just behind the disk wall, which has a spatial scale of  ~26 AU.

The second component – outside of 0.5″– comes from farther out in the disk, and has not been analysed before at this wavelength. Extended emission is seen up to 1.4″, corresponding to a physical radius of 145 AU at 103 pc. This emission probes the thermal continuum emission from the outer disk surface. The bump around 1.0″ corresponds to the first diffraction ring in the VISIR PSF, and does not represent a real structure in the disk. Because the radial profile traces the disk surface over almost an order of magnitude in radius, we can use it to put limits on the surface density of small grains, which we will do in the next section.

A.3. Limits on the surface density distribution of small grains

To test if the assumed power law index of p = 1 for the surface density profile matches our VISIR image, we compare model images of varying SDP with the observed image, keeping the dust mass and outer radius fixed, as well as the inner disk surface density. Because the pixel size is much smaller than the VISIR PSF, we convolve it with our model images before constructing the radial intensity profiles, and take into account the errors on the width of PSF (FWHM = 0.509 ± 0.008″) in the fitting procedure.

Because we solve for the hydrostatic disk structure, the shape of the radial profile depends only on the surface density distribution (Fig. A.2). The model is best fit for a surface density power law with index $\mathit{p}\mathrm{=}\mathrm{1.}{\mathrm{1}}_{-0.5}^{\mathrm{+}\mathrm{0.4}}$.

Appendix B: Continuum subtraction

To measure feature strengths and compare continuum subtracted spectra, we have to fit and subtract the continuum from both our model spectra and data. An example for model R20 is shown in Fig. B.1, for the Spitzer spectrum in Fig. B.2. We are interested in the narrow emission features of forsterite (a crystalline silicate) which are superimposed on much broader features from amorphous silicates, that occupy the same wavelength region – except for the 69 μm feature. The continuum – that includes these broad features – is therefore not easily described by a low-order polynomial.

The continuum is better fit using a spline function, which is essentially a piecewise polynomial between (arbitrary) continuum points – or knots. At every knot, the polynomials are connected such that they have a continuus derivative. The best result is obtained when knots are placed just outside the features, and the knots (“continuum points”) between features are spaced by a distance roughly equal to the feature width.

The fitted spline follows the shape of the continuum outside the features very well. Inside the features, the continuum is smooth, but has more structure than to be expected from the opacities of amorphous silicates only. Especially the bumps at 12, 27 and 32 μm are artifacts from the continuum subtraction, and not from the underlying amorphous silicate features. They arise because the opacity of forsterite between features is not always zero, especially between the big features and the little “wings” at 11, 24/28 and 33 μm. These wings are much weaker or not seen in the Spitzer observations (Fig. 4), and do appear much weaker in other forsterite opacities (Servoin & Piriou 1973).

When comparing the model to the data, we therefore want to exclude these wings, and choose our knots between feature and wings. Although the resulting continuum is not completely smooth – which systematically underestimates the integrated flux – we apply the same continuum subtraction to the observations. This cancels out the systematic error, allowing for a comparison of the shape of the feature, without a systematic error from the wings – which is not in the data. For the 33 μm feature, there is the additional complication that the Spitzer spectrum runs only up to 36.9 μm, and it is not possible to choose a continuum point outside the little wing. We therefore force the tangent at this point to zero to get a reasonable continuum, although this makes the 33 μm flux appear a little stronger than it really is.

Appendix C: Implementation of temperature-dependent opacities

When we implement temperature-dependent opacities into the radiative transfer code, we have to take into account that the opacities that set the temperature depend on that same temperature⁷. We therefore recalculate the opacities at every iteration of the hydrostatic vertical structure loop, to keep the solution selfconsistent.

Optical constants for forsterite have been measured at temperatures of T_i = 50, 100, 150, 200 and 295 K (Suto et al. 2006) for iron free forsterite (Mg₂SiO₄), but no measurements exist for crystalline olivines with small iron contaminations (Fe_xMg_2−xSiO₄). For temperatures between 50 and 295 K, we calculate the opacity from the optical constants for every temperature component separately, and combine them using a linear interpolation. The opacity at temperature T is calculated as: $\mathit{\kappa }\mathrm{\left(}\mathit{T}\mathrm{\right)}\mathrm{=}\left(\frac{{\mathit{T}}_{\mathit{i}\mathrm{+}\mathrm{1}}\mathrm{-}\mathit{T}}{{\mathit{T}}_{\mathit{i}\mathrm{+}\mathrm{1}}\mathrm{-}{\mathit{T}}_{\mathit{i}}}\right) {\mathit{\kappa }}_{\mathit{i}}\mathrm{+}\left(\frac{\mathit{T}\mathrm{-}{\mathit{T}}_{\mathit{i}}}{{\mathit{T}}_{\mathit{i}\mathrm{+}\mathrm{1}}\mathrm{-}{\mathit{T}}_{\mathit{i}}}\right) {\mathit{\kappa }}_{\mathit{i}\mathrm{+}\mathrm{1}}$(C.1)

where κ_i is the opacity at temperature T_i, and T is between T_i and T_i + 1. Below T = 50 K, we use the opacities of 50 K forsterite, and we use the 295 K opacity above a temperature of 295 K. This approach does not take into account the shift in wavelength that might occur outside the 50–295 K range. However the peak shift becomes smaller at lower temperatures (Suto et al. 2006) and is not so important below 50 K. Temperatures above 295 K are not present in the SED due to the disk gap, except in the inner disk which is does not contribute significantly at 69 μm because it is much hotter.

---

Acknowledgments

This research project is financially supported by a joint grant from the Netherlands Research School for Astronomy (NOVA) and the Netherlands Institute for Space Research (SRON). J.C.A. acknowledges CNES and PNPS for financial support

References

All Tables

All Figures