Radial Variations in Grain Sizes and Dust Scale Heights in the Protoplanetary Disk around HD 163296 Revealed by ALMA Polarization Observations

The disk of HD 163296 shows ring and gap substructures in observations with the Atacama Large Millimeter/submillimeter Array. In addition, this is the only disk where the rings and gaps are spatially resolved in millimeter-wave polarization measurements. In this paper, we conduct radiative transfer modeling that includes self-scattering polarization to constrain the grain size and its distribution. We found that the grain size and dust scale height are the key parameters for reproducing the radial and azimuthal distributions of the observed polarization signature. Radial variation is mainly determined by grain size. The polarization fraction is high if the particle size is ∼λ/2π; it is low if the particle size is larger or smaller than this. In contrast, azimuthal variation in polarization is enhanced if the dust scale height is increased. Based on detailed modeling of the polarization of HD 163296, we found the following radial variations in the grain size and dust scale height. The maximum grain size was 140 μm in the gaps and significantly larger or smaller in the rings. The dust scale height is less than one-third of the gas scale height inside the 70 au ring, and two-thirds of it outside. Furthermore, we constrained the gas turbulence to be α ≲ 1.5 × 10−3 in the 50 au gap and α ∼ 0.015–0.3 in the 90 au gap. The transition in the turbulence strength at the boundary of the 70 au ring indicates the existence of a dead zone.


