© ESO, 2010

1. Introduction

It is widely accepted that the disk is an important component of young stellar objects (YSOs) because its structure and dust properties reflect the physics of star formation, and it is a field for planet formation. In the past 20–30 years, T Tauri disks have been intensively studied (e.g., Adams & Shu 1986; Kenyon & Hartmann 1987; Bertout et al. 1988; Calvet et al. 1994; Chiang & Goldreich 1997). Their radiative transfer analyses and hydrodynamic treatments (see also Shakura & Sunyaev 1973; Pringle 1981) predict flared disk structures and can explain spectral energy distributions (SEDs) in the mid- (MIR) and far-infrared (FIR) in many objects. Beckwith & Sargent (1991) detected low dust opacity spectral indices β in several T Tauri stars from their submillimeter observations. This is the first discovery of evidence for large grains in disks of these objects.

These analytical methods have also been applied to the intermediate-mass counterparts. In terms of the infrared SED, the Herbig Ae / Be stars are classified into two groups (e.g. Meeus et al. 2001; Dullemond et al. 2001; Natta et al. 2001; Dullemond & Dominik 2004a; Isella & Natta 2005). Group I objects have a flat mid-infrared flux. The disks probably have a flared structure, similar to T Tauri stars. On the other hand, a significantly lower far-infrared flux is seen in Group II objects. In this object class, which is thought to be more evolved than the Group I objects, the outer part of the disk is deposited into the midplane and is shaded from the radiation of the central star by the puffed-up inner rim. In addition, the dust is thought to grow in size in this advanced phase. In fact, in some objects, such as CQ Tau (Testi et al. 2003) and TW Hya (Wilner et al. 2005), evidence of large grains has been obtained from millimeter flux excess and 10 μm spectra (e.g. van Boekel et al. 2003, 2005).

As these works have shown, investigations of disk properties benefit mainly from SED analyses in the thermal emission regime. On the other hand, other attempts by optical- to near-infrared (NIR) imaging polarimetry have been carried out (Bastien & Ménard 1988, 1990; Kenyon et al. 1993; Whitney & Hartman 1993; Fischer et al. 1994, 1996; Whitney et al. 1997; Lucas & Roche 1997, 1998; Wood et al. 1998). In polarization images of YSOs seen nearly edge-on, what is called a polarization disk, where the polarization vectors are aligned in the equatorial direction, is often detected at the dark lane. This feature is reproduced with models where an optically thick dust disk or torus structure exists around the central star (Bastien & Ménard 1988, hereafter BM).

Although a polarization disk does not directly indicate the real (dust) disk, i.e. a region with a Keplerian rotating matter within the centrifugal radius, detection of a polarization disk is, in fact, regarded as evidence of a disk according to the aforementioned robust works. However, the polarization analyses have been limited until recently, because the angular resolutions of the images were not high enough to resolve the disk, which causes pseudo polarization patterns. In addition, the BM scattering model is originally intended to explain the vector alignment in typical T Tauri disks. Disk properties for other object classes such as high-mass counterparts are more poorly understood. Recent high-resolution imaging techniques relieve the problem of the pseudo polarization pattern and allow us to improve our interpretations of YSO disks. In our recent work, we generalized the BM model and developed an analytical method for deriving information on the grain sizes and the disk geometry based on high-resolution imaging (Murakawa 2010, hereafter M10).

In our previous work, we obtained JHK-band polarimetric images of the Herbig Be star R Mon using the adaptive optics-equipped instruments VLT / NACO and Subaru / CIAO (Murakawa et al. 2008a, hereafter M08a). The important results are as follows. (1) The polarization disk has a butterfly shape with an  ~ 4″ extension. This size does not change much in the JHK bands. (2) The polarization in the polarization disk is as low as 1.1–1.5% in the HK bands. The polarization vectors towards the polarization disk are not clearly aligned along the equatorial plane, but are somewhat centro-symmetric. A comparison with submillimeter/millimeter observations by Fuente et al. (2003, 2006, hereafter F06) is also interesting. The extension of the polarization disk is closer to the CO gas phase disk (1500 AU  ~ 4″) and the central feature in the 2.7 mm continuum (3000 AU  ~ 8″) than is detected both in the 1.3 mm continuum and an estimated disk radius of 150 AU () by a previous modeling (Alonso-Albi et al. 2009). The disk height of the gas phase disk was estimated to be 1500 AU at the outer disk radius of 1500 AU (F06). In contrast, the appearance of the polarization disk looks thinner.

From the above results, some questions may rise: (1) what are the physical reasons for producing the low NIR polarizations and for not detecting the vector alignment; (2) what kind of disk properties are expected in R Mon from our analysis; and (3) is the dust distribution in the disk significantly different from the gas phase, i.e. do these results prove that the dust layer concentrates more towards the midplane than the gas phase molecules? To provide quantitative interpretations of R Mon’s disk, we performed radiative transfer calculations. In our modeling, we take not only our previous NIR polarimetric images into account but also the previously observed SED and CO observations (F06) to better constrain the model parameters. We also apply our generalized scattering model for the R Mon case and study the physical conditions for producing the polarization disk, as questioned above. Our approach demonstrates the power of a complementary combination with an SED and polarization analyses and a gas phase observation. In Sect. 2, we describe the details of our modeling and present the results. In Sect. 3, we show the effects of the grain size in the disk and the disk geometry to justify our model results and discuss the properties of the R Mon disk.

