HIGH-CONTRAST NEAR-INFRARED IMAGING POLARIMETRY OF THE PROTOPLANETARY DISK AROUND RY TAU^*

The modeled system consists of a pre-main-sequence star with an axisymmetric circumstellar disk, an infalling flattened envelope, and an outflow cavity. The SED is determined using full radiative transfer in the disk and envelope, i.e., including absorption, scattering and re-radiation of light by dust grains in addition to the flux directly observed from the star. The heating sources of the disks are stellar radiation and viscosity in the disk. The stellar spectrum is parameterized by the star's effective temperature and radius, via comparisons with modeled spectra for stellar photospheres (Kurucz 1993; Brott & Hauschildt 2005).

The density distribution of a standard flared accretion disk (e.g., Shakura & Sunyaev 1973; Lynden-Bell & Pringle 1974) is described in cylindrical coordinate (r,z) by

$$\begin{equation} \rho (r,z) = \rho _0 \left[1-\sqrt{\frac{R_*}{r}} \right] \left(\frac{R_*}{r} \right)^\alpha \exp \left\lbrace - \frac{1}{2} \left[\frac{z}{h} \right]^2 \right\rbrace , \end{equation} \tag{ 1 }$$

where ρ₀ is a constant to scale the density, R_* is the stellar radius, α is the radial density exponent, and h is the disk scale height. The scale height h increases with radius as h = h₀r^β, where β is the flaring power (β > 0). According to more detailed models of disk structures by Chiang & Goldreich (1997) and D'Alessio et al. (1999), the scale height exponent varies with radius in the inner disk but follows a similar power law beyond a few AU (Cotera et al. 2001). In addition to the above parameters the minimum and maximum radii of the disk are included as free parameters. The gas-to-dust mass ratio is assumed to be 100.

Throughout the simulations, Robitaille et al. (2006, 2007) assume α = β + 1. This yields the surface density distribution Σ(r)∝r⁻¹, approximately agreeing with that inferred from millimeter interferometry for disks associated with many low-mass pre-main-sequence stars (Andrews et al. 2009, 2010). This surface density power index may not be consistent with the one for the RY Tau disk seen at millimeter wavelengths by Isella et al. (2010). Even so, models with the above assumption fit the observed SEDs well, as shown in Section 4.2.

The density structure for the envelope is given by Ulrich (1976) and Terebey et al. (1984). The envelope is associated with a cavity whose shape varies as z∝r^1.5. Robitaille et al. (2006) contain further details for the envelope and outflow cavity. Although the envelope mass, cavity density, etc., will be calculated using the SED fitter, their contribution to the SEDs of optically visible pre-main-sequence stars such as RY Tau will be negligible (Section 4.2).

The dust composition assumes a mixture of astronomical silicates and graphite in solar abundance without an ice coating. The following two grain size distributions are used based on Wood et al. (2002) and Whitney et al. (2003a, 2003b): (1) the distribution in the denser regions of the disk ($m_{{\rm H}_2} > 10^{10}$ cm⁻³) where one would expect significant grain growth up to 1 mm; and (2) the distribution in the more diffuse regions ($m_{{\rm H}_2} < 10^{10}$ cm⁻³) in which the grain sizes are slightly larger than those in the diffuse interstellar matter. The grain size distribution for the former is described in Wood et al. (2002), and that for the latter is similar to Kim et al. (1994) (hereafter KMH; Whitney et al. 2003a). The former affects the thermal balance in the disk via radiative transfer, and is successful in fitting the SED of the HH 30 disk (Wood et al. 2002; Whitney et al. 2003a, 2003b). The latter grains are responsible for scattered light on the disk surface (and in the envelope) at optical and near-infrared wavelengths. The optical constants for silicate and graphite are taken from Laor & Draine (1993).

The 200,000 SEDs provided by the Robitaille et al. (2006) code include flux and polarization spectra for 250 wavelengths (from 0.01 to 5000 μm). The authors used 20,000 parameter sets and computed the results for 10 viewing angles from face-on to edge-on at equal intervals in the cosine of the inclination. The SED fitter developed by Robitaille et al. (2007) searches for best-fitting SEDs using the minimum χ² method to fit the observed fluxes at a range of wavelengths. We set a distance to the object of 140 pc, and an acceptable range for the visual extinction A_V of 2.0–2.4, based on Calvet et al. (2004). See Robitaille et al. (2007) for further details of the fitting process.

