Estimate on Dust Scale Height from the ALMA Dust Continuum Image of the HD 163296 Protoplanetary Disk

We consider a disk hosting an axisymmetric dust ring whose radial distribution of the dust surface density and the vertical distribution of the dust spatial density are described as Gaussian functions. We take a cylindrical coordinate system as r = 0 at the central star, ϕ = 0 along the major axis, and z = 0 at the midplane. The surface density of the ring at the peak (r = r₀) is Σ₀, the ring width is σ_r , and the scale height is σ_z . The dust spatial density of this object is expressed as

$$\begin{eqnarray}&&\rho (r,z)=\displaystyle \frac{{{\rm{\Sigma }}}_{0}}{\sqrt{2\pi }{\sigma }_{z}}\exp \left(-\displaystyle \frac{{\left(r-{r}_{0}\right)}^{2}}{2{\sigma }_{r}^{2}}-\displaystyle \frac{{z}^{2}}{2{\sigma }_{z}^{2}}\right).\end{eqnarray} \tag{ 1 }$$

We assume the disk inclination angle to be i, and calculate the optical depth of this object along the line of sight. Figure 1 illustrates the disk geometry. In the calculation of the optical depth, we ignore the curvature of the ring. This approximation is appropriate if the ring radius is sufficiently larger than the ring width. The optical depth at r = r₀ + Δr and ϕ is expressed as

$$\begin{eqnarray}\begin{array}{rcl}\tau ({r}_{0}+{\rm{\Delta }}r,\phi ) & = & {\displaystyle \int }_{-\infty }^{\infty }\kappa \rho ({r}_{0}+{\rm{\Delta }}r+z\tan i\sin \phi ,z)\displaystyle \frac{{dz}}{\cos i}\\ & = & \displaystyle \frac{\kappa {{\rm{\Sigma }}}_{0}}{\cos i}\displaystyle \frac{{\sigma }_{r}}{\sqrt{{\sigma }_{r}^{2}+{\sigma }_{z}^{2}{\tan }^{2}i{\sin }^{2}\phi }}\\ & & \times \exp \left(-\displaystyle \frac{1}{2}\displaystyle \frac{{\rm{\Delta }}{r}^{2}}{{\sigma }_{r}^{2}+{\sigma }_{z}^{2}{\tan }^{2}i{\sin }^{2}\phi }\right),\end{array}\end{eqnarray} \tag{ 2 }$$

where κ is the dust opacity.

Zoom In Zoom Out Reset image size

The optical depth depends on the azimuthal angle, as shown in Equation (2). We use the word "ridge" to describe the peak position in the radial direction (r = r₀) at each azimuthal angle. In particular, the optical depth along the ridge of the ring is the maximum on the major axis and the minimum on the minor axis.

We calculate the azimuthal intensity variation of both optically thin and thick disks. Here, we ignore the effects of beam smearing in interferometric observations, which we discuss in Section 2.3.

In the following discussion, we assume that the ring is isothermal in both vertical and radial directions. In reality, there are temperature gradients in both the radial and vertical directions. We assume that the dust ring width is small, so the radial temperature variation can be ignored. In the vertical direction, the temperature is higher in the upper layers that are irradiated by the central star (Chiang & Goldreich 1997; D'Alessio et al. 1998; Dartois et al. 2003; D'Alessio et al. 2005; Rosenfeld et al. 2013). For simplicity, we ignore the vertical temperature gradient. We discuss the effect of the vertical temperature gradient in Appendix A.

For optically thick disks, the intensity can be approximated as

$$\begin{eqnarray}&&{I}_{\nu }({r}_{0},\phi )={B}_{\nu }(T)(1-\exp (-\tau ))\simeq {B}_{\nu }(T),\end{eqnarray} \tag{ 3 }$$

where T is the temperature of the dust ring, and B_ν (T) is the Planck function. In other words, there is no variation in intensity if the ring is optically thick.

For optically thin disks, the intensity can be approximated as

$$\begin{eqnarray}&&{I}_{\nu }({r}_{0},\phi )={B}_{\nu }(T)(1-\exp (-\tau ))\simeq {B}_{\nu }(T)\tau ({r}_{0},\phi ).\end{eqnarray} \tag{ 4 }$$

On the major axis(ϕ = 0°), the intensity on the ridge is

$$\begin{eqnarray}&&{I}_{\nu }({r}_{0},0^\circ )\simeq {B}_{\nu }(T)\displaystyle \frac{\kappa {{\rm{\Sigma }}}_{0}}{\cos i},\end{eqnarray} \tag{ 5 }$$

and on the minor axis (ϕ = 90°), the intensity is

$$\begin{eqnarray}&&{I}_{\nu }({r}_{0},90^\circ )\simeq {B}_{\nu }(T)\displaystyle \frac{\kappa {{\rm{\Sigma }}}_{0}{\sigma }_{r}}{\sqrt{{\sigma }_{r}^{2}{\cos }^{2}i+{\sigma }_{z}^{2}{\sin }^{2}i}}.\end{eqnarray} \tag{ 6 }$$

The ratio of the intensity on the ridge on the major axis to that on the minor axis is

