Optical Properties of Non-stoichiometric Amorphous Silicates with Application to Circumstellar Dust Extinction

Non-stoichiometric amorphous samples were prepared by means of the PLD technique at the Institute of Geology, Mineralogy, and Geophysics of the Ruhr-University Bochum (RUB) (e.g., Dohmen et al. 2002). Floating-zone Si (111) crystal structure substrates with size of $10\times 10\,\mathrm{mm}$ and thickness of 740 μm (SILTRONIX) were employed since the Si wafer is highly transparent in the IR region.

An Excimer-Laser with nanosecond pulses of 193 nm wavelength at a frequency 10 Hz was used to ablate a synthesized target pellet material under high vacuum condition (10⁻⁵ mbar) in the PLD chamber. The laser fluences were between 1 and 5 J cm⁻², which typically yielded a film thickness of about 50 nm after 10 minutes. The different stoichiometry of each olivine target pellet was prepared by thoroughly mixing in a mortar three synthetic powders, which were $99.99 \%$ pure SiO₂, MgO, and Fe₂O₃ with the nominal ratio, and this mixed powder was cold-pressed into a pellet of 8 mm diameter. The fabricated target pellets were then annealed in a furnace for approximately 20 hr at temperatures of 1500, 1300, 1225, and 1100°C under an oxygen fugacity ${f}_{{{\rm{O}}}_{2}}$ of 1, 10⁻⁹, 10⁻¹⁰, 10⁻¹¹, and 10⁻¹² bar for Mg₂SiO₄, Mg_0.8Fe_1.2SiO₄, Mg_0.4Fe_1.6SiO₄, and Fe₂SiO₄, respectively. The fO₂ was controlled by the flow of gas mixture (CO and CO₂) during the annealing processes (e.g., Dohmen et al. 2002). The T-fO₂ conditions were selected for the Fe-bearing pellets such as to reduce the Fe⁺³ to Fe⁺² and to be within the stability field of the respective olivine composition according to Nitsan (1974).

In the PLD chamber, the $10\times 10\,\mathrm{mm}$ Si substrate was placed on a substrate holder that was mounted directly in front of the target material on a rotating holder. This Si substrate was first heated to 400°C under a pressure of 10⁻⁵ mbar for 15 minutes in order to remove unnecessary volatile absorbents on the surface, especially H₂O. The deposition was carried out at approximately room temperature and a background gas pressure of around 10⁻⁶ mbar. When an incoming laser beam hit the target, a plasma cloud containing a mixture of neutral and ionic atoms was created in an anterior direction of the mounted Si substrate. Hence, the particles in the plasma are deposited on the mounted Si substrate (more details in Dohmen et al. (2002)).

The chemical compositions of the deposited silicate films were determined by Rutherford backscattering spectroscopy (RBS). The RBS measurements were carried out at the RUBION facility of the RUB using the dynamitron tandem accelerator to produce a 2 MeV beam of alpha particles. The measurements were performed with a final aperture of 1 mm, a beam current typically between 20 and 50 nA, a detecting angle of 160° with a silicon particle detector at an energy resolution of about 16–20 keV, and the sample surface was tilted at about 5° relative to the beam direction to avoid channeling. The RBS spectra were simulated using the software RBX (Version 5.18 Kotai 1994) assuming a density of the amorphous silicate layers as of the respective crystalline olivine. The thin film thickness was not completely homogeneous over the sample; it was maximum for the central area. Uncertainties are ±5 nm minimum related to the energy resolution of the detector. However, the fabricated silicate films are chemically homogeneous at least down to the 10 nm scale and totally amorphous, which has been demonstrated by TEM investigations of thin films with similar compositions (Dohmen et al. 2002; Le Guillou et al. 2015). The oxidation state of Fe is mainly 2+, as shown by Le Guillou et al. (2015). Table 1 shows the physical and chemical characteristics of the deposited amorphous silicate samples analyzed by the RBS.

The IR spectroscopic transmittance measurements have been performed in order to examine the optical properties of the deposited amorphous silicates. A Bruker IFS66 v/S Fourier transform infrared (FTIR) spectrometer with deuterium triglycine sulfate (DTGS) detector was employed for the MIR transmittance measurements in the frequency range between 500 and 5000 cm⁻¹ (20–2 μm), as well as a Bruker Vertex80v with a special FIR-DTGS detector together with a Mylar beam splitter for the FIR measurements in the frequency range between 50–600 cm⁻¹ (200–16.7 μm). The sample compartment of both FTIR spectrometers was evacuated to a pressure below 4 mbar for all measurements. All transmittance spectra were normalized to the spectrum of the bare Si(111) wafer and obtained with a spectral resolution of 4 cm⁻¹. As shown in Figure 3, the transmittance spectra measured in both MIR and FIR regions overlapped each other at around 500 cm⁻¹ and were feasible to connect (more details in Wetzel et al. 2012a, 2012b). Hence, it is possible to determine the optical constants through the incorporated transmittance spectra that can cover a wide spectral range within the IR region.

Zoom In Zoom Out Reset image size

The measured transmittance spectrum of Si-film3 together with a Brendel oscillator model fit (Brendel & Bormann 1992) for the dielectric function

$$\begin{eqnarray}\varepsilon {(\omega )}_{\mathrm{IR}} & = & {\varepsilon }_{\infty }+\displaystyle \sum _{j=1}^{N}\displaystyle \frac{1}{\sqrt{2\pi }{\sigma }_{j}}\\ & & \times {\displaystyle \int }_{-\infty }^{\infty }{dz}\,{e}^{\displaystyle \frac{-{(z-{\omega }_{0,j})}^{2}}{2{\sigma }_{j}^{2}}}\,\displaystyle \frac{{\omega }_{{\rm{p}},j}^{2}}{{z}^{2}-{\omega }^{2}+{\rm{i}}{\gamma }_{j}\omega }\end{eqnarray} \tag{ 1 }$$

using five oscillators is shown in Figure 4. This model accounts for the amorphous structure of a material by assuming that the different IR modes can be represented by Lorentz-oscillators with randomly shifted resonance frequencies distributed according to Gaussian probability distributions. We performed the spectral fits in the range between 100 and 1400 cm⁻¹ using the software package SCOUT (Theiss 2011). To model the transmittance spectrum, the value of ${\varepsilon }_{\infty }$ has to be used as an input parameter. However, no values for ${\varepsilon }_{\infty }$ obtained from amorphous silicates equivalent to our non-stoichiometric samples are available.

Zoom In Zoom Out Reset image size

Fundamentally, the optical properties are sensitively influenced by chemical compositions (Koike et al. 2003; Tamanai et al. 2009), and thus the value of ${\varepsilon }_{\infty }$ also varies with different concentrations of Fe and Mg in the silicate system. Consequently, the values of the high-frequency dielectric constant (${\varepsilon }_{\infty }$) derived from the ellipsometric fitting have been applied for the transmittance modeling in order to determine the accurate dielectric function (see Section 3.3). Since the oscillator damping constant γ cannot properly be determined from the fits for vibrational spectra of disordered solids (cf. Ishikawa et al. 2000), we fixed the value of γ to the resolution of our measurement 4 cm⁻¹. Nearly the same type of procedure has previously been reported in the literature (Naiman et al. 1985; Grosse et al. 1986; Brendel & Bormann 1992; Ishikawa et al. 2000). Our parameters of the best fit for the oscillator model are given in Table 2. The Brendel oscillator model provides a reasonable fit to the experimental transmittance spectra of our amorphous silicate samples.