Introduction
Millimeter and submillimeter polarization has been detected from more than 10 protoplanetary disks thanks to the high spatial resolution and high sensitivity of the Atacama Large Millimeter/ submillimeter Array (ALMA) (e.g., Kataoka et al. 2016a;Stephens et al. 2017;Alves et al. 2018;Cox et al. 2018;Harris et al. 2018;Hull et al. 2018;Lee et al. 2018;Sadavoy et al. 2018). However, the mechanism causing the polarization is not a straightforward extension of that in star-forming regions, where the polarization is believed to be emitted from elongated dust grains aligned with magnetic fields (e.g., Lazarian 2007;Andersson et al. 2015). This polarization is used as a tracer of magnetic field directions (e.g., Rao et al. 1998;Lai et al. 2001;Girart et al. 2006;Stephens et al. 2013;Rao et al. 2014;Cox et al. 2015;Hull et al. 2017;Maury et al. 2018;Kwon et al. 2019;Le Gouellec et al. 2019;Takahashi et al. 2019). Multiple origins have been proposed to explain the millimeter polarization by protoplanetary disks, including grain alignment with magnetic fields (Cho & Lazarian 2007;Bertrang et al. 2017), radiation fields (Lazarian & Hoang 2007;Tazaki et al. 2017), gas flow , and the self-scattering of thermal dust emission (Kataoka et al. 2015;Yang et al. 2016).
The polarization mechanism depends strongly on the grain size. Tazaki et al. (2017) theoretically showed that dust grains smaller than ∼100 μm are aligned with the direction of magnetic fields whereas those larger than ∼100 μm are aligned with the direction of radiation fields; the threshold of 100 μm depends on the properties of the dust grains and the strength of the magnetic field. Note that Yang et al. (2019) and Kataoka et al. (2019) also found that grains may be aligned with gas flow, but Tazaki et al. (2017) did not report this effect. Kataoka et al. (2015) found that the self-scattering of thermal dust emission can be detected only when the grain radius is around ∼λ/2π, where λ is the observation wavelength. Considering the significant diversity in grain growth in protoplanetary disks, we expect that the mechanisms of polarization will be different for different disks and will depend on location in the disks.
Previous observations have shown that the polarization mechanism depends on location for a lopsided disk. Kataoka et al. (2016b) found a sharp change in the direction of the polarization vectors, from the radial direction to the azimuthal direction, in the northern region of the disk of HD 142527, which is consistent with self-scattering in a face-on disk. Ohashi et al. (2018) further analyzed the same target and reported that the lopsided protoplanetary disk of HD 142527 has polarization caused by self-scattering in the northern region and caused by grain alignment with toroidal magnetic fields in the southern region at a wavelength of 0.87 mm. This indicates that the grain size is smaller in the southern region and larger in the northern region, which is consistent with the scenario of azimuthal dust trapping, where dust grains are trapped at the gas pressure bump azimuthally (e.g., Fukagawa et al. 2013;van der Marel et al. 2013;Pérez et al. 2014;Casassus et al. 2015;Pinilla et al. 2015;Boehler et al. 2017).
In contrast, inclined smooth disks usually show self-scattering features. Several disks have polarization vectors parallel to the disk minor axis (e.g., Stephens et al. 2014Stephens et al. , 2017Kataoka et al. 2017;Bacciotti et al. 2018;Cox et al. 2018;Girart et al. 2018;Hull et al. 2018;Harris et al. 2018;Lee et al. 2018;Sadavoy et al. 2018;Dent et al. 2019;Harrison et al. 2019;Mori et al. 2019), which is consistent with the self-scattering model for an inclined disk (Kataoka et al. 2016a;Yang et al. 2016Yang et al. , 2017. Note that HL Tau shows a self-scattering feature in Band 7 but azimuthal polarization in Band 3, which can be interpreted as radiative alignment or gas flow alignment (Kataoka et al. 2017;Yang et al. 2019).
In particular, polarization of the protoplanetary disk of HD 163296 was observed with a spatial resolution sufficient to resolve ring and gap structures (Dent et al. 2019). It was found that the polarization vectors are basically parallel to the disk minor axis in the central region; however, there are additional azimuthal components in gaps. Furthermore, the polarization fraction has a peak in the gaps between dust rings. These radial and azimuthal variations of the polarization fraction may reflect the dust properties and distributions in the protoplanetary disk. In this study, we perform radiative transfer modeling to determine what can be constrained based on spatially resolved polarimetric observations. We then model the continuum and polarized emissions of the HD 163296 disk to extract information on dust grains.
The protoplanetary disk around the Herbig Ae star HD 163296 has been well studied not only with ALMA (e.g., de Gregorio-Monsalvo et al. 2013;Rosenfeld et al. 2013;Flaherty et al. 2015;Isella et al. 2016Isella et al. , 2018Notsu et al. 2019) but also with optical and infrared telescopes such as the Hubble Space Telescope (HST), the Very Large Telescope (VLT), and the Spectro-Polarimetric Highcontrast Exoplanet Research instrument (SPHERE) (e.g., Grady et al. 2000;Kraus et al. 2008;Muro-Arena et al. 2018). The distance was derived to be 101.5 pc by Gaia Collaboration et al. (2018) and Bailer-Jones et al. (2018). The ALMA DSHARP project  found several rings using a high angular resolution of ∼3 au and reported detailed azimuthal asymmetries in the disk Isella et al. 2018). These rings and asymmetries may be produced by planets Zhang et al. 2018). The presence of planets is also suggested by deviations from Keplerian rotational motions (Pinte et al. 2018;Teague et al. 2018). However, optical and infrared high-contrast imaging has not identified companions of planetary mass yet, but has set upper limits on the masses of giant planets of a few M J (Guidi et al. 2018;Mesa et al. 2019).
The mechanism of ring formation is still under discussion even though ring structures have been discovered in many protoplanetary disks with ALMA (e.g., ALMA Partnership et al. 2015;Andrews et al. 2016Andrews et al. , 2018Canovas et al. 2016;Isella et al. 2016;Tsukagoshi et al. 2016;Cieza et al. 2017;Kraus et al. 2017;Loomis et al. 2017;van der Plas et al. 2017;Ansdell et al. 2018;Bertrang et al. 2018;Boehler et al. 2018;Clarke et al. 2018;Dipierro et al. 2018;Dong et al. 2018;Fedele et al. 2018;Long et al. 2018;Sheehan & Eisner 2018;van Terwisga et al. 2018;van der Marel et al. 2019). The presence of planets in a gap in the disk around PDS 70 was detected based on the Hα emission of accreting gas from protoplanets (Keppler et al. 2018;Haffert et al. 2019). Furthermore, the presence of asymmetries may also indicate planets or stellar companions. However, besides planets, variations in dust opacity at the snow lines of several molecular species (Zhang et al. 2015;Okuzumi et al. 2016), secular gravitational instability (Takahashi & Inutsuka 2014;Tominaga et al. 2018Tominaga et al. , 2019, dust accumulation at the edge of a "dead zone" (Flock et al. 2015;Ruge et al. 2016), and other mechanisms (Lorén-Aguilar & Bate 2015; Béthune et al. 2016;Dullemond & Penzlin 2018) can create ring and gap structures. A recent sample study found no correlation between gap radii and disk temperature, suggesting that a snow-line model is unlikely to be a common origin of multiring systems (van der Marel et al. 2019).
The rest of this paper is organized as follows. In Section 2, we analyze the polarization data of HD 163296 and describe the observed features in the disk. Section 3 explains the method used to obtain the disk model and the radiative transfer calculations. In Section 4, we investigate the polarization patterns of self-scattering in terms of grain size and dust scale height based on a smoothed disk without any rings as a simple case. In Section 5, we apply the self-scattering model to the HD 163296 disk and investigate the distributions of grain size and dust scale height. In Section 6, we discuss the disk structure by considering the observed variations of grain size and dust scale height. We estimate the turbulence strength α and discuss the existence of a dead zone. Finally, the results and the conclusions of this work are summarized in Section 7.

Analysis of Archival Data
We used the archival data of the 0.87 mm dust polarization from the protoplanetary disk around HD 163296 (2015.1.00616.S; PI: C. Pinte). These data were originally presented by Dent et al. (2019). The reduction and calibration of the data were done with CASA version 4.5.3 (McMullin et al. 2007) in a standard manner.
Stokes images were reconstructed using the CASA task tCLEAN, with Briggs weighting with a robust parameter of 0.5. The pixel size was set to 0 02. In addition, to improve sensitivity and image fidelity, self-calibration for both phase and amplitude was performed. The beam size of the final product was 0 20×0 18, corresponding to a spatial resolution of ∼20au ×18 au at the assumed distance of 101.5 pc. Note that Dent et al. (2019) examined the calibrated data in more detail, and applied further flagging. However, we found that our simple calibration is consistent with the results of Dent et al. (2019), as shown in the following section.
Stokes Q and U components produce polarized intensity (PI). Note that we ignore the Stokes V component in this study because it has not been well characterized for ALMA. The PI value has a positive bias because it is always a positive quantity. This bias has a particularly significant effect in measurements with low signal-to-noise ratio. We thus debiased the PI map as The rms noise values of the Stokes Q and U, and the PI were derived to be σ Q = 1.7×10 −5 Jy beam −1 , σ U =3.0×10 −5 Jy beam −1 , and σ PI =4.0×10 −5 Jy beam −1 , respectively. The rms values of the Stokes Q and U emissions were calculated from an area in the images in which emission was not detected. The polarization fraction (P frac = PI/I) was derived only where detection was above the threshold 3σ PI . Figure 1 shows images of the Stokes I, Q, and U, and the polarization fraction of the protoplanetary disk around HD 163296 at a wavelength of 0.87 mm. The polarization vectors are overlaid on the image of the polarization fraction. The location of the continuum ring at a radius of 70 au is shown by a black or white ellipse. These images are consistent with those obtained by Dent et al. (2019).

Results
The polarization fraction is ∼(0.8±0.04)% in the center and increases to ∼(2.3±0.5)% in the first gap at 50 au and to ∼(3±0.7)% in the second gap at 90 au along the major axis. The polarization fraction along the minor axis is lower than that along the major axis. The PI outside the 70 au ring was not detected along the minor axis.
The polarization vectors are mainly aligned with the disk minor axis. In addition, we confirmed the azimuthal twist of ±10°in the polarization vectors reported by Dent et al. (2019). This twist is located outside the 70 au ring and can be recognized in the Stokes Q image because the Stokes Q emission represents polarization pointing along the disk major axis in this disk geometry.
In this paper, we perform radiative transfer calculations by taking into account the polarization due to only self-scattering. This assumption is consistent with the results reported by Dent et al. (2019), who indicated that the polarization is produced by the selfscattering without any additional dust alignment mechanism.

Disk Model
To perform radiative transfer calculations, we constructed an axisymmetric dust disk model with several additional rings that reproduce the continuum image. We discuss the physical parameters of the disk in the following section. We adopted the following temperature profile with a smooth power-law distribution derived by Isella et al. (2016): We adopted a power-law surface density profile with an exponential tail and added exponential bright regions to mimic the observed profile. We took the intensity profile to be on the major axis in the northwest direction to avoid the additional crescent-shaped structure reported in the opposite direction (Isella et al. 2018). The intensity can be expressed as = - Using this equation, we assume the radial optical depth profile to be t = - The surface density is derived with the assumption that τ=κ abs ×Σ dust , where κ abs is the absorption opacity and Σ dust is the dust surface density. The radial profile of the intensity is set using Equations (1)-(3). The unit of intensity is set to be Jy beam −1 for comparison with the observations. The calculation of opacity followed that by Kataoka et al. (2016a). The dust grains were assumed to be spherical and to have a power-law size distribution with an exponent of q=−3.5 (Mathis et al. 1977) and maximum grain size a max . This maximum grain size is considered to be the representative grain size in the following discussion. The opacity was calculated using Mie theory. The composition was assumed to be a mixture of silicate (50%) and water ice (50%) (Pollack et al. 1994;Kataoka et al. 2014). We used the refractive index of astronomical silicate (Weingartner & Draine 2001) and water ice (Warren 1984) and calculated their mixture based on effective medium theory using the Maxwell-Garnett rule (e.g., Bohren & Huffman 1983;Miyake & Nakagawa 1993). Note that a different abundance may lead to a different absolute value of the degree of polarization, which should be investigated in future studies.
The vertical density distribution is assumed to be Gaussian with a dust scale height h d such that ), to mimic grain settling. We performed radiative transfer calculations with RADMC-3D 3 (Dullemond et al. 2012), taking into account multiple scattering, as done by Kataoka et al. (2015). The inclination and position angle were assumed to be 47°and 133°.3, respectively, as derived by the DSHARP project (Isella et al. 2018).
The Stokes I, Q, and U images were convolved with a Gaussian beam with a full width at half maximum of 0 2 to produce the final model images. Figure 2 shows the intensity profile (Stokes I) obtained from the continuum data on the major axis overlaid with the model intensity profile. The figure confirms that the model reproduces the continuum image.
The three rings are located at 70, 104, and 160 au, and their widths are 3.8, 5.6, and 16 au, respectively. Isella et al. (2018) derived the locations of the two inner rings as 67 au (width: 6.6 au) and 101 au (width: 5.8 au), respectively. These values are similar to those for our model; the difference is a few au, which is much smaller than the spatial resolution of ∼20 au.

