3 Model analysis

In the present study we attempt to derive the dust emissivity from the observed variation in the mid-infrared spectrum. If the dust emissivity $\kappa$ is proportional to $\lambda^{-\beta}$, where $\lambda$ is the wavelength and $\beta~(>0)$ is a constant, the total flux from a dust grain $F_{\rm gr}$ is expected to vary as $F_{\rm gr} \propto \int \kappa B_\lambda(T)~{\rm d}\lambda \propto T^{4+\beta}$ with the temperature T, where $B_\lambda$ is the Planck function. Thus the variation in the shell flux is directly connected to the emissivity and to the variation in the dust temperature. In reality, the emissivity is a complicated function of the wavelength and the emergent dust shell flux is a summation of the emission from dust grains of different temperatures. To investigate the actual emissivity of the dust grains around Z Cyg we adopt a simple dust shell model. We assume that the dust shell is optically thin and spherically symmetric. The spectrum of the emergent dust shell flux $F(\lambda)$ is then given by

$$F(\lambda) = 4 \pi \int_{r_{\rm i}}^{r_{\rm m}} \sigma(\lambda) B_\lambda(T(r)) r^2 n(r) {\rm d}r / D^2,$$ (1)

where $\sigma(\lambda)$ is the emission cross-section of the dust grains, T(r) is the temperature at the distance rfrom the central star, n(r) is the number density of the dust grains, D is the distance from the earth to the star, and $r_{\rm i}$ and $r_{\rm m}$ are the inner and outer boundaries of the dust shell, respectively. We set $T(r_{\rm m}) = 30$ K and the inner dust shell temperature $T_{\rm i}=T(r_{\rm i})$ as a free parameter in the model fitting. For spherical grains $\sigma(\lambda)$ can be written as $\pi a^2 Q(\lambda)$, where $Q(\lambda)$ is the absorption efficiency factor or the dust emissivity, and a is the grain radius. The efficiency factor $Q(\lambda)$ is proportional to a for $a \ll \lambda$, which is a good approximation in the infrared range (van de Hulst 1957). As for the density distribution n(r), we take into account the acceleration of the flow near the bottom of the envelope and adopt the "third simple model'' of Habing et al. (1994):

$$n(x) = n_{\rm i} / \sqrt{(x - y_0) x^3}$$ (2)

with

$$x = r/r_{\rm i} \qquad \hbox{and} \qquad y_0 = 1 - (v_{\rm s}/v_{\rm out})^2,$$ (3)

where $v_{\rm s}$ is the local sound velocity at $r_{\rm i}$ and $v_{\rm out}$ is the terminal gas flow velocity. In this model, the dust velocity is assumed to be proportional to the gas velocity. This model provides a reasonable approximation to more accurate models in most cases (Habing et al. 1994). The parameter y₀ depends on the efficiency of the momentum transfer (e.g., Ivezic & Elitzur 1995) and thus is a function of the mass-loss rate. In this paper we simply take it as a parameter and make calculations for the cases with y₀ = 0, 0.1, and 0.4. The parameter y₀ = 0.1 represents the most likely case for low mass-loss rates (Habing et al. 1994). We find that the results do not depend sensitively on the value of y₀and thus we present results for y₀ = 0.1 in the followings unless explicitly indicated otherwise. Note that the case y₀ = 0 corresponds to a constant flow velocity model ( $\rho \propto r^{-2}$). The case y₀=0.4 provides an extreme situation whose density distribution deviates strongly from the case of constant flow velocity near the bottom of the shell.

The temperature of the grains at r is calculated from the radiative equilibrium condition:

$$\int_{0}^{\infty} Q(\lambda) F^*(\lambda) ( r_*/r)^2 {\rm d}\... ... 4\pi\int_0^\infty Q(\lambda) B_\lambda(T(r)) {\rm d}\lambda,$$ (4)