To derive the star/disk/envelope parameters for RY Tau we used their SED fitter using the photometric data tabulated in Robitaille et al. (2007) (Table 1). These parameters are used to model the PI images with the Robitaille et al. (2006) radiative transfer code. We used 10⁷ photons for each case. The resultant images for the Stokes parameters Q and U were convolved with the PSF of the reference star (Section 2) before obtaining the PI image using ${\rm PI}=\sqrt{Q^2+U^2}$. Then the PI flux is normalized to the stellar I flux I_*, and scaled to match the pixel size of Subaru/HiCIAO.

Figure 2 shows the observed and 10 best-fit SEDs. We note that the fluxes at different wavelengths were measured at different epochs (Table 1), and may be highly time variable. In particular, Herbst et al. (1994) reported a variation in the optical fluxes of a factor of ∼8 from 1961 to 1980 (see also Petrov et al. 1999 for longer monitoring observations). Even so, the fluxes in Figure 2 are well fitted by a single SED except for the 12 μm flux observed using the Infrared Astronomical Satellite. A larger excess at this wavelength may be attributed to bright silicate emission (Honda et al. 2006) in the filter coverage (8.5–15 μm).

Zoom In Zoom Out Reset image size

The parameters for the modeled SEDs are shown in Table 2 along with the χ² value and the model ID specified by Robitaille et al. (2006). The constant used to scale the scale height relation (h₀) is fixed at r = 50 AU from the star. The modeled SEDs are based on four sets of physical parameters for the star/disk/envelope with inclination angles of 57°–76°. The inclination angle for most of the SEDs and the disk mass for all the models are approximately consistent with those of Isella et al. (2010) obtained using millimeter interferometry (≳ 65° and ≳ 3 × 10⁻³ M_☉, respectively).

Table 2. Best-fit Models Using the SED Fitter

Model ID	Inclination	χ2		Stellar Parameters			Disk Parameters							Others
				R*	T*a		Rmin	Rmax	h50 AU	β	Mdisk	$\dot{M}_{\rm disk}$		Menvelope	ρcavity
	(°)			(R☉)	(K)		(Rsub)	(AU)			(M☉)	(M☉ yr−1)		(M☉)	(g cm−3)
3000949	75.5	1.3 × 102		4.9	4.7 × 103		1.0	83	2.4	1.13	1.1 × 10−1	1.2 × 10−6		1.1 × 10−4	1.7 × 10−21
	69.5	1.6 × 102
	63.2	2.2 × 102
3012376	69.5	1.6 × 102		5.3	5.1 × 103		7.8	84	2.7	1.16	8.0 × 10−2	6.1 × 10−8		3.4 × 10−4	2.0 × 10−21
	63.3	1.8 × 102
	75.5	2.0 × 102
	56.6	2.3 × 102
3007615	69.5	2.1 × 102		4.6	6.4 × 103		7.0	4.0 × 102	1.6	1.12	4.6 × 10−2	1.1 × 10−7		2.0 × 10−6	0
	75.5	2.3 × 102
3012808	56.6	2.3 × 102		5.0	5.1 × 103		4.7	1.2 × 102	2.1	1.17	2.8 × 10−2	2.4 × 10−7		3.5 × 10−8	0

Notes. R_*, stellar radius; T_*, stellar temperature; R_min, minimum disk radius; R_max, maximum disk radius; M_disk, disk mass; $\dot{M}_{\rm disk}$, disk mass accretion rate; M_envelope, envelope mass; ρ_cavity, density in the outflow cavity; R_sub, sublimation radius. The sublimation radius is empirically determined by R_sub = R_*(T_sub/T_*)^−2.1 (Robitaille et al. 2006). The inclination is measured from the polar axis. See the text for definition of χ². The masses and mass accretion rates are described for gas+dust assuming the gas-to-dust mass ratio of 100. ^aAny temperature listed here does not agree the measured spectral type of RY Tau of F8-G2 (Petrov et al. 1999; Calvet et al. 2004; Mendigutía et al. 2011), corresponding to an effective temperature of 5.8 − 6.3 × 10³ K. The different stellar temperatures listed here do not affect the modeled SEDs and PI image at 1.65 μm (Section 4).