Table 2.  Parameters for the Brendel Oscillators

Osc.	Param.	Si-film1	Si-film2	Si-film3	Si-film4
	(cm−1)
1	${\omega }_{0}$	377	372	366	307
	${\omega }_{p}$	514	625	614	391
	σ	125	139	135	94
2	${\omega }_{0}$	502	525	520	493
	${\omega }_{p}$	390	410	401	452
	σ	75	73	61	61
3	${\omega }_{0}$	739	721	708	687
	${\omega }_{p}$	223	181	191	174
	σ	70	42	56	74
4	${\omega }_{0}$	972	942	921	920
	${\omega }_{p}$	692	752	670	628
	σ	77	77	67	62
5	${\omega }_{0}$	1118	1098	1064	1062
	${\omega }_{p}$	212	299	268	396
	σ	63	74	58	79
${\varepsilon }_{\infty }$		2.43	2.67	3.15	3.69

Note. The values of ${\varepsilon }_{\infty }$ are obtained from the IRSE modeling.

Download table as:  ASCIITypeset image

To obtain the real and the imaginary part of the optical constants of the amorphous silicate samples, we carried out infrared spectroscopic ellipsometry (IRSE) as well. IRSE is a characterization technique for obtaining physical (e.g., thickness and surface roughness) and optical properties of both anisotropic and isotropic materials. Fundamentally, ellipsometry measures in the polarization state of light when it is reflected from a surface of a sample. As a result, Ψ (relative amplitude change) and Δ (relative phase change) spectra can be obtained from the measurements. In this investigation, the measurements of Ψ and Δ spectra have been carried out in the frequency range between 300 and 6000 cm⁻¹ within an angle of incidence for 60° and a resolution of 4 cm⁻¹ under dry air (IR-WASE: J.A. Woollam Col, Inc.). These Ψ and Δ parameters correspond to the Fresnel reflection coefficients R_p (p-polarized: parallel to the plane of incidence) and R_s (s-polarized: perpendicular to the plane of incidence). So the complex reflectance ratio ρ as the basic equation for ellipsometry can be described as

$$\begin{eqnarray}&&\rho =\displaystyle \frac{{R}_{{\rm{p}}}}{{R}_{{\rm{s}}}}=\tan ({\rm{\Psi }})\,{{\rm{e}}}^{-{\rm{i}}{\rm{\Delta }}},\end{eqnarray} \tag{ 2 }$$

where tan(Ψ) denotes the ratio of reflected amplitudes, and Δ is the phase shift of p- and s- linearly polarized components as a result of reflection.

A dielectric function that is derived from the ellipsometry measurements for uniform and isotropic bulk materials is described as a "pseudo-dielectric function" that can be derived from the measured values of Ψ and Δ spectra,

$$\begin{eqnarray}&&\varepsilon ={\varepsilon }_{1}+{{\rm{i}}\varepsilon }_{2}={\sin }^{2}({{\rm{\Phi }}}_{1})\left\{1+{\tan }^{2}({{\rm{\Phi }}}_{1})\left(\displaystyle \frac{1-\rho }{1+\rho }\right)\right\},\end{eqnarray} \tag{ 3 }$$

where ${\varepsilon }_{1}$ is the real part, ${\varepsilon }_{2}$ is the imaginary part of the complex dielectric function, and ${{\rm{\Phi }}}_{1}$ is the angle of the incident radiation (Fujiwara 2007).

Because of the inhomogeneous line broadening of the amorphous samples, the Gaussian oscillator model of the ellipsometry software is applied to model the IR vibrational modes of the amorphous silicates. The spectral line shape of the measured Ψ and Δ spectra is approximated by a Gaussian line profile, which can be expressed by the Kramers–Kronig relations. The real part ${\varepsilon }_{1}$ of the complex dielectric function is related to the imaginary part ${\varepsilon }_{2}$ at all frequencies and is represented as

$$\begin{eqnarray}&&{\varepsilon }_{1}(\omega )=\displaystyle \frac{2}{\pi }\,P{\int }_{0}^{\infty }\displaystyle \frac{{\rm{\Omega }}\,{\varepsilon }_{2}({\rm{\Omega }})}{{{\rm{\Omega }}}^{2}-{\omega }^{2}}d{\rm{\Omega }},\end{eqnarray} \tag{ 4 }$$

where P is the Cauchy principal part (Peiponen & Vartiainen 1991).

The imaginary part describes each Gaussian oscillator mathematically,

$$\begin{eqnarray}&&{\varepsilon }_{2}({\rm{\Omega }})={{A}{\rm{e}}}^{-\tfrac{{({\rm{\Omega }}-{\omega }_{0})}^{2}}{2{\sigma }^{2}}}-{{A}{\rm{e}}}^{-\tfrac{{({\rm{\Omega }}+{\omega }_{0})}^{2}}{2{\sigma }^{2}}},\end{eqnarray} \tag{ 5 }$$

where A denotes the amplitude, ${\omega }_{0}$ is the center/resonance energy, and σ is related to the broadening Br of a band by

$$\begin{eqnarray}&&\sigma =\displaystyle \frac{\mathrm{Br}}{2\sqrt{2\mathrm{ln}(2)}}\end{eqnarray} \tag{ 6 }$$

(De Sousa Meneses et al. 2006). All three parameters, A, σ, and ${\omega }_{0}$, are ellipsometry model fit values that are summarized in Table 3.

Table 3.  Parameters for Ellipsometric Modeling (Gaussian)

Osc.	Param.	Si-film1	Si-film2	Si-film3	Si-film4
	(cm−1)
1	${\omega }_{0}$	⋯	⋯	369	⋯
	A	⋯	⋯	4.36	⋯
	σ	⋯	⋯	139	⋯
2	${\omega }_{0}$	461	448	501	412
	A	3.83	4.51	1.66	4.61
	σ	96	106	65	134
3	${\omega }_{0}$	776	751	746	910
	A	0.72	0.44	0.59	2.57
	σ	39	26	20	54
4	${\omega }_{0}$	964	931	903	981
	A	3.62	3.61	3.43	2.09
	σ	87	88	70	123
5	${\omega }_{0}$	1157	1129	1021
	A	0.36	0.29	0.76	⋯
	σ	47	56	77	⋯
${\varepsilon }_{\infty }$		2.43	2.67	3.15	3.69

Note. "A" has the unit "1".

Download table as:  ASCIITypeset image

In this analysis, we consider a very slight composition gradient with silicate thickness by applying the "graded index layer model" of the WVASE software. A composition gradient depth like this might occur as a result of, e.g., different oxygen contents (e.g., oxidation from atmosphere) and/or surface roughness (more voids on a surface) (Synowicki 1998; Von Rottkay et al. 1998). Figure 5 shows a comparison between the Gaussian model fit with and without considering the graded effect.

Zoom In Zoom Out Reset image size

