© ESO, 2012

1. Introduction

IRAS 16342 − 3814 (hereafter I16342), also known as OH 344.1+5.8, is an oxygen-rich proto(pre)-planetary nebula (PPN) exhibiting a striking pair of bipolar lobes. He et al. (2008) detected a low ¹²CO/¹³CO intensity ratio of 1.7 and concluded that ¹²C is being converted into ¹³C via the hot bottom burning process, suggesting that the initial mass of the central star is as high as  ≳ 5 M_⊙. Another characteristic of this object is high-velocity maser components. The velocities of OH masers at 1612, 1665, and 1667 MHz are  ~65 km s^-1, which are much faster than those in many OH/IR stars (Likkel & Morris 1988). The H₂O masers are even faster  ~130−180 km s^-1 (Zuckerman & Lo 1987; Likkel & Morris 1988). The apparent dynamical ages of the masers are very short at  ≲ 100 yr (Claussen et al. 2009; Imai et al. 2012). The H₂O masers are more like jets in appearance than the typical AGB winds. Because of these maser properties, objects like I16342 are called “water fountain source (WFS)” (Likkel & Morris 1988; Imai et al. 2007). So far, 14 objects are classified into this group (Imai et al. 2007; Suárez et al. 2009; Walsh et al. 2009; Gómez et al. 2011). Many PNs exhibit jets, which cannot be explained with only the generalised interacting stellar wind (GISW) model (Sahai & Trauger 1998). The WFSs are thought to be progenitors that provide crucial clues to understanding the shaping mechanisms of these PNs.

Because of the peculiar maser properties in WFS, many works have studied H₂O and OH masers. In I16342, the best-studied object, the circumstellar dust shell (CDS) has been studied. The optical and near-infrared (NIR) images show a distinct bipolar appearance with an extension of  ~5″ lying at a position angle of 73° (Sahai et al. 1999). In the optical, the eastern lobe is fainter by a factor of  ~0.1 with respect to the western one. Sahai et al. (1999) estimated the viewing angle to be  ~40° from the 1612 MHz OH maser distribution, i.e. the western polar axis is tilted by  ~40° toward the observer. The Keck II L′-band images show a corkscrew-like feature in the bipolar lobes (Sahai et al. 2005). Although a clear correlation with the distribution of the H₂O maser jets is not found, the feature is probably a kind of precessing jet that carved out the inner part of the bipolar lobe. In the mid-infrared (MIR), an elliptic flux structure was detected (Meixner et al. 1999; Dijkstra et al. 2003). However, the high-resolution VISIR/VLT data resolved a pair of bipolar lobes with a separation of the flux peaks of 0$\stackrel{\mathrm{″}}{.}$92 at 11.85 μm (Verhoelst et al. 2009; Lagadec et al. 2011). These results suggest an optically and geometrically thick torus in the equatorial plane of this object.

In this present work, we have modeled the dust shell of I16342 by means of two-dimensional radiative transfer calculations. The spectral energy distribution (SED) collected from various sources and our new NIR polarimetric data obtained using the VLT/NACO instrument were used to constrain the physical parameters of the CDS. In Sects. 2 and 3, our VLT/NACO observations and radiative transfer calculations are presented. In Sect. 4, we discuss the structure and the formation of the inner most region of the dust shell.

2. VLT/NACO imaging polarimetry

2.1. Observations and data reduction

We obtained standard (non-polarimetric) images in the J-band and polarimetric images in the H and K_S band using the CONICA camera with the NAOS adaptive optics (AO) system mounted on the Very Large Telescope (VLT). The S27 camera with a pixel scale of 27.15 mas pix^-1 was operated in the Double_RdRstRd read-out mode. For the wavefront sensing, the N20C80 dichroic filter was used. Because the target is extended in the optical and near-infrared and the central star feature is invisible, a single star 2MASS J16374030–3820087 (m_K = 8.0 mag) at a 9$\stackrel{\mathrm{″}}{.}$7 separation from the target was monitored for the AO reference.

The J-band imaging was carried out on August 4, 2009. The natural seeing was 0$\stackrel{\mathrm{″}}{.}$5 − 0$\stackrel{\mathrm{″}}{.}$6 according to the ESO observatories ambient conditions data base. The exposure time was 30 s. per frame and six images were taken at each four jitter positions, yielding a total integration time of 12 min. The raw data were reduced by subtracting the dark frame, dividing by the flat frame, and combining the all science data. We used the photometric data of GSPC S875-C (m_J = 11.085 mag) for the flux calibration, which was offered in the ESO standard calibration plan. The signal-to-noise ratios per beam is 40 − 240 in regions with the surface brightnesses of 2 × 10^-15 W m^-2   μm^-1 arcsec^-2 or brighter.