Download table as:  ASCIITypeset image

Despite the similarity in the shape of SEDs, some parameters are significantly different between models. The outer disk radius is ∼80 AU for two out of four models, but ∼120 and ∼400 AU for the remaining. The inner disk radius calculated with regard to the sublimation radius varies by a factor of ∼8 between models; the disk mass by a factor of ∼4; disk accretion rate by a factor of ∼20; and the envelope mass by a factor of ∼10⁴. The gas+dust density in the outflow cavity is ∼2 × 10⁻²¹ g cm⁻³ in two out of the four models, but zero for the remaining.

Figure 3 shows the simulated PI images for the four physical parameter sets, with different inclination angles. Despite significant differences in the parameters described above, the four parameter sets result in PI distributions that are strikingly similar to each other, except for the fact that model 3007615 produces a slightly fainter distribution than the others, due to the small flaring of the disk (h_50 AU = 1.6 AU, β = 1.12). This similarity between models, both for the SEDs and PI distribution, can be attributed to the fact that both the infrared radiation and the scattered light at the disk surface are determined by the surface geometry of the disk. The PI distribution shown in Figure 3 is dominated by scattering on the disk surface, and the contribution from the envelope and outflow cavity is significantly smaller. To investigate this, we removed the envelope from the models with a relatively massive envelope and diffuse dust layer in the cavity (models 3000949 and 3012376 according to the SED fitter), and ran the simulations again. We find that the PI distribution is almost identical to the model with the envelope and cavity, in particular within ∼100 AU of the star.

Zoom In Zoom Out Reset image size

The models at inclination angles i = 57° and 76° show two separated bright regions similar to the observations. Furthermore, the modeled PI flux normalized to the stellar I flux (PI/I_*) is of the same order as the HiCIAO observations. However, the offset from the major axis of the disk is significantly smaller for all of the modeled PI images compared to the observations. These offsets for the models and observations are clearly shown in the bottom right plots of Figure 3. While that of the observed image is clearly offset from the major axis, the modeled images are more symmetric about the major axis of the disk. This indicates that the disks inferred from the SEDs are geometrically thinner than that producing the observations,²⁷ if the shape of the emission is attributed to the disk.

We follow the method described in Fischer et al. (1994). We place a central unpolarized light source equivalent to the star as the starting point for calculating the scattering of photons from the disk. The light path for the next scattering position is calculated for an opacity distribution based on the disk in Equation (1) and the dust opacity described below. The scattering angle and Stokes parameters after scattering are calculated based on Mie theory. The Stokes parameters for each photon are initially set to (I₀,Q₀,U₀,V₀) = (1,0,0,0) and normalized to $I_{\rm out} = {\rm albedo} \cdot I_{\rm in}$ after each scattering.

We use 10⁶ photons for each simulation. The photons escaping from the disk are collected in imaginary detectors at different viewing angles. For our purpose unscattered stellar photons are not collected. In order to normalize the PI flux to the stellar I flux, we separately calculate the expected number of photons for each viewing angle based on the incident number of photons and extinction.