A clear deviation can be seen at longer ($\gtrsim 20\ \mu {\rm{m}}$) and shorter ($\lesssim 2\ \mu {\rm{m}}$) wavelength regions in Figures 5(a) and (c). Especially, the discrepancy at shorter wavelength region influences the accurate determination of the real part value of the dielectric function at high photon energies (${\varepsilon }_{\infty }$). In ellipsometry measurements, the value of the high-frequency dielectric constant (${\varepsilon }_{\infty }$) is determined as a fit parameter. The value of ${\varepsilon }_{\infty }$ increases with increasing concentration of Fe in the amorphous silicate system (Table 3). When the "graded index layer model" is not applied for the model fitting, the ${\varepsilon }_{\infty }$ value of Si-film3 is 3.09, which is 0.06 lower than the "graded layer" value. The maximum difference in the ${\varepsilon }_{\infty }$ value between the graded layers and single-layer fitting is 0.12, which is obtained for Si-film4. Furthermore, in the case of silicates, the high-accuracy fitting in the longer wavelength region is necessary to pin down the exact location of the Si-O-bending vibration band at around 20 μm. Therefore, the graded index effect is taken into account to improve the model fit of silicate samples, as shown in Figures 5(b) and (d).

There are several discrepancies between the oscillator data of Tables 2 and 3. Whereas transmittance at normal incidence measures only vibrational dipoles parallel to the layer, ellipsometry is also sensitive to those perpendicular to the layer. In our fit we assumed isotropic and smooth layers, and hence, any small anisotropy and inhomogeneity is not taken into account in the ellipsometric analysis. Nevertheless, ellipsometry delivered accurate ${\epsilon }_{\infty }$ values, the use of which improved the transmittance analysis.

In Figure 6 the measured transmittance spectra have been converted into absorbance spectra by the equation

$$\begin{eqnarray*}&&\mathrm{absorbance}=-\mathrm{log}(T)\end{eqnarray*}$$

(T is the transmittance). Figure 6 shows these absorbance spectra in the wavelength range between 8 and 25 μm for four various stoichiometries of amorphous silicate samples. Table 4 lists the positions of the two strongest absorbance features. They mainly include Si-O-Si stretching and bending vibrations, respectively. As obvious in Figure 6, the Si-O asymmetric stretching vibration band undergoes a redshift as the Fe concentration in the amorphous silicate system increases. However, this tendency is not systematic up to Si-film4 because of the different behavior of the two modes that contribute to this peak, see Table 2. A weak Si-O-Si symmetric stretching vibration band around 13–14 μm has been observed for Si-film1, Si-film2, and Si-film3.

Zoom In Zoom Out Reset image size

As mentioned in Sections 3.2 and 3.3, we have derived the optical constants of each amorphous silicate film by two different methods. We inspect whether these methods deliver corresponding optical constants and how we can combine the results for further astronomical simulations in this section.

A comparison of the real and imaginary parts derived from modeling the ellipsometric and transmittance spectra of the Si-film3 sample is shown in Figures 7(a) and (b), respectively. Since we employed the value of ${\varepsilon }_{\infty }$ derived from the ellipsometric model for the transmittance, the spectral fit could be improved significantly. Concerning the imaginary part, the peak positions of the resonance frequency around 10 μm as well as the peak around 20 μm are located at almost the same positions. However, an obvious difference is observed around 12.5 μm: absorption due to the Si-O symmetric vibration clearly appears as a shoulder in the ${\varepsilon }_{2}$ spectrum derived from the ellipsometric modeling, but in the transmittance fit, it appears only as a plateau between two strong features. The fit parameters for this resonance position are different for the two methods, see Tables 2 and 3. The same result has been obtained for the Si-film1 and Si-film2 samples. These deviations are related to the very weak IR activity of this mode, which together with the inhomogeneous line broadening complicates the spectral analysis. The strong disparity at the longer wavelength side ($\gtrsim 20\ \mu {\rm{m}}$) is related to the limited measurement range of the IRSE, which is 1.7–33 μm. When we approach the limit of the measurement range ($\gtrsim 20\ \mu {\rm{m}}$), the signal-to-noise ratio decreases, see this spectral range in Figure 5. In this noisy range, the ellipsometric data are no longer used for further simulations.

Zoom In Zoom Out Reset image size

The two derived data sets were combined at reasonable points in order to cover a broad spectral range between 2 and 200 μm for astronomical simulations. In Figure 7 the short dashed line exhibits the combined dielectric function. For both the real and imaginary parts, there were some points of contact where it was easy to connect these data sets without critically losing accuracy. Figure 7(a) shows the real part of Si-film3. The transmittance data were used at the short wavelength between 1.25 and 4.17 μm. The ellipsometric model data were integrated at 4.17 μm up to 14.3 μm because the Si-O-Si symmetric stretching vibration band appears more clearly in this region. The ellipsometric data were again combined with the transmittance data at the point of 15.4 μm. Because there was no exact intersection point between those two data sets, we here interpolated the ellipsometric data between 14.3 and 15.4 μm where the ellipsometric data line was descended along the transmittance line curve up to 15.4 μm. Then the transmittance data were selected from 15.4 to 200 μm. As mentioned in Section 3.3, the Δ spectra are highly sensitive to thin films; therefore, we selected the ellipsometric data in the MIR region where the Si-O-Si asymmetric and symmetric vibration bands appear clearly. A similar modification of the complex dielectric function data sets was performed for Si-film1, Si-film2, and Si-film4.

Dust associated with highly evolved cool stars is freshly formed in the massive outflows that evolve as the stars ascend the red giant branch in the Hertzsprung–Russell diagram (HR) of stellar evolution toward very high luminosities. Depending on whether the oxygen-to-carbon abundance ratio remains $\lesssim 1$ or exceeds this value, the stars form mineral dust from the abundant rock-forming elements (Mg, Si, Al, Ca, and Fe) or they form soot and carbides, respectively.

The low- and intermediate-mass stars (with initial masses $\lesssim 8\,{M}_{\odot }\,{{\rm{a}}}^{-1}$) are in their dust-forming phase on their second giant branch, the AGB, where they have a degenerated carbon–oxygen core and alternatively burn hydrogen and helium in two shell sources. Some of them become carbon stars (with C/O $\gtrsim 1$ at their surface) if some of the freshly produced carbon is dredged up to the stellar surface. The massive supergiants (with initial masses $\gtrsim 8\,{M}_{\odot }\,{{\rm{a}}}^{-1}$ and $\lesssim 40\ {M}_{\odot }\,{{\rm{a}}}^{-1}$) are on their first giant branch, where they burn hydrogen in a shell source on top of their helium core. Owing to their different mass and internal structures, the properties of AGB stars and supergiants are somewhat different during their final evolution. The most important differences in our context are that supergiants (i) retain their initial abundances of refractory dust forming elements during their whole evolution until they explode as a supernova, and (ii) they have much higher luminosities and slightly higher surface temperatures than AGB stars, such that they do not cross one of the instability strips in the HR-diagram where stars become large-amplitude pulsators. These differences makes supergiants better suited as test cases for investigations of dust condensation in stars than AGB stars.

First, the critical factor ruling the chemistry in the outflow, the C/O abundance ratio, does not increase by mixing substantial amounts of carbon from the core region to the surface as in AGB stars, but instead, it slightly decreases by mixing CNO-cycle equilibrated material in which most C is converted into ¹⁴N to the surface. This guarantees that there is always a sufficient excess of oxygen over carbon that all the refractory rock-forming elements can be oxidized. This leaves no doubt about the chemical nature and composition of the dust mixture that could be formed.