$$\begin{eqnarray}&&\displaystyle \frac{{I}_{\nu }({r}_{0},90^\circ )}{{I}_{\nu }({r}_{0},0^\circ )}=\displaystyle \frac{{\sigma }_{r}}{\sqrt{{\sigma }_{r}^{2}+{\sigma }_{z}^{2}{\tan }^{2}i}}.\end{eqnarray} \tag{ 7 }$$

Thus, there is a difference in intensity between the major and minor axes.

For the optically thin case, we further discuss the difference of the intensity between the major and minor axes. As shown in Equation (7), the projected dust scale height, ${\sigma }_{z}\sin i$, relative to the projected ring width, ${\sigma }_{r}\cos i$, is observed as the difference in the intensity between the major and minor axes. If ${\sigma }_{r}\cos i\gg {\sigma }_{z}\sin i$, there is no difference in the intensity between the major and minor axes. On the other hand, if the projected ring width, ${\sigma }_{r}\cos i$, is not much wider than the projected dust scale height, ${\sigma }_{z}\sin i$, the difference in the intensity between the major and minor axes depends on the dust scale height.

Figure 2 shows images of radiative transfer simulations of protoplanetary disk models composed of a dust ring with different optical and geometric thicknesses. We assume that all rings have the same width (σ = 5 au) and distance from the center (r₀ = 100 au). We assume that the rings are isothermal, and we set the temperature such that the brightness temperature at the peak position on the major axis is 10 K. In the left panel, the optical depth at the ridge of the ring is κΣ₀ = 5, and the temperature is 15.0 K. In the middle and right panels, the optical depth at the ridge of the ring is κΣ₀ = 0.01 and the temperature is 717.8K. The scale height in the middle panel is 5/16 au, and that in the left and right panels is 5 au. Figure 3 plots the intensity along the ridges of the simulated images shown in Figure 2 and Equation (12) from the analytical calculations. We set the inclination to be 45^◦. We can see that for the optically thick (the left panel) or geometrically thin (the middle panel) cases, the intensity does not depend on the azimuthal direction, but for the optically thin and geometrically thick case (the right panel), the intensity on the minor axis is smaller than that on the major axis. We confirm that Equation (12) describes the simulation results well. From the above, we conclude that the intensity along the ridges depends on the azimuthal angle if the ring is optically thin, the dust scale height is not too wide compared to the dust ring width, and the inclination is not too small.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Recent studies have pointed out that optically thick emission can be fainter than the Planck function if the albedo is high (Miyake & Nakagawa 1993; Liu 2019; Zhu et al. 2019; Ueda et al. 2020; Lin et al. 2020). If the scattering works, one can imagine that the azimuthal variation may disappear due to the scattering-induced intensity reduction. This does not change our conclusion if the ring is optically thin (τ_abs ≪ 1) because the scattering effects are not significant. Moreover, this does not change if the ring is extremely optically thick (τ_abs ≫ 1) because we do not expect any azimuthal variation even if we do not consider the scattering effects. We should be careful if the ring is marginally optically thick (τ_abs ∼ 1). In this case, if we do not consider the effects of scattering, we expect an azimuthal intensity variation to some extent. However, if scattering plays a role the intensity on the major axis decreases more than that on the minor axis. As a result, the azimuthal variation is suppressed compared with the nonscattering case, and it would be observed as if it were extremely optically thick. Therefore, the scattering effects may diminish the azimuthal variation of marginally optically thick rings. Note that the scattering may not significantly contribute to the emission in the case of HD 163296 because τ_ext = τ_abs + τ_scat < 1 from the observation of CO backside emission (Isella et al. 2018).

We discuss the conditions in which we can constrain the dust scale heights, including the effects of the beam smearing of the interferometric observations. We approximate the beam smearing as a Gaussian convolution of an ellipsoidal shape. Based on the previous section, we only consider the optically thin case because we are not able to constrain the dust scale height if the dust ring is optically thick. From Equation (2), the intensity at the location (r + Δr, ϕ) can be expressed as

$$\begin{eqnarray}\begin{array}{rcl}{I}_{\nu }({r}_{0}+{\rm{\Delta }}r,\phi ) & = & {B}_{\nu }(T)\displaystyle \frac{\kappa {{\rm{\Sigma }}}_{0}}{\cos i}\displaystyle \frac{{\sigma }_{r}}{\sqrt{{\sigma }_{r}^{2}+{\sigma }_{z}^{2}{\tan }^{2}i{\sin }^{2}\phi }}\\ & & \times \exp \left(-\displaystyle \frac{1}{2}\displaystyle \frac{{\rm{\Delta }}{r}^{2}}{{\sigma }_{r}^{2}+{\sigma }_{z}^{2}{\tan }^{2}i{\sin }^{2}\phi }\right).\end{array}\end{eqnarray} \tag{ 8 }$$

We convolve this equation with a Gaussian beam. We deproject the observed Gaussian beam for the inclination and denote the standard deviation of the deprojected beam along the ϕ direction as σ_beam(ϕ). Here we ignore the curvature of the ring and convolve the radial profile and the beam smearing. The convolved intensity is