Simple Disk Model
In this section, we first review the basic morphology of polarization vectors in inclined disks with a schematic illustration. Then, we perform radiative transfer calculations, focusing on the effects of grain size and dust scale height distribution, which are the key parameters in the following modeling of the HD 163296 disk. Figure 3 shows schematic illustrations of the polarization vectors for face-on, inclined, and edge-on disks (e.g., Kataoka et al. 2016a;Yang et al. 2016). The intensity profile is assumed to be smooth and axisymmetric. The polarization vectors are in the azimuthal directions for the face-on disk. This is because the radiation gradients are in the radial directions. In contrast, the edgeon disk has polarization vectors parallel to the disk minor axis because the dominant radiation is parallel to the disk midplane. The inclined disk has a combination of the polarization patterns of the face-on and edge-on disks. In the central region, the polarization vectors are mostly parallel to the minor axis; those in the azimuthal directions are around the edge of the disk. With the limited sensitivity of polarization observations, the polarized flux is mainly detected from the central region of inclined disks. Therefore, polarization vectors parallel to a disk minor axis have been observed in the central part of inclined disks (e.g., Kataoka et al. 2016a;Stephens et al. 2017;Bacciotti et al. 2018;Hull et al. 2018). In particular, Bacciotti et al. (2018) detected the polarization vectors parallel to the disk minor axis in the center and those in the azimuthal directions outside of the disk of DG Tau using highsensitivity observations. These studies observationally demonstrate that the polarization pattern changes radially in an inclined disk.

Dust Grain Size and Dust Scale Height Distribution
In this subsection, we discuss how dust settling can change the polarization pattern. To investigate the dependence of the polarization pattern and the polarization fraction on grain size and dust scale height, we performed radiative transfer calculations. For simplicity, we considered a smooth inclined disk. The density and optical depth profile were assumed to be The profile is the same as that in Figure 2 but without the ring structures.
To investigate effects of grain size distribution, we considered two different grain populations. One consists of dust grains with a power-law size distribution and a maximum grain size of 140 μm, and the other consists of those with a maximum grain size of 1 mm ( Figure 4). As shown in Figure 4, the number of dust grains is dominated by the smallest grains. However, the minimum grain size is small enough not to affect the results. The two populations of dust grains are representative of those that contribute to polarization and those that do not: the 140 μm population produces the most efficient polarization at an observation wavelength of 0.87 mm whereas the 1 mm population does not contribute to the polarization, but does to the Stokes I emission. To investigate the effects of dust population, we changed the mass ratio of the 140 μm and 1 mm populations while keeping the total optical depth fixed. With this setting, we expect that the 1 mm population will contribute to reducing the polarization fraction without changing the Stokes I intensity.
We also investigated the effects of dust scale height distribution. We performed radiative transfer calculations for various f set values to investigate the effects of dust settling. Figure 5 shows the results obtained for various dust scale heights and dust grain populations. The left, middle, and right columns show the cases for the 140 μm population, and combinations of the 140 μm and 1 mm populations with mass ratios of 7:3 and 5.5:4.5, respectively. Figure 5 clearly shows the dependence of the polarization fraction and vectors on dust grain population and dust scale height.   The ratio of the 140 μm and 1 mm populations changes the polarization fraction but does not change the polarization vectors. With increasing dust scale height, the azimuthal component is enhanced, especially in the outer region. As a result, the polarization fraction increases along the minor axis in the outer region.
The enhancement of the azimuthal components with larger dust scale height can be understood by considering the anisotropies of the radiation fields. Let us consider an increase in dust scale height with the total optical depth and column density fixed. As shown in Figure 6, although the optical depth along the line of sight does not change when the dust scale height increases, the optical depth perpendicular to the line of sight decreases because the dust spatial density decreases. Therefore, radiation from greater distances also contributes to the scattering. Thus, the incoming flux to be scattered by dust grains is more in the radial direction and becomes more anisotropic with increasing dust scale height. This explains the higher polarization fraction with azimuthal vectors for the thicker disk.
This enhancement of the azimuthal components in the disk with larger dust scale height results in the azimuthal variation of the polarization fraction. Along the minor axis, we identified a drop of the polarization fraction in the region where the polarization vectors change from the disk minor axis directions to the azimuthal directions. This drop is caused by the canceling-out of these two components, whose polarization directions are orthogonal to each other. In contrast, the polarization fraction along the major axis at a given orbital radius has a peak because the polarization components are in the same direction. Outside the regions of the peaks, the polarization fraction decreases with increasing radius because the radiation field becomes more isotropic along the major axis.
The other interesting feature is that the polarization vectors do not change when the dust components are changed, but the polarization fraction does. This is because the flux anisotropies determine the vector pattern, but the grain size does not change the flux anisotropies. We can thus conclude that the grain size only affects the polarization fraction; it does not affect the polarization vectors in self-scattering.
In summary, the azimuthal variation is mainly determined by dust scale height, and the absolute value of the polarization fraction can be tuned by changing the ratio of the two populations. Based on these qualitative findings, we discuss the detailed modeling of HD 163296 in the next section.