Second, although the supergiants are variable to some extent as all stars in the upper right part of the HR-diagram, they do not show the high-amplitude visual magnitude variations characteristic of the pulsationally unstable stars found on or close to the AGB (Arroyo-Torres et al. 2015). For this reason, the hydrodynamic structure of the outflow from supergiants is expected to be much simpler than for their lower-mass relatives on the AGB because there are no strong shock waves running through the atmosphere and the inner part of the dust envelope, as in the case of Miras or semi-regular variables. A simple steady-outflow model represents the density structure of the dust shell sufficiently accurately for our purpose.

We restrict our considerations to objects with optically thin dust shells. This prevents radiative transfer effects and scattering from being important for the formation of the IR-emission spectrum, which complicates the interpretation of the spectra.

Dust in the environment of supergiants has been studied by Speck et al. (2000) and Verhoelst et al. (2009) for a large sample of stars for which ISO-spectra are available. We select from these studies a small number of objects that satisfy the following requirements:

• 1.  
  The dust shells are clearly optically thin in the IR region. For this, we consider only objects of the category 2.SEc as defined by Kraemer et al. (2002).
• 2.  
  The spectra show marked but not very strong emission bands around 10 and 18 μm that are typical for amorphous silicate dust. This serves to ensure that the optical thickness in the bands is $\ll 1$.
• 3.  
  The spectra show no obvious indication of the presence of crystalline silicates to avoid a significant contribution of such material to the bands around 10 and 18 μm.
• 4.  
  The spectra show low noise in the spectral region $\lambda \lesssim 25\ \mu {\rm{m}}$.

The objects that fulfill these requirements such that they are suited for our purpose belong to the objects of the category 2.SEc as defined by Kraemer et al. (2002), which is essentially defined as satisfying the properties of items 1–3. Our set of objects selected from the database of Sloan et al. (2003) is shown in Table 5 together with some basic parameters of the star and the outflow. Completely reduced ISO-spectra for these four stars, as used in Verhoelst et al. (2009), were kindly provided to us by S. Hony.

Table 5.  The Sample of Supergiants and Their Basic Data

	Spectrala	Variablea	Perioda	Teffb	Lb	${\dot{M}}_{\mathrm{gas}}$ c	vexpc	Dd	AVb
Name	Type	Type	(day)	(K)	(L☉)	$({M}_{\odot }\,{{\rm{a}}}^{-1})$	(km s−1)	(pc)
μ Cep	M2Ia	SRc	730	3700	3.4 × 105	5.0 × 10−6	20	870	2.01
RW Cyg	M3Iab	SRc	550	3600	1.45 × 105	3.2 × 10−6	23	1320	4.49
W Per	M3Iab	SRc	485	3550	5.5 × 104	2.1 × 10−6	16	1900	4.03
RS Per	M4Iab	SRc	245	3550	1.4 × 105	2.0 × 10−6	20	2400	2.63

Notes.

^aTaken from SIMBAD database. ^bLevesque et al. (2005). ^cMauron & Josselin (2011). ^dJones et al. (2012).

Download table as:  ASCIITypeset image

We fit the MIR spectral data of these four test objects with calculated spectra from a model for an optically thin circumstellar dust shell. The properties of the dust shell including the dust mixture are determined by an optimization procedure where we vary the parameters of the problem such as to obtain the best possible fit between observed and synthetic spectrum.

We do not model in this study the formation of dust in the stellar outflow and the acceleration of the outflow by radiation pressure on the dust particles and subsequent momentum transfer from dust to gas by frictional coupling. Instead, we follow the practice of most other laboratory studies that compare their results with observed stellar spectra and assume a spherically symmetric dust density distribution around the star and instantaneous condensation of dust when the dust temperature drops below a given condensation temperature. For optically thin dust shells, this simple model is not as unrealistic as it looks at a first glance, either, because the details of the spatial distribution of the dust are of minor importance (see the Appendix).

The main purpose is to check how well the amorphous silicates for which we determined optical constants can reproduce the observed MIR emission spectrum of our test stars. The amorphous silicate dust generally dominates the IR emission from circumstellar dust in the 9–20 μm spectral region where the diagnostically important resonances from stretching vibrations of the Si-O bond ($\sim 10\ \mu {\rm{m}}$) and O-Si-O bending ($\sim 18\ \mu {\rm{m}}$) vibrations are located as well. To fit the synthetic model, we therefore use the spectral region between 9 and 25 μm. Furthermore, we consider in the fitting procedure the spectral region between 3 and 4.5 μm because (i) in this range the shape of the observed spectrum shows no suspiciously strong molecular absorption bands that render it difficult to find the level of stellar continuum emission, and (ii) the absorption of silicate dust in this range is very low such that we observes the unattenuated stellar radiation. This NIR spectral range helps to fix the baseline of zero dust emission.

Figure 8 shows the mass-absorption coefficients calculated for the four different amorphous silicate materials discussed in Section 4.1. To average over particle sizes, a MRN distribution between a minimum diameter of 10 nm and a maximum of 500 nm in 10 bins is used for all species. In applications to circumstellar dust shells, it usually turns out that the assumption of spherical grains is not suited to reproduce observed band profiles. These are much better reproduced by using, e.g., a continuous distribution of ellipsoids (CDE) (Bohren & Huffman 1983), which is therefore used in our calculations for all species, except for iron, for which we use the Mie theory. We compare the absorption coefficient calculated for our amorphous silicates with absorption coefficients using data for optical constants of amorphous silicates from Dorschner et al. (1995) and Jäger et al. (2003), which are synthesized by methods different from ours.

Zoom In Zoom Out Reset image size

To calculate dust temperatures, we need to also know the optical properties in the range λ < 2 μm, where we have no measured data for our silicate materials. For the purpose of calculating the dust temperature, we use in the 0.2 $\leqslant$ λ < 4 μm wavelength range the x = 0.5 data from Dorschner et al. (1995). The cut at λ = 4 μm is chosen such that it is located in a wavelength range where the silicate opacity is very low, such that it contributes neither to dust heating nor to emission from dust. This choice is somewhat arbitrary, but it influences only the calculated radius of the inner edge of the dust shell for a given condensation temperature, but not the infrared model spectrum at given condensation temperature and dust mass.

In principle, it would be desirable to use in the fitting process only the amorphous silicate dust material studied in this paper since we are interested in the first place whether it is possible to already obtain reasonable fits with these materials alone. However, we also find in most oxygen-rich objects a non-negligible contribution to emission from an additional dust component, which is identified with amorphous corundum (Onaka et al. 1989a, 1989b; Little-Marenin & Little 1990; Jones et al. 2014). Unfortunately, the emission from corundum contributes to the spectral range between 10 and 15 μm is such that it interferes with the silicate dust emission. In particular, it fills up the deep absorption minimum between the two silicate features. For this reason, we include amorphous corundum as a dust component in the modeling.

Optical constants for the corundum for $\lambda \geqslant 7.8\ \mu {\rm{m}}$ are taken from Begemann et al. (1997) (the constants for porous amorphous aluminium oxide). For shorter wavelengths they are augmented by data from Palik (1985).