where r_* is the stellar radius and $F^*(\lambda)$ is the photospheric flux of the star at r_*. The continuum around 4 $\mu$m is assumed to be due to photospheric emission in the present analysis. The underlying photosphere should include the contribution from the extended atmosphere or the MOLsphere, which has been detected in several late-type stars (e.g., Tsuji et al. 1997; Tsuji 2000a, 2000b; Yamamura et al. 1999a, 1999b) and theoretically investigated in pulsating stars (Fleischer et al. 1992; Bessel et al. 1996; Höfner et al. 1998; Höfner 1999). This component is expected to show a complicated dependence on the visual variability and it is difficult to predict detailed spectra of the outer atmosphere at present (Lobel et al. 2000). Investigations of the near-infrared spectra suggest that the depression of the continuum due to absorption of water vapor in the outer atmosphere may be significant (up to 30%) even at 4 $\mu$m (Matsuura et al. 2002). In this paper we simply assume that the photospheric emission has a functional form with a stellar effective temperature T^* given by Engelke (1992), which is scaled to the observed flux around 4 $\mu$m, and that most of the emission at $\lambda > 7$ $\mu$m comes from the dust shell. The contribution from the photosphere is estimated to be less than 20% in the 10 $\mu$m region even at minimum and the uncertainties in the subtraction procedure do not affect the present results significantly. In Eq. (4) the attenuation due to absorption by grains is neglected since in the present case the optical depth of the grains is still small in the near-infrared, where most of the stellar radiation energy resides, and since the effect is small compared to the uncertainty in the spectrum of the photospheric emission. For T^* we adopt values of 2300 K, 2200 K, 2700 K, 2500 K, 2300 K, 2200 K, and 2700 K at phases $\phi =$ 0.55, 0.79, 0.97, 1.20, 1.42, 1.63. 1.93, respectively, according to the results for o Cet by Danchi et al. (1994). These temperatures may not directly correspond to the photospheric temperatures of the star. The effect of uncertainties in the stellar temperatures is found to be insignificant in deriving the dust emission and it does not affect the present results.

Since the photospheric contribution becomes dominant for wavelengths shorter than 7 $\mu$m, the dust emissivity cannot be estimated from the observations accurately. The optical properties of the dust grains shorter than 7 $\mu$m were calculated for astronomical silicate (Draine & Lee 1984), assuming spherical dust grains with a radius of 0.1 $\mu$m. The absorptivity in the near-infrared determines the physical size of the dust shell, but does not affect the emergent spectrum because of the assumption of an optically thin dust shell. For the optical properties of dust grains at wavelengths longward of 7 $\mu$m, we created a set of dust absorption efficiency factors from the spectrum observed at the first maximum ( $\phi = 0.97$) by changing the assumed inner dust shell temperature. From Eq. (1) the relative spectral shape of the absorption efficiency factor $Q(\lambda)$ is derived straightforwardly from the observed spectrum owing to the nature of the optically thin dust shell. The derived absorption efficiency factor was scaled such that Q/a = 1.35 $\mu$m^-1 at the peak around 9.7 $\mu$m to agree with that of the astronomical silicate. This scaling was needed to make the optical properties connect smoothly at 7 $\mu$m. The absolute values of the absorption coefficient cannot be determined by the present analysis.

By adopting inner dust shell temperature of 400 K, 500 K, 600 K, 700 K, 800 K, 900 K, and 1000 K at the first observation maximum ( $\phi = 0.97$), we have finally obtained a set of 7 dust absorption efficiency factors. We denoted them as Q_i (i=4,...10). To remove the noise originating from the observed spectrum, we fitted Q_i to the spectrum by using an analytical function. For spherical grains of $a \ll \lambda$, the absorption efficiency factor is given by

$$Q(\lambda) = \frac{8 \pi a}{\lambda} \hbox{Im}\left\{ \frac{\epsilon(\lambda)-1} {\epsilon(\lambda) +2} \right\},$$ (5)

where $\epsilon(\lambda)$ is the complex dielectric constant of the dust (van de Hulst 1957). It can be approximated by a summation of oscillators which satisfy the dispersion relation (e.g., Koike et al. 1989):