The H and K_S band imaging polarimetry was carried out on 2 May, 2009. The natural seeing varied between 1$\stackrel{\mathrm{″}}{.}$0 and 1$\stackrel{\mathrm{″}}{.}$5 during the observations. The turnable half-waveplate (HWP) and the Wollaston prism with a slit width of 3″ were used to measure linear polarization. The slit position is aligned along the long axis of the target at a position angle of 73°. Four image sets were taken at position angles of the HWP of 0°, 45°, 22.5°, and 67.5°. Seven and ten jitter positions were performed and exposure times were 60 s. and 30 s. per frame, in the H and K_S band, respectively. The total integration times were 28 min. in the H band and 20 min. in the K_S band. After the subtraction of the dark frames, the sky background levels and flat fielding, the Stokes IQU parameter images were obtained by combining the ordinary and extra-ordinary ray frames. The data of a single star HD 147283 were used for polarization calibration (P_H = 0.91    ±    0.05% and P_K = 0.48    ±    0.04%; Whittet et al. 1992). Our uncalibrated measurements are P_H = 3.6% and P_K = 1.7%. These discrepancies include the instrumental polarization and some other error sources. In observations that require a polarization accuracy better than 1%, such as the line spectropolarimetry, a careful calibration is essential (see Witzel et al. 2011). On the other hand, an  ~ 2% accuracy is acceptable in most cases of optical and infrared imaging polarimetry of bipolar reflection nebulae. The difference of these tolerances also results from the different interests of the observations. In imaging polarimetry, the spatial information is important, but is affected by the photon noise and the temporal variations in the shape of the point-spread function (PSF) because of the use of small “aperture” sizes, e.g. beam sizes. It is opposite in the aperture polarimetry or line spectropolarimetry. We assumed that the discrepancies between our observations and the previous aperture polarimetry (Whittet et al. 1992) are caused by the instrumental polarization and the offset Stokes Q/I and U/I values. The estimated signal-to-noise ratios per beam of the Stokes I images in regions with surface brightness of 2    ×    10^-15 W m^-2 μm^-1 arcsec^-2 or brighter are 60 − 120 in the H band and 20−120 in the K_S band. The errors of linear polarization are P_H ~ 3% and P_K ~ 2%.

Fig. 1 Radial plot of the PSFs in the J (red), H (green), and KS (blue) bands. The surface brightnesses are normalized. The observations are compared with the two one-dimensional PSF models of Moffat and Gaussian functions.

Fig. 2 J, H, and KS band images of IRAS 16342–3814 obtained using NACO on the VLT. The polarization vector lines are overlaid on the H and KS band images. The field of view is $\mathrm{2}\stackrel{\mathrm{″}}{.}\mathrm{7}\hspace{0.17em}\mathrm{\times }\hspace{0.17em}\mathrm{5}″$. The images are rotated counterclockwise by 107°. The color scale bar indicates the surface brightness in W m-2   μm-1 arcsec-2.

We evaluated the shapes and sizes of the beam. In the J band, a single star detected in the science frame was examined. In the H and K_S band the standard star was used because no single stars are detected in the science frames. We found that most of flux from the reference stars is detected within $\mathrm{1}\stackrel{\mathrm{″}}{.}\mathrm{5}$ from the central star, and the star feature is nearly circular. Therefore, the star feature was tangentially averaged and was examined as a one-dimensional PSF. Figure 1 shows the radial plot of the PSFs in the J, H and K_S bands. The solid curves are for the observations. The measured full width at half maximum (FWHM) are $\mathrm{0}\stackrel{\mathrm{″}}{.}\mathrm{11}$, $\mathrm{0}\stackrel{\mathrm{″}}{.}\mathrm{12}$, and $\mathrm{0}\stackrel{\mathrm{″}}{.}\mathrm{07}$ in the J, H and K_S bands, respectively. The real PSF sizes of the object would be somewhat larger than them because an offset AO referencing was applied. We examined a Gaussian model (e.g., Bendinelli et al. 1990) and a Moffat model (Moffat 1969). The model PSF forms are given by ${\mathit{F}}_{\mathrm{Gauss}}{}^{\mathrm{\right(}}\mathit{r}^{\mathrm{\left)}}\mathrm{=}\mathit{A}\exp \mathrm{[}\mathrm{-}\mathrm{1}\mathit{/}\mathrm{2}{{}^{\mathrm{(}}\mathit{r}\mathit{/}{\mathit{\sigma }}^{\mathrm{)}}}^{\mathrm{2}}\mathrm{]}$ and ${\mathit{F}}_{\mathrm{Moffat}}{}^{\mathrm{\right(}}\mathit{r}^{\mathrm{\left)}}\mathrm{=}\mathit{A}\mathit{/}\mathrm{[}\mathrm{1}\mathrm{+}{{}^{\mathrm{(}}\mathit{r}\mathit{/}{\mathit{\sigma }}^{\mathrm{)}}}^{\mathrm{2}}{\mathrm{]}}^{\mathit{\beta }}$, where σ is the beam radius and β determines the overall shape of the PSF. The resulting parameters and the radial plot of the modeled PSF are listed in Table 1 and Fig. 1, respectively. The Gaussian function is only able to fit around the central star and can provide a good estimate of the beam size. On the other hand, the Moffat function fits better throughout the entire PSF. Although the resulting model of this simple analysis would be inaccurate compared to more sophisticated image reconstruction methods such as the Richardson-Lucy deconvolution, it is useful for our purpose of comparing images produced by radiative transfer modeling with the observations.