$$\begin{eqnarray}&&\begin{array}{l}{I}_{\nu }({r}_{0}+{\rm{\Delta }}r,\phi )={B}_{\nu }(T)\displaystyle \frac{\kappa {{\rm{\Sigma }}}_{0}}{\cos i}\\ \displaystyle \frac{{\sigma }_{r}}{\sqrt{{\sigma }_{r}^{2}+{\sigma }_{z}^{2}{\tan }^{2}i{\sin }^{2}\phi +{\sigma }_{\mathrm{beam}}{\left(\phi \right)}^{2}}}\\ \times \exp \left(-\displaystyle \frac{1}{2}\displaystyle \frac{{\rm{\Delta }}{r}^{2}}{{\sigma }_{r}^{2}+{\sigma }_{z}^{2}{\tan }^{2}i{\sin }^{2}\phi +{\sigma }_{\mathrm{beam}}{\left(\phi \right)}^{2}}\right).\end{array}\end{eqnarray} \tag{ 9 }$$

Particularly, the intensities along the major axis and the minor axis are expressed as

$$\begin{eqnarray}\begin{array}{rcl}{I}_{\nu }({r}_{0}+{\rm{\Delta }}r,0^\circ ) & = & {B}_{\nu }(T)\displaystyle \frac{\kappa {{\rm{\Sigma }}}_{0}}{\cos i}\displaystyle \frac{{\sigma }_{r}}{\sqrt{{\sigma }_{r}^{2}+{\sigma }_{\mathrm{beam}}{\left(0^\circ \right)}^{2}}}\\ & & \times \exp \left(-\displaystyle \frac{1}{2}\displaystyle \frac{{\rm{\Delta }}{r}^{2}}{{\sigma }_{r}^{2}+{\sigma }_{\mathrm{beam}}{\left(0^\circ \right)}^{2}}\right),\end{array}\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}\begin{array}{rcl}{I}_{\nu }({r}_{0}+{\rm{\Delta }}r,90^\circ ) & = & {B}_{\nu }(T)\displaystyle \frac{\kappa {{\rm{\Sigma }}}_{0}}{\cos i}\\ & & \displaystyle \frac{{\sigma }_{r}}{\sqrt{{\sigma }_{r}^{2}+{\sigma }_{z}^{2}{\tan }^{2}i+{\sigma }_{\mathrm{beam}}{\left(90^\circ \right)}^{2}}}\\ & & \times \exp \left(-\displaystyle \frac{1}{2}\displaystyle \frac{{\rm{\Delta }}{r}^{2}}{{\sigma }_{r}^{2}+{\sigma }_{z}^{2}{\tan }^{2}i+{\sigma }_{\mathrm{beam}}{\left(90^\circ \right)}^{2}}\right).\end{array}\end{eqnarray} \tag{ 11 }$$

Along the ridge of the ring, the intensity is expressed as

$$\begin{eqnarray}\begin{array}{rcl}{I}_{\nu }({r}_{0},\phi ) & = & {B}_{\nu }(T)\displaystyle \frac{\kappa {{\rm{\Sigma }}}_{0}}{\cos i}\\ & & \displaystyle \frac{{\sigma }_{r}}{\sqrt{{\sigma }_{r}^{2}+{\sigma }_{z}^{2}{\tan }^{2}i{\sin }^{2}\phi +{\sigma }_{\mathrm{beam}}{\left(\phi \right)}^{2}}}.\end{array}\end{eqnarray} \tag{ 12 }$$

From Equation (10), the intensity along the major axis is independent of the dust scale height. Therefore, we can construct a disk model that is independent of the dust scale height from the observed intensity along the major axis. On the other hand, from Equations (11) and (12), the intensity along the ridge except for the major axis or along the minor axis depends on the dust scale height. Therefore, we can constrain the dust scale height by the consistency of the observed intensity along the ridge or the minor axis.

We discuss the conditions for constraining the dust scale height from the observed images. The dust scale height contributes to the intensity as the term ${\sigma }_{z}^{2}{\tan }^{2}i$ in the denominator of Equations (11) and (12). If ${\sigma }_{z}^{2}{\tan }^{2}i$ is not too small compared to the other terms (${\sigma }_{r}^{2}$, ${\sigma }_{\mathrm{beam}}^{2}$), the effect of the dust scale height appears in the observed intensity. Then, we can constrain the dust scale height from observations. We further discuss the effect of the beam smearing in Appendix B.

To summarize the discussion above, we can constrain the dust scale height in the following steps. First, we construct a disk surface density model from the observed intensity along the major axis. Then, we search for the dust scale height that reproduces the observed intensity along the ridge or the minor axis. We require the following four conditions to constrain dust scale heights: 1. The dust ring is optically thin, 2. the dust ring width is not much wider than the dust scale height, 3. the inclination is not too small, and 4. the dust ring is spatially resolved.

In this study, we use the fits data of the image of the HD 163296 disk taken in the DSHARP campaign (2016.1.00484.L), which is one of the ALMA large programs (Andrews et al. 2018). We show the image in Figure 4. This image was taken with ALMA Band 6, which corresponds to the wavelengths of λ = 1.25 mm. This image was made by combining the data of project 2013.1.00366.S13 (Flaherty et al. 2015) and project 2013.1.00601.S4 (Isella et al. 2016) to improve the sensitivity and uv-coverage on shorter baselines. The details of this image are shown in Andrews et al. (2018) and Isella et al. (2018). We assume that the distance to this object is 101.5 pc from Gaia DR2 parallax measurements (Gaia Collaboration et al. 2018). We assume that the central star mass M_star = 2.02M_sun is a dynamical mass based on the line observation from Teague et al. (2019).