Interferometric observations have shown that the corundum dust seems to be formed much closer to the star than the silicate dust (e.g., Danchi et al. 1994; Karovicova et al. 2013; Wittkowski et al. 2015, and references therein), as could be expected because of its much higher thermal stability (Salpeter 1977). The grains may then either be heterogeneously composed of corundum grains that condense earlier than silicates in a stellar outflow, which later become overgrown by silicate material as the temperature drops in the outflow, or they are the result of agglomeration of corundum grains and silicate grains. To model the opacity, we have to include corundum grains either as a separate dust component, or as inclusions within silicate grains, or in both variants. Practically, we found that including the corundum only in the form of a separate dust component resulted in slightly better fits than also including silicate particles with corundum inclusions. The difference is, however, only marginal such that we cannot clearly distinguish between the possible cases on the basis of spectral modeling alone. We consider in our modeling corundum and amorphous silicate as separate dust components.

Furthermore, metallic iron particles seem to be present as inclusions in the silicate grains. They are required to explain the NIR opacity of real circumstellar dust (Ossenkopf et al. 1992) and the elemental composition of presolar silicate grains (see Section 2). The contribution of such inclusions to the opacity of the composed particles is low in the diagnostically important spectral region at wavelengths longer than 9 μm, and they cause only an almost constant upward shift of dust opacity in this region that cannot be separated from the effect of a slight increase in silicate dust abundance. The main effect of these iron inclusions is their strong contribution to the optical-to-NIR dust opacity, which rules the dust temperature. This effect of the iron inclusions becomes important only for an iron-free silicate matrix, however. If the silicate material already contains a significant iron content, the influence of possible iron metal inclusions becomes moderate. We found that we cannot clearly distinguish from the quality of fit between models with and without iron grains because of a degeneracy between the content of iron inclusions and the amount of dust that is present. We therefore neglect iron inclusions.

Another possibility is that iron forms as a separate dust component. This follows from considerations on thermal stability (e.g., Salpeter 1977; Gail & Sedlmayr 1999) and has been found to be a likely contributor to opacity in circumstellar dust shell in the $\lambda \lt 8\ \mu {\rm{m}}$ spectral region by Kemper et al. (2002) and Verhoelst et al. (2009). This case differs from the case of iron inclusions in so far as inclusions are thermally coupled to the silicate matrix, while separate iron particles have their own temperature distribution. In effect, this results in significant differences in the calculated model spectra compared to the case of iron inclusions. We also consider iron grains in our modeling. The optical constants are taken from Ordal et al. (1988).

We fit the MIR spectral data of the four test objects with calculated spectra from circumstellar dust shell models (also called synthetic spectra) by an optimization procedure described in Appendix.

The stellar luminosities and effective temperatures of the objects are taken from the literature; the corresponding values and references are given in Table 5. No dereddening was finally performed for the ISO-spectra because it turned out that the corresponding correction is very small in the wavelength range of interest for all our objects.

The dust temperature in the models drops below 100 K at a distance of about 10³ stellar radii. Regions farther out do not contribute significantly to the infrared spectral region λ < 30 μm we used in our comparison of observed with synthetic model spectra. The outer radius of the dust shell therefore cannot be pinned down by our model optimization; it was set to a fixed value of 10⁴ stellar radii in all cases.

5.3.1. Fits with Two Dust Components

In the following models, we consider two dust species, amorphous silicate particles and corundum particles, as opacity sources.

W Per: A fit to W Per was tried first using the optical constants of Si-film1, Si-film2, and Si-film3 for the silicate component. The results are shown in Figure 9. The iron-free Si-film1 can be excluded because the resonance at $\sim 10\ \mu {\rm{m}}$ does not fit the data. The optical data of Si-film2 result in a good fit for the resonance at $\sim 18\ \mu {\rm{m}}$ and a reasonable fit to the peak and the long-wavelength flank of the $\sim 10\ \mu {\rm{m}}$ resonance. The short-wavelength flank is too extended, however. The dust model Si-film3 gives a low-quality fit, but it deviates from the ISO data just in the opposite direction than the Si-film2 model. This means that an amorphous silicate of the type studied in this paper with an intermediate iron content between Si-film2 and Si-film3 would fit the observed spectral energy distribution.

Zoom In Zoom Out Reset image size

A linear interpolation of the data for Si-film2 and Si-film3 using a fraction of 0.6 of Si-film2 and a fraction of 0.4 of Si-film3 results in an rather good fit between the ISO data and the corresponding synthetic spectrum. The corresponding model is shown in Figure 10. The correspondence between observed spectrum and optimized synthetic spectrum appears almost perfect, except in the range between about 6–8 μm. In particular, the two silicate bands are well reproduced by the model.

Zoom In Zoom Out Reset image size

We discuss the missing flux shortward of 8 μm below.

RW Cyg: In the case of RW Cyg, the observed stellar spectrum is embraced by the optimized models using the data of Si-film2 or Si-film3. A model using a weighted mean of the opacity data with a fraction of 0.4 from Si-film2 and 0.6 from Si-film3 results in a model that fits the observational data well. This model is shown in Figure 10.

RS Per: Although the observed stellar spectrum is again embraced by the optimized models using the data of Si-film2 or Si-film3, it is not possible to obtain a simultaneous close fit of the 10 μm and the 18 μm silicate bands for RS Per. The 10 μm band is fitted reasonably well with a weighted mean of the opacity data using a fraction of 0.4 from Si-film2 and 0.6 from Si-film3. The 18 μm band is best fit with Si-film3 alone. A compromise is to a use a weighted mean of the opacity data using a fraction of 0.7 from Si-film3 and 0.3 from Si-film2. The resulting model spectrum is shown in Figure 10. The peak position of the 18 μm band in the model appears at a position slightly shortward of the peak in the observed spectrum, and there is a noticeable discrepancy in the model flux at wavelengths longward of 20 μm.

μ Cep: A similar situation is encountered for μ Cep. A compromise here is to a use a weighted mean of the opacity data using a fraction of 0.8 from Si-film2 and 0.2 from Si-film3. The resulting model spectrum is shown in Figure 10. Again, the peak position of the 18 μm band in the model appears at a position slightly shortward of the peak in the observations, and here there is also a noticeable discrepancy in the model flux at wavelengths longward of 20 μm.

5.3.2. Fits with Three Dust Components

In the following models, metallic iron as a separate dust component is added to the amorphous silicate and corundum particles in order to provide opacity in the $\lambda \lt 8\ \mu {\rm{m}}$ spectral region where the silicate and corundum dust components are transparent. For the optimization we also included the wavelength range between 5.5 and 7.5 μm in the calculation of ${\chi }^{2}$ because iron dust is the sole source of the emission from the dust shell in this region.

W Per: The addition of a third dust component changes the relative abundances of amorphous silicate and corundum, and additionally, the iron content of the amorphous silicate dust required to obtain a good fit. An interpolated dielectric function using a fraction of 0.5 for Si-film2 and Si-film3 results in a close fit between observed and synthetic spectrum. The resulting model spectrum is shown in Figure 11. Now the flux deficit in the 5–8 μm spectral region of the two-component models seen in Figure 10 has completely disappeared. The remaining slight flux excesses of the synthetic spectrum over the observed spectrum at around 4.8, 6.8, and 8 μm are obviously due to molecular absorption bands of CO, H₂O, and SiO, respectively. The stretching and bending modes at 10 and 18 μm of the amorphous silicate are well reproduced by the model spectrum.