2.2. Results

Figure 2 shows the intensity (Stokes I) images in the J, H, and K_S bands. The polarization vector lines are overlaid on the H- and K_S band images. The intensity images show a peanut-like pair of bipolar lobes with approximate extensions of $\mathrm{2}\stackrel{\mathrm{″}}{.}\mathrm{5}\mathrm{\times }\mathrm{4}″$ and the appearance is nearly constant in all bands. The western (upper) lobe is about three times brighter than the eastern (lower) one. These results are consistent with Keck II images presented by Sahai et al. (2005). Our images show uneven flux structures in the bipolar lobes but do not clearly detect the corkscrew-like structures that are seen in the Keck II L_p-band image.

The degree of polarization reaches P_H ~ 25% and P_K ~ 20% in the brighter lobe and P_H ~ 38% and P_K ~ 38% in the fainter lobe. The polarization vectors are aligned along the equatorial plane. The physical reason is the PSF smoothing effect, namely, the highly scattered polarized light in the bipolar lobes spreads with an extended PSF halo component (see Fig. 1). This cancels the polarization component perpendicular to the equatorial direction and the parallel component remains. Indeed, our model results convolved with the modeled PSF reproduce the vector alignment as seen in Fig. 5.

One of the aims of our NIR imaging polarimetry was to detect the three-dimensional dusty structure in I16342’s CDS. If the corkscrew-like pattern, as detected in Sahai’s L_p band image, is formed by a jet, the degree of polarization would gradually vary along the spiral because the scattering angle also varies along it. Another is the hollow structure in the bipolar lobe, which could be detected along the rim of the bipolar lobe as a polarization enhancement for the optically thin case or the opposite for the optically thick case. The hollowness is most likely the result of an an interaction between the fast post-AGB wind and the slow AGB wind, as discussed in the (G)ISW model (e.g., Kwok 1982; Balick et al. 1987; Icke 1988; Soker & Livio 1989). Although an uneven distribution of polarization is detected in our data, we are not able to provide suitable interpretations that explain these factors. Therefore, we used the polarization data only to estimate the particle sizes in the bipolar lobes, but not to discuss the dust shell structure in the other sections.

3. Dust shell model

We have performed radiative transfer calculations of I16342’s CDS using our Monte Carlo code (Murakawa et al. 2008a). The SED, and intensity and polarization images are considered simultaneously, as done before (Murakawa et al. 2008b, 2010a). We first tried some test models to obtain approximated solutions and to understand which parameters affect which results. Then, many parameter sets (several thousands in total) were examined by SED fit, which takes less computation time. From this result, a dozen good parameter sets were chosen and their images were evaluated. We finally selected one that reproduced the shape of the SED, bipolar appearance, and obscuration of the central star in the optical to NIR, and had nearly constant NIR polarizations in the bipolar lobe.

Sections 3.1–3.3 describe the model assumptions and we compare the model results with the observations in Sect. 3.4.

3.1. Stellar parameters

The stellar temperature T_eff was estimated to be  ~15 000 K from the optical photometry with a 15″ aperture (van der Veen et al. 1989), whereas a strong 2.3 μm CO absorption was detected with a 5″ aperture, indicating T_eff ≲ 3500   K. Sahai et al. (1999) pointed out that this strong discrepancy is due to a contamination by the fluxes of two background stars (m_V = 14.9 and 13.9 mag) in the 15″ aperture, which are brighter than I16342 (m_V = 15.6 mag). Because the central star is invisible in the optical and infrared, no precise investigation of the stellar properties have been reported. We assumed a T_eff of 3000 K. The stellar luminosity and the distance are also not well constrained. Therefore, we adopted a luminosity of 6000 L_⊙, typical for AGB and post-AGB stars, and a distance of 2 kpc, as considered by Sahai et al. (1999, 2005). The central star, as a radiation source, was assumed to have a blackbody spectrum. We find that an effective temperature of 3000 K fits better than 15 000 K.

3.2. Model geometry