For dust grains we use homogeneous spherical particles as commonly used in other studies, including Whitney & Hartmann (1992), Whitney et al. (2003a, 2003b), Cotera et al. (2001), Wood et al. (2002), and Robitaille et al. (2006). We use the grain compositions of Robitaille et al. (2006), Robitaille et al. (2007), Cotera et al. (2001), and Wood et al. (2002), i.e., astronomical silicate and carbon dust without an ice coating. For the size distribution, we use (1) that of interstellar dust measured by KMH (R_V = 3.1); (2) the larger size distributions used by Cotera et al. (2001) and Wood et al. (2002) to reproduce the scattered light observed in the HH 30 disk (hereafter C01); and (3) same as C01 but with the grain size scaled by a factor of 15, preserving the total dust mass (hereafter C01 × 15) (Figure 4). Note that the C01 × 15 distribution would be too large for the disk surface of RY Tau. Honda et al. (2006) conducted mid-infrared spectroscopy and analysis of silicate emission toward a number of pre-main-sequence stars. These authors fit the spectrum of individual stars using models with two dust sizes (0.1 and 1.5 μm), and showed that the disk surface of RY Tau has one of the smallest ratios of large-size grains to small-size grains of their sample (34%). This suggests that the grain size distribution in RY Tau is fairly close to that of the HH 30 disk, which has a central star that is significantly younger than most optically visible pre-main-sequence stars (e.g., Burrows et al. 1996; Watson & Stapelfeldt 2007). However, we still show the results obtained with C01 × 15 to investigate how the results are affected by the assumed size distribution.

Zoom In Zoom Out Reset image size

Different authors use different types of carbon dust, either graphite or amorphous carbon (Cotera et al. 2001; Wood et al. 2002). While graphite has been extensively used (e.g., Draine & Lee 1984; Laor & Draine 1993; Kim et al. 1994; Whitney et al. 2003b; Robitaille et al. 2006; Dong et al. 2012a, 2012b), far-infrared SEDs of young stellar objects and evolved stars suggest the absence of graphite and presence of amorphous carbon in circumstellar dust (Jäger et al. 1998 and references therein). We still use graphite for the KHM distribution for consistency with the authors as their size distribution is determined assuming graphite for the carbon dust. We use amorphous carbon for C01 and C01 × 15 for the above reason and following Cotera et al. (2001) and Wood et al. (2002). The use of different carbon dusts does not significantly affect the modeled PI distribution, and it does not affect the conclusions of the paper (Appendix A).

Calculations for Mie scattering are made using the code developed by Wiscombe (1996). The optical constants for astronomical silicate and graphite are obtained from Draine & Lee (1984): they are identical in the near-infrared with the Laor & Draine (1993) values used for the simulations in Section 4. For amorphous carbon we use the optical constants of Jäger et al. (1998) with a pyrolysis temperature of 600°C. As for the use of different carbon materials (amorphous carbon or graphite), amorphous carbon with different pyrolysis temperatures does not significantly affect the modeled PI distribution and conclusions of the paper (Appendix A).

Table 3 shows the physical properties (opacity, albedo, forward throwing parameter, and the maximum degree of polarization) for the above three dust models. Figure 5 shows the scattering properties of these dust models at 1.65 μm, i.e., I, PI, and the degree of polarization as a function of scattering angle, in the case where the incident light is not polarized. The I and PI fluxes for each scattering angle are derived by normalizing the scattering matrix elements S₁₁ and −S₁₂ by a constant I₀ to match the results to the "weighted photon method" described above (i.e., ∫I/I₀dΩ = albedo). The forward scattering is more significant for large grain models, but different dust models show a relatively similar distribution of scattered intensity I/I₀ at scattering angles of 40°–180°. The polarization shows a similar dependency on scattering angle, peaking at ∼90°, but the absolute value for the C01 and C01 × 15 models are 30%–40% lower than for the KMH. The PI flux shows a maximum at scattering angles of θ = 75°, 60°, and 28° for the KMH, C01, and C01 × 15 models, respectively, and decreases toward both sides for all these dust models.

Zoom In Zoom Out Reset image size

We assume an outer radius for the disk of 80 AU, based on our HiCIAO observations and millimeter interferometry by Isella et al. (2010). The inner radius is set to 1 AU. We find that our results are not significantly affected by the choice of inner radius as long as it is within the coronagraphic mask (r ∼ 30 AU at a distance to the target of 140 pc) and the disk surface has a flared geometry. Note that the millimeter interferometry of Isella et al. (2010) suggests the presence of a disk hole with a radius of ∼15 AU, while the presence of near-infrared emission suggests the presence of warm dust within 1 AU (Akeson et al. 2005; Pott et al. 2010).