2. Radiative transfer modeling

The stellar system of R Mon has been extensively studied. We adopted a distance of 800 pc (Jones & Herbig 1982), a stellar mass of 10.4 M_⊙, and an effective temperature of 30 000 K (Hillenbrand et al. 1992). Cohen et al. (1984) estimated a luminosity of 1400 L_⊙. Close et al. (1997) discovered a companion with a luminosity of  ~50 L_⊙. Since this is much fainter than the primary, we assume that the spectrum of the radiation source has a single blackbody function with the aforementioned effective temperature and the luminosity of only the primary.

R Mon exhibits a complex monoconical lobe towards the north (Walsh & Malin 1985; Lightfoot 1989; Close et al. 1997; Murakawa et al. 2008a). However, since our purpose is to study the dust disk and not the detailed appearance of the outer envelope, we apply an axi-symmetric geometry with an inner disk and infalling envelope structure.

For the disk, we examine two different forms of the mass density distribution (Shakura & Sunyaev 1973; Lazareff et al. 1990; Fischer et al. 1996), which are given with the two-dimensional cylindrical coordinate by $\begin{array}{ccc}{\mathit{\rho }}_{\mathrm{disk}}^{\mathrm{F}\mathrm{1}}\mathrm{(}\mathit{r,}{\mathit{z}}^{\mathrm{)}}& \mathrm{=}& \\ \mathit{h}\mathrm{(}{\mathit{r}}^{\mathrm{)}}& \mathrm{=}& {\mathit{R}}_{\mathrm{disk}}{\mathit{H}}_{\mathrm{0}}{\mathrm{(}\mathit{r}\mathit{/}{\mathit{R}}_{\mathrm{disk}}\mathrm{)}}^{\mathrm{9}\mathit{/}\mathrm{8}}\mathit{,}\\ {\mathit{\rho }}_{\mathrm{disk}}^{\mathrm{F}\mathrm{2}}\mathrm{(}\mathit{r,}{\mathit{z}}^{\mathrm{)}}& \mathrm{=}& \\ \mathit{h}\mathrm{(}{\mathit{r}}^{\mathrm{)}}& \mathrm{=}& {\mathit{R}}_{\mathrm{disk}}{\mathit{H}}_{\mathrm{0}}{\mathrm{(}\mathit{r}\mathit{/}{\mathit{R}}_{\mathrm{disk}}\mathrm{)}}^{\mathit{\beta }}\mathit{,}\end{array}$where R_disk is the outer disk radius and H₀ the ratio of the disk height to the disk radius. The F1 form is most widely used for YSO disk models; however, this form does not guarantee a better fit to all observations, and an alternative form of F2 is also exceptionally considered (e.g. Lazareff et al. 1990; Fischer et al. 1996; Murakawa et al. 2008b). In this form, the values of α and β are free parameters. The difference between these models is the concentration of the matter particularly toward the midplane. While the F1 form explicitly indicates a flared structure, the F2 form can be a flat (non-flared) structure in the case of β = 1. The mass density coefficient ρ_d is determined with the specified disk mass M_disk, assuming a gas-to-dust mass ratio of 100. We set the inner disk radius R_in to be 5 AU, which is approximated as the dust sublimation radius for the assumed luminosity. Under hydrostatic equilibrium, the disk height is correlated with the pressure, the mass density, and the viscosity (Shakura & Sunyaev 1973; Pringle 1981; Chiang & Goldreich 1997). However, in our analysis, such physical quantities are either not clearly determined or not given from observations. Thus, the H₀ is parameterized and is determined by fits of SED and intensity and polarization images.

For the envelope, we use the model of the slowly rotating, infalling envelope (Ulrich 1976; Cassen & Moosman 1981; Terebey et al. 1984). The mass density distribution is given by $\begin{array}{ccc}{\mathit{\rho }}_{\mathrm{env}}\left(\mathit{R,\mu }\right)\mathrm{=}\frac{\mathit{Ṁ}\mathrm{env}}{\mathrm{4}\mathit{\pi }\sqrt{\mathit{G}{\mathit{M}}_{\mathit{\star }}{\mathit{R}}^{\mathrm{3}}}}{\left(\mathrm{1}\mathrm{+}\frac{\mathit{\mu }}{{\mathit{\mu }}_{\mathrm{0}}}\right)}^{\mathrm{-}\mathrm{1}\mathit{/}\mathrm{2}}{\left(\frac{\mathit{\mu }}{{\mathit{\mu }}_{\mathrm{0}}}\mathrm{+}\frac{\mathrm{2}{\mathit{\mu }}_{\mathrm{0}}^{\mathrm{2}}{\mathit{R}}_{\mathrm{c}}}{\mathit{R}}\right)}^{-1}\mathit{,}& & \end{array}$(3)where R, G, Ṁ_env, and R_c are the radial distance from the central star , the gravitational constant, the envelope mass infall rate, and the centrifugal radius, respectively. The quantity μ is the cosine of the polar angle μ = cosθ, and μ₀ is the angle of a stream line as related through the equation . The polar cavity is the region for μ > μ₀. The density inside the cavity is reduced to 0.01 compared to the surrounding envelope. From the large-scale optical images, the full opening angle of the NGC 2261 nebula was measured as  ~80° (Lightfoot 1989). We set the μ₀ to be cos⁽80°/2⁾. The centrifugal radius R_c is set to be R_disk; that is a region inside R_disk is regarded as the disk and outside as the envelope. Although clear definitions of the disk and envelope exist in theory, these are in fact difficult to distinguish. We estimate R_disk from our radiative transfer modeling. The extension of the nebula is 4′ (Lightfoot 1989), which corresponds to  ~0.93 pc. We apply this value for the outer radius R_out.