Zoom In Zoom Out Reset image size

RW Cyg: In this case, a linear interpolation of the optical constants of Si-film2 and Si-film3 with weights 0.7 and 0.3, respectively, is required for a best fit. The resulting model spectrum is shown in Figure 11. The correspondence between observed spectrum and model spectrum is very close.

RS Per: In this case, a linear interpolation of the optical constants of Si-film2 and Si-film3 with weights 0.2 and 0.8, respectively, is required for a best fit. The resulting model spectrum is shown in Figure 11. Again, the correspondence between observed spectrum an model spectrum is good in the range between about 9–20 μm. At longer wavelengths, the model flux is slightly lower than the stellar flux.

μ Cep: In this case, a linear interpolation of the optical constants of Si-film2 and Si-film3 with weights 0.6 and 0.4, respectively, is required for a best fit. The resulting model spectrum is shown in Figure 11. The correspondence between observed spectrum and model spectrum now is almost perfect, also in the range $\lambda \gt 20\ \mu {\rm{m}}$.

We also tried to reproduce the observed spectra without corundum, but retaining iron to fill up the silicate absorption trough between the two silicate bands. This combination gave no good fit because we can fit either the wavelength range between the two silicate features or the wavelength range between 5.5 and λ = 8 μm, but not both simultaneously.

5.3.3. The 10 μm Silicate Band

The 10 μm band often attracts particular attention in discussions of the quality of the fit between radiative transfer models and observed spectra, therefore we consider this feature in more detail.

Figure 12 shows details of the fits between observed and synthetic spectra of our comparison stars for the wavelength range around the 10 μm band corresponding to the Si-O stretching vibrations. The models using the optical data of the set of amorphous silicates synthesized and investigated by us result in fits that are able to reproduce the width, peak position, and general shape of this silicate band with considerable accuracy when they are combined with the dust materials corundum and iron. They also reproduce the change in slope of the long-wavelength flank at ∼12 μm, which seems to be present in most supergiants (Speck et al. 2000). In our model this particular detail results from a the weak Si-O-Si symmetric stretching band in the amorphous silicate absorption in this region and the broad maximum of the corundum absorption. This weak silicate feature is lacking in other data sets for silicate dust absorption because the corresponding vibration mode cannot be easily separated from the background because its absorption strength is much weaker than that of the asymmetric band (see Sections 3.3 and 4.2).

Zoom In Zoom Out Reset image size

The fits are not perfect, however, and the observed spectra show more structure than the calculated spectra, which may be associated with the presence of some minor additional dust components. It is not necessary, however, to include other minerals than the components used by us to explain the main structure of the observed 10 μm band. In particular, it is not necessary to invoke substantial amounts of Ca-Al-silicates (melilite) to fit the observed shape of the profile, as in Verhoelst et al. (2009). Such a compound is unlikely to exist as wide-spread dust component because if it were, it would be found in meteorites as presolar dust grains like other aluminum oxide particles, which is not the case.⁴ The overall fit is, however, substantially improved when an iron dust component is included, as has been found in Verhoelst et al. (2009).

The quality of the fit of the 10 μm feature is nearly the same for the models with or without iron dust particles. The properties of the model with respect to dust content or iron content of the silicate material are slightly changed by including or omitting an iron dust component (cf. Table 6), but in both cases, a parameter combination can be found that reproduces the 10 μm feature. The presence of iron dust is inconsequential in this respect.

Table 6.  Optimized Dust Shell Parameters

	Silicate		Corundum		Iron
Object	Mdu	Ti	Mdu	Ti	Mdu	Ti	${\tau }_{10\mu }$	g
W Per	33.1	554	16.4	1 400			.061	.6
	25.0	623	12.4	1 400	53.2	379	.050	.5
RW Cyg	23.9	709	17.7	1 400			.058	.7
	21.7	753	16.1	1 400	21.7	485	.053	.4
RS Per	3.6	917	3.8	1 400			.027	.7
	3.1	1 087	2.9	1 400	8.6	579	.028	.2
μ Cep	10.0	854	7.0	1 400			.030	.8
	7.1	959	4.0	1 400	28.2	484	.022	.6

Note. M_du is mass of the dust species in units ${10}^{-5}\,{M}_{\odot }$, T_i is temperature of the dust species at the inner boundary (in K), ${\tau }_{10\mu }$ means optical depth of the dust shell at 10 μm, and g and $1-g$ are the weighting fractions of the optical constants for Si-film2 and Si-film3, respectively.

Download table as:  ASCIITypeset image

5.3.4. The 18 μm Silicate Band

Figure 13 shows details of the fits between observed and synthetic spectra of our comparison stars for the wavelength range around the 18 μm band corresponding to the Si-O stretching vibrations. The quality of the data for the observed flux is not as good for the feature at 10 μm, but the band profiles are clearly marked, and this second important emission band of silicate minerals is also reasonably well fit by the three-component mixture of our silicate materials augmented with corundum and iron, although slightly worse than the 10 μm band.

Zoom In Zoom Out Reset image size

One problem is, however, that the model flux in the spectral range $\lambda \gt 20\,\mu {\rm{m}}$ is slightly but systematically below that of the stellar spectrum (cf. Figures 10 and 11) when most weight is given to the peaks by the optimizing procedure. One reason for this deficit may be that a value of p = 2 as assumed in our models for the density variation $\varrho \propto {r}^{-p}$ is unrealistic because of, e.g., variable mass-loss rate or outflow velocity. Including p as an additional variable in the optimization procedure, however, returns a value of $p\approx 2.02$ and only slightly improves the situation. The flux deficit seems to be related to a deficit of the extinction coefficient of the dust mixture at long wavelength for at least one of the three dust components.

5.3.5. Discussion

From these results we conclude that it is possible to completely explain the silicate dust emission of supergiants by the presence of some type of amorphous silicate material with a non-standard stoichiometry corresponding to a composition between olivine and pyroxene. We only have to add corundum and iron dust particles in order to supply the required opacity in the range between the two silicate bands and for $\lambda \lt 8\ \mu {\rm{m}}$ to account for the observed emission from the dust shell in these regions, which cannot result from the silicate dust. Experiments with additional dust components (e.g., spinel or iron oxide) resulted in low abundances that were assigned to such additional components by the optimization algorithm.

The formation of a non-stoichiometric amorphous silicate material can be expected under the conditions of an outflow from a supergiant. Because of the high stellar atmospheric temperature, the silicate dust condensation commences at considerable distances from the star when the equilibrium temperatures of the dust drop sufficiently low for the dust to become stable against vaporization. In fact, the temperature of the inner edge of the dust shell, T_i, in our radiative transfer model is substantially lower than an estimated vaporization temperature of 1100 K, except for one case (RS Per, the model with three dust components, see Table 6), where it is close to the limit. This means that the dust is formed under conditions that resemble the vapor deposition on a cold substrate to some extent, where growth species attached to the surface almost immediately become immobile at the surface of the condensed phase.