The remaining free parameters for the density distribution of the disk are ρ₀, h_50 AU, β, and α. The first three parameters all have a large effect on the offset of the PI distribution from the major axis of the disk. This makes it difficult to search for the parameter sets which best fit the results. To overcome this problem, we set the total optical thickness between the star and the edge of the disk to 1 at a direction θ from the midplane of the disk. This angle will be adjusted to fit the degree of the offset of the PI distribution from the major axis of the disk (see Figure 6 for an example). This determines the parameter ρ₀, and as a result, our parameter searches will be made with the remaining two free parameters, i.e., h_50 AU and β. For the simulation in Section 5.2 we adopt α = β + 1 as we did in Section 4. The use of an independent α produces almost identical results.

Zoom In Zoom Out Reset image size

In addition to ρ₀, h_50 AU, β, and α, Equation (1) also includes the stellar radius R_*, but its uncertainty does not affect the results of simulations. The stellar radius is significantly smaller than the radii for the density distribution we use (i.e., >1 AU), and therefore the term $\sqrt{R_*/r}$ produces only a very minor contribution to Equation (1). The equation also includes the term (R_*/r)^α, but the stellar radius here can be replaced by any radius, depending on the radius at which we define ρ₀. Again, this parameter is not explicitly provided but scaled using the constraint of the optical thickness described above.

Throughout the simulations we assume that re-radiation from the disk is negligible, or that it occurs in very close proximity to the star so that the effect of the extended structure in the emission region far from the star is negligible. To prove the validity of this assumption, we ran this simulation code with the disk parameters in Table 2, and found that the code reproduces the PI distribution in Figure 3 well.

Figure 7 shows examples for the PI distribution using the C01 grain size distribution (i.e., that of the HH 30 disk). The combination of the angle for τ = 1 (θ = 30°) and the viewing angle (49° from the face-on view) are selected to approximately fit the observed morphology. As for the observations and models in Sections 3 and 4, the PI flux displayed is normalized to the stellar I flux. The image convolution is made in the same manner as Section 4. The contrast for each figure is adjusted in each image for the best morphological comparison for the bright part.

Zoom In Zoom Out Reset image size

In Figure 7, all except the case of h_50 AU = 25 AU, β = 1.3 show a butterfly-like morphology similar to the observed PI distribution. The peak PI flux in these models, however, is three to six times larger than the observed value. Increasing β produces a brighter PI distribution. In the case of β = 2.0 and 2.7, the modeled PI/I_* is larger for h_50 AU = 25 AU than h_50 AU = 5 AU by a factor of ∼1.5. This difference is attributed to the fact that, for the latter model, the star suffers extinction from the edge of the disk which decreases I_*. In the case of β = 2.0 and 2.7, increasing h_50 AU also produces a larger PI flux at the other side of the disk. The PI distribution is fairly centrosymmetric for h_50 AU = 25 AU and β = 1.3 for any contrast.

In the same figure we also show the density profiles for the individual disks. The dashed, solid, and dotted curves show the positions at which τ = 0.5, 1, and 2, respectively, measured from the star. These curves indicate that the scattering layer is geometrically thin at the disk surface. Increasing β produces a large flaring angle in the outer disk, which explains the brighter PI flux as the surface is more easily illuminated. The curvature of the surface remains similar for the same β; however, the disk is not flared for β = 1.3 and h_50 AU. In this case, the scattering is dominated by the inner part of the disk, producing a centrosymmetric PI distribution as described above.

The viewing angle selected for Figure 7 (49°) is smaller than that determined from millimeter emission by Isella et al. (2010) (>65°). A larger viewing angle does not reproduce the observed PI morphology. Figure 8 shows examples for a viewing angle of 63°. The disk parameters are identical to those for Figure 7, but the optical thickness along the radial direction has been set to τ = 1 at θ = 25° to fit the offset of the bright PI emission from the major disk axis. The results are similar to those of Figure 7, but the bright part shows a thinner distribution in the vertical direction than the observations. Furthermore, a large β results in a larger spatial extension at the bright side. To clearly show the former discrepancy, we extract the one-dimensional (1D) profiles at the positions indicated in Figures 7 and 8 and show them in Figure 9. The peak PI fluxes are comparable to or larger than those in Figure 7, and therefore significantly larger than the observations.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 10 shows the PI images for the same disk parameters as Figures 7 and 8 but with different dust models. Changing the dust model results in little variation in the PI images. This can be explained by the fact that the PI flux at the bright side is dominated by a single scattering with scattering angles 60°–135°, and the PI flux for this range of angles is similar for the three different dust models, as shown in Figure 5. In particular, the PI fluxes shown in Figure 5 for these scattering angles explain the fact that the PI flux for the KMH model is larger than the others in Figure 10 by a factor of ∼2. Table 4 shows the peak PI/I_* for individual models. These are larger than the observations by a factor of 2–12 at the peak.