The absorption and scattering properties depend on various properties of the dust, such as the chemical composition, grain size distribution, and surface micro-structure. However, our goal is to investigate the grain size. We adopted spherical grains with the DL-chemistry (Draine & Lee 1984) to make it simpler. The dust was assumed to have an MRN-like size distribution (Mathis et al. 1977); 0.005 μm  ≤ a ≤ a_max   μm and . We applied multiple dust models; i.e. grains with different a_max for the disk and envelope. The a_max in the envelope was found to be 0.25 μm from a one-dimensional single scattering model based on our previous Subaru / CIAO data (M08a). However, this model produces a slightly higher polarization than the VLT / NACO data in the K_S band. Thus, we chose 0.5 μm, which results in an intermediate value. The disk was divided into three mass density domains with the thresholds of mass density ρ_th, and a different a_max was applied for each domain. In our modeling, we determine the a_max in the disk and the ρ_th values.

2.1. SED and polarimetric data

The photometric data in the UBVRIJHKLMNQ bands are from Lorenzetti et al. (1983), Cohen (1973), Hillenbrand et al. (1992), and Pezzuto et al. (1997). The MIR to FIR photometric data are from the Midcourse Space Experiment (MSX) point source catalog and the Infrared Astronomical Satellite (IRAS) point source catalog. The 10–20 μm spectra are from the IRAS low-resolution spectra (LRS). The submillimeter and millimeter photometric data are from Mannings (1994). The typical aperture sizes of the above photometric data are 6″–12″ in the infrared and 15″–22″ otherwise. Because the extension of the nebulosity (~4′) is much greater than these aperture sizes, the measured fluxes are lower than the total fluxes. We have confirmed this aperture size effect and find that factors of up to 2–3 lower fluxes in the optical, submillimeter, and longer wavelength ranges. On the other hand, these are less than a factor of 2 in the NIR to FIR, because a large fraction of fluxes are localized towards the central region. Since our interest is the disk instead of the outer envelope, we ignored the aperture size effect in the resulting SED. However, the interstellar extinction should be taken into account because R Mon is located towards the galactic plane (). We corrected the flux assuming the standard galactic extinction of 1.6 mag kpc^-1 (Allen 1973).

For the polarimetric data, we use our Subaru / CIAO data in the J and H bands and VLT / NACO data in the K_S band (M08a). The estimated point-spread function (PSF) sizes of these data are 0$\mathrm{″̣}$54, 0$\mathrm{″̣}$35, and 0$\mathrm{″̣}$066 in the J, H, and K_S bands, respectively. The model images are convolved with the Gaussian functions of the equivalent FWHM sizes to compare with the observations.

2.2. Results

To perform radiative transfer calculations, we used our own Monte Carlo STSH code, which can compute the SED, the dust temperature, and the Stokes IQUV images of models with arbitrary geometries (up to three-dimension) and multiple dust models (Murakawa et al. 2008b). We first tried some parameter sets to find approximate solutions and appropriate parameter grids and ranges to scan. The parameters and their ranges follow. For α and β of the F2 disk form, we chose (α,β) = (1.5,1.0), which fits the nebulous feature towards the central star and the appearance of the polarization disk better than (15 / 8,9 / 8). If an outer disk radius of 1500 AU, which was estimated from a gas phase observation by F06, is applied, the apparent size of the polarization disk becomes slightly smaller than the observed size. We find that 3000 AU is better. This size is, in fact, the same as an extension of the central feature detected at 2.7 mm (Fuente et al. 2003), although they considered that this feature may include the envelope. The inclination angle θ_inc of the disk was estimated to be 70° ± 10° measured from the polar axis (Close et al. 1997). We found that this is quite reasonable. Thus, we adopted θ_inc as 70°. For grain sizes, we tried a_max,1 of 100.0, 1000.0, and 10 000.0 μm in the highest density region of ρ ≥ ρ_th,1 and a_max,2 of 1.0 and 10.0 μm in the intermediate density region of ρ_th,2 ≤ ρ ≤ ρ_th,1. The examined threshold densities were ρ_th,1 of 5 × 10^-21 g cm^-3 and 1 × 10^-20 g cm^-3 and ρ_th,2 of 1 × 10^-21 g cm^-3 and 3 × 10^-21 g cm^-3. For others, we examined the ratio of the height to the radius of the disk H₀ of 0.1, 0.2, 0.3, 0.4, 0.5, and 1.0, the disk mass of 0.01, 0.02, 0.03, and 0.05 M_⊙, and the mass infalling rate Ṁ_env of 1 × 10^-6 to 1 × 10^-5   M_⊙ yr^-1.