Application to HD 163296
In the previous section, we showed that the polarization fraction of the azimuthal component increases with increasing dust scale height, enhancing the azimuthal variation. The polarization observations of the HD 163296 disk show such azimuthal variation, especially beyond the 70 au ring, which may indicate that dust grains do not fully settle to the midplane of the disk. Note that the azimuthal component of the polarization can be recognized by the Stokes Q image in Figure 1 because the Stokes Q emission represents polarization pointing along the disk major axis in this disk geometry. Therefore, the detection of the Stokes Q emission outside the 70 au ring indicates the polarization of the azimuthal component. Here, we describe possible distributions of dust grains and dust scale heights in HD 163296. First, we show the case for a single size distribution, then we consider two dust populations, and finally we discuss the variations of dust scale height.

Single Grain Population
As a simple case, we consider a grain size population with a single power-law distribution with n(a)∝a −3.5 and a min = 0.1 μm. The maximum grain size a max is the main parameter.
We fix a max =140 μm, which produces polarization with the maximum efficiency. We vary the dust scale height of the disk and investigate whether our model matches the observations.
First, we investigate the effects of dust scale height on the total intensity image. Figure 7 shows the total intensity maps with f set =1 and 10. The total intensity images for f set =1 and 10 are almost identical and the projected ring width is different between the major and minor axes. This is because the convolved beam size (0 2) is insufficiently small to observe differences in the total intensity image.
Second, we investigate the effects of dust scale height on the scattering-induced polarization. Figure 8 shows the polarization fraction overlaid with the polarization vectors for various dust scale heights. In all cases, the polarization fraction is higher at locations in the gaps and lower at locations in the rings. In addition, the polarization vectors are parallel to the disk minor axis in the center and in the azimuthal directions in the outer regions. Different polarization patterns can be clearly seen, unlike the case for the total intensity. The polarization vectors in all cases are the combination of those parallel to the minor axis in the center and those in the azimuthal directions in the outer regions. However, the contribution of the azimuthal pattern to the polarization decreases as the dust scale height decreases. In addition, the variation of the polarization fraction in the azimuthal direction decreases with decreasing dust scale height. To confirm the tendency of the azimuthal patterns of the polarization vectors to become enhanced with increasing dust scale height, we plot histograms of the polarization vectors for each model in Figure 9. The number of polarization vectors is normalized so that the maximum number is one. As shown in the figure, the polarization vectors are more widely distributed in position angle with increasing dust scale height. This is consistent with the above results.
Here, we discuss why the polarization fraction is high in the gaps and low in the rings. To investigate the isotropy of the radiation in the rings and gaps, we derive the physical length (l) where the scattering optical depth becomes unity. The scattering optical depth is defined as  where κ scat is the scattering opacity, ρ dust is the spatial dust density, and Σ dust is the dust surface density.
where τ abs =Σ dust ×κ abs . Therefore, l decreases with f set assuming that the surface density is fixed. The scattering is contributed to only by nearby radiation if dust grains are settled. We calculated the scattering and absorption opacities using Mie theory and derived them to be ∼4.5 g cm −2 and ∼0.85 g cm −2 , respectively. Here, we assumed the grain size of a max =140 μm and the observing wavelength of λ=870 μm. At the 70 au ring of HD 163296, l<H gas because the optical depth at the ring is τ abs 1 and κ scat /κ abs ∼5. Therefore, the radiation field that contributes to scattering at the ring is almost isotropic because the width of the rings is similar to the gas scale height . In contrast, the optical thickness at the gaps is   Figure 8 for models with various f set values. The number of polarization vectors is normalized so that the maximum number is one. τ abs 10 −2 . Therefore, l>H gas , which means that the largerscale radiation gradient contributes to the scattering and thus the radial radiation gradients produce radial anisotropy and azimuthal polarization. This explains why the polarization fraction is higher in the gaps and lower in the rings. Figure 10 shows radial plots of the polarization fraction along the major and minor axes. The colors represent differences in dust scale height. The error bars indicate ±1σ of the observed polarization fraction. As shown in the figure, along the major axis, the polarization fraction increases in the gaps and decreases in the rings with increasing dust scale height. This opposite behavior is explained by canceling-out of the two orthogonal components of the scattering in the ring due to beam convolution.
In Figure 11, to understand this complicated behavior, we show the Q and U maps for f set =1 and f set =10 before beam convolution. The black circle indicates the 70 au ring. In this setup, the Stokes U emission basically corresponds to the polarization vectors parallel to the minor axis, and the butterfly-like pattern in the Stokes Q image corresponds to the azimuthal polarization vectors. The Stokes Q image for f set =1 shows emission at the 70 au ring, which indicates that the azimuthal component is produced. In contrast, the Stokes Q image for f set =10 shows almost no emission, which indicates that the polarization vectors are almost parallel to the disk minor axis. The Stokes U image for f set =1 shows negative values only in the 70 au ring of the major axis; the other regions have positive values. Therefore, for f set =1, the Stokes U intensity is decreased at the ring of the major axis due to beam dilution after beam convolution. The polarization for f set =1 in Figure 8 shows that the polarization vectors have different directions on the scale of the beam size and that the polarization fraction drops in the ring of the major axis.
We now explain why the Stokes U image for f set =1 produces negative values in the ring. The negative emission of the Stokes U in the ring of the major axis can be understood by regarding the ring as a thin disk, where the height is much smaller than the width, or as a thick tube, where the height is comparable to the width. For the thin disk, the scatteringinduced polarization has polarization vectors parallel to the disk minor axis, as described in the previous section. In contrast, as Kataoka et al. (2015) showed, the tube-like structure of the disk produces scattering-induced polarization pointing in radial directions on the ridge of the tube. This is because the net flux in the azimuthal direction is larger than that in the radial direction.
As shown in Figure 10, the polarization fraction along the minor axis at r100 au decreases with increasing dust scale height. This is because of the canceling-out of the different orthogonal polarization directions (the azimuthal and disk minor axis directions), as described in Section 4.2. Furthermore, the polarization fraction is slightly higher in the northeast direction than in the southwest direction, which indicates that the polarization fraction of the near (northeast) side is higher than that of the far (southwest) side (Flaherty et al. 2017). This is because the disk emission is marginally optically thick and the local inclination angles for the near and far sides are different (see also Figure 6 of Yang et al. 2017).
Here, we compare the observations and models. As shown in Figure 10, along the major axis, the polarization fraction for f set =2 is well matched. However, none of the f set values match the observations along the minor axes for both near and far sides. These models show a higher polarization fraction (∼2%-3%) than the observation (∼0.5%-1.5%) on the minor axis within r100 au. The low polarization fraction indicates that the grains are larger or smaller than 140 μm. Therefore, we performed the same calculations for a different grain size. Figure 12 shows the same plot as that in Figure 10 but for a grain size of 160 μm for the polarization fraction to fit along the minor axis. As shown in the previous section, the polarization vectors are the same as those in Figure 10, and the polarization fraction becomes lower when the grain size is changed. Even though the polarization fraction becomes lower than that in the previous models, the radial distributions of the polarization fraction do not match the observations for all f set . These models indicate that a model with a single grain population cannot reproduce the observations.