This object satisfies the requirements for constraining the dust scale height described in Section 2.3. The observation's beam size is 0020 × 0016 with a standard deviation that corresponds to 2.1 × 1.6 au, and the beam position angle is 817 (Isella et al. 2018). This source has narrow rings; the width is 4.0 au for the inner ring and 3.9 au for the outer ring with a standard deviation, which we discuss in Section 3.2. The disk inclination is 467 ± 01, and the disk position angle is 1331 ± 01 (Isella et al. 2018). The gas scale heights of this source are 3.8 au in the inner ring and 6.7 au in the outer ring, using the temperature model discussed in Section 3.2. As described above, this object is sufficiently inclined, and the ring width and beam are not too large compared to the gas scale height, so we think that we can constrain the dust scale height of this source.

We describe our simulation setup to compare with the observation. We perform radiative transfer simulations with RADMC-3D (Dullemond et al. 2012). In our simulation, we use spherical coordinates (r, θ, ϕ), and we assume axial symmetry in ϕ direction and plane symmetry with respect to the midplane. The radial grid is linearly spaced between r_in = 10 au and r_out = 200 au with 512 grids. The θ grid is linearly spaced with 512 grids in the range 0 ≤ θ ≤ 0.2, where θ = 0 is the midplane, and θ = π/2 is the north pole.

We adopt the following temperature profile with the smooth power-law distribution derived by Dullemond et al. (2020),

$$\begin{eqnarray}&&T(r)=18.7\ [{\rm{K}}]{\left(\displaystyle \frac{r}{400[\mathrm{au}]}\right)}^{-0.14}.\end{eqnarray} \tag{ 13 }$$

Dullemond et al. (2020) discussed the midplane temperature based on the CO line emission. We assume that the disk is isothermal in the vertical direction. In reality, the disk surface is hotter than the midplane, but the vertical distribution of the gas is approximated well by hydrodynamic equilibrium at midplane temperature, so using the midplane temperature is a good approximation when we consider the vertical distribution of the gas. We discuss the effect of the vertical temperature gradient on the intensity in Appendix A.

We use the dust opacity model of Birnstiel et al. (2018) with κ_abs = 0.484 cm² g⁻¹ at a wavelength of λ = 1.25 mm. Here, we consider only the absorption of the dust and do not include dust scattering (see Miyake & Nakagawa 1993; Kataoka et al. 2014; Liu 2019; Zhu et al. 2019; Ueda et al. 2020; Lin et al. 2020). We need the dust surface density and dust opacity for the input parameters of RADMC-3D. From the observed intensity and temperature distribution model of Equation (13), we can derive the optical depth but cannot solve the degeneracy of the dust surface density and dust opacity, so we need a dust opacity model to derive the dust surface density. However, output images depend on the optical depth regardless of the combination of the dust opacity and surface density. Therefore, our result does not depend on the dust opacity model.

We make a surface density profile based on the observational radial intensity profile along the major axis. This is possible because the intensity along the major axis does not depend on the dust scale height, as shown in Equation (10). We use the radial intensity profile along the major axis in the northwest direction to avoid the additional crescent-shaped structure reported in the southeast direction (Isella et al. 2018). We perform numerical simulations to search for parameters that reproduce the observational intensity along the major axis. We assume that the dust surface density model is composed of an inner disk and two Gaussian rings as

$$\begin{eqnarray}\begin{array}{rcl}{{\rm{\Sigma }}}_{\mathrm{dust}}(r) & = & {{\rm{\Sigma }}}_{\mathrm{disk}}\exp [-{\left(r/{r}_{\mathrm{disk}}\right)}^{5}]\\ & & +{{\rm{\Sigma }}}_{\mathrm{ring}1}\exp [-{\left(r-{r}_{0,\mathrm{ring}1}\right)}^{2}/2{w}_{\mathrm{ring}1}^{2}]\\ & & +{{\rm{\Sigma }}}_{\mathrm{ring}2}\exp [-{\left(r-{r}_{0,\mathrm{ring}2}\right)}^{2}/2{w}_{\mathrm{ring}2}^{2}].\end{array}\end{eqnarray} \tag{ 14 }$$

We use the parameters Σ_disk = 20.66 g cm⁻², r_disk = 25.0 au, Σ_ring1 = 1.26 g cm⁻², r_0,ring1 = 67.9 au, w_ring1 = 4.0 au, Σ_ring2 = 0.58 g cm⁻², r_0,ring2 = 100.5 au, and w_ring2 = 3.9 au. In Section 3.4 we discuss the consistency of this model and observation. We summarize the parameters related to the rings on which we focus in this study in Table 1.