We examined 2880 parameter sets. Determining the bestfit-model is an important process. One can apply an objective method such as χ² test. This kind of analysis is preferred, for example, for fitting spectral data, such as the 10 μm silicate feature. However, when different types of data are simultaneously taken into account, as in our model, the situation is different. In general, the best SED-fit model does not guarantee fitting the best in the intensity images or the polarization data, and vise versa. Even if we attempt to apply χ² test, it is not clear how to determine the weights for SED, intensity distribution, and polarization data without a subjectivity. In addition, we did not intend to reproduce all aspects of the observations. Thus, we fit by eye and chose a good solution. We believe that our selected model can be somewhat different the best model without changing our conclusion. The parameters of the selected model have a disk mass of 0.02 M_⊙, a disk scale height of 0.4, a mass infalling rate of 5 × 10^-6   M_⊙ yr^-1, an a_max,1 of 1000.0 μm, a ρ_th,1 of 1 × 10^-20 g cm^-3, an a_max,2 of 10.0 μm, and a ρ_th,2 of 1 × 10^-21 g cm^-3. The optical depth of the disk is 4.2 at V, which is calculated by integrating radially outward through the entire disk R_disk in the midplane (z = 0). Table 1 summarizes all parameters and the uncertainty of the corresponding parameters. In later discussions, we justify and present the uncertainties of important disk properties: a large radius of 3000 AU, an intermediate disk height of H₀ = 0.4, extremely large grains of a_max = 1000.0   μm, and an intermediate optical depth of τ = 4.2 in the disk midplane. In the next section, we discuss them in detail and briefly describe how the uncertainties of some parameters listed in Table 1 are derived.

Figure 1a shows the absorptive opacity κ, the scattering opacity σ, and the total opacity of the selected dust models as functions of the wavelength. In the dust models of a_max = 0.5   μm and 10.0 μm, the opacities suddenly change in the NIR and MIR, respectively, where the size parameter x = 2πa/λ becomes  ~1. On the other hand, a nearly flat spectrum is seen in the a_max = 1000.0   μm model. Figure 1c presents the central region of the mass density distribution of the selected model. The two contour lines denote the density thresholds of ρ_th,1 and ρ_th,2. Because of our assuming mass density domains, millimeter-sized grains exist even up to  ~100 AU above the midplane. According to theories of grain growth and dust settling (e.g. Miyake & Nakagawa 1993; Dullemond & Dominik 2004b; Schräpler & Henning 2004; Dullemond & Dominik 2005; Tanaka et al. 2005), such extremely large grains could sediment closer to the midplane. In a future work, one can examine the distribution of large grains in disks with a large radius and an intermediate optical depth. The dust temperature is  ≳ 60 K in the whole disk and has a nearly cylindrical distribution (see Fig. 1d).

Because of an intermediate optical depth of 4.2 at V-band, the entire disk is moderately heated. Figure 1b compares the model SED with the observations. The overall shape of the model SED agrees with the observations. The scattered light dominates in the wavelengths up to  ~ 2 μm. A smooth emission bump at 10 μm, but not a sharp triangle shape, in the observation is quantitatively reproduced in our model, which is also due to large grains in the disk (e.g. van Boekel et al. 2003, 2005). In submillimeter to millimeter wavelength ranges, we see higher fluxes by a factor up to 2–3 in our results, in which a large fraction of the flux comes from the outer envelope. The results fit better if a 20″ aperture is applied, as mentioned in Sect. 2.1. The uncertainty in the mass infalling rate of the envelope listed in Table 1 is estimated from the above flux error. Figure 2 shows the J, H, and K_S band intensity and polarization images of our observations (upper panels) and the selected model results (lower panels). Except for the lobe shape, a nearly spherical nebulosity close to the central star is reproduced in the model images. Table 2 shows the result of aperture polarimetry.

Because an arc-like feature extending towards the northwest is seen in the polarization disk and this probably indicates the scattered light from the disk surface, we used a small aperture of 0$\mathrm{″̣}$1 in circular diameter at the central star to reduce a dilution by other polarized components than from the polarization disk. The polarization values in the selected model are 1.7%, 1.6%, and 1.6% in the JHK bands, respectively, and are in good agreement with observations, excluding the J band. We believe that this discrepancy comes from contamination by highly polarized scattered flux from the envelope. The butterfly-shaped polarization disk and a weak wavelength dependency of the extension are also reproduced in our model results. In the next section, we justify our estimates of the grain size and the disk geometry and show the complementarity of these analyses.

