SURFACE GEOMETRY OF PROTOPLANETARY DISKS INFERRED FROM NEAR-INFRARED IMAGING POLARIMETRY^*

• 1.  
  The analysis below applies to disks for which the surface emission can be observed, with the viewing the angle i not close to 90° (the edge-on view).
• 2.  
  The disk is optically thick in the near-infrared, at least for the wavelength range encompassed by our observational data (see Watson et al. 2007, for a review). This fact has been observationally confirmed in some young edge-on disks, where the disk appears as a dark lane (e.g., Padgett et al. 1999), and in disks silhouetted against bright nebular backgrounds (see McCaughrean et al. 2000, for a review). The optically thick nature of the disk is also consistent with some simulations that model near-infrared scattered light from protoplanetary disks (e.g., Silber et al. 2000; Duchêne et al. 2004; Dong et al. 2012a).
• 3.  
  A disk surface is defined by the locus of points, between the star and outer disk radius, with an optical depth of approximate unity (Watson et al. 2007; Jang-Condell 2007, 2008; Perrin et al. 2009; Muto 2011). The width of the scattering layer on the disk is significantly smaller than the geometrical thickness of the disk. This approximation has been used for theoretical determinations of the radiative heating of the disk by the central star (Kenyon & Hartmann 1987; Chiang & Goldreich 1997). Takami et al. (2013) numerically examined this issue for some disk models, finding that the positions of the τ = 0.5, 1 and 2 curves approximately match, agreeing with the above approximation. Observationally, this approximation is corroborated by optical and near-infrared observations of some edge-on disks, which show a clear boundary between the shadowed flared disk and the reflection nebulae (e.g., Padgett et al. 1999; Watson et al. 2007).We note that the above definition of the disk surface does not allow the exploration of regions of the disk which do not receive any photons, or in other words, regions of the disk where we cannot define a position corresponding to an optical thickness τ ∼ 1. These cases include (a) a disk surface parallel to the light path from the star; (b) areas outside the disk where dust does not exist, or regions shadowed by the outer edge of the disk; or (c) regions shadowed by an inner disk structure (e.g., a puffed-up inner rim due to heating from stellar radiation (e.g., Dullemond et al. 2001) or by the protruding inner edge of a disk gap (e.g., Jang-Condell & Turner 2012, 2013, see also Section 3).We also note that we cannot accurately determine the disk surface if an optically thin layer above the disk (e.g., a remnant envelope) contributes to the optical thickness between the star and the disk surface (e.g., Stapelfeldt et al. 2003; Follette et al. 2013; Takami et al. 2013).
• 4.  
  The disk is geometrically thin (z/r ≪ 1, where z and r are the height from the midplane and the distance from the star in the midplane, respectively). This assumption is reasonable for most protoplanetary disks (see Watson et al. 2007; Dullemond et al. 2007, for reviews). We note that a "geometrically thin disk" does not always imply that the model disk is flat (e.g., Whitney & Hartmann 1992). For instance, a variation in z or z/r, including a concave-up or concave-down geometry, is still possible with the constraint z/r ≪ 1 (see Sections 3 and 4 for examples).
• 5.  
  As shown in Figure 1, α is the angle of the surface of the disk at the position of scattering with respect to the midplane. θ is the elevation angle of this position from the midplane. We assume that α is not significantly greater than θ (i.e., α − θ ≪ 1). This is not a standard assumption for protoplanetary disks but we confirm its validity for our sources in Sections 3 and 4. We note that α is always equal to or larger than θ due to the definition of the surface described above.
• 6.  
  The star is treated as a point source of illumination. In practice, a pre-main sequence star has a stellar radius of a few solar radii (Stahler & Palla 2005). As shown below, this approximation does not cause significant errors.
• 7.  
  The contribution of photons, with multiple scattering, to the PI flux is negligible. This approximation is supported by Takami et al. (2013). For the scattering geometries used in that study, the contribution of multiply scattered photons to the PI flux was within 10% of that for singly scattered photons. This is because, in the case of small grains, photons are relatively isotropically scattered and the scattered photons have polarizations with a variety of position angles (P.A.s), thus canceling each other out. In the case of large grains, most of the photons are scattered forward, resulting in a significant reduction in individual photon polarization after the first scattering.

In the case where the contribution of multiple scatterings is negligible (Approximation 7), the scattered PI intensity from the disk surface is described as follows.

$$\begin{equation} I_{\nu ; \rm out} (r,\phi) {dA}_{\rm out} \,d \Omega _{\rm out} = \frac{L_\nu }{\rm 4 \pi R^{\rm 2} (r,\phi)} \left(\frac{ {\rm PI}}{I_{\rm 0}} \right) {dA}_* d \Omega _{\rm out}, \end{equation} \tag{ 1 }$$

where $I_{\nu ; \rm out}(r,\phi)$ is the observed PI intensity at the distance r from the star and the azimuthal angle ϕ in the plane of the disk, dA_out is the differential area of the disk surface (dA) projected onto the plane of the sky, dΩ_out is the solid angle from the disk surface toward the observer, and L_ν is the stellar luminosity per unit frequency. R(r, ϕ) is the distance from the star to the surface of the disk, i.e., the position of a scattering event. This parameter is equal to r/cos θ(r, ϕ), where θ(r, ϕ) is the elevation angle of the position of a scattering event from the midplane. (PI/I₀) is the fraction of the scattered PI intensity relative to the incident I₀ flux (sr⁻¹). dA_* is the differential area dA projected on the plane perpendicular to the light path from the star to the disk.

The right side of Equation (1) is a product of the incident light on a unit differential area at the disk surface (L_νdA_*/4πR²) and the fraction of the energy radiated in the given direction ((PI/I₀)dΩ_out). The parameter (PI/I₀) in Equation (1) depends on the dust particles and the scattering angle (Figure 5 of Takami et al. 2013). For nonpolarimetric disk imaging at optical and near-infrared wavelengths, this ratio would be replaced by the "phase function."

The differential area of the disk surface projected onto the plane of the sky dA_out is

$$\begin{equation} {dA}_{\rm out} = {dA} \cos (i+\alpha), \end{equation} \tag{ 2 }$$

where i is the disk inclination or the viewing angle.

The parameter dA_*, the differential area projected on the plane perpendicular to the light path from the star to the disk, is

$$\begin{equation} {dA}_* = {dA} \sin (\alpha - \theta). \end{equation} \tag{ 3 }$$

Substituting Equations (2) and (3) into Equation (1), we derive

$$\begin{equation} \left(\cos i - \sin i \tan \alpha \right) {I_{\nu ; \rm out} }=\frac{L_\nu }{4\pi R^2} \left(\frac{\rm PI}{I_{\rm 0}} \right) (\tan \alpha \cos \theta - \sin \theta). \end{equation} \tag{ 4 }$$

Since z/r = tan θ and dz/dr = tan α, where z is the height of the disk surface from the midplane,

$$\begin{equation} \tan \alpha =dz/dr = \tan \theta + (r/\cos ^2 \theta) d\theta /dr. \end{equation} \tag{ 5 }$$

Substituting Equation (5) into Equation (4) we derive

$$\begin{equation} {d\theta \over dr}={\cos i - \sin i \tan \theta \over C \sin i + \cos \theta }{4\pi r I_{\nu,{\rm out}} \over L_\nu ({\rm PI}/I_0)}, \end{equation} \tag{ 6 }$$

where

$$\begin{equation} C \equiv \frac{4\pi R^2 I_{\nu ; \rm out}}{L_\nu ({\rm PI}/I_0)} = \frac{{dA}_*}{{dA}_{\rm out}} = \frac{\sin (\alpha - \theta)}{\cos (i+\alpha)} \end{equation} \tag{ 7 }$$

As described above, 0 ⩽ α − θ ≪ 1 (Approximation/Assumption 5) and the disk is not close to edge-on (i.e., i ≫ α, Approximation/Assumption 1), therefore C ≪ 1. Furthermore, the disk is geometrically thin (tan θ ≪ 1, Approximation/Assumption 4) and the disk is not close to edge-on (i.e., not cos i ≪ sin i, Approximation/Assumption 1), therefore cos i ≫ sin itan θ in Equation (6). Therefore Equation (6) is approximated by

$$\begin{eqnarray} {dS} & \sim & \frac{ \rm 4 \pi r I_{\nu ; \rm out }(r,\phi) \rm cos\ i }{L_\nu ({\rm PI}/I_0)} dr. \end{eqnarray} \tag{ 8 }$$

where S = sin θ.

Using the following respective equations, we can replace $I_{\nu ; \rm out} (r,\phi)$ and L_ν with the observed number of photons per second (n_{ν; PI}, n_{ν; *}):

$$\begin{eqnarray} n_{\nu ;{\rm PI}} (r,\phi)& = &\frac{I_{\nu ; \rm out} (r,\phi)}{h\nu } A_{\rm pix} \Omega _{\rm Tel}, \nonumber \\ n_{\nu ;*} & = & \frac{1}{h \nu } {\frac{L_\nu }{\rm 4 \pi }} \Omega _{\rm Tel}, \nonumber \\ \frac{n_{\nu ;{\rm PI}} (r,\phi)}{n_{\nu ;*}} &=& \frac{\rm 4 \pi I_{\nu ; \rm out} (r,\phi) A_{\rm pix}}{L_\nu } \end{eqnarray} \tag{ 9 }$$

where hν is the photon energy; A_pix is the area corresponding to the HiCIAO pixel scale (9.5 mas) at the distance to the target; Ω_Tel is the solid angle corresponding to the telescope mirror ($\Omega _{\rm Tel} = A_{\rm Tel}/4 \pi d^{\rm 2}$, where A_Tel and d are the area of the telescope mirror and the distance to the target, respectively). Substituting Equation (9) in Equation (8), we derive:

$$\begin{eqnarray} {dS} & \sim & r \frac{n_{\nu ;{\rm PI}}(r,\phi)}{n_{\nu ;*}} \frac{\rm cos\ i }{A_{{\rm pix}}} \left(\frac{{\rm PI}}{I_{\rm 0}} \right)^{-1} dr. \end{eqnarray} \tag{ 10 }$$

The surface geometry z(r, ϕ) = r tan [θ(r, ϕ)] can then be derived via numerical integration of dS, using Equation (10), along a given azimuthal angle in the disk coordinate (ϕ) or the observed P.A. The input parameters are the observed PI distribution normalized to the stellar I flux (n_PI/n_*), the inclination i, and an assumed dust property (PI/I₀). For a given ϕ or P.A., the scattering angle should be approximately constant over any position if the disk is geometrically thin. Thus, (PI/I₀) is approximately constant if the dust properties do not vary along with P.A. For simplicity, we show all the results in later sections in terms of P.A.s of the observations, rather than azimuthal angles ϕ in the disk coordinate.

We use the extended Newton–Cotes formula (Press et al. 1992) to numerically integrate Equation (10). We tested this algorithm as follows. First, we provided the surface geometry using an analytic function and analytically derived the PI distribution. Then we performed the numerical integration and investigated if the calculated surface geometry matches the original geometry. We found that the surface geometry is reproduced well.

The initial value of the numerical integration S₀ (=sin θ₀, where θ₀ is the elevation angle of the position of a scattering event from the midplane) can be arbitrarily given at radius r₀. This convention corresponds to an arbitrary inclination of the disk surface with respect to the midplane. The observed flux is independent of S₀ or θ₀ because an arbitrary inclination will not change the variation in the illumination of the disk.

As described above, we have used the assumption of 0 ⩽ α − θ ≪ 1 (where α is the angle of the surface with respect to the midplane and θ is the elevation angle of the position of a scattering event from the midplane), which may not be applicable for all cases. This can be verified by investigating whether C ≪ 1 using Equation (7). Substituting Equations (9) to Equation (7) we derive

$$\begin{eqnarray*} &&C = \frac{n_{\nu ;{\rm PI}}(r,\phi)}{n_{\nu ;*}} \frac{R^2}{A_{\rm pix}} \left(\frac{{\rm PI}}{I_{\rm 0}} \right)^{-1}, \nonumber \end{eqnarray*}$$

or

$$\begin{equation} C \sim \frac{n_{\nu ;{\rm PI}}(r,\phi)}{n_{\nu ;*}} \frac{r^2}{A_{\rm pix}} \left(\frac{{\rm PI}}{I_{\rm 0}} \right)^{-1}. \end{equation} \tag{ 11 }$$

if the disk is geometrically thin (i.e., r ∼ R). Therefore, we can determine whether C ≪ 1 using the observed parameters (n_{ν; PI}/n_{ν; *}, r), the instrument parameter (A_pix) and the scattering properties based on the dust model (PI/I₀).

To derive Equation (10), we regard the star as the point source of illumination. In practice, a pre-main sequence star has a stellar radius of a few solar radii (Stahler & Palla 2005). Given the considered angular radius of the star from the disk surface, the point-source approximation yields an error on order of 10⁻⁴ for θ, S, and z/r, at 100 AU. This error is significantly smaller than the disk structure scales we discuss in later sections.

Table 6 shows the list of objects we have analyzed. This includes disks with spiral structures (SAO 206462, MWC 758), ring-like morphologies (PDS 70, 2MASS J16042165–2130284), and a relatively uniform distribution (MWC 480). The initial SEEDS results for these objects have been published by Muto et al. (2012), Hashimoto et al. (2012), Kusakabe et al. (2012), Mayama et al. (2012), and Grady et al. (2013). All except PDS 70 show bright millimeter/submillimeter emission (Hamidouche et al. 2006; Brown et al. 2009; Isella et al. 2010; Andrews et al. 2011; Mathews et al. 2012), indicating the youth of these disks, and therefore suggesting an optically thick nature in the near-infrared (Section 2). A 03 diameter circular occulting mask was used to suppress the bright stellar flux for SAO 206462, MWC 758 and MWC 480. The occulting mask was not used for the observations of PDS 70 or 2MASS J1604. See the above papers for the details of data reduction.

We need to scale the observed PI flux by the spatially integrated stellar I flux to apply Equation (10). The stellar I flux for PDS 70 was measured using the same data set as those used for the PI observations. For the other objects we used short exposures with a neutral density filter without a coronagraphic mask. These images were taken on the same date as the coronagraphic images. See Takami et al. (2013) for technical details of the measurement.

We use the same dust model used in Section 3. As shown in Table 3 this model yields (PI/I₀) = 0.012 sr⁻¹ with a scattering angle of 90° (Takami et al. 2013). This is the approximate scattering angle expected along the major axis of the PI flux distribution, or the entire disk for a face-on view. Dust grains may be smaller than those defined by the KMH model in some disks (see the Appendix). As seen in Equation (10), the use of different dust models (and the corresponding PI/I₀) results in the linear scaling of the surface geometry derived. This scaling does not affect the main conclusions of the paper.

To derive the surface geometry several P.A.s were initially chosen for SAO 206462, MWC 758, and 2MASS J1604, which clearly show nonaxisymmetric distributions³⁶ (Figures 4–6). These disks have inclination angles of 11°, 21°, and 10°, respectively, implying that the scattering angles are slightly different at different P.A.s. Thus, the true (PI/I₀) should be slightly smaller in practice than that determined at 90°, depending on the position. These differences are 6%–7% for SAO 206462 and 2MASS J1604, and 17% for MWC 758 along the minor axis, where we expect the largest deviation, assuming the same scattering properties of dust grains measured for PDS 70 and MWC 480 (see the Appendix). For simplicity we use the (PI/I₀) derived for scattering angles of 90° over the entire disk. The variations in the corrections for different scattering angles are significantly smaller than the calculated variation in the disk surface structure or illumination shown below, therefore this approximation does not affect our conclusions. The PDS 70 and MWC 480 disks have fairly large inclination angles, therefore we derive the surface geometries along the major axis only (see Section 3). We calculate the constant C defined in Equations (7) and (11) to be ∼0.3 or smaller for all the objects and P.A.s described above. This agrees with the approximation we have used to derive Equation (10) (C ≪ 1).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The numerical integrations for 2MASS J1604 and PDS 70 were made starting at 025 (∼35 AU) from the star, approximately corresponding to the minimum in the radial PI flux in the ring structure. The integrations for the other objects were made starting from 02 (28 AU for SAO 206462 and MWC 480; 40 AU for MWC 758) from the star, i.e., slightly outside the coronagraphic mask and the innermost radius at which we have reliable measurements of the PI flux distribution.

Figures 4–8 show the observed PI images, the position at which the PI intensity profile was extracted, and the calculated surface geometry. For each object a few values are arbitrarily chosen for S₀ (=sin θ₀, where θ₀ is the elevation angle of the position of a scattering event from the midplane at the radius r₀) to show how the location of the disk surface varies with this parameter. The outermost regions with constant or nearly constant S are shown with fainter curves, since we cannot determine if disk material is present or the region is shadowed (Section 2). We note that the diffuse emission at these radii may also suffer from an artificial halo, a result of the faint halo of the point-spread function smearing the emission in the bright inner regions. In all the plots the height of the disk surface from the midplane is significantly smaller than the corresponding radius, consistent with the first approximation of Section 2.

In Figures 4–8 all the objects except MWC 480 show local maxima in the plotted one-dimensional PI distributions. However, none of them show local maxima in the disk height z in these figures. Instead, the parameter S (therefore θ) shows the largest increases near the bright PI maxima, slightly offset toward smaller r. These positions are: r ∼ 50 AU, P.A. = 300° in SAO 206462; r ∼ 90 AU, P.A. = 160° in MWC 758; and the position of the ring in the PI distribution in 2MASS J1604 and PDS 70. For the rest of the paper we describe such geometries as "local concave-up structures" at the disk surface. The presence of local concave-up structures is independent of the assumed S₀ (=sin θ₀).

As shown in Figures 6 and 7 and explained in Section 2, a strong deficit in the PI flux in the ring structure does not imply the absence of an inner disk. This can be explained with either a geometrically thin disk whose surface is nearly parallel to the light path of the star (i.e., a nearly constant S and θ), or self-shadowing by the inner disk (see Section 2). The disk geometries shown in Figures 6 and 7 do not show the presence of a sharp vertical inner wall. This contrasts with discussions or assumptions in many observational studies of disks with a hole or gap (see Espaillat et al. 2014, for a review). We will further explore this issue in Section 4.4 with a monochromatic radiative transfer simulation.

Zoom In Zoom Out Reset image size

We plot the surface geometry for individual objects for S₀ = 0 in the left column of Figure 9. All the disks except MWC 480 show large nonaxisymmetries. This contrast to the axisymmetric geometry which is assumed for many disk models (see, e.g., Dullemond et al. 2007, for a review). We investigate below if a more axisymmetric geometry is possible assuming either: (A) different values of S₀ (=sin θ₀, where θ₀ is the elevation angle of the position of a scattering event from the midplane at distance r₀) at different azimuthal angles in the disk coordinate ϕ (or the observed P.A.); or (B) a nonaxisymmetric illumination of the disk (see Section 5.2 for possible origins of such illumination). In this subsection we will also discuss other behavior in the concave-up and concave-down structures not clearly seen in Figures 4–8.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

For Case (A) the parameter S₀ is obtained as follows. From Equation (10) we derive

$$\begin{eqnarray} &&S (r_{1}, \phi) = S_0 (\phi) + \frac{\rm cos\ i }{n_{\nu ;*} A_{{\rm pix}}} \left(\frac{{\rm PI}}{I_{\rm 0}} \right)^{-1} \int _{r_0}^{r_{1}} r n_{\nu ;{\rm PI}}(r,\phi)\, dr, \nonumber\\ \end{eqnarray} \tag{ 13 }$$

where r₁ is the outer radius; ϕ is the azimuthal angle in the disk plane; i is the viewing angle of the disk; n_* is the observed stellar I flux; A_pix is the area corresponding to the pixel scale of the instrument; (PI/I₀) is the PI flux normalized to the incident I flux determined by the dust properties; r₀ is the inner radius where we define S₀; n_PI is the observed PI flux distribution from the disk. If the disk height is constant at the radius r₁, S(r₁, ϕ) = S₁ (i.e., a constant). Substituting this into Equation (13) we derive S₀(ϕ) and a disk height z(ϕ) of

$$\begin{eqnarray} S_0 (\phi) & = & S_1 - \frac{\rm cos\ i }{n_{\nu ;*} A_{{\rm pix}}} \left(\frac{{\rm PI}}{I_{\rm 0}} \right)^{-1} \int _{r_0}^{r_{1}} r n_{\nu ;{\rm PI}}(r,\phi)\, dr \nonumber \\ & = & S_1 - S(r_{1},\phi, S_0=0). \end{eqnarray} \tag{ 14 }$$

For Case (B) Equation (10) becomes

$$\begin{equation} {dS} = r \frac{n_{\nu ;{\rm PI}}(r, \phi)}{f(\phi) n_{\nu ;*} } \frac{\rm cos\, i }{A_{{\rm pix}}} \left(\frac{{\rm PI}}{I_{\rm 0}} \right)^{-1} dr, \end{equation} \tag{ 15 }$$

where f(ϕ) is the fraction of stellar flux passing through the dust screen between the star and the disk surface. Thus,

$$\begin{equation} S(r_{1}, \phi) = \frac{\rm cos\, i }{f(\phi) n_{\nu ;*} A_{{\rm pix}}} \left(\frac{{\rm PI}}{I_{\rm 0}} \right)^{-1} \int _{r_0}^{r_{1}} r n_{\nu ;{\rm PI}}(r,\phi)\, dr, \end{equation} \tag{ 16 }$$

if S₀ = sin θ₀ = 0. Substituting S(r₁, ϕ) = S₁ we derive

$$\begin{eqnarray} f(\phi) & = & \left[ \frac{\rm cos\ i }{n_{\nu ;*} A_{{\rm pix}}} \left(\frac{{\rm PI}}{I_{\rm 0}} \right)^{-1} \int _{r_0}^{r_{1}} r n_{\nu ;{\rm PI}}(r,\phi) dr \right] \left/\right. S_1 \nonumber \\ & = & \frac{S (r_{1}, \phi, f=1, S_0=0)}{S_1}. \end{eqnarray} \tag{ 17 }$$

In Figure 9 we plot the revised locations of the disk surface for Cases (A) (i.e., different S₀ for different P.A.s) and (B) (nonaxisymmetric illumination, S₀ = 0). Table 7 shows the parameters used for Cases (A) and (B) in Figure 9. The parameter S₁ at radius r₁ is chosen to be an approximately minimum possible value for each object, which results in the minimum S₀ and the maximum  f  being ∼0 and ∼1, respectively. The use of different S₁ values does not significantly change the results described below.

In Figure 9 the surface geometries at different P.A.s are similar, in particular for Case (B) (i.e., with nonaxisymmetric illumination) with the exception of the MWC 758 disk, P.A. = 160°. The MWC758 disk shows a clear concave-up geometry at r = 60–100 AU at this position angle. This concave-up geometry is associated with a faint arm in the southeast, seen in the PI image in Figure 5.

Figure 9 shows that Equations (8) and (11) produce the same local concave-up and concave-down structures. The surface geometry is close to a straight line for all the other P.A.s for SAO 206462 and MWC 758 which are associated with spiral structures. In other words, most of the spiral structures in these objects are attributed to surface structures whose spatial variation in the z-direction is significantly smaller than the thickness of the disk. As shown in Figures 6–8, the disk surfaces in 2MASS J1604 and PDS 70 show a distinct concave-up geometry at small radii, and a marginal concave-down geometry at large radii. The surface in MWC 480 shows a marginal concave-up geometry.

What is the relationship between the spiral structures and the local concave-up/down surface geometries in the SAO 206462 and MWC 758 disks? To clearly demonstrate this, we plot the disk height in these objects with P.A.s at 30° intervals for Case (B) (nonaxisymmetric illumination, S₀ = 0) in Figure 10. To easily separate concave-up/down structures from the effects of the surface inclination of the disk, we also subtract the linear component of the surface geometry characterized by the positions (r₀,0) and (r₁, z₁), i.e., z = z₀(r − r₀)/(r₁ − r₀), and plot them in the same figure.

Zoom In Zoom Out Reset image size

In Figure 10 all the local concave-up structures are attributed to a few prominent spiral arms identified in the PI distribution. There may also be a concave-up structure (thereby a spiral arm) at P.A. = 30°–90°, r = 100–120 AU in the MWC 758 disk, but a higher signal-to-noise in the PI flux is required to confirm this.

As discussed in the above subsections, we demonstrate that our approach is very useful for determining the surface geometry from the observed PI distribution. In order to reinforce our approach and results, we will now reproduce the observed PI flux distribution for PDS 70 and SAO 206462 using the disk geometry computed by Equation (10) and Monte-Carlo scattering simulations as described below.

As for Section 3 a central unpolarized point source (star) serves as the starting point for calculating the scattering of photons from the disk surface. We use the same interstellar dust model as in Section 3. For each photon reaching the disk surface, we calculate the scattering angle and Stokes parameters 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. The photons escaping from the disk are collected in a imaginary detector at a viewing angle of the observations. As outlined in Section 2, multiple scatterings are not included, and photons scattered inward to the disk are discarded.

For PDS 70 we use an axisymmetric surface geometry parameterized by S (=sin θ, where θ is the elevation angle of the position of a scattering event from the midplane) as for Section 4.2. The centroid of the observed ring structure is offset by ∼6 AU along the minor axis due to its thickness in the vertical direction (Hashimoto et al. 2012). The parameter S₀ is set to 0.1 at r₀ = 35 AU to reproduce this offset. We derive the height of the surface z by averaging the values at two P.A.s along the major axis (159° and 339°). This surface geometry produces disk heights of 3.5 and 21 AU (57 and 85 from the midplane) at r = 35 and 140 AU, respectively. The viewing angle of the simulation is set to be 50°, matching the observed image, integrated over a range of ±10°. We do not include illumination below 57 from the midplane, which may cause back scattering at the inner edge of the ring. To produce a sufficient signal-to-noise, we used 10⁷ photons for the model. In order to normalize the PI flux to the stellar I flux, we separately calculate the expected number of photons based on the incident number of photons.

Figure 11 shows the comparison between the observed and modeled PI distributions. The figure shows that the modeled PI distribution explains the observed PI distribution fairly well, typically within 30% in flux. The observed and modeled PI distributions match particularly well along the minor axis, where one might expect back scattering from the vertical wall at the far side of the ring. Figure 11(c) shows the residuals of the modeled image subtracted from the observed image. The observed fluxes are marginally larger than the model fluxes in the southwest of the ring, and marginally lower in the northeast of the ring. These may be due to either of the following: (1) the size distribution (and thereby the scattering properties) of dust grains deviating from the dust model we have used (see the Appendix), or (2) local deviation from symmetry of the surface geometry.

Zoom In Zoom Out Reset image size

SAO 206462 shows a nonaxisymmetric PI distribution a close to face-on view. We therefore derived the disk surface for 36 P.A.s at 10° intervals from the observed PI image for S₀ = 0. The surface positions between these P.A.s are linearly interpolated. The imaginary detector is set to collect photons at a P.A. of 155° integrated over a range of ±10° and a viewing angle of 11° integrated over a range of ±5° in order to match the observed image. We have used 10⁹ photons for the simulated image. Figure 12 shows the observed and simulated images, which match well.

Zoom In Zoom Out Reset image size

Young planets may tidally interact with protoplanetary disks, resulting in spiral structures in the disk (e.g., Papaloizou et al. 2007, for a review). For example, the spiral structures observed in the SAO 206462 and MWC 758 disks have been discussed as a possible signpost of ongoing planet formation (Muto et al. 2012; Grady et al. 2013).

Tidal interaction with planets may cause a variety of waves in protoplanetary disks. Lubow & Ogilvie (1998) categorized as follows based on analogy with the theory of stellar oscillations: (1) p-modes, i.e., acoustic modes responsible for global spiral density waves; (2) f-modes, which are surface gravity waves; (3) g-modes, gravity modes relying on buoyancy forces resulting from an entropy gradient in the vertical direction; and (4) r-modes, which propagate by inertial forces. Spiral structures in the disk have been investigated mainly in the surface density distribution (e.g., Papaloizou et al. 2007, for a review), where they are incited by p-modes. In contrast, f- and g-modes may also create near-infrared observable spiral structures at the disk surface, but they do not significantly affect the surface density distribution of the disk. In particular, Lubow & Ogilvie (1998) suggest that the f-modes carry most of the torque exerted at the resonance. These modes may induce spiral shock waves observable at the surface of the disk (Boley & Durisen 2006). Zhu et al. (2012b) discussed the importance of g-modes in disks based on their simulations which show three-dimensional structures of spiral waves induced by tidal interaction with a young planet.

Comparisons between near-infrared and millimeter/submillimeter images would give useful constraints for understanding which modes are responsible for spiral arms, and thereby understanding the internal structures in the disk in a manner similar to stellar seismology. The millimeter/submillimeter dust continuum emission from the disk is optically thin (e.g., Williams & Cieza 2011; Espaillat et al. 2014, for reviews), therefore the emission would trace the surface density distribution. Submillimeter images for SAO 206462 and MWC 758 have been obtained by Brown et al. (2009), Isella et al. (2010), and Andrews et al. (2011) with angular resolutions of 03–08 (corresponding to 42–110 AU and 60–160 AU, respectively). These images do not show the presence of spiral structures expected for p-modes. However, these modes cannot be fully ruled out without millimeter/submillimeter images at angular resolutions similar to our near-infrared images (005–01).

Disk holes have been directly observed in dozens of protoplanetary disks, including those in the 2MASS J1604 and PDS 70 rings, with near-infrared and millimeter/submillimeter imaging techniques. Such structures have also been inferred from deficits at certain wavelengths in infrared SEDs (Williams & Cieza 2011; Espaillat et al. 2014, for reviews). As with spiral structures, these holes have been discussed as potential signatures of ongoing planet planet formation and tidal interaction between the disk and planet(s). Such disks, which have been termed "transitional disks," are often assumed to be associated with a vertical wall at the edge of the hole or the gap (Espaillat et al. 2014, for a review). Scattered light at optical and near-infrared wavelengths agrees with the presence of such a wall in GG Tau (Silber et al. 2000; Duchêne et al. 2004) and LkCa 15 (Thalmann et al. 2010). Such a geometry is consistent with three-dimensional simulations by Zhu et al. (2012a) for disk clearing by planetary-mass companions without dust filteration (i.e., without dust and gas accreting through the gap).

In contrast, our analysis of the 2MASS J1604 and PDS 70 rings shows that the observed near-infrared PI distributions can be explained without back scattering at the inner wall. The surface geometries in Figures 6, 7, and 9 do not show a discontinuous boundary at the inner edge of the ring. The disk may still exist at small radii with very low PI fluxes. Millimeter/submillimeter interferometry with high angular resolutions, which would allow us to sufficiently resolve the ring, is necessary to further investigate the nature of ring structures (e.g., Espaillat et al. 2014 for a review).

The 2MASS J1604 ring shows a remarkable deficit in the PI flux distribution at P.A. ∼ 90° and a marginal deficit at P.A. = 180–270° (Mayama et al. 2012, Figure 6). PDS 70 and MWC 480 show asymmetries in the PI distribution with respect to the minor axis (Hashimoto et al. 2012; Dong et al. 2012a; Kusakabe et al. 2012, Figures 7 and 8). Possible explanations for such asymmetries include differences in the scattering geometry, dust mass distribution and dust properties (Mayama et al. 2012; Dong et al. 2012a). Mayama et al. (2012) also discussed the possibility that the deficit at P.A. ∼ 90° ("Dip D" according to Mayama et al. 2012) is associated with an embedded protoplanet.

Some results of our analysis in Section 4.2 are consistent with a nearly axisymmetric surface geometry with nonaxisymmetric illumination. Such illumination has been discussed for several disks, including those associated with HH 30 (e.g., Watson & Stapelfeldt 2007), GG Tau (Silber et al. 2000; Itoh et al. 2002; Krist et al. 2005), HD 163296 (Wisniewski et al. 2008), and RY Tau (Takami et al. 2013). Possible mechanisms for nonaxisymmetric illumination include: (1) hot (or cool) spots on the star; (2) shadowing by a nonaxisymmetric inner disk, including a puffed-up inner rim; (3) obscuration by a companion star; and (4) obscuration by a disk associated with a companion star (Watson & Stapelfeldt 2007; Wisniewski et al. 2008; Takami et al. 2013; Dullemond et al. 2003). In addition to these mechanisms, clumpy accretion of dust and gas onto the star or an inner wind associated with the star may also explain the observed asymmetry, as has been proposed for UX-Ori type variables (see Herbst et al. 1994; Grady et al. 2000 and references therein).

Dip D in the 2MASS J1604 ring is narrow in P.A., and may be explained by obscuration by a nonaxisymmetric distribution of dust in the inner disk or clumpy mass accretion. The same interpretation was made by Silber et al. (2000), Itoh et al. (2002), and Krist et al. (2005) for a similar structure in the scattered light of the GG Tau ring. In contrast, the asymmetric distribution in 2MASS J1604 and MWC 480 may be explained by either hot/cool spots or by the presence of a companion, both of which provide an asymmetry at a wide range of P.A.s. Alternatively, these asymmetried may be caused by self-shadowing by a warped or misaligned inner disk (Whitney et al. 2013).

Synoptic observations of the PI distribution would be extremely useful to discriminate between the above mechanisms. For example, a rotation period of up to 10–20 days would be expected if nonaxisymmetric illumination is due to the non-uniform surface brightness of the star (Herbst et al. 2007, for a review). In contrast, a significantly longer period (e.g., 10 yr or more) would be expected if the nonaxisymmetric illumination is due to shadowing by a nonaxisymmetric inner disk (Krist et al. 2005; Watson & Stapelfeldt 2007; Wisniewski et al. 2008; Takami et al. 2013). The rotational period of the disk depends on radius, therefore measurements of the speed of the nonaxisymmetric patterns would be useful for identifying the radius of the obscuration. Synoptic observations would also allow us to further investigate whether the rotation speed of the pattern agrees with predictions of the spiral density wave theory (Muto et al. 2012).