The absorption optical depth (${\tau }_{\mathrm{abs}}={\kappa }_{\mathrm{abs}}{{\rm{\Sigma }}}_{d}/\cos i$) of this model at the rings is 0.89 for the inner ring and 0.41 for the outer ring. Isella et al. (2018) estimated the extinction optical depth, τ_ext, of this source from the extinction of CO backside emission at the rings, and constrained that τ_ext = 0.64 ± 0.05 for the inner ring and τ_ext = 0.74 ± 0.05 for the outer ring. Since κ_ext = κ_abs + κ_scat, τ_ext must be higher than τ_abs. However, τ_abs in this study and the constraint of τ_ext in Isella et al. (2018) do not satisfy this relation for the inner ring. A possible explanation for this inconsistency is that the dust surface density is overestimated because the temperature model used in this study underestimates the temperature, or τ_ext in Isella et al. (2018) is underestimated due to beam smearing or contamination of frontside emission.

Now we are ready to determine the dust scale height. To estimate the scale height that reproduces the observation, we change the dust settling parameter f_set. We define the dust settling parameter such that f_set = h_gas/h_dust, where h_dust is dust scale height and h_gas is the gas scale height. The gas pressure scale height is h_gas = c_s /Ω_kep from the hydrostatic equilibrium, where c_s is the isothermal sound speed and Ω_kep is the Keplerian orbital velocity. We run multiple simulations in the range of 1/2 ≤ f_set ≤ 16 with every 2^1/16 steps, which corresponds to 81 models. We convolve the simulated images with the beam and compare them with the observed image.

The analytical expression of the azimuthal variation described as Equation (12) shows that the intensity along the ridge depends on the dust scale height. We constrain the dust scale height by comparing the azimuthal intensity variation along the ridge of the observation with that of the simulations.

To demonstrate the azimuthal intensity variation, we show in Figure 5 images of the observation and simulations in the case of f_set = 1 (no settling) and 8 (settled). We use a different color range for each row to emphasize the azimuthal variations of the inner and outer rings. Figure 6 shows the azimuthal variation along the ridges of the observation and simulations in the case of f_set = 1, 2, and 8. As shown in the simulation images of f_set = 1, where the dust is flared, the intensity along the ridge is higher on the major axis (ϕ = 0°, 180°) than on the minor axis (ϕ = 90°, 270°). In contrast, as shown in the simulation images of f_set = 8, where the dust is settled, the intensity along the ridge does not depend on the azimuthal angle. The observational images show that the intensity along the ridge of the inner ring is lower on the minor axis (ϕ = 90°, 270°) than on the major axis (ϕ = 0°, 180°). In contrast, the intensity along the ridge of the outer ring is uniform for the azimuthal direction. This means that the dust scale height is large for the inner ring and small for the outer ring.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We constrain f_set using the χ². We run simulations varying f_set from 1/2 to 16 by a factor of 2^1/16, and calculate χ². To estimate the parameter f_set, we calculate the χ² between the observation and the simulations along the ridge. To exclude the effect of the crescent-shaped substructure in the southeast direction, we excluded the region between − 45° and 45^◦ in our calculations of the χ². Figure 7 shows the χ² at each f_set. The best-fit parameters are f_set = 1.1 for the inner ring and f_set = 16 for the outer ring. We search the range where ${\chi }^{2}\lt {\chi }_{{\rm{\min }}}^{2}+{3}^{2}$. We obtain ${f}_{{\rm{set}}}={1.1}_{-0.1}^{+0.1}$ for the inner ring and f_set > 9.5 for the outer ring with 3σ uncertainty. With a 3σ uncertainty, the inner ring is more flared than the outer ring. Figure 7 also shows that our method has high sensitivity where ${f}_{{\rm{set}}}$ is small while it has low sensitivity where ${f}_{{\rm{set}}}$ is large. Note that ${f}_{{\rm{set}}}$ cannot be smaller than unity because the dust scale height cannot be larger than the gas scale height, while our results satisfy this condition.

Zoom In Zoom Out Reset image size

We discuss whether the derived dust scale height from the azimuthal intensity variation is consistent with the observed radial intensity profile along the minor axis. From Equations (10) and (11), the intensity along the major axis does not depend on the dust scale height, while that along the minor axis depends on the dust scale height. Along the minor axis, the higher the dust scale height, the wider the observed dust ring width and the lower the observed intensity at the peak. We plot the intensity of both simulations and the observation along the major and minor axes to check the consistency of f_set, which we estimate in Section 3.3.

Figure 8 shows the radial intensity profiles of the observation along the major and minor axes. Around the inner ring at 68 au, the intensity along the major axis is higher than that along the minor axis. On the other hand, around the outer ring at 100 au, the intensity along the major and minor axes is almost the same.

Zoom In Zoom Out Reset image size

Figure 9 shows the radial intensity profiles along the major (top) and minor (bottom) axes of the observation and simulations in case of f_set = 1 and 8. We can see that the simulations reproduces the observation well along the major axis from the upper panel of Figure 9. The bottom panel of Figure 9 shows that the radial profile along the minor axis around the inner ring is well reproduced by f_set = 1, while that around the outer ring is well reproduced by f_set = 8. This result is consistent with the result in Section 3.3 that the inner ring is flared and the outer ring is settled.

Zoom In Zoom Out Reset image size