Fig. 1 Results of the selected model. a) The dust opacities as functions of the wavelength. The red, blue, and green curves denote the dust models with amax values of 0.5 μm, 10.0 μm, and 1000.0 μm, respectively. The dotted, dashed, and solid lines denote the scattering opacity σ, the absorptive opacity κ, and the total opacity κ + σ, respectively. b) Comparisons of the model SEDs with the previous observations. The solid, dashed, dashed-dotted, and dotted curves denote the total flux, the attenuated stellar flux, the thermal emission, and the scattered flux. The red, green, orange, light blue, violet, and yellow dots are data from Hillenbrand et al. (1992), the MSX point source catalog (PSC), the IRAS LRS spectra of short and long wavelengths, IRAS PSC, and Mannings (1994). c, d) The mass density and dust temperature distribution maps in the inner 5000 × 5000 AU region. The two contour lines indicate the threshold densities of 1 × 10-20 g cm-3 and 1 × 10-21 g cm-3 to change the grain sizes. From the inner, higher density region, the grain models with amax = 1000.0   μm, 10.0 μm, and 0.5 μm are applied.

Fig. 2 The intensity and polarization images in the J, H, and KS bands. The polarization vector lines are overlaid on the polarization images. The upper two rows show our previous observational results (Murakawa et al. 2008a). The JH-band data and the KS-band data are the results of the Subaru / CIAO and VLT / NACO observations, respectively. The bottom two rows are our selected model results. In these results, Gaussian functions that have the corresponding PSF sizes are convolved.

3. Discussion

3.1. Grain size

Evidence of large grains in R Mon’s disk was already discovered from the detection of a low opacity spectral index β (Fuente et al. 2003). Thus, in this section, we investigate the consistency with other methods such as SED and polarization analyses.

We evaluate two models with different grain sizes in the disk: (1) the grains with a_max = 10.0   μm for ρ_th ≥ 1 × 10^-21 g cm^-3, i.e. a_max = 1000.0   μm grains are replaced by the a_max = 10.0   μm grains in the selected model; and (2) only a_max = 0.5   μm grains for the entire disk, hereafter models G1 and G2, respectively. Figures 3 and 4 show the SED and the model images in the K_S band. The result of aperture polarimetry is also presented in Table 2. In the optical wavelength range, the flux is too low because smaller grains result in higher extinction in this wavelength range. The same reason stands in the NIR, although the difference in flux is much less. A slightly better fit in the G1 model is found. The real R Mon disk possibly has a slightly higher optical depth than the selected model. It is worth considering whether an optical depth of the disk in the NIR is kept constant for different dust models. Such a model SED probably satisfies the observations, as intended. However, we should also take the disk mass into account (some relative discussions, Chiang & Goldreich 1997; D’Alessio et al. 2001; Wolf et al. 2003; Murakawa et al. 2010). Because the optical depth τ is approximately proportional to ⁽σ + κ⁾M_disk, changing the dust size has a strong impact on the resulting disk mass. With the selected dust model of a_max = 1000.0   μm, a disk mass of 0.02 M_⊙ is required, which roughly agrees with an estimated value of  ~0.01 M_⊙ from an observation of gas phase CO emission lines (F06). In contrast, the G1 and G2 models have  ~100 and  ~10 times higher dust opacities in the optical, causing  ~100 and  ~10 times lower disk masses. Thus, small grain models are inadequate.

Fig. 3 Comparison of model SEDs with different grain sizes in the disk. The black, red, and blue lines denote the results of the selected model, model G1 (amax = 10.0   μm grains for ρ ≥ 1 × 10-21 g cm-3), and model G2 (amax = 0.5   μm grains for the entire disk and envelope), respectively. All other parameters are the same as the selected model.

Detection of a flux excess in the submillimeter and millimeter wavelength ranges is regarded as evidence of large grains (e.g. Beckwith & Sargent 1991), because large grains have a low opacity spectral index β (Draine 2006, and references therein). It is known that the β is 2 for the interstellar population and  <1 for extremely large grains. In our dust models, the β values for the a_max = 1000.0   μm and a_max = 10000.0   μm are 1.3 and 1.0, respectively, at a wavelength of 1.3 mm, which are reasonable values compared to similar ones calculated by Draine (2006, see Fig. 3). In contrast, F06 find a low value of 0.3. It is very likely that a compact feature detected at 1.3 mm towards the central star (~150 AU) is the only region to contribute the millimeter flux. Therefore, the β should be characterized by large dust grains existing in this region. However, the dust opacity with this low β is significantly lower than the a_max = 1000.0   μm model. As a result, a much higher disk mass (≫0.01 M_⊙) is required to fit the optical to NIR fluxes, but this is unrealistic. This low β issue still remains puzzling in R Mon. The most plausible reason for this discrepancy lies in the different assumptions. For example, F06 assume a disk temperature of 215 K and a mean opacity of  cm² g^-1, but these are not in our modeling. In addition, when the β is estimated from the submillimeter / millimeter flux slope α, a compensation term Δ should be taken into account, i.e. α ≈ 2 + β / (1 + Δ), as discussed in Sect. 4.1 in a paper by Beckwith & Sargent (1991). It is obvious that a more self-consistent model is preferable for more qualitative study.