Zoom In Zoom Out Reset image size

In summary, the modeled morphology in the PI distribution does not match the observations with the viewing angle inferred from the millimeter observations (>65°). As shown in Figures 7, 8, and 10, the modeled PI distribution from the bright side of the disk does not vary significantly with the value of h_50 AU, β, or the choice of dust model once we determine the outer radius of the disk and the optical thickness in a given direction. Some models may also be excluded due to the degree of polarization, the inferred disk mass or extinction toward the star. However, these do not provide constraints as clearly as those from the PI flux. We discuss these constraints in Appendix B.

In Sections 4 and 5, we found inconsistencies between the models obtained with conventional disk and dust models and the observed PI distribution. The disk geometry obtained from the SED corresponds to a PI distribution with a significantly smaller offset from the major axis of the disk (Section 4). On the other hand, the disk geometries which reproduce this offset require a smaller viewing angle (i.e., close to face-on) than that determined through millimeter observations (Section 5). We have not explored the full parameter space, but once we fix the parameters to reproduce the angular scale along the major axis and the offset of the bright PI emission on the far side of the disk, varying disk and dust parameters has little effect on the resulting PI distribution.

The results of the two sets of simulations (full radiative transfer and monochromatic scattering) may imply that the system consists of (1) a geometrically thin disk which is partially responsible for the infrared SED but does not contribute to the PI flux in the near-infrared and (2) an optically thin and geometrically thick upper layer which is responsible for the observed PI distribution in the near-infrared and the remaining mid-to-far infrared flux. A similar geometry has been proposed for SEDs observed for some Herbig AeBe stars (Meeus et al. 2001, Group I). Follette et al. (2013) also discuss such a geometry for the SR 21 disk in which the presence of a inner cavity is indicated by submillimeter observations but not seen in the scattered light in near-infrared.

Figure 11 shows a schematic view for the disk, an optically thin layer, stellar radiation, scattered light, and infrared re-radiation. A rough analytic estimate for the optical thickness, density, and mass of such a layer is described in Appendix C. The observed PI image does not show evidence for emission along the major axis without an offset as modeled in Section 4. This can be explained if the thickness of the disk is significantly smaller than those used in Section 4, which provide a PI flux level comparable to the observations.

Zoom In Zoom Out Reset image size

Such a scattering geometry can easily explain the observed PI flux and offset from the major axis of the disk, by adjusting the column density and vertical distribution of the optically thin layer, respectively. Here we qualitatively demonstrate this using the Sprout code with an optically thin upper layer using Equation (1) and an optically thick disk with a geometrical thickness of zero. To focus on the scattered light outside the coronagraphic mask, we set the density to be zero for disk radii within the mask (r < 28 AU).

Figure 12 shows an example of simulations using the KMH dust model, with a viewing angle of 70°. The butterfly-like bright PI distribution approximately matches the observations. We measure the modeled polarization degree to be 25% at a position corresponding to the PI peak of the observations. This is also consistent with the observations of ≫9% measured in Section 3. The faint extended emission at the other side of the disk is brighter than the upper limit of the observations by a factor of ∼4. The modeled flux in this region could be adjusted, e.g., if we included extinction by a remnant envelope, or if we define an exponential cutoff for the outer radii of the optically thin layers like that used to reproduce the millimeter flux distribution of the disk (e.g., Andrews et al. 2009; Isella et al. 2010).

Zoom In Zoom Out Reset image size