Figure 10 shows the radial intensity profile of the best model along the minor axis. Based on the result in Section 3.3, we set f_set = 1.1 for the inner ring and f_set = 16 for the outer ring. This result reproduces the onservation well for the entire region.

Zoom In Zoom Out Reset image size

We discuss the dust scale height for the entire region of the protoplanetary disk. Figure 11 shows a schematic illustration of the dust scale height of HD 163296. Ohashi & Kataoka (2019) determine the dust scale height at the gaps by using the polarized intensity of this object. They concluded that the dust scale height is less than one-third of the gas scale height at the gap at 48 au, and the dust scale height is two-thirds of the gas scale height at the gap at 86 au. Taken together with the results of this study, we can draw the entire picture of the radial distribution of the dust scale height except for the central disk. From the inner side to the outer side, the scale height is unknown at the central disk, thin at the 48 au gap, thick at the 68 au ring and the 87 au gap, and thin at the 100 au ring.

Zoom In Zoom Out Reset image size

We discuss whether the radial variation of the dust scale height can be explained by a simple dead-zone model (Gammie 1996; Okuzumi et al. 2016; Ueda et al. 2019). The dead zone exists in the inner part of the protoplanetary disk because the low-ionization ratio in the inner part of the disk suppresses the magnetorotational instability. The dust is settled in the dead zone due to weak gas turbulence. Therefore, the dust scale height is smaller in the dead zone and larger outside the dead zone. As shown in Figure 11, however, the dust scale height of HD 163296 varies in a complex way in the radial direction and cannot be explained by a simple disk model with a dead zone. This suggests that other factors cause the change in dust scale height, e.g., the radial change of the dust size or the turbulence.

We compare our dust scale height model with an infrared polarization observation. Muro-Arena et al. (2018) made a 3D disk model of HD 163296 using millimeter dust continuum observation using ALMA and infrared polarized scattered light using SPHERE. They reported that the polarized infrared light is observed only at the inner ring at 68 au and is not observed at the outer ring at 100 au. This result can be interpreted to mean that the outer ring is in the shadow of the inner ring, which is consistent with our results of the dust scale height.

We estimate α/St from the dust scale height. By assuming a balance between vertical settling and turbulent diffusion, the dust scale height is written as (e.g., Dubrulle et al. 1995; Youdin & Lithwick 2007)

$$\begin{eqnarray}&&{h}_{d}={\left(1+\displaystyle \frac{\mathrm{St}}{\alpha }\displaystyle \frac{1+2\mathrm{St}}{1+\mathrm{St}}\right)}^{-1/2}{h}_{g},\end{eqnarray} \tag{ 15 }$$

where h is the scale height, and the subscripts g and d indicate gas and dust, respectively. This equation can be approximated as

$$\begin{eqnarray}&&{h}_{d}={\left(1+\displaystyle \frac{\mathrm{St}}{\alpha }\right)}^{-1/2}{h}_{g},\end{eqnarray} \tag{ 16 }$$

if St ≪ 1. Using the f_set constrained in Section 3.3, we derive α/St > 2.4 for the inner ring and α/St < $1.1\times {10}^{-2}$ for the outer ring with 3σ uncertainties. We note that Equation (16) may not be appropriate if ${h}_{d}\approx {h}_{g}$.

We discuss the consistency with measurements in other studies. Dullemond et al. (2018) and Rosotti et al. (2020) estimated α/St by the radial diffusion of dust due to the gas turbulence. Dullemond et al. (2018) assumed ring formation by dust traps at gas bumps and derived the following relationship between dust ring width w_d and gas ring width w_g :

$$\begin{eqnarray}&&{w}_{d}=\sqrt{\displaystyle \frac{\alpha }{\mathrm{St}}}{w}_{g}.\end{eqnarray} \tag{ 17 }$$

Dullemond et al. (2018) estimated the dust ring width from the dust continuum emission, and Rosotti et al. (2020) estimated the gas ring width from the deviation of the gas rotation velocity. Their results are α/St = 0.23 for the inner ring and α/St = 0.04 for the outer ring. We note that Rosotti et al. (2020) and our results show the same tendency that the inner ring is more diffuse than the outer ring. However, their $\alpha /{\rm{St}}$ is somewhat different from that of us. There are two possible reasons for this disagreement. One is that the turbulence is not isotropic, i.e., the turbulence is different in radial and vertical directions. The other is that the ring formation mechanism is different from that assumed by Rosotti et al. (2020). Rosotti et al. (2020) assume that ring formation is due to gas bumps. In other scenarios, however, the dust ring width and gas ring width do not necessarily satisfy the relationship of Equation (17). We discuss other ring formation mechanisms in Section 4.5.

In this section, we constrain the dust size based on the α/St. In Section 4.2 we obtained α/St based on f_set. The Stokes number is expressed as

$$\begin{eqnarray}&&\mathrm{St}=\displaystyle \frac{\pi {\rho }_{\mathrm{mat}}{a}_{\mathrm{dust}}}{2{{\rm{\Sigma }}}_{\mathrm{gas}}},\end{eqnarray} \tag{ 18 }$$