For the type of dust shell models calculated by us it remains open, however, whether a low temperature at the inner edge, T_i, well below the stability limit of silicates as found in most of our models results (i) from slow particle growth in the continuously diluting and cooling gas, or (ii) from an onset of silicate dust condensation at these low temperatures. In the first case, silicate dust formation commences at higher temperature, but dust grows only slowly in the outflowing gas, such that the dust density peaks at large distances, i.e., low temperature, when sufficient dust has formed for a rapid acceleration of the dusty gas. This then results in a rapid density decrease with increasing distance from the star. In the second case, we possibly have an onset of dust formation at rather low temperatures (e.g., Gail et al. 2013). More detailed models considering dust particle growth are required to choose one of these possibilities, which is beyond the scope of this paper, however.

For all our sample stars a considerable iron content of the silicate material is required to obtain a satisfactory spectrum fit. The iron content is determined by interpolation between the Fe/(Fe+Mg) ratio of our four amorphous silicate samples. The fit quality of the peak position of the two silicate bands and of the short-wavelengths flank of the 10 μm feature allows us to fix the relative contribution of the two components between which we have to interpolate sufficiently well. It appears desirable, however, to also determine optical data for iron contents intermediate between that of our four samples. This would enable a more accurate determination of the iron content than from our admittedly rather coarse grid.

A fit of the spectra with the iron-free silicate Si-film1 and also with the iron-free silicate of Jäger et al. (2003) was not possible. It remains to check whether this is a peculiarity of our selected sample or if this looks different when more objects are studied.

Table 6 shows the dust masses required for a best fit in the range from the inner edge R_i where a dust species condenses to the outer radius R_a at ${10}^{4}{R}_{* }$ assumed in the model calculation. It would be more instructive to convert this into a dust mass-loss rate or into the fraction of the dust-forming elements condensed into dust. Unfortunately, this requires knowledge of the velocity structure of the stellar wind. The assumed variation ∝r⁻² of the density for the model calculation requires a constant velocity, which is not realistic since the radiation pressure on dust increasingly accelerates the outflow with progressing dust condensation. The simplification of a constant velocity is acceptable to calculate the dust mass in optically thin shells, but it may be a problem when we determine dust mass-loss rates. This is particularly a problem for the radius range where corundum but no other dust species exists, because in this range the velocity may be significantly lower than in the range where the dominant dust species exist, i.e., silicates and possibly iron, and where they accelerate the dusty gas. We cannot simply use the observed outflow velocity as listed in Table 5 either because this refers to the asymptotic value of the velocity at large distances. In order to give a crude estimate despite the uncertainty of this approach, we assume an average velocity in the dust condensation zone of v = 10 km s⁻¹ for all stars and species (typically, one half of the final velocity) and calculate a dust mass-loss rate from ${\dot{M}}_{\mathrm{du}}={M}_{\mathrm{du}}\,v/({R}_{{\rm{a}}}-{R}_{{\rm{i}}})$. The result is shown in Table 7. This can be compared to the gas mass-loss rate, which shows that as expected, the dust mass-loss rate is a fraction of about 10⁻² to $\lt {10}^{-3}$ of the gas mass-loss rate. The accuracy of the estimated dust mass-loss rate is low, however, because the value of ${\dot{M}}_{\mathrm{gas}}$ is only known for μ Cep with some accuracy; in the other cases, the values are only rough estimates.

Table 7.  Estimated Dust Mass-loss Rates

	Silicate		Corundum		Iron
Object	${\dot{M}}_{\mathrm{du}}$	f	${\dot{M}}_{\mathrm{du}}$	f	${\dot{M}}_{\mathrm{du}}$	f	${\dot{M}}_{\mathrm{gas}}$
W Per	18.3	1.6	9.0	14.3	40.3	10.9	2100
RW Cyg	10.0	0.57	7.4	7.7	10.2	1.8	3200
RS Per	1.4	0.12	1.3	2.2	4.0	1.1	2000
μ Cep	2.3	0.083	1.3	0.84	9.1	1.03	5000

Note. ${\dot{M}}_{\mathrm{du}}$ is the mass-loss rate of the dust species in units of ${10}^{-9}\,{M}_{\odot }\,{{\rm{a}}}^{-1}$, f is the estimated fraction of the key element for dust formation (Si for silicates, Al for corundum, Fe for iron) that is condensed into dust.

Download table as:  ASCIITypeset image

We can even proceed and calculate from $f={\dot{M}}_{\mathrm{du}}/(A\,\epsilon {\dot{M}}_{\mathrm{gas}})$ a condensation fraction f of the key elements (Si for silicates, Al for corundum, Fe for iron) for condensation into dust, where A is the atomic weight of the chemical formula unit of the solid and is the element abundance of the key element. The results are also shown in Table 7. In principle, the fraction f should satisfy $f\lt 1$. This is satisfied for μ Cep, where the value of ${\dot{M}}_{\mathrm{gas}}$ is accurately known. In the other cases, we partly obtain distinctely higher values. For corundum dust, this may result from a lower-than-assumed value of v in the inner range of the dust shell where only corundum dust exists, as mentioned above. For iron, it may result from the fact that we used optical data for pure iron, while iron in space is always a nickel-iron alloy that may have different extinction properties. Finally, the estimated dust mass-loss rates as derived by Mauron & Josselin (2011) from the formula of Jura & Kleinmann (1990) may be systematically too low. A more accurate modeling of the dust shell considering the growth of grains and the outflow dynamics is required to determine the condensation degrees for the dust species.

The synthetic spectrum of a dust shell is calculated by taking advantage of the assumption of an optically thin dust shell. This allows a very rapid calculation of a model spectrum, and this in turn enables computing large sets of such models to obtain optimized fits between model spectra and observations. The model is calculated as follows.

• 1.  
  For the emission from the central star, we assume a blackbody spectral energy distribution with effective temperature T_eff. Alternatively, a read-in spectral energy distribution can be used.
• 2.  
  The dust temperature T_d is determined from radiative equilibrium between absorption and emission of a dust grain

  $$\begin{eqnarray}&&{\int }_{0}^{\infty }{\kappa }_{\nu }^{(\mathrm{abs})}{B}_{\nu }({T}_{\mathrm{eff}})W(r){{\rm{e}}}^{-{\tau }_{\nu }}\,d\nu \\ &&\quad ={\int }_{0}^{\infty }{\kappa }_{\nu }^{(\mathrm{abs})}{B}_{\nu }({T}_{{\rm{d}}})d\nu ,\end{eqnarray} \tag{ 7 }$$

  where B_ν is the Kirchhoff–Planck function, ${\kappa }_{\nu }^{(\mathrm{abs})}$ is the mass-absorption coefficient of the dust material, and W(r) is the geometric dilution factor of the stellar radiation field

  $$\begin{eqnarray}&&W=\displaystyle \frac{1}{2}\left(1-\sqrt{1-\displaystyle \frac{{R}_{* }^{2}}{{r}^{2}}}\right),\end{eqnarray} \tag{ 8 }$$

  where R_* is the stellar radius. Each dust species has its own dust temperature T_d. The quantity

  $$\begin{eqnarray}&&{\tau }_{\nu }=\displaystyle \sum _{j}{\int }_{{R}_{+}}^{r}{\varrho }_{j}({\kappa }_{j,\nu }^{(\mathrm{abs})}+{\kappa }_{j,\nu }^{(\mathrm{sca})}){dr}^{\prime} \end{eqnarray} \tag{ 9 }$$

  is the optical depth between the stellar surface and the radial distance r of the dust grain from the stellar center. The density ϱ_j is the mass density of the dust species j, and the summation is over all dust species. Our assumption of an optical thin dust shell mainly refers to the mid- to FIR wavelength for which we compare the emission from the dust shell with observations. The heating of the dust grains is due to absorption of stellar radiation which peaks around $1...1.5\ \mu {\rm{m}}$. Since the extinction by dust grains strongly increases from the near-infrared to the ultraviolet, the dust shell may be optically thin in the mid- to FIR wavelength even if its optically thickness is not small in the optical spectral region. For this reason, we retain the factor $\exp (-{\tau }_{\nu })$ despite our assumption of an optically thin dust shell.