The use of different dust models (C01, C01 × 15) provides similar results to Figure 12, but the polarization degree is significantly smaller: 13%–14% at the position corresponding to the PI peak of the observations. Such a polarization may be too close to the lower limit of the observations of ∼9%. It is beyond the scope of the paper to identify the best dust models using the modeled degree of polarization. See Section 6.2 for more complex and realistic dust models than used in our study.

Use of a geometrically thinner disk than Section 4 implies that the disk gets less stellar flux, thereby resulting in less infrared excess than that shown in Figure 2. This issue could be solved as follows. As suggested by Chiang & Goldreich (1997), an optically thin layer on the disk surface scatters or reemits directly to space about half the radiation, while the other half is scattered/emitted inward. The extra illumination from the optically thinner layer would warm up the disk and enhance the observed infrared flux. Although such an optically thin layer is usually assumed to be geometrically thin (e.g., Chiang & Goldreich 1997; Dullemond et al. 2007), the same physical mechanism should work even in the case where the layer is geometrically thick. If the scattered light is dominated by the optically thin layer, this implies that the disk surface would be illuminated more significantly from this layer than the direct stellar flux. Therefore, the presence of this layer may significantly affect the temperature structure of the disk.

In addition, the thermal radiation from the optically thin layer may directly contribute to the infrared SED (e.g., Chiang & Goldreich 1997). However, the fraction compared to the disk flux would be significantly smaller than that of the scattered light in the near-infrared due to significantly smaller dust opacities in the mid-to-far infrared (e.g., Wood et al. 2002; Whitney et al. 2003a, 2003b; Dong et al. 2012b).

While homogeneous spherical particles are widely used for the dust in disks (Whitney & Hartmann 1992; Fischer et al. 1994; Whitney et al. 2003a, 2003b; Robitaille et al. 2006, 2007), dust particles in practice appear to be aggregated (e.g., Henning & Stognienko 1996; Dominik 2009), and may also be coated with ice (Malfait et al. 1999; Meeus et al. 2001; Honda et al. 2009). These facts would produce different optical properties. Min et al. (2012) conducted simulations for scattered light from the disk using both spherical and aggregated particles, and show that a lower PI flux is expected for the latter due to the extremely large forward scattering. This would allow us to better match the PI flux we simulated in Section 5 for the disk surface, as our values were significantly larger than the observations for many cases.