The grain-size effect also appears in the intensity and polarization images. The nebulosities are more asymmetric in both G1 and G2 models compared to the selected model because of the higher opacities in the smaller grains. The grain size effect is more prominent in the polarization data. In the G1 model, the appearance of the polarization disk is acceptable, but the polarizations of  ~10% in the JHK-bands are too high. In the G2 model, a strong X-shaped feature with high polarization is seen, and the result of an aperture polarimetry shows low (14% at J) to high (34% at K) polarizations, which is the opposite wavelength dependence to typical T Tauri stars. Strictly speaking, the grain shape (i.e. asphericity and fluffiness) and the alignment efficiency should be taken into account for the aperture polarimetry. The dichroism (or deferential opacity ΔQ ≡ Q_// − Q_⊥) appears when the dust size parameter x = 2πa / λ is around unity, i.e. a submicron size in the NIR. However, a spherical grain approximation would be enough if large grains (x ≫ 1) are expected. The polarization analysis also rules out small grain models.

3.2. Disk geometry

We examined three models with (1) an F1 disk form of Eq. (1), (2) a small disk radius of 100 AU (= 0$\mathrm{″̣}$125 at 800 pc), and (3) a high disk height (H₀ = 1.0), hereafter the D1, D2, and D3 models. For the D1 model, we set an H₀ of 0.7, a disk mass of 0.01   M_⊙, and threshold mass densities of ρ_th,1 = 3 × 10^-21 g cm^-3 and ρ_th,2 = 3 × 10^-22 g cm^-3 to fit the SED. For the D2 model, we set a disk scale height of 0.1 and a disk mass of 0.01 M_⊙. We applied only the a_max = 1000.0   μm dust model in the entire disk. This disk geometry is typical of T Tauri and Herbig Ae stars and is similar to the ones estimated from the observations at 1.3 mm continuum (F06) and a previous radiative transfer model (Alonso-Albi et al. 2009). For the D3 model, the disk height is from previous CO emission line observations and modeling (F06). The other parameters are the same as the selected model. Figures 5 and 6 show the model results of the SED and the images in the K_S band.

Fig. 4 The model results with different grain sizes in the disk. The left and right columns correspond to the models G1 and G2, respectively. The wavelength is 2.20 μm.

The D1 model result is actually fairly good. However, the nebulosity around the central star is slightly more asymmetric and the MIR flux too low compared to the selected model and the observations. Such a slight asymmetric appearance of the nebulosity is caused by a higher optical depth due to a higher α value of 15 / 8 compared to the selected model (1.5). A similar discussion is found in Fischer et al. (1996). We believe that a lower α value (e.g. 1.5) is better for the R Mon disk. For the other aspects, such as whether a flaring (β > 1) or flat (β = 1) structure or whether the inner rim is puffed up or not, we do not find any clear clue in our modeling.

In the D2 model, the intensity image has a cup-like appearance, which is often seen in T Tauri stars, but not in the case of R Mon. Such a feature with a dark lane is produced by a high optical depth in the disk midplane. In this case, the matter of a 0.01 M_⊙ is condensed in a 100 AU radius. This high optical depth model provides another problem for fluxes that are too low in the optical and NIR. To solve this problem, one can decrease the disk mass. However, as already discussed in the last section, this contradicts an estimation of the disk mass with gas phase data. A model by Alonso-Albi et al. (2009) applies a disk radius and a dust disk mass of 150 AU and 1.4 × 10^-4   M_⊙, which corresponds to a total mass of 0.014 M_⊙ if a typical gas-to-dust mass ratio of 100 is assumed. Their model SED fits the observation closely. Unfortunately, the intensity and polarization images do not satisfy the observation because their model has a dust-free region between the outer disk boundary at 150 AU and the innermost part of the outer envelope at 700 AU. In the intensity and polarization images, a feature of this gap should appear. We interpret that their modeled disk is the inner part of the disk, where large grains exist, instead of the entire disk. The polarization image provides another important clue to distinguish it from compact disks. With such a disk, the polarization disk looks like the hands of a clock at 8:20 rather than a disk, and a strong X-shaped feature is produced, which is not seen in R Mon.

The D3 model fits the SED and the intensity image well. However, the appearance of the polarization disk looks spherical instead of butterfly-shaped as seen in the observations. Thus, the dust disk is expected to be geometrically thinner than the gas disk. Also considering that R Mon’s disk has a somewhat thick geometry, it is likely that the dust layer is concentrated towards the midplane more than the gas phase molecules.

In this section, we have presented comparison models with different disk geometries. We find that even though the resulted model SEDs satisfy the observations, the intensity and polarization images do not guarantee a fit to the observations. The opposite should also be true, although we did not show it. We expect that R Mon’s disk has a large radius and possibly a thinner geometry than the gas phase disk; i.e. the dust settling is ongoing.

Fig. 5 Comparison of model SEDs with different disk geometry parameters. The black, red, blue, and green lines denote the results of the selected model, the F1 disk model (D1), the Rdisk = 100 AU model (D2), and the H0 = 1.0 model (D3), respectively.

3.3. The polarization and dust disks of R Mon

From these results, we find that R Mon’s disk is unique and contrasts with typical T Tauri and Herbig Ae disks. The polarization disk in our high-resolution data traces the real disk closely and the polarization vector alignment is weak. The point is that the optical depth of R Mon’s disk (τ = 4.2) is not expected to be as high as those of typical T Tauri disks (τ ~ 10⁵–10⁶, e.g. Whitney et al. 2004). Such large differences in the optical depth come from large differences in the disk volumes but comparable disk masses. R Mon’s disk has an outer radius of 3000 AU and a disk height H₀ of 0.4, while those of typical T Tauri stars are  ~100 AU and 0.01–0.1, respectively; that is R Mon’s disk is 1000 to 10 000 times bigger. In contrast, the disk masses are comparable (~0.01   M_⊙). Because of this difference in optical depth of the disk, an alternative interpretation of the polarization disk other than the BM mechanism, should be considered.

In the BM mechanism, the direct star light and single-scattered light vanish, but the double-scattered light can reach the observer and causes a vector alignment. Because photons behave as if they avoid the optically thick midplane, we call this phenomenon the roundabout effect (M10). In addition, the second scattering event can happen outside the disk, and the apparent size of the polarization disk can be larger than the real disk (M10). The scenario is different in R Mon’s disk and is explained with our new model (M10). The roundabout effect is weak, and the single-scattered light still reaches the observer. In this condition, the polarization status somewhat shows the single-scattering properties of dust in the disk, such as the grain size and the disk height, which determines the scattering angle. It is worth mentioning that, although the scattering properties do not change for large grains (x ≫ 1) (see Fig. 1 in our previous paper M10), the result of aperture polarimetry is different between our selected model and the G1 model. In both cases, the photon penetrate the disk. However, the roundabout effect is stronger in the G1 model because the optical depth of the disk is higher for a higher dust opacity model. The flux towards the central star contains more highly polarized, scattered light from the disk surface, where small grains exist than does our selected model. As a result, a higher polarization is obtained in the G1 model. Furthermore, the projection of the real disk appears as the polarization disk, and the vector pattern looks more centro-symmetric than linearly aligned.

Fig. 6 Model results with different disk geometries. From the left, the results in the KS band of the F1 disk form (D1), a small disk radius of 100 AU (=0$\mathrm{″̣}$125) (D2), and an H0 of 1.0 (D3) are presented. The disk scale height of the D3 model is the same as the one estimated from the gas phase disk (F06).

Here, we would like to recall the definition of a polarization disk. BM argue that the polarization disk is a region where the vectors are aligned along the equatorial plane. The polarization in this region is lower than in the outer lobe or envelope in general. At the same time, another low polarization region with nearly centro-symmetric vector alignment is often detected around the polarization disk; for example, a feature of the hands of a clock is seen in the D2 model. These features are sometimes difficult to distinguish in the observations. However, we have pointed out that the polarization disk contains important information on the grain size and / or grain shape. In contrast, the other is simply reproduced by depolarization due to multiple scattering and only indicates the existence of a high-density region near the central star (M10). Thus, we call the latter region a polarization shadow (or polarization dark lane) and distinguish it from the polarization disk.

We stress that the R Mon is a special case. According to this definition, the butterfly-shaped low-polarization feature with  ~4″ extension towards the central star appears to be a polarization shadow because the vectors are not horizontally aligned. Nevertheless, we still regard this feature as the polarization disk because it turns out from our modeling that this feature traces well the projected angular size of the real dust disk. A low NIR polarization and not detecting the vector alignment indicate the presence of large grains and the intermediate NIR optical depth of the disk.

The disk size is also different from typical T Tauri disks and more similar to ones in some massive protostars such as IRAS 20126+4104 (1000 AU, Sridharan et al. 2005), AFGL 490 (in a transition stage to Herbig Be, 1500 AU, Schreyer et al. 2006), M 17 (3000 AU, Steinacker et al. 2006), and CRL 2136 (2000 AU, Murakawa et al. 2008b). Theoretical calculations predict that a thousand AU size can form round a 10 M_⊙ star (e.g. Yorke & Bodenheimer 1999). On the other hand, the expected large dust grains in R Mon’s disk agrees with T Tauri disks. For massive protostars, estimation of the grain sizes is less precise or is done with comparably fewer samples, because it is difficult from only an SED fit and from a detection of millimeter flux excesses. However, in some cases, high NIR polarizations are detected (e.g. Whittet et al. 2008; Murakawa et al. 2008b). In particular, for an object where the central feature is visible, such as CRL 2136, a high NIR polarization can be reproduced with only submicron-sized grains in the disks despite the grain shape. Therefore, in terms of the disk geometry and grain sizes, R Mon’s disk shows an intermediate property between T Tauri stars and massive protostars. Furthermore, it is also interesting that ongoing dust settling is expected in the R Mon disk, as mentioned in the previous section. Some authors have investigated the disk evolution, both dust settling and dust growth (e.g. Miyake & Nakagawa 1993; Dullemond & Dominik 2004b; Schräpler & Henning 2004; Dullemond & Dominik 2005; Tanaka et al. 2005). The grain growth time t_grow, in which a 0.1 μm dust grain grows to micron size, is expected to be shorter than the formation time of the dust layer t_settle, (e.g. t_grow ~ 0.1t_settle, Tanaka et al. 2005). Although these studies are based on T Tauri disks, a similar conclusion is still applicable in R Mon’s disk.

4. Conclusion

We have modeled the R Mon dust disk by means of radiative transfer calculations. In our approach, we combined the polarization disk analysis with an SED fit. In addition, we used previous results from CO emission line observations in the radio frequencies and modeling to constrain the dust size and to compare them with the dust disk geometry. The estimated disk parameters are a disk of 3000 AU, a disk height of H₀ = 0.4, an optical depth in the midplane of 4.2 at V, a disk mass of 0.02 M_⊙, and a grain size of a_max = 1000.0   μm. We justified these parameters with some alternative models and showed how different analytical methods constrain them. For the grain size, an SED fit, a detection of the millimeter flux excess, and a polarization analysis suggest that large grains should exist in the disk. Taking practical factors into account such as observation precisions, the millimeter flux excess, and polarization analyses, can provide a quantitative evaluation, and the SED fit combined with the dust and gas phase models allows a qualitative estimation in any size range. For the disk geometry, even if it is possible to satisfy the SED with a small radius of  ~ 100–150 AU, the NIR intensity and polarimetric images do not permit this and require a large radius. We believe that such a small radius, which is detected at a radio continuum of 2.7 mm (F06), is the inner part of the disk where large grains exist.

We attempted to apply our generalized scattering model to interpreting R Mon’s polarization disk. Because of a large outer disk radius and an intermediate disk height, it turned out that the polarization disk roughly traces the projected real dust disk. Not detecting a vector alignment and a nearly constant, low NIR polarizations (P ~ 1–2%) are explained similarly, where the roundabout effect is weak. From these results, we found that R Mon’s disk is unique and contrasts typical T Tauri and Herbig Ae disks. The large disk radius is fairly similar to those of some massive (proto)stars. In R Mon’s disk, the dust grains grow in size, and the dust condensation towards the midplane is probably currently in progress because an intermediate disk height is expected. This agrees with a theoretical prediction that the dust settling time is longer than the dust growth time.

We have demonstrated that our complementary combination method (polarization analysis, SED fit including the dust and gas phase modeling) is a powerful technique for deriving important properties of YSO disks. A multiple analytical methods based on different physics will be important in the coming ELT, HERSCHEL, JWST, and ALMA era.

Acknowledgments

Our NIR polarimetric images were obtained at the European Southern Observatory (proposal ID: 078.C-0554) and the Subaru telescope. We would like to thank an anonymous referee for his / her meaningful suggestions and comments.

References

All Tables

All Figures

Fig. 1 Results of the selected model. a) The dust opacities as functions of the wavelength. The red, blue, and green curves denote the dust models with amax values of 0.5 μm, 10.0 μm, and 1000.0 μm, respectively. The dotted, dashed, and solid lines denote the scattering opacity σ, the absorptive opacity κ, and the total opacity κ + σ, respectively. b) Comparisons of the model SEDs with the previous observations. The solid, dashed, dashed-dotted, and dotted curves denote the total flux, the attenuated stellar flux, the thermal emission, and the scattered flux. The red, green, orange, light blue, violet, and yellow dots are data from Hillenbrand et al. (1992), the MSX point source catalog (PSC), the IRAS LRS spectra of short and long wavelengths, IRAS PSC, and Mannings (1994). c, d) The mass density and dust temperature distribution maps in the inner 5000 × 5000 AU region. The two contour lines indicate the threshold densities of 1 × 10-20 g cm-3 and 1 × 10-21 g cm-3 to change the grain sizes. From the inner, higher density region, the grain models with amax = 1000.0   μm, 10.0 μm, and 0.5 μm are applied.