• 3.  
  For each dust species it is assumed that the dust is spherically symmetrically distributed around a central star and extends from some inner radius R_i to an outer radius R_a. The radial distribution of the dust density ϱ is assumed to follow a power law ${r}^{-p}$, where p is generally assumed to equal p = 2, which corresponds to a stellar wind with constant outflow velocity. The inner radius, R_i, is determined such that the dust temperature of the species takes a prescribed value T_c at this radius. The outer radius, R_a, is prescribed.
• 4.  
  The observed spectral energy flux from the object observed at Earth is

  $$\begin{eqnarray}{F}_{\nu } & = & {e}^{-{\tau }_{\nu }^{\mathrm{ISM}}}\displaystyle \frac{4\pi }{{D}^{2}}\ \left\{\displaystyle \frac{1}{4}{B}_{\nu }({T}_{\mathrm{eff}}){{\rm{e}}}^{-{\tau }_{\nu }}{R}_{* }^{2}\right.\\ & & +\left.\displaystyle \sum _{j}{\displaystyle \int }_{{R}_{{\rm{i}}}}^{{R}_{{\rm{a}}}}{B}_{\nu }({T}_{{\rm{d}},j}){\varrho }_{j}^{}{\kappa }_{j,\nu }^{(\mathrm{abs})}\,{r}^{2}\,{dr}\right\},\end{eqnarray} \tag{ 10 }$$

  where D is the distance of the object from Earth and ${\tau }_{\nu }^{\mathrm{ISM}}$ the optical depth of the ISM along the sightline. The first term in curly brackets is the contribution from the star and the second term is the contribution of the dust shell.For optically thin dust shells, the dust temperature distribution calculated from Equation (7) is spherically symmetric even if the dust distribution is not. Therefore we introduced in the emission term spherical coordinates. More generally, the corresponding integral is a volume integral

  $$\begin{eqnarray*}&&\int {B}_{\nu }({T}_{{\rm{d}},j}){\varrho }_{j}^{}{\kappa }_{j,\nu }^{(\mathrm{abs})}\,{dV},\end{eqnarray*}$$

  where we may decompose the integration domain into spherical shells and perform the angular integrations over each shell separately. Then we can write the volume integral equally well in terms of the average mass densities of the shells. Hence, we have to interpret ϱ in Equation (10) as the mass density averaged over concentric shells around the star. In this sense, Equation (10) also holds when the real dust distribution is not strictly spherically symmetric.
• 5.  
  The model allows us to consider an arbitrary number of dust species. The mass absorption and scattering coefficients are calculated from the complex index of refraction. For inhomogeneously composed particles, an effective complex index of refraction according to the Bruggeman mixing rule (cf. Berryman 1995) can be calculated. The calculation can be made for either spherical particles (Mie theory) or for small ellipsoidal particles (Bohren & Huffman 1983) with a numerically specified distribution of axis ratios and radii, for a continuous distribution of ellipsoids (CDE, see Bohren & Huffman 1983), or for small cubes (Fuchs 1975).

In order to calculate a model for the spectral energy distribution of a given object with observed infrared spectrum, we have to specify the luminosity L and effective temperature T_eff of the star. The stellar radius follows from the standard relation $4\pi {R}_{* }^{2}\sigma {T}_{\mathrm{eff}}^{4}=L$. Furthermore, we have to specify the dust species considered in the model, the corresponding optical properties, and the particle shape and size distribution. For all species, ${\kappa }_{j,\nu }^{(\mathrm{abs})}$ is calculated at the specified wavelength grid. For the dust distributions of all species j, we fix the power p, the outer radius of the dust shell, R_a, and the total dust masses, M_j, and condensation temperatures, ${T}_{j,{\rm{c}}}$, of the dust species.

For a given set of parameters, the radial variation of dust temperature for each species is calculated according to Equation (7) for a small set of assumed values for R_i, and the value of R_i corresponding to the specified value of ${T}_{j,{\rm{c}}}$ is determined by backward interpolation. Then the dust temperatures for $r\gt {R}_{i}$ are calculated, and the spectral energy distribution is calculated according to Equation (10).

For given dust properties, a model depends on the two parameters L and T_eff and on the $2J$ parameters M_j and ${T}_{j,{\rm{c}}}$, where J is the number of dust species. The parameters L and T_eff are derived from astronomical observations and are given quantities for an object.⁵ The parameters M_j and ${T}_{j,{\rm{c}}}$ of the dust shell are unknown and have to be determined by fitting the observed spectrum with the model spectrum.

The goodness-of-fit between the synthetic and the observed spectrum is checked by calculating the quality function

$$\begin{eqnarray}&&{\chi }^{2}=\displaystyle \sum _{l}{\left(\displaystyle \frac{{ \mathcal N }{F}_{\nu }^{\mathrm{obs}}-{F}_{\nu }^{\mathrm{calc}}}{{F}_{\nu }^{\mathrm{calc}}}\right)}^{2}.\end{eqnarray} \tag{ 11 }$$

To calculate ${\chi }^{2}$, not all frequency points used to calculate the spectrum are taken into account, but only those where no strong molecular features from the underlying star are seen in the observed spectrum. The quantity ${ \mathcal N }$ is a normalization factor that accounts for the circumstance that the distance D is usually not accurately known, such that to enable a comparison of the two spectra, we have to apply a shift to one of them to place them on top of each other. This factor is determined such as to minimize ${\chi }^{2}$ as the solution of

$$\begin{eqnarray}&&\displaystyle \frac{\partial \,{\chi }^{2}}{\,\partial { \mathcal N }}=0\end{eqnarray} \tag{ 12 }$$

for ${ \mathcal N }$. The value of ${\chi }^{2}$ calculated with this normalization is taken as the quality function. Then we vary the $2J$ parameter dust masses M_j and condensation temperatures ${T}_{j,{\rm{c}}}$ and search for the parameter combination for which ${\chi }^{2}$ takes the lowest value. This parameter set is then considered as the best-fit model for the dust shell.

Such nonlinear minimization problems are notoriously difficult to solve. Here the minimization of ${\chi }^{2}$ was made by applying a genetic algorithm that usually allows an efficient minimum search even if many parameters are involved. We used the special variant of such an algorithm described in Charbonneau (1995) because this worked quite well in other projects where we used this method.

The genetic algorithm operates with discrete values of the model parameters. The optimized parameter set resulting from this method is therefore slightly off from the true minimum. Instead of improving the accuracy of the genetic algorithm by using a finer resolution of the variable space, we perform a final refinement step by minimization with a gradient method. This converges rapidly because we start already from a position in variable space that is very close to the searched-for minimum.