where ρ_mat is the dust material density, a_dust is the dust radius, and Σ_gas is the gas surface density. Therefore, we can constrain the dust size from α/St obtained in Section 4.2 and α, Σ_gas, and ρ_mat. We assume that the dust is compact and icy (ρ_mat = 1.0 g cm⁻³), and the gas surface density is Σ_gas = 21.1 g cm⁻² for the inner ring and Σ_gas = 10.0 g cm⁻² for the outer ring, following Rab et al. (2020) based on molecular line observations. Flaherty et al. (2017) placed an upper limit of α < 3 × 10⁻³, so we calculate the dust radius for the typical three turbulences as α = 3 × 10⁻³, 1 × 10⁻³, and 1 × 10⁻⁴. We show the results in Table 2. If we assume α = 1 × 10⁻³, we obtain a_dust < 5.6 × 10⁻³ cm for the inner ring and a_dust > 5.7 × 10⁻¹ cm for the outer ring.

Table 2. Estimation of the Dust Size at Typical Gas Turbulence

Ring	fset	α/St	Σgas	α	adust
			(g cm−2 )		(cm)
(1)	(2)	(3)	(4)	(5)	(6)
				3 × 10−3	<1.7 × 10−2
inner ring	<${1.1}_{-0.1}^{+0.1}$	>2.4	21.0	1 × 10−3	<5.6 × 10−3
				1 × 10−4	<5.6 × 10−4
				3 × 10−3	>1.7
outer ring	>9.5	<1.1 × 10−2	10.0	1 × 10−3	>5.7 × 10−1
				1 × 10−4	>5.7 × 10−2

Note. Column (2): Settling parameters that we constrain in Section 3.3. Column (3): α/St that we constrain in Section 4.2. Column (4): Gas surface density at the peak position of the ring (Rab et al. 2020). Column (5): Gas turbulence. We assume three typical parameters based on the upper limit from Flaherty et al. (2017). Columns (6): Dust size calculated from the left parameters.

Download table as:  ASCIITypeset image

Since Flaherty et al. (2017) placed an upper limit of α < 3 × 10⁻³, we can constrain the dust size of the inner ring without any assumption of turbulence. As shown in Table 2, the dust size in the inner ring is a_dust < 1.7 × 10⁻²cm. The result rules out models in which the midplane dust consists of centimeter-sized pebbles. Note that Flaherty et al. (2017) assumed that α is constant throughout the disk, but this may not be the case.

If the turbulence is the same in the two rings, we find that the dust size is larger in the outer ring than in the inner ring. Guidi et al. (2016) found a trend of the dust size from the spectral index. They concluded that the dust size decreases toward the outside from the increasing trend of the spectral index, while the spatial resolution is not high enough to resolve the rings. Dent et al. (2019) obtained the specrtal index, α_SED, between Band 6 and Band 7 (0.87 mm) from high-resolution observations that spatially resolve the ring. Dent et al. (2019) obtained that α_SED in the rings is smaller than that in the gaps, and α_SED in the inner ring (α_SED = 2.2) is larger than that in the outer ring (α_SED = 1.8). Since the systematic uncertainty is large, α_SED for the total disk flux differs among Pinilla et al. (2014; 2.73 ± 0.44), Guidi et al. (2016) and Notsu et al. (2019; 2.7), and Dent et al. (2019; 2.1 ± 0.3). The value of α_SED itself may not be accurate, but we can consider that the radial trend of α_SED is correct. If we assume that this disk is optically thin, the dust is smaller in the inner ring than in the outer ring. This result is consistent with our result assuming the same turbulence between the two rings. We do not rule out the possibility that the gas turbulence is different between the inner and outer rings. We note that if the outer ring is optically thick, α_SED becomes smaller in the outer ring, and the azimuthal variation in intensity is small.

Ohashi & Kataoka (2019) estimated the maximum dust size at the dust gap to be 140 μm based on the polarization fraction of ALMA Band 7 (0.87 mm) dust continuum. This result is smaller than the dust size of the outer ring shown in Table 2. The dust size estimated from the polarization observation is just the dust size in the gap, and we cannot directly compare it with the dust size in the ring. As shown in Table 2, however, it is consistent to assume that the dust size in the inner ring is the same as that in the gaps.

In this section, we discuss the possible range of α and St. Since our method cannot solve for the degeneracy of α and St from α/St, we estimate α under the following two assumptions: one assumption is that the dust size is limited by the collisional fragmentation due to the turbulence, and the other assumption is that the dust size is fixed to 1 mm.

First, we consider the case where the dust size is limited by the collisional fragmentation of dust grains induced by gas turbulence. In this case, α can be written as (Birnstiel et al. 2012; Rosotti et al. 2020)

$$\begin{eqnarray}&&{\alpha }_{\mathrm{frag}}=\sqrt{\displaystyle \frac{1}{3}\displaystyle \frac{{v}_{\mathrm{frag}}^{2}}{{c}_{{\rm{s}}}^{2}}\left(\displaystyle \frac{\alpha }{\mathrm{St}}\right)},\end{eqnarray} \tag{ 19 }$$