Model with Two Dust Grain Populations
Next, we investigate a model with two populations with different grain sizes. Here, we assume that these populations are those shown in Figure 4. The maximum grain sizes are set to 140 μm and 1 mm, respectively. The 140 μm population contributes to both polarization and total intensity, whereas the 1 mm population contributes to only total intensity.
In the previous section, the observed polarization fraction in the gaps was well explained by the model with 140 μm dust grains. However, for the central part and the ring regions, there was inconsistency between the observations and models. Therefore, we changed the dust grain populations in the central part and the rings while keeping the 140 μm population in the gaps.
We simply assume that the rings have only a 1 mm dust grain population because dust grains with a size of the order of 1 mm or larger have been suggested in studies of the spectral index β (e.g., Isella et al. 2007;Ricci et al. 2010;Testi et al. 2014). We assume that the 1 mm dust grain population is distributed within ±5 au from the centers of the rings. This is consistent with ring sizes of 3.8 and 5.6 au derived from highresolution observations (Isella et al. 2018).
In the central part, which is within r40 au, we use a mixture of the two dust grain populations (140 μm and 1 mm). We set a ratio of surface density (and density) of 5.5:4.5 (corresponding to a ratio of 11/9) between the 140 μm and 1 mm dust grain populations to reproduce the observed polarization fraction. Note that it is also possible to explain the lower polarization fraction by assuming a slightly larger (∼170 μm) or smaller (∼80 μm) dust grain size than 140 μm instead of the mixture of 140 μm and 1 mm dust grains.  Figure 13 shows images of the polarization fraction overlaid with the polarization vectors of the models of the two populations of dust grains. The panels show different dust scale heights. These models show that the polarization fraction increases on the major axis and decreases on the minor axis with increasing dust scale height. This trend is the same as that in Figure 8 and can be explained by the combination or canceling-out of the two components of the azimuthal and disk minor axis directions (see Section 4.2). Figure 14 shows radial plots of the polarization fraction along the major and minor axes, as done in Figures 10 and 12. A comparison of Figures 10 and 14 indicates that the models with two populations of dust grains match the observations better than do the models with a single grain population. For example, as shown in Figure 10, the polarization fraction along the major axis depends strongly on dust scale height for the single grain population. In contrast, the models with two populations of dust grains are less sensitive to dust scale height and thus better reproduce the observations. This is because the beam dilution does not affect the models with two populations of dust grains. As discussed in the previous section, the drop is caused by the beam dilution of the different polarization directions on the scale of the beam size on the rings. However, the 1 mm dust grains produce no polarized emission by scattering on the rings.
Although the models with f set =1-2 look similar to the observations in terms of the radial plots of the polarization fraction along the major and minor axes, these models do not exactly match the observations. For example, the radial profile of the minor axis of the near side shows that the polarization fraction of the models is slightly lower than the observations at 50au r100 au.
To compare the models directly with the observations, we investigate the Stokes Q and U images. Differences were found in the Stokes Q and U images between the models with two grain populations with f set =1-2 and observations, as shown in Figure 15. The observations show that the Stokes Q emission is only detected outside the 70 au ring with an intensity range of −8.7×10 −5 to 7.8×10 −5 Jy beam −1 . In contrast, the models produce the Stokes Q emission inside the 70 au ring as well as outside it. We investigated each model by comparing the absolute intensity value of the Stokes Q emission to fit the observations. For f set =1, the Stokes Q emission ranges from −1.1×10 −4 to 7.8×10 −5 Jy beam −1 , which is a slightly wider intensity range than the observations both inside and outside the 70 au ring. For f set =1.5, the Stokes Q emission ranges from −8.5×10 −5 to 7.6×10 −5 Jy beam −1 , which is similar to the observations outside the 70 au ring but is higher than the observations inside it. To reproduce the Stokes Q emission inside the ring, the dust scale height needs to be further reduced. However, for f set =2, the Stokes Q emission ranges from −8.1×10 −5 to 6.7×10 −5 Jy beam −1 , which is a narrower intensity range than the observations outside the 70 au ring. Therefore, we conclude that the radially constant dust settling parameter f set cannot reproduce the observations.
To decrease the Stokes Q emission, the maximum grain size can be made slightly larger (or smaller) than 140 μm or the population of the 1 mm dust grains can be increased. However, these conditions would also decrease the Stokes U, resulting in a lower polarization fraction. The non-detection of the Stokes Q requires dust grains larger than 200 μm or the mass ratio of 140 μm grains to 1 mm grains to be 0.4:0.6 with the assumption that f set =1.5. These dust grains give a polarization fraction of 0.4% in the central part, which is lower than the observations. Therefore, the only way to reduce the Stokes Q emission while keeping the Stokes U emission is to change the dust scale height. Therefore, the low Stokes Q emission cannot be explained by larger dust grains or a different mass ratio between the two grain sizes.