$$\epsilon = \epsilon_\infty + \sum_j \frac{2 P_j \Gamma_j \lam... ... \frac{1}{(\lambda^2 - \lambda_j^2) - i~\Gamma_j \lambda_j },$$ (6)

where $\epsilon_\infty$ is the dielectric constant in the spectral range much shorter than the wavelengths in question and P_j, $\Gamma_j$, and $\lambda_j$ are the strength, width and position of each oscillator. We find that the emissivity derived from the observed spectrum can be fitted quite well with a set of five oscillators. The "best fit'' parameters of the five oscillators are determined by least square fitting method. The efficiency factors used in the model fitting for y₀=0.1 are shown in Fig. 2. For other values of y₀, we derived similar Q_i's. The ratio of the 18 $\mu$m to 10 $\mu$m bands increases from Q₄ to Q₁₀ as expected.

Figure 2: Trial absorption efficiency factors for the models with y0=0.1. Each curve labeled by Qi ( i=4,... 10) is derived with the assumed inner dust shell temperature at the first maximum ( $\phi = 0.97$) of 400 K, 500 K, 600 K, 700 K, 800 K, 900 K, and 1000 K, respectively (see text).

Figure 3: Variations in the inner dust shell temperature $T_{\rm i}$ (a) and the normalized inner dust density $n_{\rm i}$ (b) for the "best fit'' models with y0 = 0.1. The solid lines represent the parameters for models with Q7, the dotted lines those with Q4, and the dot-dashed lines those with Q10. The inner dust density is normalized to the visual maximum ( $\phi = 0.97$) and is derived based on the assumption that the inner dust shell radius $r_{\rm i}$ is constant.

With each Q_i we fit all the observed spectra with the model, in which the inner dust shell temperatures $T_{\rm i}$ and the inner dust density $n_{\rm i}$ are the only free parameters. The best fit spectra with y₀=0.1 and Q₇ for each observed spectrum are superposed in Fig. 1. Figure 3 shows the parameters of the best fit models for Q₄, Q₇, and Q₁₀, where $n_{\rm i}$ is normalized to the maximum at $\phi = 0.97$ for each Q_i and $r_{\rm i}$ is assumed to not vary from phase to phase. The density $n_{\rm i}$ scales as $r_{\rm i}^{-3}$ if $r_{\rm i}$ changes (see Eq. (1)). Although we obtained different values of $T_{\rm i}$ and $n_{\rm i}$ for the best fit models with each Q_i, the best fit spectra appear to not differ from each other. All of them provide similarly good fits to the observations. Some examples are shown in Fig. 4 for y₀ = 0.1. The fitted spectra are almost identical to each other. Similar results are also obtained for the cases with y₀=0 and y₀=0.4. Therefore, we cannot distinguish between models with the dust emissivities and density distributions examined here solely from the quality of the fit.

Figure 4: Examples of the model fits to the spectra at $\phi = 0.79$ (a) and 1.20 (b). The "best fit'' models with Q4, Q7, and Q10 are plotted by smooth thin lines together with the observed spectra , but they are hardly distinguishable from each other.

Dust formation, if it occurs, adds new dust grains to the shell and could lead to a variation in the observed spectrum. However, dust formation is believed to occur mostly at minimum (e.g., Winters et al. 1994) and it is difficult to imagine that dust formation will explain the observed increase in the flux at maximum. The dust velocity near the bottom of the shell is several km s^-1for optically thin shells (Habing et al. 1994) and during the interval of the observations the dust grains travel only several 10¹⁰ m, which causes a decrease of only a few degrees in the dust temperature even at the bottom of the shell. This hardly affects the emergent spectrum. Furthermore, the motion of the dust is so slow that the decrease in the dust temperature from maximum to minimum cannot be due to removal of hot dust grains even if the observed change at maximum is attributed to the production of newly-formed hot dust grains near the bottom of the dust shell. Therefore we assume that the observed variations originate in the luminosity variation of the central star and that the density distribution does not change over the present observation period. Under these assumptions we investigate whether any of the seven dust emissivities can provide the flux variation compatible with the observations.

We directly use the observed dust shell flux to estimate the variation in the emergent flux. If the dust grains stay at the same position and the amount does not change, the variation in the dust shell flux at each variability phase can be predicted from the derived $T_{\rm i}$ and compared to the observations. We use the integrated flux of the observed dust shell from 5 to 45 $\mu$m as the shell luminosity. There are some differences in the absolute flux levels between the SWS and PHT data (see Appendix A). In the present analysis we use the SWS flux to have consistency in the comparison of the model with the observations.

Figure 5: The variation in the integrated dust shell flux in the range 5-45 $\mu$m. The solid line shows the observed variation. The open triangles are the results from the best fit models with Q10, the open circles those with Q7, and the open squares those with Q4. The error bars indicate the uncertainties estimated from the SWS spectra and do not include the absolute calibration uncertainty.

Figure 5 shows the comparison between the observations and the models with Q₄, Q₇, and Q₁₀ for y₀ = 0.1. In the models the amount of the dust grains in the shell is fixed for each Q_i to provide the best match with all the observed spectra. We plot the uncertainties estimated from the quoted flux errors of the SWS spectra. The differences between the models are greatest at the maxima. If we take a high temperature dust model (Q₁₀), the model predicts too large a variation from minimum to maximum because it requires a large change in the dust temperature to explain the observed variation in the infrared spectral shape. In contrast the model of low dust temperature (Q₄) shows too small a variation between minimum and maximum since a small change in $T_{\rm i}$ with Q₄ can explain the observed spectral variations.

Figure 6: Reduced $\chi ^2$ values for the best fit models (see text for more details). The solid lines indicated the results for y0=0.1, the dotted line for y0=0.4, and the dot-dashed line for y0=0.

To make a quantitative comparison we plot reduced $\chi ^2$ values ( $=\sum_{i=1}^{7} (f_{{\rm model}, i} - f_{{\rm obs}, i})^2/\sigma_{i}^2$/6, where $f_{{\rm model}, i}$and $f_{{\rm obs}, i}$ are the model and observed fluxes, respectively, $\sigma_{i}$ is the uncertainty in the flux, and 6 is the number of degrees of freedom of the fit, the number of the observation epochs (7) minus the number of the dust amount scaling factor (1)) for each Q_i model in Fig. 6. The emissivity Q₇ gives the smallest $\chi ^2$ values for y₀=0, 0.1, and 0.4. The increase in $\chi ^2$ in the high temperature models depends on the assumed dust optical properties in the 5-7 $\mu$m region. We find that $\chi ^2$ is minimum for the emissivity Q₇ unless the real dust properties differ by more than a factor of 3 from those of the assumed astronomical silicate. Even if we change the 5-7 $\mu$m absorptivity by a factor of 10, the minimum is shifted to either Q₆ or Q₈. Thus the present results are not sensitive to the absorptivity in 5-7 $\mu$m. The models with the emissivity Q₇ provide the most consistent results with the observed variations in the spectral shape ($T_{\rm i}$) and the flux ($n_{\rm i}$). This can also be inferred from Fig. 3, where Q₇ gives the least variation in $n_{\rm i}$ over the variability phase. The model with Q₇ predicts a low integrated flux at the first minimum ( $\phi = 0.55$), which marginally agrees with the observed value, but it agrees with the observed integrated flux at the rest of the phases within the uncertainties. The discrepancy between $\phi = 0.55$ and 0.79 will be discussed in terms of possible dust formation near minimum in the next section.

Models with Q₆ and Q₈ differ insignificantly from the Q₇ model in $\chi ^2$ values. Comparison with the PHT data suggests that the SWS spectrum of the last observation ( $\phi = 1.93$) may be systematically brighter by approximately 15-20% relative to the observations at the other phases (Appendix A). If we reduce the integrated flux of the last observation by 20%, then we find that the Q₆ model gives a slightly smaller $\chi ^2$ value in the integrated flux than the Q₇ model. In the spectrum fit we also find that a change in $T_{\rm i}$ by 10% increases $\chi ^2$ by a factor of 2, which corresponds to about 20% change in the predicted flux level. This changes the minimum in $\chi ^2$ either to Q₆ or Q₈. Therefore, we estimate that the uncertainty in the determination of the most likely emissivity is about one division in the emissivity set (100 K) and conclude that the models with $T_{\rm i}=700 \pm 100$ K at visual maximum provide the best fit to both the spectral shape and the flux level. The derived inner dust temperature at maximum is somewhat low compared to theoretical predictions, but is not significantly different from the range 800-1000 K suggested by the latest investigation (Gail & Sedlmayr 1999).

The best fit for Q₇ implies a mass-loss rate $\dot{M} = 7 \times 10^{-8}$ (r_*/3 $\times$  $10^{13}~\hbox{cm})$  $(D/490~\hbox{pc})^2$  $(v_{\rm i}/5~\hbox{km s}^{-1})$ $(\delta/200)$ $M_{\odot}$, where $v_{\rm i}$ is the dust flow velocity at the inner shell boundary and $\delta$ is the gas to dust ratio. If we take $r_* = 3 \times 10^{13}$ cm, $\delta = 200$, $v_{\rm i} = 5$ km s^-1, and D=490 pc (Young 1995), the estimated mass-loss rate is in fairly good agreement with the results of the CO observations ( $4 \times 10^{-8}~M_{\odot}$, Young 1995). The optical depth at 9.7 $\mu$m is derived to be 0.08. It is larger than the previous estimate (Onaka et al. 1989b) because the present study derives a higher $T(r_{\rm i})$, but confirms that the optically thin approximation is still valid.