where v_frag is the fragmentation velocity and c_s is the sound speed isothermal. We obtain α_frag > 3.0 × 10⁻² for the inner ring and α_frag < 2.1 × 10⁻³ for the outer ring when we assume the fragmentation velocity to be v_frag = 10 m/s. The dust size may not be determined by the fragmentation limit by the turbulence, but if the turbulence is higher than α_frag, the dust is destroyed and becomes smaller. Therefore, α_frag is the upper limit of the turbulence.

Next, we consider the turbulence when the dust size is 1 mm. As we discuss in Section 4.3, the dust size may not be the same for the two rings, but as an example, we calculate the turbulence α_{1 mm} when the dust size a_dust = 1 mm for the two rings. The St is determined as Equation (18), and we assume that the gas surface density and dust material density are the same as those in Section 4.3. Under these assumptions, we obtain α_{1 mm} > 1.8 × 10⁻² for the inner ring and α_{1 mm} < 1.8 × 10⁻⁴ for the outer ring.

We summarize the results in Table 3 and Figure 12. The constraint of α in this study at the inner ring is larger than the observational upper limit of α in Flaherty et al. (2017). There are three possible explanations for this inconsistency. First, there may be a mechanism that keeps the dust size small in the inner ring. For example, the fragmentation velocity may be lower than the assumption in this study due to dust sintering (Okuzumi et al. 2016). Second, the turbulence may vary in the radial direction. Flaherty et al. (2017) constrained α throughout the disk from gas emission line observations, and this value corresponds to the intensity-weighted average of α. Since the gas disk is larger than the dust ring, low turbulence in the outer side results in a small α throughout the disk. Third, the α of the dust diffusion we estimate may be different from the α of the gas motion measured by Flaherty et al. (2017). Pinilla et al. (2021) discuss the case where α for gas evolution, turbulent velocity, vertical diffusion, and radial diffusion takes different values.

Zoom In Zoom Out Reset image size

Liu et al. (2018) constrained α from hydrodynamic simulations assuming gap opening by planets. They showed that if α is small ( ≤ 10⁻⁴) in the inner disk and large ( ≥ 7.5 × 10⁻³) in the outer disk, they can reproduce both dust and gas observations. However, our study shows that α/St is larger at the inner ring than at the outer ring, which is different from Liu et al. (2018). This disagreement indicates that it is difficult to reproduce the gas and dust profiles with the dust trap alone. The ring formation mechanism is discussed in Section 4.5.

In this section, we discuss the ring formation process based on the dust scale height that we estimate in this study. Previously, several ring formation mechanisms have been proposed, and we consider the following four mechanisms: 1. Dust accumulation by gas gaps induced by planets, 2. enhanced fragmentation of dust grains by sintering at snowlines, 3. dust accumulation at the outer edge of the dead zone, and 4. secular gravitational instability (SGI). This study shows that the settling conditions of the two rings are different. Therefore, we consider ring formation mechanisms of each ring independently.

A possible mechanism is dust accumulation by the pressure gaps induced by planets (Lin & Papaloizou 1979; Goldreich & Tremaine 1980; Duffell & MacFadyen 2013; Fung et al. 2014; Kanagawa et al. 2015). Zhu et al. (2012) show that large dust accumulates on the outside of the pressure gap, while small dust couples with gas and flows into the inside of the pressure gap. This process is called dust filtration. Since the rings formed by this process consist of large dust, the dust is settled in the midplane. Compared to the results of our study, the outer ring is consistent with this process because the dust is settled, while the inner ring is inconsistent with this process because the dust is flared.

Another possible mechanism is dust sintering (Okuzumi et al. 2016; Zhang et al. 2015). The sintering theory predicts that the rings consist of small dust fragments. The inner ring may be explained by dust sintering because the dust is flared, but the outer ring cannot be explained because the dust is settled, which suggests large dust grains.

Another possible mechanism is dust accumulation at the outer edge of the dead zone (Flock et al. 2015; Pinilla et al. 2016; Mori & Okuzumi 2016; Ueda et al. 2019). The dead zone is the inner region of the disk where the MRI is stabilized and the gas turbulence is weak. However, since the dust in the inner ring is not settled, the turbulence is not weak, and the dead zone is not the formation mechanism of the ring. The outer ring is far from the central star when we consider the dead zone. We do not rule out the possibility that the central disk is due to the dead zone.

The last possible mechanism we consider is SGI (Ward 2000; Youdin 2011; Michikoshi et al. 2012; Takahashi & Inutsuka 2014; Tominaga et al. 2019). The SGI is stabilized by the turbulence, so the turbulence must be weak for the instability to develop. Therefore, the dust is settled in the dust ring that is formed by the SGI. From the perspective of the dust scale height, the outer ring is consistent, and the inner ring is inconsistent.

To summarize the above, from the perspective of the dust scale height, the formation mechanism of the inner ring is consistent with dust sintering. In contrast, the formation mechanism of the outer ring is consistent with dust accumulation by gas gaps induced by planets and secular gravitational instability. We do not rule out other ring-forming mechanisms, such as dust replenishment from a planet around the substructure.

We assume the temperature to be isothermal in the vertical direction and a simple power law in the radial direction, but this may be different from the actual temperature.

In this study, we focus on only one object, HD 163296. This method can be applied to other objects to estimate the dust scale height. By using high spatial resolution observations, it is possible to limit the dust scale height of other objects.