Several geometry forms have been proposed for the CDSs around AGB and post-AGB stars (e.g., Kahn & West 1985; Meixner et al. 1999; Ueta & Meixner 2003; Oppenheimer et al. 2005). Those proposed by Kahn & West (1985), Meixner et al. (1999), and Ueta & Meixner (2003) have a latitudinal density gradient with a higher density at the equatorial region. Although these structures reproduce fan-shape morphologies for high equator-to-polar mass density ratios, e.g.  ≳ 10, the interaction of the stellar winds and an optically thick disk are not explicitly expressed. Oppenheimer et al. (2005) used one including an inner disk, a pair of bipolar lobes, and an outer spherical AGB shell. We employed a similar form for our previous work of a PPN M 1–92 where the bipolar appearance was reproduced well with this geometry (Murakawa et al. 2010a).

We began with the latter form for I16342. However, we encountered a problem. The SED of this object has moderate NIR and strong MIR fluxes. At first glance, this (semi-) double-peaked flux structure can be explained with the standard classification of the post-AGB objects (Hrivnak et al. 1989; Kwok 1993). In this evolutionary phase, the intensive mass loss stops and the dust shell detaches from the central star. This geometry with a large inner boundary causes less emission from hot dust at 3 to 10 μm and increases the scattered light in the optical-to-NIR. On the other hand, I16342 is expected to have a torus that is optically thick even in the MIR (Verhoelst et al. 2009; Lagadec et al. 2011). With the aforementioned geometry form, we could not find a good solution to fit the SED and to reproduce a bipolar appearance in the NIR-to-MIR simultaneously. We attempted a geometry that includes an additional disk inside the inner disk, i.e. the inner region consists of an inner disk and outer torus. The mass density form of the entire dust shell is given by $\begin{array}{ccc}\mathit{\rho }& \mathrm{=}& {\mathit{\rho }}_{\mathrm{disk}}\mathrm{+}{\mathit{\rho }}_{\mathrm{torus}}\mathrm{+}{\mathit{\rho }}_{\mathrm{lobe}}\mathrm{+}{\mathit{\rho }}_{\mathrm{AGB}}\mathit{,}\\ {\mathit{\rho }}_{\mathrm{disk}}\mathrm{(}\mathit{r,}{\mathit{z}}^{\mathrm{)}}& \mathrm{=}& {\mathit{\rho }}_{\mathrm{d}}{\mathrm{(}\mathit{r}\mathit{/}{\mathit{R}}_{\mathrm{disk}}\mathrm{)}}^{-2.25}\exp \left(\mathrm{-}\frac{{\mathit{z}}^{\mathrm{2}}}{\mathrm{2}{\mathit{h}}^{\mathrm{2}}\mathrm{(}{\mathit{r}}^{\mathrm{)}}}\right)\mathit{,}\\ \mathit{h}\mathrm{(}{\mathit{r}}^{\mathrm{)}}& \mathrm{=}& \mathrm{0.1}{\mathit{R}}_{\mathrm{disk}}{\mathrm{(}\mathit{r}\mathit{/}{\mathit{R}}_{\mathrm{disk}}\mathrm{)}}^{\mathrm{1.25}}\mathit{,}\\ {\mathit{\rho }}_{\mathrm{torus}}\mathrm{(}{\mathit{R}}^{\mathrm{)}}& \mathrm{=}& {\mathit{\rho }}_{\mathrm{t}}{\mathrm{(}\mathit{r}\mathit{/}{\mathit{R}}_{\mathrm{disk}}\mathrm{)}}^{-1.8}\mathit{,}\\ {\mathit{\rho }}_{\mathrm{lobe}}\mathrm{(}\mathit{R,}{\mathit{\theta }}^{\mathrm{)}}& \mathrm{=}& {\mathit{\rho }}_{\mathrm{e}}{\mathrm{(}\mathit{R}\mathit{/}{\mathit{R}}_{\mathrm{lobe}}\mathrm{)}}^{-2}\mathrm{\times }\{\begin{array}{c}\\ \\ \end{array}\\ {\mathit{\rho }}_{\mathrm{AGB}}\mathrm{(}\mathit{R,}{\mathit{\theta }}^{\mathrm{)}}& \mathrm{=}& {\mathit{\rho }}_{\mathrm{e}}{\mathrm{(}\mathit{R}\mathit{/}{\mathit{R}}_{\mathrm{lobe}}\mathrm{)}}^{-2}\mathit{,}\end{array}$where the coordinates ${\mathrm{(}}^{}\mathit{R,}{\mathit{\theta }}^{\mathrm{)}}$ and ⁽r,z⁾ = (Rsinθ,Rcosθ) are the two-dimensional spherical and cylindrical coordinates. The disk and torus are placed in the regions where R_in ≤ R ≤ R_disk and R_disk ≤ R ≤ R_torus, respectively. The most plausible definition of a disk should be a region where the matter is gravitationally bounded by the central star and has a Keplerian rotating motion. Although dusty disk structures have been spatially resolved in several evolved stars (e.g., Matsuura et al. 2006; Lykou et al. 2011), the dynamics are not well known or much less understood than accretion disks of young stellar objects. Therefore, we regarded the disk simply as a geometrically thin, optically thick dust structure. In Sect. 3.4 we show that such an inner disk is required in I16342.

The outer envelope consists of the bipolar lobe and spherical AGB shell. For the bipolar lobe, the regions ${\mathit{R}}_{\mathrm{torus}}\mathrm{\le }\mathit{R}\mathrm{\le }{\mathit{R}}_{\mathrm{lobe}}\left(\mathrm{\Gamma }{}^{\mathrm{\right(}}\mathit{\theta }^{\mathrm{\left)}}\mathrm{-}\mathit{\gamma }\right)$ and ${\mathit{R}}_{\mathrm{lobe}}\left(\mathrm{\Gamma }{}^{\mathrm{\right(}}\mathit{\theta }^{\mathrm{\left)}}\mathrm{-}\mathit{\gamma }\right)\mathrm{\le }\mathit{R}\mathrm{\le }{\mathit{R}}_{\mathrm{lobe}}$ are the inner cavity and the rim of the lobe, respectively, where $\mathrm{\Gamma }{}^{\mathrm{\right(}}\mathit{\theta }^{\mathrm{\left)}}\mathrm{=}{\left|\mathit{\theta }\mathit{/}\mathit{\pi }\mathrm{-}\mathrm{1}\mathit{/}\mathrm{2}\right|}^{\mathit{\beta }}$ determines the lobe shape (see Oppenheimer et al. 2005). The AGB shell is placed outside the bipolar lobe with an outer boundary of R_out. The mass density coefficients are determined by the masses of the corresponding components. The dust mass is converted by assuming a gas-to-dust mass ratio of 200.

Fig. 3 Mass density map showing the inner 3000 AU region in a plane cutting through the symmetry axis.

Figure 3 shows the mass density map in the inner 3000  ×  3000 AU region. Our modeled geometry is complicated and it is in essence impossible to determine all parameters accurately from our observational data. Therefore, some parameters are set to be fixed values. Important parameters such as the inner disk radius R_disk, the disk mass M_disk, the torus radius R_torus, and the torus mass M_torus are estimated from our radiative transfer calculations.

Fig. 4 Comparison of SEDs. Black lines are the result of the model with an inner disk. The solid, dotted, and dash-dot lines denote the total flux, the scattered light, and the thermal emission, respectively. Red and blue lines denote the results of the same geometry form but the disk mass of 0 and a model with only a torus with an inner radius of 500 AU. Colored dots are collected from previous observations: US Naval Observatory (USNO), the HST Guide Star Catalog (GSC)-II, 2MASS all-sky catalog of point source, Catalog of Infrared Observations, Edition 5 (Gezari et al. 1999), ISO photometry, SWS and LWS, MSX infrared point source catalog, IRAS point source catalog, and Akari/IRC mid-IR all-sky Survey (Ishihara et al. 2010).

3.3. Dust particles

Dijkstra et al. (2003) discussed the mineralogy of I16342 based on their ISO spectra obtained using the Short Wavelength Spectrometer (SWS). Their data show a variety of absorption and emission features attributed to amorphous and crystalline silicate, frostetite, crystalline water ice, etc. Investigating the dust chemistry is one of the interesting and important tasks in dust formation of evolved stars. However, a thorough treatment of so many different chemical compositions and identification of the spatial distributions of individual species are very complicated tasks and are, therefore, beyond the scope of this paper. In our modeling, we adopted simple dust models. We used the optical constants for oxygen-deficient silicate (Ossenkopf et al. 1992).

We assumed different particle sizes for the outer envelope, the torus, and the inner disk, but the same sizes for the bipolar lobes and AGB shell. The particle size in the outer envelope was estimated with the polarization data. Since the polarizations are high and do not change much in the H- and K_S-bands, the fraction of small particles is expected to be higher than that of the interstellar population. In our separated paper (Murakawa et al. 2012), we found a similar result in PPN I18276 and that a size distribution function with a steep power index provides a good estimate: 0.05 μm  ≤ a and $\mathit{n}{}^{\mathrm{\right(}}\mathit{a}^{\mathrm{\left)}}\mathrm{\propto }{\mathit{a}}^{-5.5}\exp {}^{\mathrm{(}}\mathrm{-}\mathit{a}\mathit{/}{\mathit{a}}_{\mathrm{c}}{}^{\mathrm{)}}$, where a_c is the cut-off size (Kim et al. 1994). In I16342, a_c = 5.0   μm fits the H and K_S band polarizations.

Although large grains are expected in the disk and torus, we do not have (sub-)millimeter flux data or their gas masses and therefore we cannot constrain the sizes. Therefore, we provisionally adopted a size distribution of 0.005 μm  ≤ a ≤ a_max μm and $\mathit{n}{}^{\mathrm{\right(}}\mathit{a}^{\mathrm{\left)}}\mathrm{\propto }{\mathit{a}}^{-3.5}$ (Mathis et al. 1977). a_max for the disk and torus are set to be 100 μm and 10 μm, respectively. We will not discuss dust growth in this object with our results.

3.4. Results

Table 2 lists the parameters of the selected model. Figure 4 shows the SEDs. The black lines are the result of the selected model (model 1). The viewing angle was chosen to be 30° measured from the equatorial plane. Although this result is slightly different compared to that (40°) estimated by Sahai et al. (1999), this discrepancy is not important in the following discussion. The model result fits the observation well from the optical to FIR, except for the 10 μm silicate absorption feature.

We justify the presence of a geometrically thin, inner disk. First, we compared with a model without an inner disk (model 2), in which M_disk is set to be 0 and the other parameters are identical. The difference is obvious in the wavelength range between 3 μm and 10 μm. Although the disk mass in model 1 is only 0.01 M_⊙, the MIR flux is governed by the thermal emission from the hot dust in the inner disk. The other consideration is a model without an inner disk where we tried to fit the observations as best we can (model 3). This is our first attempt in the modeling (see Sect. 3.2). We find that models with small inner radii (a few tens AU) show too strong thermal emission at 3 μm to 10 μm and that a model with a large radius of R_in = 500 AU and M_torus = 2.0   M_⊙ fits best. However, as Fig. 4 shows, the flux between 3 μm to 6 μm is too low and model 1 fits better. The interpretation of I16342’s inner region is summarized as follows. An optically and geometrically thick torus should exist, as previous observations have shown. The torus is expected to have a large inner radius (a few hundred AU) and a high mass (1 M_⊙). An additional inner geometrically thin disk is required, which is responsible for the thermal emission in the 3 μm to 10 μm.

Fig. 5 Results of our two-dimensional radiative transfer calculations. The wavelengths are the 1.25 μm (J-band), 1.65 μm (H-band), and 2.05 μm (KS-band) from the left panels. The color scale bar indicate the surface brightness in W m-2   μm-1 arcsec-2. The Moffat PSF model images are convolved (see Sect. 2.1).

Figure 5 shows the monochromatic model images at 1.25 μm (J-band), 1.65 μm (H-band), and 2.05 μm (K_S-band) from left to right, which are simulated by scattering and absorption only. In these images, the Moffat PSF model images (see Sect. 2.1) are convolved. We confirmed that the Gaussian PSF convolution produces a centro-symmetric polarization pattern in the bipolar lobe, but the vectors are aligned with the Moffat PSF model. This is quite understandable because the surface brightness of a Gaussian function decays much faster than a Moffat function. The vector alignment in the entire nebula is reproduced by the PSF smoothing effect. Compared with the observations, the striking bipolar appearance and the presence of a dark lane in the equatorial plane are well reproduced. Although we do not present the results, we also confirmed that our model image in the 11.85 μm shows a bipolar appearance, as seen in the results of Verhoelst et al. (2009) and Lagadec et al. (2011).

4. Discussion

The mechanisms of PN shaping have been discussed in a number of publications before. To form a striking bipolarity with an equatorial narrow waist, the mass loss material should be highly concentrated in the equatorial plane. Morris (1987) modeled mass loss in a binary system from the primary red giant. Depending on the physical conditions, e.g. the mass and the separation of the secondary, a part of the mass loss material from the primary can be transferred onto the secondary, which forms a circum-companion disk. A part of accreting matter can be blown along the polar direction as a jet. Many bipolar PNs and PPNs, presumably I16342 too, are thought to form by the binary interaction based on this idea. Both theories and observations have developed the binary interaction hypothesis.

Silicate carbon stars are objects in which oxygen-rich species exist around carbon stars. According to the current interpretation of the dual chemistry, the silicate dust that was ejected when the central star was an M giant is stored in a disk around the hypothetical binary companion and then the central star is transformed into a carbon star through the third dredge-up (Lloyd-Evans 1990). V778 Cyg is an object of this class, which shows the 10 μm silicate emission feature, see Yamamura et al. (2000). These authors found that the dust has sub-micron sizes and temperatures of 300 − 600 K and is located at about 24 AU from the primary. The present-day mass-loss rate is estimated to be about 10^-8 M_⊙ yr^-1. They interpret that the mass loss material can be trapped in the circum-companion disk, but leaves soon after within a few orbital periods because the disk is optically thin and the dust is blown away by the strong radiation pressure from the primary. With these low mass-loss rates, the disks are unlikely to be responsible for a bipolar appearance. It is of interest to mention the results of the smoothed particle hydrodynamic simulations for high mass-loss rate cases of  ≳ 10^-5 M_⊙ yr^-1 by Mastrodemos & Morris (1998, 1999,hereafter M98 and M99). Their simulations assumed a mass-losing AGB primary and a main-sequence companion (M_S = 0.25 − 2   M_⊙) with orbital separations of 3.6 − 50 AU. In M 1, 2, 10, 12, and 15 out of their 18 models, the density ratios of the equator to the pole ρ_eq/ρ_pol exceeds  ~ 10 and the nebulae are classified as bipolar. M 1 and 2 are extreme cases, where the separations of the companion is small (6.3 AU in their models) and ratios of the accretion rate to the mass-loss rate are as high as  ~ 0.1, yielding Ṁ_acc ≳ 10^-6   M_⊙ yr^-1. ρ_eq/ρ_pol exceeds 100. These systems are likely to be progenitors of bipolar PNs with narrow waists. The thick winds sufficiently shield the strong radiation from the primary and the mass loss matter is stabilized there for a long time (Yamamura et al. 2000; Soker 2000), as expected in the Red Rectangle (Jura et al. 1995; van Winckel et al. 1998), AFGL 2688 (e.g., Bieging & Nguyen-Q-Rieu 1996), and probably I16342.

We now consider a simple model of AGB mass loss to see the competition of the radiation and the shielding based on a similar concept to “equatorial back-flow” introduced by Soker (2000). We assumed that a dust particle is directly irradiated by the radiation from the central star. The ratio of the radiation pressure to the gravitational attraction β_dust is given by β_dust = L _⋆ C_rp/4πGM _⋆ c, where C_rp is the cross section of the radiation pressure per unit mass, G the gravity constant and c the velocity of light. Applying the stellar parameters of L _⋆  = 6000   L_⊙ and M _⋆  = 1   M_⊙, a dust model that is derived from the outer envelope in our modeling, C_rp = 6.73    ×    10³ cm² g^-1, where the stellar flux is assumed to dominate at λ = 1   μm, β_dust is found to be  ~ 3    ×    10³. This high value is quite reasonable in AGB winds, where the dust is blown away by the radiation pressure. Next, the shielding by the mass loss ejecta is considered. We assumed a spherically symmetric mass loss at a constant Ṁ and constant radial expansion velocity V_r. The mass density has a power law distribution $\mathit{\rho }{}^{\mathrm{\right(}}\mathit{r}^{\mathrm{\left)}}\mathrm{\propto }{\mathit{r}}^{-2}$, where r is the distance from the central star. The optical depth between the inner boundary (r = R_in) and a position at r = R is given by $\begin{array}{ccc}\mathit{\tau }\mathrm{=}\frac{{\mathit{C}}_{\mathrm{ext}}\mathit{Ṁ}}{\mathrm{4}\mathit{\pi }{\mathrm{V}}_{\mathit{r}}\mathit{\gamma }}\left(\frac{\mathrm{1}}{{\mathit{R}}_{\mathrm{in}}}\mathrm{-}\frac{\mathrm{1}}{\mathit{R}}\right)\mathit{,}& & \end{array}$(6)where C_ext is the extinction cross section per unit mass of the dust and γ the gas-to-dust mass ratio of 200. If the mass-loss rate reaches 10^-5 M_⊙ yr^-1 at an expansion velocity of 15 km s^-1, and R ≫ R_in = 5 AU is assumed, we derive an optical depth $\mathit{\tau }\mathrm{\approx }\mathrm{20.3}\left(\frac{\mathrm{15} \mathrm{km}\hspace{0.17em}{\mathrm{s}}^{-1}}{{\mathit{V}}_{\mathrm{e}}}\right)\left(\frac{\mathit{Ṁ}}{{\mathrm{10}}^{-5} {\mathit{M}}_{\mathrm{\odot }}\hspace{0.17em}{\mathrm{yr}}^{-1}}\right)\left(\frac{\mathrm{5} \mathrm{AU}}{{\mathit{R}}_{\mathrm{in}}}\right)\mathit{,}$(7)where C_ext = 9.10    ×    10³ cm² g^-1. The dust shell is optically thick. When dust particles move away, the radiation is shielded by the dust that is formed by the subsequent mass loss event. The optical depth attains  ~8.0   ( = lnβ_dust) at R ~ 8.3 AU and the attenuated β_dust values becomes  ~1. The radiation pressure force cannot govern the motion of the particles. If the system has a close binary, the gravitational attraction of the secondary can affect the motion of the mass loss ejecta and can concentrate the flow toward the orbital plane. This leads to a flatter geometry, i.e. disk.

We have a rough scenario of the inner part of I16342’s CDS as follows. The inner disk of I16342 in our model, which is introduced to be responsible for the MIR flux, is a circum-companion or circum-binary disk that is formed by this mechanism. Because the optical depth of the disk is so high at  ~ 400 in the V-band that the effect of the radiation is negligible, the disk can be gravitationally bounded by the central star(s). On the other hand, the torus is formed in a somewhat different way. If we assumed that I16342 has a binary companion, an angular momentum that is supplied by a close companion (say M_S ~ 1   M_⊙) is insufficient for the entire torus to have a keplerian rotating motion. Therefore, the torus is probably a structure that is accumulated from dust that is escaping from the central star(s).

5. Conclusion

We modeled the CDS of the WFS IRAS 16342–3814 with two-dimensional radiative transfer calculations. The model parameters were constrained by fitting the SED collected from various sources and the polarization images obtained using the VLT/NACO instrument. I16342 shows a striking bipolar appearance from the optical to the MIR wavelength ranges. Taking into account a large viewing angle of 30° − 40°, an optically thick dust torus structure is expected in the equatorial plane. We first tried to fit a model consisting of a geometrically thick torus and an outer envelope with bipolar lobes and a spherical AGB shell. However, these models cannot reproduce the SED and the bipolar appearance simultaneously. We found that an additional component of an optically thick, geometrically thin disk is required inside the torus. In our modeling, different dust sizes were assumed for the disk, the torus, and the outer envelope. The particle sizes in the outer envelope were estimated from the polarization data. The polarization reaches 25 − 38% in the H band and 20 − 38% in the K_S band in our observations. This weak wavelength dependence is explained with a steeper power index of  − 5.5 of the size distribution function than that of  − 3.5 for the ISM. Similar results are found in some other PPNs. Because of the submillimeter and millimeter flux information and the estimation of the disk are missing, the particle sizes of these component cannot be constrained well. We adopted dust models of a_max = 100   μm for the disk and a_max = 10   μm for the torus. With this assumption, the masses of these components were determined to be 0.01 M_⊙ and 1 M_⊙, respectively, but have large uncertainties.

We discussed a possible formation scenario for the disk and torus. The important factors are the mass-loss rate and the binary interaction. For low- to intermediate-rates, the system is likely to have an optically thin circum-companion disk. To form an optically thick disk, the mass-loss rate would need to be so high that the thick wind blocks the radiation from the primary. This allows a part of the mass loss ejecta to be condensed in the equatorial region and to stay there for a long time. The remains, which are leaving from the central star(s), form the torus.

Acknowledgments

The NIR polarimetric images were obtained at the European Southern Observatory (proposal ID: 383.D-0197). The photometric and spectroscopic data are data produces from the Two Micron All Sky Survey, which is a joint project of the University of Massachusetts and the Infrared Processing and Analysis Center/California Institute of Technology, funded by the National Aeronautics and Space Administration and the National Science Foundation, AKARI, a JAXA project with the participation of ESA, Infrared Space Observatory, and Infrared Astronomical Satellite. We thank Luke T. Maud for his encouraging discussions and proofreading of the paper.

References

All Tables

All Figures

	Fig. 1 Radial plot of the PSFs in the J (red), H (green), and KS (blue) bands. The surface brightnesses are normalized. The observations are compared with the two one-dimensional PSF models of Moffat and Gaussian functions.
In the text

Fig. 2 J, H, and KS band images of IRAS 16342–3814 obtained using NACO on the VLT. The polarization vector lines are overlaid on the H and KS band images. The field of view is $\mathrm{2}\stackrel{\mathrm{″}}{.}\mathrm{7}\hspace{0.17em}\mathrm{\times }\hspace{0.17em}\mathrm{5}″$. The images are rotated counterclockwise by 107°. The color scale bar indicates the surface brightness in W m-2   μm-1 arcsec-2.

In the text

	Fig. 3 Mass density map showing the inner 3000 AU region in a plane cutting through the symmetry axis.
In the text

Fig. 4 Comparison of SEDs. Black lines are the result of the model with an inner disk. The solid, dotted, and dash-dot lines denote the total flux, the scattered light, and the thermal emission, respectively. Red and blue lines denote the results of the same geometry form but the disk mass of 0 and a model with only a torus with an inner radius of 500 AU. Colored dots are collected from previous observations: US Naval Observatory (USNO), the HST Guide Star Catalog (GSC)-II, 2MASS all-sky catalog of point source, Catalog of Infrared Observations, Edition 5 (Gezari et al. 1999), ISO photometry, SWS and LWS, MSX infrared point source catalog, IRAS point source catalog, and Akari/IRC mid-IR all-sky Survey (Ishihara et al. 2010).

In the text

	Fig. 5 Results of our two-dimensional radiative transfer calculations. The wavelengths are the 1.25 μm (J-band), 1.65 μm (H-band), and 2.05 μm (KS-band) from the left panels. The color scale bar indicate the surface brightness in W m-2   μm-1 arcsec-2. The Moffat PSF model images are convolved (see Sect. 2.1).
In the text