In the text

Fig. 2 The intensity and polarization images in the J, H, and KS bands. The polarization vector lines are overlaid on the polarization images. The upper two rows show our previous observational results (Murakawa et al. 2008a). The JH-band data and the KS-band data are the results of the Subaru / CIAO and VLT / NACO observations, respectively. The bottom two rows are our selected model results. In these results, Gaussian functions that have the corresponding PSF sizes are convolved.

In the text

	Fig. 3 Comparison of model SEDs with different grain sizes in the disk. The black, red, and blue lines denote the results of the selected model, model G1 (amax = 10.0   μm grains for ρ ≥ 1 × 10-21 g cm-3), and model G2 (amax = 0.5   μm grains for the entire disk and envelope), respectively. All other parameters are the same as the selected model.
In the text

	Fig. 4 The model results with different grain sizes in the disk. The left and right columns correspond to the models G1 and G2, respectively. The wavelength is 2.20 μm.
In the text

	Fig. 5 Comparison of model SEDs with different disk geometry parameters. The black, red, blue, and green lines denote the results of the selected model, the F1 disk model (D1), the Rdisk = 100 AU model (D2), and the H0 = 1.0 model (D3), respectively.
In the text

	Fig. 6 Model results with different disk geometries. From the left, the results in the KS band of the F1 disk form (D1), a small disk radius of 100 AU (=0$\mathrm{″̣}$125) (D2), and an H0 of 1.0 (D3) are presented. The disk scale height of the D3 model is the same as the one estimated from the gas phase disk (F06).
In the text