However, the disks used in Section 5 would produce significantly different SEDs than those shown in Figure 2. A large flaring allows the disk to receive more stellar photons, which warms up the disk more efficiently, thereby producing a larger excess than the thinner disks we used in Section 4. Furthermore, some disks used in Section 5 have a large extinction toward the star (Appendix B). This would result in a double peaked SED in the optical to near-infrared, and far-infrared, respectively (e.g., D'Alessio et al. 1999; Robitaille et al. 2007). The former problem could be overcome if the dust particles had a smaller emissivity in mid-to-far infrared, or large albedo at UV-to-IR in order to impede the absorption of the stellar photons. To investigate further, one would need a detailed understanding of the optical properties of the aggregates over a wide range of wavelengths, plus simulations with full radiative transfer.

Alternatively, observations of the PI distribution at different wavelengths would give useful constraints for the nature of the dust grains. We expect that the PI flux is dominated by the photons with a single scattering, even in the case of optically thick disks (such as those used in Sections 4 and 5) in which multiple scattering occurs as well as single scattering. The contribution from multiple scattering would be negligible for the reasons below. In the case of small grains (e.g., KMH), where photons are fairly isotropically scattered, the scattered photons have polarizations with a variety of P.A., canceling each other out. In the case of large grains (e.g., C01 and C01 × 15), in which most of the photons are scattered forward, the polarization of individual photons is significantly reduced after the first scattering.

To demonstrate this, we re-ran the simulations for the disks in Section 5, removing all the photons which experience multiple scatterings. We find that the results are almost identical for the PI flux distribution, with a flux difference within 10% at the peak. This is in contrast with the I flux for which multiple scatterings are responsible for 20%–40% of the flux. Throughout, studies of the PI flux allow a more accurate investigation of the dust properties than studies of the I flux, in particular for optically thick disks.

As shown in Section 3, the observed PI distribution is asymmetric about the rotation axis of the disk: the southwest side is brighter than the northern side by a factor of ∼2. A similar asymmetry in the distribution of the scattered light has been extensively studied for the HH 30 disk, i.e., a low-mass protostar at a similar evolutionary stage, using the Hubble Space Telescope (e.g., Burrows et al. 1996; Stapelfeldt et al. 1999; Cotera et al. 2001; Watson & Stapelfeldt 2007). These authors attributed such a distribution to uneven illumination by the central object. These observations also show that the distribution of the scattered light is time variable. Watson & Stapelfeldt (2007) summarize possible mechanisms to explain the asymmetric distribution and its time variation. These include (1) hot (or cool) spots on the star, (2) shadowing by a non-axisymmetric inner disk, (3) obscuration by a companion star, and (4) obscuration by a disk associated with a companion star. The time variability in the scattered light is also observed in the more evolved disk system HD 163296 (Wisniewski et al. 2008). The authors attributed this to the asymmetric shadowing of the inner disk.

Of the explanations described above, shadowing by the inner disk may be the most likely for the asymmetric PI distribution in RY Tau. While uneven illumination by hot (or cool) spots is one of the favored explanations for the HH 30 disk (Wood & Whitney 1998; Stapelfeldt et al. 1999; Cotera et al. 2001; Watson & Stapelfeldt 2007), it is not likely for RY Tau: the optical continuum associated with this star does not show clear evidence for time variability due to such spots (e.g., Petrov et al. 1999; Chou et al. 2013). Obscuration by a secondary star/disk may also be excluded, since there is no clear evidence for the presence of a close companion associated with RY Tau. Near-infrared interferometric observations have ruled out the presence of a companion at 0.35–4 AU from the star and down to a stellar flux ratio of 0.05 (Pott et al. 2010).

The shadowed-disk scenario is corroborated by the UXOR-type variability observed toward this star, i.e., time-variable obscuration by circumstellar dust (e.g., Herbst et al. 1994; Petrov et al. 1999). While many authors attribute this obscuration to the disk, its nature is not clear. It could be clumpy accretion of dust+gas onto the star (see Herbst et al. 1994; Grady et al. 2000, and references therein), the puffed-up inner rim of the disk (Dullemond et al. 2003), or the outer edge of the disk (Thé 1994). Herbst et al. (1994) also state the possibility that such a dust screen is associated with a wind, not a disk or disk accretion. Of the above explanations, the outer edge of the disk may be excluded if the disk associated with RY Tau is geometrically thin, and we are not observing this target close to an edge-on view (Sections 4 and 6.1).

RY Tau showed two abrupt brightening events at optical wavelengths in 1983/1984 and 1996/1997 (e.g., Herbst et al. 1994; Petrov et al. 1999). If we attribute this to a single orbital period around a 2 M_☉ star, the corresponding radius of the disk is 7 AU. This radius is comparable to the radius where Isella et al. (2010) identified emission peaks at two sides of the disk using millimeter interferometry. If obscuration occurs in the inner disk, it would require time variation of the disk structures including warping and precession. This may be possible via tidal interaction with a very low mass companion such as a protoplanet (Hughes et al. 2009).

Synoptic observations of the scattered light of RY Tau would allow us to identify the origin of the asymmetric PI distribution, determine a typical radius where the obscuration occurs, and constrain the disk geometry (and its time variation) within the coronagraphic mask. Bastien (1982) made monitoring observations of optical polarization integrated over the object, revealing variation of the position angle of polarization between −25° and 45° over a few years.²⁸ Spatially resolved observations like ours would have significant advantages for searching for a periodic variability toward this active and complicated object. We also note that the coronagraphic observations for the optical I image of RY Tau were made in 2007–2009, and these show a brighter lobe in the north than the southwest (J. McCleary et al. 2007, private communication) in the outer envelope, i.e., the opposite side as we observed in the PI distribution. Although these observations show scattered light in significantly more distant outer regions, this is consistent with the idea of time variation of scattering light associated with this object.