Radial Variation of Dust Scale Height
In the previous subsection, we showed that the models with two populations of dust grains match the observations better than do the models with a single grain population. However, there are still inconsistencies between the models and observations in terms of the radial profile and the Stokes Q and U distributions. Here, we introduce the radial variation of dust scale height to improve the models.
Outside the 70 au ring, the model with f set =1.5 is consistent with the observations. However, that model shows higher Stokes Q intensity than the observations inside the 70 au ring. As shown before, the Stokes Q emission is reduced by decreasing the dust scale height. Therefore, we need to reduce the dust scale height inside the 70 au ring. We investigate the dust settling parameter ( f set ) inside the 70 au ring to explain the non-detection of the Stokes Q emission. Figure 16 shows the absolute value of the peak intensity of the Stokes Q emission inside the 70 au ring for each model with f set =1-10. The absolute value is used because the Stokes Q emission has both positive and negative values. As shown in the figure, the Stokes Q emission is reduced by decreasing the dust scale height. The Stokes Q emission may saturate at f set ∼6 because it is also produced by the component of the disk minor axis. The position angle of the disk is set to be 133°.3. Therefore, the Stokes Q emission is slightly produced by the polarization of the disk minor axis even if the azimuthal component is not produced. By taking into account the non-detection of the Stokes Q emission inside the 70 au ring, we can set the upper limit of the dust scale height as H dust (1/3)H gas from the 3σ upper limit of the Stokes Q observations.
To reproduce the observations, we apply f set =10 within r68 au and f set =1.5 for r>68 au. Figure 17 shows a comparison between our model and the observations. We compare the images not only for the polarization fraction but also for the Stokes Q and U distributions. The model matches the observations well not only for the polarization fraction but also for the Stokes Q and U images. Furthermore, the twist of the polarization vectors pointed out by Dent et al. (2019) is reproduced. Therefore, the polarization vectors in our model are consistent with the observations. The original images of the Stokes I, Q, and U and the PI before convolution by the Gaussian beam are shown in Appendix B. We show radial plots of the polarization fraction for the best model and observations in Figure 18. The figure indicates that the model matches the observations well. In particular, the radial profile along the minor axis of the near side is improved compared to that in Figure 14. This may be because the radiation becomes more anisotropic due to the gradient of the dust scale height. Figure 19 shows a schematic view of the grain size distribution and the dust scale height for our best model. The best model is as follows. The 140 μm dust grains are distributed in the gaps, and dust grains contributing no polarization are distributed in the rings. In the rings, we simply assume 1 mm dust grains as a representative population that does not produce polarization at the observation wavelengths. However, 1 mm dust grains are just one possible example. The constraint from the polarization observations is the exclusion of grains with a size of around 140 μm. There are two possibilities for the central part because the observed polarization fraction is less than the prediction with 140 μm but still detected. One possibility is two populations of dust grains with grain sizes of 140 μm and 1 mm and a ratio of surface density of 5.5:4.5. The other possibility is a grain size distribution that can be expressed as a single power law; the maximum grain size would be ∼170 μm or ∼80 μm rather than 140 μm. Furthermore, the dust scale height inside the 70 au ring is at least three times smaller than the gas scale height, but it is two-thirds of the gas scale height outside the 70 au ring. Note that we cannot constrain the dust scale height in the rings because the polarization due to self-scattering is not produced in our model. An additional test of the single grain population with changing dust scale height is shown in Appendix A. As shown in the Appendix A, the observation results cannot be reproduced if the radially constant particle size is assumed even if the dust scale height is changed.

Discussion
In this section, we discuss the implication of these constraints and their physical origin.

Distributions of Dust Grains
We discuss the consistency in grain sizes derived from analyses of polarization and spectral index. Dent et al. (2019) derived the spectral index, α, using the ALMA 0.87 and 1.3 mm dust continuum observations and showed that it is significantly less than 2.0 in the central region,~2.0 in the 70 and 104 au rings, and slightly higher in the gaps than in the rings. In general, a low spectral index can be explained by optically thick emission or large dust grains. However, a spectral index of less than 2 can be explained by the effect of dust scattering. The continuum thermal dust emission becomes fainter than the blackbody radiation if the scattering is efficient (e.g., Liu 2019; Zhu et al. 2019). This high efficiency of dust scattering at a wavelength of 0.87 mm is consistent with the detection of self-scattering polarization. Therefore, a consistent grain size distribution between the spectral index and the polarization is as follows. The central region inside 40 au, where self-scattering is efficient, but not maximally efficient, and the spectral index is less than 2, is likely to be dominated by 80-170 μm dust grains. The rings, which have a low spectral index and polarization consistent with 1 mm grains, can be explained in two ways. (i) The rings are optically thick and the grain size is much larger or smaller than 140 μm; (ii) the rings are optically thin and the grain size is much larger than 140 μm. The case with the larger grains distributed in the rings is consistent with the scenario of dust trapping (e.g., . In contrast, the case with the smaller grains distributed in the rings is consistent with the grain growth scenario of the non-sticky ice such as CO 2 . The gaps, where the spectral index is high Figure 14. Same as Figure 10 but the models have two populations of dust grains, with grain sizes of 140 μm and 1 mm, respectively. The central part (r 40 au) is a combination of 140 μm and 1 mm dust grains at a ratio of surface density (and density) of 5.5:4.5; the rings within ±5 au have 1 mm dust grains, and the gaps have 140 μm dust grains. The radii of the three intensity rings are indicated by the dotted vertical lines.
Although the radial grain size distribution with large grains in rings and small grains in gaps is generally consistent with the scenario of radial dust trapping, it is still not fully consistent with grain growth and migration theory. The synthetic observations of dust trapping and a simulation of grain growth/migration with a Jupiter-mass planet show that the polarization fraction is only 1% at the peaks around the gaps (Pohl et al. 2016). This is because the grain size distribution is not radially constant in gaps. In contrast, the 0.87 mm polarization of HD 163296 shows that the peak polarization fraction is ∼4% in the gap center. The polarization results are more consistent with 140 μm dust grains widely distributed in the gaps. For example, shallow bumps in the gas for a less massive planet would produce a shallow grain size distribution.

Estimate of Turbulence Strength
The dust scale height is determined by the balance between dust settling and gas turbulence. In previous sections, we constrained the dust scale height and dust grain size based on the polarization measurements. Here, we constrain the gas turbulence. The dust scale height is constrained to be H dust /H gas 1/3 for the r70 au region, and H dust /H gas ∼ 2/3 for the r70 au region. 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)  To find the radial difference in turbulence, we discuss the turbulence strength in the two gaps at 50 and 90 au. We conclude that these gaps are filled with grains with a size of 140 μm. The gas surface density has been constrained by   (Isella et al. 2016). Using these values and assuming that the density of the dust material is ρ∼1.67 g cm −3 , the turbulence strength α is derived to be α1.5×10 −3 in the 50 au gap and α;0.015-0.3 in the 90 au gap.
The low value of α1.5×10 −3 inside the 70 au ring is consistent with the upper limit derived from molecular line observations where α<10 −3 (Flaherty et al. , 2017. Magnetohydrodynamic (MHD) simulations have shown that there is a floor value of turbulence strength of α  10 −4 as long as the surface layers are magnetorotational instability (MRI)active (Suzuki et al. 2010;Simon et al. 2015). Even without magnetic fields, there are hydrodynamic instabilities that make the disk turbulent, such as vertical shear instability, which has a turbulence strength of α∼10 −4 (Nelson et al. 2013). Lin (2019) studies the dust settling against the vertical instability by taking into account the dust back-reaction onto the gas, and shows that the dust can settle down to 0.1H gas . Therefore, the upper limit of α1.5×10 −3 is consistent with the simulations of low-turbulence disks.
The turbulence strength of α;0.015-0.3 outside the 70 au ring corresponds to v turb /c s ;0.12-0.55 by assuming a~v c s turb 2 ( ) , which is higher than the upper limit of v turb /c s 0.04-0.05 toward the midplane derived from the molecular line observations (Flaherty et al. 2017). One possible reason for these different estimations of turbulence strength is that the molecular lines and polarization are observed in different regions. The molecular line observations mainly trace the ring structures where the dust and gas surface densities are much higher than those in the gaps. The intensities of continuum and molecular lines are also high. Therefore, the line profile obtained with a larger beam size is dominated by the emission from the ring structures. In contrast, the polarization observations allow us to derive the turbulence strength only in the gaps because the dust grains in the gaps dominate the polarization. Therefore, the turbulence strength estimated from the molecular line and polarization observations may be for radially different regions of the disk. Even if the dust settling parameter ( f set ) is the same for the rings and gaps, the turbulence strength α in these regions could be different because the gas surface densities are different.
It should be noted that the gas surface density was estimated from the CO observations with the assumption of a cosmic molecular abundance ([ 12 CO]/[H 2 ] = 5 × 10 −5 ). If the CO abundance is lower or the CO intensity is reduced but the gas surface density does not drop in the dust gaps, the gas surface density will be higher than that used here. Then, the Stokes number and turbulence α will be lower. Booth et al. (2019) detected the rarest stable CO isotopologue, C O 13 17 , in the disk and derived the gas mass with a factor of 2-6 larger than previous gas mass estimates using C O 18 . In this case, the turbulence strength may be decreased up to a´- 2.5 10 4 in the 50 au gap and a~0.002 in the 90 au gap, respectively. Yet, the difference of an order of magnitude in turbulence strength between the inner and outer regions is robust. One of the consistent scenarios for low turbulence in the inner region and high turbulence in the outer region is a dead zone of MRI (Gammie 1996). The dead zone is considered as a region where the MRI is stabilized by ohmic dissipation (Sano & Miyama 1999). Furthermore, the MRI is also suppressed by the other non-ideal MHD effects of ambipolar diffusion (Desch 2004;Bai & Stone 2011;Dzyurkevich et al. 2013) and the Hall effect (Wardle 1999;Balbus & Terquem 2001;Bai 2014). Theoretical calculations suggest that the outer edge of the dead zone is located around 10-20 au (Sano et al. 2000;Semenov et al. 2004;Bai & Goodman 2009;Dzyurkevich et al. 2013). In contrast, the polarization observations suggest that the boundary between turbulence-active and -inactive zones is located at the 70 au ring, which is much larger than 20 au.
The radial transition from low to high turbulence strength may be explained by the electron heating (e-heating) zone rather than the dead zone (Okuzumi & Inutsuka 2015). In the e-heating zone, electrical conductivity is suppressed by heating of electrons due to turbulence. As a result, the turbulence strength is significantly reduced. Therefore, the e-heating zone may explain the low turbulence in the region where MRI is expected to be active. Mori & Okuzumi (2016) showed that the e-heating zone can extend to ∼100 au and that the dead zone can reach <50 au depending on the magnetic field, grain size, and gas density. Therefore, the boundary between the turbulence-active and -inactive zones at the 70 au ring may indicate the transition from the e-heating zone to the active zone.
Note that in terms of the ring formation process, the dead zone scenario is qualitatively consistent with dust accumulation at the outer edge of the dead zone, where dust grains are accumulated at the gas pressure bump created by the outer edge of the dead zone (Flock et al. 2015;Pinilla et al. 2016).

Comparison with SPHERE and HST Observations
Muro-Arena et al. (2018) observed the disk surface of HD 163296 in scattered light in the H and J bands with the SPHERE instrument. Based on the lack of scattered light emission detected in the rings at 104 and 160 au, they proposed two models. In the first model, small dust grains (a few microns) settle into the midplane in the outer disk and thus the outer regions are in the shadow cast by the 70 au ring. In the second model, there is a depletion of dust grains smaller than ∼3 μm in the outer disk. In this case, the scattered light has insufficient optical depth to be detected at the infrared wavelengths. According to our study with millimeter polarization, the dust scale height is constrained to be higher outside the 70 au ring than inside it. This is consistent with both models. For the first model, because there is no constraint on the dust scale height of the 70 au ring itself, it may be as large as  the gas scale height. The dust scale height was derived to be twothirds of the gas scale height outside the 70 au ring, and thus the outer regions will be in the shadow cast by that ring. The second model is also consistent with our model of dust scale height because there is no constraint on the dust scale height of the small dust grains from the SPHERE observations. It would be possible to distinguish the two models by performing millimeter-wave polarization observations at wavelengths longer than 0.87 mm to detect the self-scattering emission from the rings to constrain the dust scale height of the 70 au ring.
Observations of the disk were also carried out with HST (Grady et al. 2000;Wisniewski et al. 2008). Dust in the disk was detected in scattered light out to ∼500 au, which indicates that the outermost regions are flaring and small dust grains are directly illuminated by the central star. Wisniewski et al. (2008) found that the disk color (V -I) becomes redder at a radial distance of 3 5 corresponding to 355 au. They suggested that the red disk color shows evolution in the grain size distribution such as grain growth. This may be consistent with the depletion of the smallest dust grains in the second model proposed by Muro-Arena et al. (2018).

Caveats of Models
In Sections 5.2 and 5.3, we used two populations of dust grains to reproduce the observations. The maximum grain sizes were set to 140 μm and 1 mm, respectively. Even though we set the 1 mm dust grains in the rings, it is possible that dust grains much larger or smaller than 140 μm are distributed in the rings rather than the 1 mm dust grain population. To further constrain the dust size in the rings, polarization observations at longer or shorter wavelengths are needed.
Another caveat is that the polarization is assumed to be produced only by self-scattering. We show that the polarization of this disk at a wavelength of 0.87 mm can be understood only by self-scattering. However, polarization can also be produced by grain alignment mechanisms such as magnetic fields, radiation gradients, and gas flows. The observed polarization may also be contributed by the thermal emission of aligned dust grains. In this case, the azimuthal component may be explained by the dust alignment rather than the dust scale height of self-scattering. To investigate the contribution to polarization from dust alignment, observations with higher spatial resolution and at different wavelength are needed.

Conclusion
The protoplanetary disk around HD 163296 shows polarized emission at a wavelength of 0.87 mm (Dent et al. 2019), where ring and gap structures are spatially resolved. Motivated by the spatially resolved polarization of the ring-and-gap disk, we performed radiative transfer calculations to determine the key parameters to calculate the polarization vectors and fraction in a disk with rings and gaps. Furthermore, we also developed a model to reproduce the polarization of HD 163296 to constrain the disk parameters. The main findings are summarized as follows.
1. The key parameters that determine the polarization vectors and fraction due to self-scattering in an inclined disk are the dust scale height and the dust grain size or population, as summarized in Figure 5. As the dust scale height increases, the polarization fraction increases more along the minor axis than along the major axis. As a result, the azimuthal variation of the polarization fraction is enhanced. To further investigate the effects of grain size, we used two populations of dust grains, with maximum grain sizes of 140 μm and 1 mm, respectively. As the proportion of the 1 mm grain population increases, the polarization fraction decreases. This is because the 140 μm population produces polarization whereas the 1 mm population does not. 2. From the radiative transfer modeling of HD 163296 polarization observations, we found that the gaps require maximum efficiency of the polarization whereas the rings and the central region do not. In terms of grain size, the gaps are likely to be filled with the 140 μm grain population and the central region and rings are likely to have a significant contribution from a dust grain population with a grain size of much larger or smaller than 140 μm. The rings are consistent with no contribution from the 140 μm population. 3. We compared the size constraints with the spectral index analysis of HD 163296. The low spectral index (α ∼ 1.5-2.0) in the central part of the disk may indicate optically thick emission that is further reduced by scattering from dust with a size of ∼100 μm. The low spectral index (α ∼ 2.0) in the 70 au ring may be caused by the large optical depth and/or effects of grain growth because the 70 au ring is optically thick (τ ∼ 1.1). 4. The model indicates that the dust scale height is less than (1/3)H gas inside the 70 au ring and as high as (2/3)H gas outside it to explain the Stokes Q image. This suggests that the turbulence strength differs between the regions inside and outside the 70 au ring. We estimated the gas turbulence parameter at the two gaps after constraining the grain size and the dust scale height. The turbulence α is α1.5×10 −3 in the gap inside the 70 au ring and α∼0.015-0.3 in the gap outside that ring. A dead zone or e-heating zone may explain this transition. The dead zone scenario is qualitatively consistent with dust accumulation at the outer edge of the dead zone and ring formation.
investigate a model with a single grain population that takes into account the radial variation of the dust scale height. We compared the observations with the model with a single grain population by changing the dust scale height near the 70 au ring. We used a single grain population with a grain size of 140 μm, which is the same model as that described in Section 5.1. We applied f set =10 within r68 au and f set =1.5 for r>68 au. We found that the Stokes Q emission is produced around the 70 au ring, as shown in Figure 20, because the dust grains that contribute to the polarization are located on the ring. However, the observations show that the Stokes Q emission is not detected on the 70 au ring. Therefore, a model with a single grain population is unlikely to explain the observations even if the dust scale height is changed.

Appendix B Original Images of the Best Model and Future Prospects
Here, we discuss possible future observations with high spatial resolution by showing the original image before convolution with the beam size. As shown in Figures 17 and 18, our best model matches the observations well.
The original images of the Stokes I, Q, and U and the PI of the model are shown in Figure 21. These images are not smoothed by the Gaussian beam. We predict that the emission of the Stokes Q and U and the PI cannot be detected in the rings because dust grains larger or smaller than 140 μm cannot produce polarization due to self-scattering at a wavelength of 0.87 mm. The peak emission of the Stokes U and the PI could be detected inside and outside the 70 au ring. Observations with a large beam size thus show peak emission of the Stokes U and the PI around the 70 au ring. The Stokes Q emission could be detected only outside the 70 au ring because it is not produced inside that ring due to the low dust scale height (H dust 0.33H gas ). We can constrain the dust scale height inside the 70 au ring if we detect the Stokes Q emission in the future with higher-sensitively observations (see Section 5.3).