Thermal Emission and Scattering by Aligned Grains: Plane-Parallel Model and Application to Multiwavelength Polarization of the HL Tau Disk

Consider a slab that is infinite along the $x$- and $y$-axis direction in a Cartesian coordinate system, but finite in $z$. The slab is put between $z=0$ and $z=\Delta z$ with arbitrary units since the radiative transfer depends only on the optical depth (defined in Section 2.2) and not the physical depth. The slab has a uniform density $\rho$ and is isothermal with a temperature $T$. As shown in Fig. 1, $\theta$ is the angle from the $z$-axis and $\phi$ is the azimuthal angle from the $x$-axis in the $xy$-plane. For convenience, we define $\mu\equiv\cos\theta$. The slab consists of aligned grains with the alignment axis along the $x$-axis.

The unit vectors, $\hat{\theta}$ and $\hat{\phi}$, are used to denote the direction in increasing $\theta$ and $\phi$ respectively as shown in Fig. 1. The Stokes parameter vector, $\mathbf{I}\equiv(I,Q,U,V)^{T}$, as viewed by an observer in the direction $\hat{n}$ are defined in the plane formed by $\hat{\theta}$ and $\hat{\phi}$, i.e., the sky or image plane. Following the definitions from Mishchenko et al. (2000), positive Stokes $Q$ is polarization parallel to $\hat{\theta}$ and positive Stokes $U$ is polarization parallel to $\hat{\theta}-\hat{\phi}$. The linear polarization is

and the polarization angle is

where the function atan2 computes the angle in the appropriate quadrant based on the signs of Stokes $Q$ and $U$ (see Mishchenko et al. 2000). As shown in Fig. 1, the polarization angle $\zeta$ is measured from the direction of $\hat{\phi}$ in the image plane and $\zeta$ increases in the clockwise direction as seen by the observer. Polarization parallel to $\hat{\phi}$ has $\zeta=0^{\circ}$ and polarization parallel to $\hat{\theta}$ has $\zeta=90^{\circ}$. Note that the convention in Mishchenko et al. (2000) is different from the IAU 1973 convention in which the position angle is defined from the North and increases in the counterclockwise direction. For Sections 2 and 3, we use the Mishchenko et al. (2000) convention for consistency with the T-matrix code described below in Section 2.2.

For the observer located at the direction ($\theta,\phi$), the radiation transfer equation for a ray is

where $B_{\nu}$ is the black body radiation at some frequency $\nu$. $\overline{\mathbf{K}}_{\nu}$, $\mathbf{A}_{\nu}$, and $\overline{\mathbf{Z}}_{\nu}$ are the extinction matrix, absorption vector, and scattering matrix for the grain, and they depend on the grain orientation and material properties. For the remainder of the paper, we ignore the subscript $\nu$ for brevity, but note that all quantities related to intensity and opacity depend on frequency. $\overline{\mathbf{K}}$ is a 4-by-4 matrix and $\mathbf{A}$ is a 4 element vector both of which depend on where the ray is headed ($\theta,\phi$). The scattering matrix $\overline{\mathbf{Z}}$ is a 4-by-4 matrix which depends on where the scattered ray is headed ($\theta,\phi$) and also where the incoming ray was headed before scattering ($\theta^{\prime},\phi^{\prime}$). The opacity matrices depend on the dust model, which we will show in more detail below in Section 2.2. The three terms on the right-hand-side of Eq. (2.1) are called the extinction, thermal emission, and scattering terms.

We are interested in solving the emergent intensity, which is the Stokes parameter vector that the observer sees from this slab (see Fig. 2 for a schematic illustration). For an observer at $\mu>0$ (or $\theta\in[0^{\circ},90^{\circ})$), the emergent intensity is $\mathbf{I}(\Delta z,\theta,\phi)$. We assume that there is no impinging radiation from outside the slab, i.e., $\mathbf{I}(0,\theta,\phi)=0$ for $\mu>0$ and $\mathbf{I}(\Delta z,\theta,\phi)=0$ for $\mu<0$. For convenience, we define the Stokes $Q$ and $U$ normalized by Stokes $I$ respectively as

At some depth $z$, a grain scatters the incoming radiation $\mathbf{I}(z,\theta^{\prime},\phi^{\prime})$ to the observer at $(\theta,\phi)$ with the Stokes parameter $\overline{\mathbf{Z}}\mathbf{I}d\Omega^{\prime}$. Including scattering (the last term) in the radiative transfer equation greatly complicates the problem, which would otherwise be a simple integration. The amount of radiation towards the observer at each location no longer depends on just the local (non-radiation) quantities such as the temperature and density, but also the incoming radiation (the scattering term). However, knowing the incoming radiation means knowing the full radiation field. That includes the radiation towards the observer which is what we seek to solve in the first place. This is the “chicken or egg” effect that makes scattering problems difficult to solve.

To solve the “chicken or egg” effect, we divide the radiation field into discrete grid points and iterate for a converging solution. This procedure is commonly called “lambda iteration” (Mihalas, 1978). We use the Gauss-Legendre quadrature to integrate the $\mu$ direction by dividing the $\mu$-grid into an even number of sampling points $N_{\mu}$. We consider a linear $\phi$-grid with $N_{\phi}$ points and a linear $z$-grid with $N_{z}$ points. Before each iteration, we use the T-matrix method to produce the opacities at the discrete angular points (details are described in the Section 2.2). For all slab calculations below, we use $N_{\mu}=32$, $N_{\phi}=32$, and $N_{z}=128$. Higher number of sampling points did not change the results much which we demonstrate in Appendix C.

We begin the iteration by first assuming the radiation field is zero. In other words, we set $\mathbf{I}(z,\theta^{\prime},\phi^{\prime})=0$. In this case, the scattering term is zero, and we know the complete source term. We solve the differential equation numerically using the DELO-linear method (Rees et al., 1989; Janett et al., 2017) to obtain a first approximation to the radiation field $\mathbf{I}(z,\theta,\phi)$. We then simply replace the incoming radiation in the scattering term with the approximated radiative field and solve the differential equation again for a new estimate of the radiation field. This procedure repeats for several iterations until the radiation field converges. The radiation field is converged if the maximum relative difference between the new and the prior iterations is less than $0.1\%$. For small grains with low albedo, it takes less than 10 iterations to converge. For large optical depth and grains with large albedo, it can take a few hundred iterations to converge. Appendix A shows an example of how the emergent intensity converges and the converged solution of Stokes $I$ is tested against the analytical solution assuming isotropic scattering.

Once a converged solution to the radiation field is obtained, we can post-process to find the emergent intensity at any ($\theta$, $\phi$) which do not have to lie on the discrete grid. We use the dust model to evaluate the opacity matrices at the new ($\theta$, $\phi$) and use the discrete radiation field to evaluate the scattering term. A simple integration along the line of sight with the same boundary conditions will give us the emergent intensity.

The T-matrix method is a numerical modeling technique for light scattering by nonspherical particles of arbitrary size (see Mishchenko et al. 1996b and Mishchenko et al. 2000 for a review). We use the python package PyTMatrix¹¹1Leinonen, J., Python code for T-matrix scattering calculations. Available at https://github.com/jleinonen/pytmatrix/. to calculate the amplitude matrix and scattering matrix for prolate grains (Leinonen, 2014)²²2See Mishchenko et al. (2000) for details on calculating the extinction and absorption matrices from the amplitude and scattering matrices.. As described in Section 1, the choice of prolate grains is motivated by the work of Yang et al. (2019) and Mori & Kataoka (2021) who demonstrated the potential for prolate grains to explain the azimuthal polarization pattern observed in HL Tau in Band 3. PyTMatrix is a python interface for a T-matrix FORTRAN code (Mishchenko & Travis, 1994; Mishchenko et al., 1996a; Wielaard et al., 1997)³³3The FORTRAN code is available at https://www.giss.nasa.gov/staff/mmishchenko/t_matrix.html..

The prolate grain geometry is characterized by its aspect ratio $s$ and size $a$. The ratio between the short axis to the axis of symmetry (i.e., the long axis) is defined as $s$. A prolate grain has $s<1$, while a spherical grain has $s=1$. Its size is defined by the radius of the volume-equivalent sphere, i.e., the volume of the prolate grain is equal to the volume of a sphere with radius $a$.

We adopt the Disk Substructures at High Angular Resolution Project (DSHARP) mixture, which consists of 20% water ice, 33 % astronomical silicates, 7% troilite, and 40% refractory organics by mass (Birnstiel et al., 2018). The optical constants for water ice are from Warren & Brandt (2008), astronomical silicates from Draine (2003), and troilite and refractory organics from Henning & Stognienko (1996). We assume the MRN size distribution of $n(a)\propto a^{-3.5}$ where $a$ is the grain size (Mathis et al., 1977). The distribution cuts off at a minimum grain size $a_{\text{min}}$ and a maximum grain size $a_{\text{max}}$. Since the results are not sensitive to $a_{\text{min}}$, we set $a_{\text{min}}=0.01\mu$m (e.g Ricci et al., 2010; Kataoka et al., 2015).

Since radiation transfer relies on the optical depth, it is convenient to define the average extinction opacity of the non-spherical grain for non-polarized light:

We use this to define the vertical maximum optical depth of the slab as $\tau_{m}=\rho\kappa_{\text{ext}}\Delta z$.

Inclination of a slab with only spherical grains induces polarization through scattering that is parallel the direction of $\hat{\theta}$ in Fig. 1 with positive Stokes $Q$ (Yang et al., 2016a). When considering aligned oblate grains in the single-scattering limit, Yang et al. (2016b) showed that the thermal polarization superposes with the inclination-induced polarization due to scattering. We first look at how scattering of aligned grains produces polarization at different inclinations with the inclusion of multiple scattering.

We fix the vertical optical depth $\tau_{m}=1$ and show $q$ (the Stokes $Q$ normalized by Stokes $I$; Eq. (4)) as a function of inclination for $\phi=0^{\circ}$ and $90^{\circ}$ in Fig. 3. Results for other azimuthal angles are presented in Section 3.3. The prolate grains, as seen by the observer, are projected vertically for $\phi=0^{\circ}$ (see Fig. 2 for an illustration) and horizontally for $90^{\circ}$. At these two azimuths, Stokes $U$ is zero because of the symmetry of the system. Thus, linear polarization is completely described by $q$. To help understand the polarization curve from scattering aligned grains, we also plot two limiting cases: the polarization for a slab of scattering spherical grains of the same effective size and the polarization for a slab of aligned grains without scattering.

The limiting case of polarization for spherical grains helps gauge the degree of the polarization induced by inclination itself. It is computed by setting $s=1.0$ while keeping the rest of the parameters the same. The polarization curve is near zero when $\theta$ is close to zero (face-on). When the inclination increases, $q$ is positive and increases up to some peak at around $\theta\sim 75^{\circ}$. The curve is similar to Fig. 2 of Yang et al. (2017) which considered a semi-infinite slab with single scattering. In Appendix A, we illustrate with spherical grains how multiple scattering gradually adds Stokes $I$ and $q$ to the thermal emission until convergence.

For the non-scattering aligned grains ($s=0.975$), we calculate the emergent intensity through Eq. (2.1) without including the scattering term (though the scattering contribution to the extinction matrix remains). In other words, the grains are emitting polarized thermal radiation which undergoes dichroic extinction (photons are absorbed or scattered away), but none of the radiation is scattered into the line of sight of the observer. At $\phi=0^{\circ}$, $q$ of the non-scattering aligned grains is positive, because the long axis of the inclined grain looks vertical to the observer which produces positive Stokes $Q$. As the inclination increases, the observer views the prolate grain more pole-on (illustrated in Fig. 2), with a reduced degree of elongation which leads to less polarization. Additionally, dichroic extinction along the line of sight increases due to an increase in path length, which decreases the polarization even further. At $\phi=90^{\circ}$, $q$ is negative, because the long axis of the prolate grain looks horizontal (parallel to $\hat{\phi}$ in Fig. 1) to the observer which produces negative Stokes $Q$. As the inclination increases, the polarization just from the shape of the grain does not change. Instead, the decrease in $q$ is entirely due to the increase in path length along the sight line and thus the dichroic extinction. For both viewing azimuths, the magnitude of polarization from direct thermal emission decreases with increasing inclination which is opposite from scattering polarization of spherical grains (at least when $\theta<75^{\circ}$).

For aligned grains including scattering, the polarization remains completely positive at any inclination for $\phi=0^{\circ}$. When viewed at $\phi=90^{\circ}$, the polarization changes direction from being perpendicular to $\hat{\theta}$ (negative $q$) to being parallel to $\hat{\theta}$ (positive $q$). We can understand the inclination effects for scattering of aligned grains from the limiting cases of scattering from spherical grains and of non-scattering aligned prolate grains.

At low inclination, polarization from scattering aligned grains is similar to the polarization from non-scattering aligned grains (essentially thermal polarization). At the same time, polarization from scattering spherical grains reaches $0$, which means scattering is inefficient at producing the inclination-induced polarization. However, there is a slight difference between scattering aligned grains and non-scattering aligned grains at $i=0^{\circ}$, which arises from scattering due to the differential cross section of the grain instead of inclination explained by the following. In the case of spherical grains, radiation coming from the $x$- and $y$-axis in Fig. 1 scatter by $90^{\circ}$ to the observer at $i=0^{\circ}$ (or along the $z$-axis), which makes the radiation maximally polarized with a direction perpendicular to the respective scattering plane. Since the scattering cross section is the same for both, the two polarization cancel out, and the resulting polarization that reaches the observer is 0. The scenario is different for the prolate grain. The radiation coming along the pole of the grain (along the $x$-axis) sees a smaller scattering cross section than the cross section seen by radiation traveling perpendicular to the grain long axis (along the $y$-axis). Even though the scattered radiation is maximally polarized for incoming photons from both directions, a larger fraction of the photons coming along the $y$-axis is scattered than that along the $x$-axis. The resulting scattered polarization has the same orientation as that of the thermal polarization of the grain, which serves as a boost to the polarization (see Yang et al. 2016b for a description of the oblate case).

At higher inclinations, the thermal polarization decreases due to increasing dichroic extinction, while the polarization of scattering aligned grains is similar to the polarization with only spherical grains at both $\phi=0^{\circ}$ and $90^{\circ}$. This means polarization by scattering becomes dominant when the inclination is large enough. For $\phi=90^{\circ}$, the change in the sign of $q$ (and consequently the $90^{\circ}$ change in the polarization direction) is due to the changing balance between the two polarization mechanisms as the inclination angle and the associated optical depth along the sight line changes.

To isolate the optical depth effects, we will fix the inclination angle to $\theta=45^{\circ}$ and vary the total (vertical) optical depth of the slab $\tau_{m}$ in this section. The variation in optical depth can come from either the spatial (i.e., radial) variation of the dust distribution at a single wavelength or a change of observing wavelength at the same disk location. The results are illustrated in Fig. 4, where we plot the linear polarization fraction $q$ as a function of the total optical depth of the slab along the line of sight $\tau_{m}/\mu$ (where $\mu=\sqrt{2}/2$ for $\theta=45^{\circ}$) at $\phi=0^{\circ}$ and $90^{\circ}$. Similar to Section 3.1, we also show $q$ from non-scattering aligned grains and from scattering spherical grains. To better view the optically thin and optically thick limits, we show the same polarization results, but plot the optical depth in linear and log scales.

For the case of non-scattering aligned grains at both $\phi=0^{\circ}$ and $90^{\circ}$, the polarization is non-zero at low optical depth and determined simply by its projected shape (Draine & Lee, 1984). As the optical depth increases, dichroic extinction reduces the polarization from direct thermal emission monotonically (Hildebrand et al., 2000).

For scattering spherical grains, we see that at low $\tau_{m}/\mu$, the polarization approaches zero because the photons hardly scatter. The polarization peaks at an optical depth of $\sim 1$ and asymptotes to a constant value at large optical depths. This result is qualitatively similar to Fig. 3 in Yang et al. (2017) which considered only single-scattering.

The polarization of scattering aligned grains is a mix of these two limiting cases. At low optical depth, polarization from scattering aligned grains follows the thermal polarization because there is a lack of scattering for inclination effects to produce significant polarization. As the optical depth increases, scattering of photons is more likely to occur and produces polarization induced by inclination. At the same time, thermal radiation undergoes more dichroic extinction as it travels through the optically thick slab and thermal polarization decreases. As a result, inclination-induced polarization dominates at large optical depth. For $\phi=0^{\circ}$, polarization stays positive entirely, because the thermal polarization and the inclination-induced polarization both produce positive Stokes $Q$. For $\phi=90^{\circ}$, polarization is initially negative at low $\tau_{m}/\mu$, but it becomes positive at higher optical depth. There is a cross-over point (or polarization reversal) at $\tau_{m}/\mu\sim 0.7$. It makes physical sense because the thermal polarization (negative $q$) is opposite to the inclination-induced polarization (positive $q$) and the former dominates at low optical depth, while the latter dominates at high optical depth. The cross-over point is not limited to $\tau_{m}/\mu$ of order unit, but instead depends on the aspect ratio or the grain. For example, a spherical grain cannot directly emit any $q$ and the cross-over happens effectively at $\tau_{m}/\mu=0$, while a slab of highly elongated prolate grains will require much larger optical depth to take out the large intrinsic polarization.

At high optical depths, there exists a slight deviation in polarization between the spherical grains and that of scattering aligned grains. The polarization at $\phi=0^{\circ}$ is slightly less than the polarization of the spherical grains slab, while the polarization at $\phi=90^{\circ}$ is slightly higher. The deviation is on the order of the deviation of polarization of non-scattering aligned grains from zero. As we demonstrate below, this is a result of the incomplete cancellation between the polarization produced by dichroic extinction and that by thermal emission.

Consider the non-scattering aligned grains case of an observer along $\phi=0^{\circ}$ or $90^{\circ}$. Since the grain does not emit Stokes $U$ or $V$ in this frame (i.e., $A_{3}=A_{4}=0$), we only have to solve for Stokes $I$ and $Q$. The thermal polarization from a grain is simply $A_{2}/A_{1}$. Also, in this frame, only two elements of the extinction matrix $\overline{\mathbf{K}}$ are independent which we define: $K_{1}\equiv K_{11}=K_{22}$ and $K_{2}\equiv K_{12}=K_{21}$. The resulting radiation transfer equation along path $s$ becomes

which is a set of first order nonhomogeneous differential equations.

Since scattering is not involved, we can easily obtain an analytical solution for $[I,Q]^{T}$. The eigenvalues for the extinction matrix are $K_{\pm}\equiv K_{1}\pm K_{2}$. The complete solution assuming no external radiation as a boundary condition is

where $A_{\pm}\equiv A_{1}\pm A_{2}$. In the optically thin limit, we retrieve $Q/I=A_{2}/A_{1}$ as expected from a single grain. In the optically thick limit, the polarization is

where the approximation on the right-hand-side applies because $A_{2}/A_{1}$ is only a few percent for grains with small aspect ratio and the two quantities, $A_{2}/A_{1}$ and $K_{2}/K_{1}$, are comparable. Note that $K_{2}=0$ and $A_{2}=0$ for spherical grains. Eq. (8) simply means the prolate grains emit thermal polarization, $A_{2}/A_{1}$, but the same grains attenuate the polarization through their dichroic extinction, $K_{2}/K_{1}$. In the limit of small grains when the scattering opacity is negligible, $A_{i}$ are equal $K_{i}$, which makes the polarization in the optically thick limit zero (e.g., Hildebrand et al. 2000). However, if the dichroic extinction does not perfectly attenuate the thermal polarization, then we get residual polarization even in the optically thick limit⁴⁴4Even if $A_{i}=K_{i}$, residual polarization can also be achieved when there is a temperature gradient (Yang et al., 2017; Lin et al., 2020a). For the chosen grains of $x=0.4$, the dichroic extinction takes out more polarization than the grains themselves can emit which leaves the polarization perpendicular to its thermal polarization. Note that the polarization in this limit, as expressed by Eq. (8), is typically small.

Understanding the azimuthal variation of polarization requires us to consider $\phi$ other than the two specific azimuths discussed thus far. We use the same slab above, but fix the inclination at $45^{\circ}$. Since the system is no longer mirror symmetric in general, we need both Stokes $Q$ and $U$ to describe the linear polarization instead of just Stokes $Q$. Each row in Fig. 5 shows how the $q$ and $u$ vary as a function of azimuth along with $p_{l}$ and $\zeta$. Each column corresponds to a different $\tau_{m}/\mu$ to see how the azimuthal profile changes with different optical depth. We plot only from $\phi=-90^{\circ}$ to $90^{\circ}$ since the system is symmetric under $180^{\circ}$ rotation for our setup.

We also plot the case of scattering spherical grains and non-scattering aligned grains. For scattering spherical grains, the emergent intensity does not depend on $\phi$. Thus, all Stokes parameters are constant across $\phi$. $q$ (and consequently Stokes $Q$) is positive due to inclination-induced polarization. On the other hand, $u$ is entirely 0, because the spherical grain cannot emit any polarized emission directly and scattering cannot produce Stokes $U$ due to of the symmetry of the system. As such, the polarization angle is always $90^{\circ}$.

It is helpful to understand the non-scattering aligned grains starting with the most optically thin case $\tau_{m}/\mu=0.05$. As seen from the $\phi=0^{\circ}$ and $90^{\circ}$ shown in Sections 3.1 and 3.2, the sign of $q$ depends on the orientation of the grain. Shown in Fig. 5, the switch in the sign occurs roughly at $\phi\sim 35^{\circ}$ as the projected grain looks slightly more horizontal than vertical. Also at $\phi\sim 45^{\circ}$, $u$ is maximized since the projected long axis of the grain looks diagonal with $\zeta\sim 135^{\circ}$ (Fig. 5b). The minimum $u$ happens at $\phi\sim-45^{\circ}$ when the grain also looks diagonal but with $\zeta\sim 45^{\circ}$.⁵⁵5Considering a dipole in the optically thin limit, the switch of the sign of $q$ is at $\phi=\pm 35^{\circ}$ when $\theta=45^{\circ}$. The maximum and minimum for $u$ is at $\phi=45^{\circ}$ and $-45^{\circ}$ respectively regardless of $\theta$. At other $\tau_{m}$, the azimuthal variations of $q$ and $u$ from the thermal emission by non-scattering prolate grains are similar except that the magnitudes of $q$ and $u$ are decreased because of an increased dichroic extinction at a higher optical depth.

For scattering aligned grains, in the optically thin case of $\tau_{m}/\mu=0.05$, the $q$ and $u$ (Fig. 5a and b) are similar to those of the non-scattering aligned grains except there is a positive offset of $q$ relative to the non-scattering aligned grain. The amount of offset is similar to $q$ of scattering spherical grains. The $u$ in this case is completed dominated by the direct emission from the aligned prolate grains with no contribution from scattering. In Fig. 5c, including scattering to aligned grains decreases the degree of azimuthal variation of $p_{l}$. This is due to the addition and canceling effects between inclination-induced polarization from scattering and thermal polarization at different azimuth. The inclination adds to the positive $q$ thermal polarization at $\phi=0^{\circ}$ and decreases $q$ at $\phi=90^{\circ}$ as shown in Section 3.1. At $\phi=\pm 45^{\circ}$, $u$ largely determines $p_{l}$ meaning that adding scattering does not alter $p_{l}$ since $u$ from scattering is 0. The slight contribution from $q$ due to scattering creates the deviation of polarization angle $\zeta$ from the non-scattering aligned grains case (Fig. 5d). Specifically, the polarization angle deviates towards $90^{\circ}$ (or parallel to $\hat{\theta}$; Fig. 1).

For the unity optical depth case of $\tau_{m}/\mu=1$, $q$ increased relative to $\tau_{m}/\mu=0.05$ as expected since scattering becomes more prominent and contributes positive $q$ due to inclination (Fig. 5e). Due to this increase in $q$, $p_{l}$ now peaks at $\phi=0^{\circ}$ and reaches a minimum at $\phi=\pm 90^{\circ}$ (Fig. 5g). The location of the peak and trough is opposite to the $\tau_{m}/\mu=0.05$ case (Fig. 5c). It is likely that there exists an optical depth $\tau_{m}/\mu$ between 0.05 and 1 that makes the polarization fraction nearly independent of the azimuthal angle $\phi$. Also, $p_{l}$ can become roughly flat across azimuth when the level of scattering is just enough at some optical depth. The polarization angle $\zeta$ shown in Fig. 5h evolves towards $90^{\circ}$ as the optical depth increases.

When the optical depth reaches $\tau_{m}/\mu=5$, $q$ of the scattering aligned grains varies around the level of $q$ of scattering spherical grains in Fig. 5i. In contrast, $q$ from non-scattering aligned grains is low. Thus, we see that the inclination-induced polarization due to scattering dominates the azimuthal profile. The thermal polarization provides a deviation to $q$, but there is an additional deviation particularly at $\phi\sim 45^{\circ}$. In Fig. 5j, $u$ from scattering aligned grains does not resemble that from non-scattering aligned grains as was the case for lower optical depths (Fig. 5b,f). Since $u$ is an order of magnitude smaller than $q$, the effect on $p_{l}$ is small. As a result, $p_{l}$ mainly follows $q$ and the polarization angle is mainly constant at $90^{\circ}$ (Fig. 5k,l).

From Fig. 3, 4, and 5, we can observe that $q$ of scattering spherical grains and $q$ of non-scattering aligned grains roughly add together linearly to produce $q$ of scattering aligned grains. For direct comparison, we add the $q$ and $u$ quantities of the non-scattering aligned grains with those of scattering spherical grains in Fig. 5 and find close agreement with the proper scattering aligned grains case for $\tau_{m}/\mu=0.05$ and $1$. The linearity is because the scattering term and thermal emission terms add linearly as the source term. Furthermore, the prolate grain of $s=0.975$ is nearly spherical which makes the scattering matrix not too different from the spherical case (e.g. Kirchschlager & Bertrang, 2020) and that is why inclination produces similar polarization behavior for both. Scattering of a non-spherical grain can cause deviations (for example, at $i=0^{\circ}$ in Fig. 3) though that is small compared to what inclination can produce (however, see Appendix D for a case with highly elongated prolate grains).

The linearity breaks down when the optical depth is large (Fig. 5i): scattering changes the radiation field in the slab drastically and has a strong impact on the scattering source term. This impact is non-local and is not important at low optical depth, but becomes obvious at high optical depths. Thus, to a good approximation, particularly in the optically thin and moderately optically thick regimes, the bulk of the polarization is the superposition of the thermal polarization from the non-spherical grain and the inclination-induced polarization approximated by spherical grains (see Appendix E for a case with larger size parameter).

We use the HL Tau multiwavelength images shown in Stephens et al. (2017), but we rotate the image such that the disk minor axis is along the vertical direction of the image (Fig. 6). For convenience, we will call this the “principle frame” because the major and minor axes form the horizontal and vertical axes of the image for an inclined axisymmetric flat disk. The disk inclination $i$ is $46.7^{\circ}$ ($i=0$ means face-on), and the position angle of the disk minor axis $\eta$ is $48.2^{\circ}$ East-of-North (ALMA Partnership et al., 2015). The spatial coordinates of the image in the principle frame $(x^{\prime},y^{\prime})$ is related to the coordinates of the original image by

where $x$ and $y$ are the RA and DEC relative to the disk center of the original image. We adopt the distance of 147 pc (Galli et al., 2018).

The coordinate rotation changes the reference frame for the Stokes parameters as well. We use $\mathbf{I}$ to represent the original Stokes parameters as observed by ALMA and use primes, $\mathbf{I}^{\prime}$, to denote the Stokes parameters in the principle frame. The Stokes parameters are related by

while Stokes $I$ equal Stokes $I^{\prime}$. Since inclination-induced polarization of a disk produces polarization parallel to the disk minor axis, such polarization will exhibit positive Stokes $Q^{\prime}$. The deviation of the polarization angle away from the disk minor axis will show up as Stokes $U^{\prime}$, which makes it easy to identify in the rotated frame. Fig. 6 shows the resulting linear polarized intensity, Stokes $Q^{\prime}$, and Stokes $U^{\prime}$ all in mJy beam^-1 (note this is not $q^{\prime}$ and $u^{\prime}$) and the polarization fraction image with the polarization vectors to denote the polarization angle.

Assuming the dust disk is geometrically thin (e.g. Pinte et al., 2016), we can deproject the disk with a known inclination $i$ and obtain the azimuthal profiles by relating the spatial coordinate to the radius and azimuth:

where $r$ is the radius. $\phi_{d}$ is the azimuth of the disk with $\phi_{d}=0^{\circ}$ starting along the right major axis going counterclockwise.

Before fitting the data quantitatively, we will benefit from visualizing the slab models in the form of an image and discuss the main features. Since the disk is geometrically thin, we can approximate each point in the disk image as an independent slab, and we piece together multiple slab calculations to form the image.

As described in Section 1, we assume the prolate grains are aligned azimuthally around the disk. As such, we let the $x$-axis of the slab follow the direction opposite to increasing $\phi_{d}$ and let the $y$-axis of the slab point away from the disk center. The azimuth of the slab is related to the azimuth of the disk by

The inclination of the disk corresponds to the polar angle of the slab $i=\theta$. Due to the different definitions of the Stokes parameters, Stokes $U$ (as defined by ALMA according to the IAU 1973 convention) differs from the Stokes $U$ (as defined by Mishchenko et al. 2000) of the plane-parallel slab by a negative sign. For direct comparisons, we convert the Stokes parameters of the slab (Mishchenko & Travis, 1994) to the Stokes parameters of ALMA.

For multiwavelength polarization images like that of HL Tau, the optical depth decreases at longer wavelength, because the opacity of the grains decreases. At the same time though, the albedo varies and the refractive index also has a wavelength dependence. To isolate the effects just from optical depth, we will consider an image at $\lambda=1$mm with grains of a maximum size parameter of $x=0.4$, which keeps the albedo and refractive index fixed, but we scale the surface density up and down to change the optical depth. We adopt a fairly simplistic but representative surface density profile (Lynden-Bell & Pringle, 1974) and parameterize the dust surface density by

where $r$ is the radius in au and $\tau_{0}$ is a characteristic optical depth. The temperature profile is 20 K at 100 au and goes as $r^{-0.5}$.

Fig. 7 shows the linear polarized intensity, Stokes $Q^{\prime}$, Stokes $U^{\prime}$, and $p_{l}$ along with the polarization vectors at different $\tau_{0}$. The different $\tau_{0}$ values are chosen to illustrate particular changes of features in the polarization images.

In the most optically thin case, $\tau_{0}=0.01$, the polarization fraction is largest along the minor axis and smallest along the major axis of the disk (Fig. 7m). Also, the polarization vectors are elliptical. At this limit, scattering is negligible, and thus the polarization pattern is a simple result of the viewing geometry of the azimuthally aligned prolate grains. Along the minor axis of the disk, the symmetry axis of the grain is directed horizontally and seen from the longest side which produces the largest thermal polarization. The Stokes $Q^{\prime}$ is negative and the Stokes $U^{\prime}$ is zero because the polarization is completely horizontal. This corresponds to the $\phi=90^{\circ}$ case in Section 2. The grains along the major axis are seen more pole-on, and thus the projection decreases the thermal polarization. Since the polarization is completely vertical, the Stokes $Q^{\prime}$ is positive, and Stokes $U^{\prime}$ is zero. In Fig. 7, the polarized intensity is the greatest near the center simply because the temperature is highest. Additionally, it is slightly larger along the minor axis than along the major axis for a given radius because the polarization fraction is higher (this is more obvious when beam convolution is considered in Section 4.2.1).

With a slight increase in optical depth, $\tau_{0}=0.05$, the outer regions of the disk remains largely similar to the most optically thin case across all four polarization properties. However, the inner regions in the polarization fraction (Fig. 7n) begin to see a horizontal bar of higher polarization fraction. This is due to the emergence of a positive Stokes $Q^{\prime}$ as we can see in Fig. 7f which is a result of scattering.

The trend continues to the $\tau_{0}=0.5$ case, where the positive Stokes $Q^{\prime}$ becomes even more obvious. The polarized intensity in Fig. 7c, the Stokes $Q^{\prime}$ in Fig. 7g, and $p_{l}$ in Fig. 7p all have two prominent lobes across the major axis. The lobes are formed because the thermal polarization and inclination-induced polarization add with each other. In addition, $p_{l}$ begins to develop two polarization “holes” along the minor axis above and below center. This corresponds to the location where Stokes $Q^{\prime}$ changes from a positive value to a negative value going away from the center, in which case the polarization vector shifts by $90^{\circ}$ in Fig. 7p (see also $\phi=90^{\circ}$ of Fig. 4). The pattern of the polarization fraction in the outer regions of the $\tau_{0}=0.01,0.05$ cases (Fig. 7m, n) disappears as thermal polarization gives way to scattering polarization due to higher optical depth for $\tau_{0}=0.5$ and $3$.

There are two symmetries that we can exploit when assuming azimuthally aligned prolate grains and a geometrically thin disk. For any inclined axisymmetric disk, the image is mirror symmetric across the disk minor axis (i.e., the left half is mirror symmetric to the right half in the principle frame).⁶⁶6Geometrically thick disks or highly inclined disks can feature near/far-side asymmetry in the Stokes $I^{\prime}$ and polarization properties (Yang et al., 2017), but the image will remain mirror symmetric across the minor axis. For Stokes $U^{\prime}$, the symmetry involves a change in sign. Since we further assumed that the disk is geometrically thin, the image is also mirror symmetric across the disk major axis (the top half to the bottom half) with a change in sign for Stokes $U^{\prime}$ again. The models in Fig. 7 clearly demonstrate the two symmetries.

The expected symmetries for an axisymmetric disk motivates us to “fold” the observed images twice: we first average $\mathbf{I}^{\prime}$ across the disk minor axis and then average across the disk major axis. The benefits of folding is to average out the asymmetries and to improve the signal-to-noise ratio of the data. The asymmetries in the observed data (Fig. 6) cannot be explained by the axisymmetric disk model adopted in this paper and in most of the disk polarization models in the literature to date. Thus, we opt to remove the asymmetries for comparison with axisymmetric disk models as a first step towards a more comprehensive model. Averaging across the two symmetries increases the signal-to-noise by a factor of 2. Table 1 lists the improvement of the signal-to-noise for each Stokes parameter at each wavelength (before the image was corrected by the primary beam). We adopt a factor of 2 improvement for simplicity in the following analysis.

We should caution that an elongated observing beam can break the two symmetries described above. If the beam is elongated and not along the disk major or minor axis, then the image does not have any symmetry. The beams of the Bands 3 and 7 images are not along the disk major or minor axis, but they are at least roughly circular. The Band 6 beam happens to be elongated roughly along the disk major axis, which should alleviate the impact of the beam orientation. Utilizing the symmetries in the visibility domain may produce more robust results, but it is beyond the scope of this paper.

Fig. 9 shows the folded images for all three wavelengths. After folding the image twice, we copy the quadrant back into the four quarters for visualization. We take the geometric mean of the beam major and minor axis as a rough estimate of the resolution. One benefit from folding is that the azimuthal variation of the $p_{l}$ at Band 3 becomes clearer. In particular, the level of $p_{l}$ is greater along the minor axis than along the major axis which is what we expect from azimuthally aligned prolate grains (see Fig. 7m, n and Fig. 8m, n) if the image is optically thin⁷⁷7In fact, Carrasco-González et al. (2019) inferred the Band 4 optical depth to be $\sim 1$ from $\sim 60$ to $100$ au which suggests that Band 3 is optically thin over most of the disk where polarization is detected..

In Section 2, we treated self-scattering of aligned grains consistently. From the assumed aspect ratio of the grain, we can observe that, to a good approximation, the thermal polarization of non-spherical grains adds linearly with scattering polarization of spherical grains to produce the polarization of scattering aligned grains (see Fig. 5). Furthermore, assuming azimuthally aligned grains, the thermal polarization varies with azimuth, while scattering polarization of spherical grains is constant across azimuth. Thermal polarization adds with scattering polarization along the major axis of the disk. On the other hand, along the minor axis of the disk, thermal polarization cancels with scattering polarization. The superposition breaks down for optically thick cases (which is more likely for Band 7), but the error is regulated since we have the least contribution from thermal polarization. Building on these two approximations, we can roughly differentiate the level of polarization from thermal emission and that from scattering from a single image.

Assuming the grain is in the dipole limit, the thermal polarization depends on the projected shape of the grain. For a prolate grain, we have (Lee & Draine, 1985; Yang et al., 2016b):

where $\theta_{g}$ is the viewing angle from the axis of symmetry of the prolate grain. If $\theta_{g}=0^{\circ}$, the grain is seen pole-on, while $\theta_{g}=90^{\circ}$ means the grain is seen from the side (perpendicular to the axis of symmetry). We call $p_{0}$ the intrinsic polarization, since it is the maximum polarization of the prolate grain determined just from its shape. The approximation in the right-hand-side applies because $p_{0}$ is much less than unity (recall Band 3 polarization is only $\sim 2\%$).

The polarization is attenuated by dichroic extinction as optical depth increases (Hildebrand et al., 2000). However, since $p_{0}$ is small, the optical depth attenuates the polarization by the same factor regardless of $\theta_{g}$. Thus, we define the $p_{0}$ attenuated by optical depth as $T_{0}$; the “$T$” stands for the “thermal” polarization component (not to be confused with the temperature $T$). The observed thermal polarization after attenuation as a function of $\theta_{g}$ is

Along the minor axis of the disk, the azimuthally aligned prolate grain produces a negative Stokes $Q^{\prime}$ since the projected elongation is parallel to the disk major axis. The magnitude of the polarization is just $T_{0}$, since the grain is seen from the edge (or $\theta_{g}=90^{\circ}$). Thus, thermal polarization cancels with the polarization due to scattering which makes:

where $S$ is the polarization fraction due to scattering approximated by spherical grains. Along the major axis of the disk, $\theta_{g}=90^{\circ}-i$ and the thermal polarization has the same sign as scattering polarization. Thus, the polarization due to scattering and thermal polarization add together:

From the two algebraic equations, we can easily solve for the two unknowns $S$ and $T_{0}$.⁸⁸8We should note that scattering aligned grains can produce Stokes V which is not captured in the simple linear combination. Given that HL Tau currently does not have a firm detection of Stokes V, the available data does not contradict the linear combination approximation.

In Fig. 10, we show the observed $q^{\prime}_{\text{major}}$ and $q^{\prime}_{\text{minor}}$ across wavelength at a radius of 100 au. The radius is chosen to maximize the number of independent beams across azimuth, but also to avoid low signal-to-noise. The uncertainties are based on statistical error from the noise levels of the folded Stokes $Q^{\prime}$ and Stokes $I^{\prime}$. For the rest of the paper, we only consider statistical noise, but for noise levels that are less than $0.1\%$, we use $0.1\%$ since the ALMA instrumental error of polarization fraction is expected to be about $0.1\%$. From just the observational data, it is easy to see why scattering polarization alone does not work and thermal polarization is necessary. Based on $q^{\prime}_{\text{major}}$ (Fig. 10a), the measured polarization fraction increases with increasing wavelength. This is opposite to what is expected for scattering alone, which should be a near monotonic decrease (Lin et al., 2020b).

Furthermore, we also plot the disentangled polarization due to scattering $S$ and that due to thermal polarization, which is $T_{0}\cos^{2}i$ along the major axis and $-T_{0}$ along the minor axis. The uncertainties are based on error propagation. We can see that $S$ is smaller at Band 3 than at Band 6 and 7, while the $S$ between the Bands 6 and 7 are comparable. The behavior is expected from scattering. As seen in Fig. 4, the scattering polarization of spherical grains increases with increasing optical depth until the optical depth is $\sim 1$. Just from the observed spectrum, we can infer that Band 3 is optically thin, while Band 6 and 7 have optical depths of unity or larger.

The contribution from thermal polarization is the least at Band 7 and gradually increases to Band 3. The monotonic increase in thermal polarization to longer wavelength is also expected because the optical depth decreases which leads to less dichroic extinction as depicted by the non-scattering aligned grains case in Fig. 4.

Once we solved for $T_{0}$, we can predict the thermal polarization at any azimuth based on geometrical arguments. Recall that $\theta_{g}$ is the viewing angle of the grain or the angle from the axis of symmetry ($x$-axis) to the observer. Since we assume the grain is azimuthally aligned in the midplane of the disk, $\theta_{g}$ simply depends on the inclination and azimuth of the disk. We can obtain the relation from geometrical arguments:

which gives the polarization in the grain frame through Eq. (16). Rotating into the lab frame will give us $q^{\prime}$ and $u^{\prime}$ from both scattering and thermal emission:

(see Appendix B for details).

In Fig. 11, we compare the predicted azimuthal profiles for $q^{\prime}$ and $u^{\prime}$ based on the values of $S$ and $T_{0}$ determined from the polarization data on the major and minor axes (and shown in Fig. 10) with the observed azimuthal variations. Although there are some discrepancies (e.g., around $\phi_{d}\sim 80^{\circ}$ for $q^{\prime}$ in Band 3 and $\sim 45^{\circ}$ for $u^{\prime}$ in Bands 3 and 6), the predicted profiles match the observed ones, especially the folded data, remarkable well overall, which adds support to our empirical method of decomposing the observed polarization into scattering and thermal components.

In Section 4.4, we obtained an empirical measurement of scattering polarization $S$ and thermal polarization $T_{0}$ for one radius. In principle, one can measure $S$ and $T_{0}$ at each point across radius and infer the property of grains as a function of radius as the resolution may allow. However, for the HL Tau data, there is roughly only one beam across the disk minor axis, which limits this empirical technique to just one independent data point at $\sim 100$au.

Nevertheless, as a consistency check, we solve for $S$ and $T_{0}$ at different radii in a single image under the influence of the finite beam size. We choose an inner radius of $50$au which gives $\sim 4$ beams around the azimuth of the Band 3 image and the outer radius is limited by sensitivity. In Fig. 12 we show the radial profile at each wavelength. As a comparison, we show the beam sizes projected along the disk minor axis, FWHM$/\cos i$, since we are mainly limited by the resolution along the minor axis for independent radial points.

We find that $S$ is roughly constant of radius for each wavelength. $T_{0}$ decreases towards the inner radius likely due to beam convolution which can average out the azimuthal variation (see Fig. 8). Similar to Fig. 10, one can see immediately that the thermal component $T_{0}$ dominates the scattering component $S$ everywhere for Band 3 (the green solid and dashed lines in Fig. 12), especially at large radii where effects of beam averaging is less. The opposite is true for Band 7 with larger $S$ than $T_{0}$ and at Band 6, the two components are comparable.

From Section 4.4, we predicted the azimuthal profiles of $q^{\prime}$ and $u^{\prime}$ based on $S$ and $T_{0}$. Since we can solve the azimuthal profile at each radius, we can produce an image of $q^{\prime}$ and $u^{\prime}$. To retrieve Stokes $Q^{\prime}$ and $U^{\prime}$, we simply multiply the solved $q^{\prime}$ and $u^{\prime}$ by the observed Stokes $I^{\prime}$ (Eq. (4)). In Fig. 13, we show the reconstructed polarization images plotted in a similar manner as the folded images in Fig. 9. The hole in the reconstructed images are simply because of our inner cut of radius. Indeed, the polarization images are quantitatively similar to the folded data.

One notable difference is the polarized intensity at Band 3. In Fig. 9a, the polarized intensity has four peaks instead of expected two peaks along the minor axis in Fig. 13a. This is likely because of the intrinsic asymmetry at Band 3. Shown in Fig. 6a, there is a very strong asymmetric peak of the polarized intensity that is not along the minor axis. Even folding the data could not cancel out the asymmetry which leaves a footprint in the folded image as a bright spot slightly offset from the minor axis. Since the folded image is symmetric across the disk major and minor axes, the bright spot is imprinted across the four quadrants. The cause of the intrinsic asymmetry is unclear.

Despite the above difference, the simple analytic decomposition of the observed polarization into the scattering and thermal components (Fig. 12) and the predicted multiwavelength polarization distributions based on the decomposition and simple geometric effects match the folded observational data remarkable well (compare Fig. 13 and Fig. 9). In particular, as the observing frequency increases, there is a transition from a more azimuthal pattern to a more pattern along the minor axis in the polarization orientation and from an azimuthal distribution of polarization fraction that is higher along the minor axis than along the major axis to the opposite (i.e., the dumbbell-shaped distribution along the major axis) in both the observation data and the model prediction. Fig. 12 quantified the degree to which the Band 3 polarization is dominated by thermal emission and the Band 7 polarization is dominated by scattering. It shows that the crossover occurs near Band 6, where the thermal and scattering components are comparable, which is consistent with the conclusion of Stephens et al. (2017) based on a simple mix of a uniform and an azimuthal polarization patterns.

The analytic decomposition based purely on the observational data and geometry without prior knowledge of the disk or dust properties (except for the assumptions of axisymmetry and the geometrically thin limit) demonstrated in Section 4.4 can be a unique tool for analyzing resolved polarization images. Plane-parallel slab calculations have been used to fit the radial properties of disks across wavelength without assuming any power-law if the multiwavelength images are resolved (e.g. Carrasco-González et al., 2019; Yen & Gu, 2020; Sierra et al., 2021). If polarization images can be resolved across wavelength, the empirical decomposition method or plane-parallel slab calculations can be used to determine even more properties of grains, such as the degree of alignment or scattering properties of the grain.

Based on Fig. 10, we can make rough predictions for polarization at other wavelengths. The remaining ALMA bands that can provide continuum polarization are Bands 1, 4, and 5. For Bands 4 and 5 (or $\lambda=1.5$ and $2.1$ mm respectively), we expect that the the scattering polarization $S$ should be in between $0.5\%$ to $0.7\%$ unless either of the bands happen to form a peak in polarization when the optical depth is of order unity (Fig. 4). We can expect that the thermal polarization will be greater than that at Band 6 and less than that of Band 3. Furthermore, we expect the polarization images of Band 4 and 5 to appear similar to a transition between Bands 3 and 6. Preliminary Bands 4 and 5 data suggests that the predictions seem consistent with the observations (private communication).

With the incorporation of Band 1 ($\lambda=6$ to $8.5$mm) on ALMA, the longer wavelength can potentially probe the HL Tau disk in the optically thin limit thereby avoiding the effects of scattering altogether. Based on previous works (Yang et al., 2019; Mori & Kataoka, 2021) and our results, one can conclude that the grains responsible for the thermal polarization are most likely azimuthally aligned prolate grains in the HL Tau disk. However, it is still important to determine the grain alignment directions with as little contamination from scattering as possible. The polarization image at an optically thinner wavelength should better reveal the intrinsic alignment field which can serve as a more robust test for azimuthally aligned prolate grains.

Since the optical depth at Band 1 should be lower than that at Band 3, the scattering polarization $S$ should be less than the $\sim 0.5\%$ at Band 3 and the thermal polarization should increase. Since the Band 3 image is similar to the case of $\tau_{0}=0.05$ in Fig. 8, we expect that the polarization image should appear similar to the case of $\tau_{0}=0.01$ with a more obvious contrast between the higher polarization fraction region along the minor axis and the lower polarization fraction region along the major axis (Fig. 8m) and, more strikingly, a vertical bar-like morphology for the polarized intensity after beam convolution (Fig. 8a).

From Section 4, we have showed that scattering mainly provides a positive Stokes $Q^{\prime}$ in the principle frame. The Stokes $U^{\prime}$ on the other hand, is not heavily affected by scattering (Fig. 11) and retains the underlying thermal polarization (see also Section 3.3). A natural result is that along the major and minor axis, we should have Stokes $U^{\prime}=0$ if the prolate grains are perfectly aligned in the toroidal direction.

Fig. 15 shows the Stokes $U^{\prime}$ for Bands 3, 6, and 7 without folding. We further plot the lines, $x^{\prime}=0$ and $y^{\prime}=0$, which is where Stokes $U^{\prime}$ should equal zero for an axisymmetric disk. Along the minor axis of the Band 3 Stokes $U^{\prime}$, the path where Stokes $U^{\prime}=0$ has a slight clock-wise tilt from $x^{\prime}=0$ and less of a tilt from $y^{\prime}=0$. For Band 6, Stokes $U^{\prime}=0$ also has a slight clock-wise tilt from $y^{\prime}=0$. These deviations can be due to an elongated beam that is not parallel to the major or minor axis.⁹⁹9The offset could also be due to uncertainty of the position angle of the disk major axis in general, but it is unlikely the case for HL Tau given the tight constraint from available high angular resolution images (ALMA Partnership et al., 2015). However, we demonstrate that it can also be due to a spiral alignment.

For illustration, we use the reference model in Section 4. We let

where $\delta$ is the tilt between the $y$-axis of the slab and the radial direction of the disk in the counterclockwise direction when the disk is viewed face-on.

In the optically thin and face-on case (left column of Fig. 16), the polarization angles directly trace the long axis of the grain which is tilted away from the perfect azimuthal pattern when $\delta=15^{\circ}$ (Fig. 16m). The polarization fraction in the center is low because of beam averaging. The Stokes $Q^{\prime}$ and $U^{\prime}$ is rotated in the clock-wise direction as seen in the plane-of-sky (Fig. 16e,i). The rotation can be easily understood by the following. Consider a grain at some positive $x^{\prime}$ (left of the vertical axis as defined in Section 4 or $\phi_{d}=0^{\circ}$) that is azimuthally aligned. The Stokes $Q^{\prime}$ is positive and Stokes $U^{\prime}$ should be $0$ because the projected long axis of the grain is in the vertical direction. With the extra spiral angle $\delta=15^{\circ}$, the grain rotates counter-clockwise and results in a positive Stokes $U^{\prime}$. The rotation applies to each location in the disk and thus Stokes $Q^{\prime}$ and $U^{\prime}$ are effectively rotated.

When the optically thin disk is inclined (second column of Fig. 16), the azimuthal variation in $p_{l}$ appears. Similar to the azimuthally aligned case (left column of Fig. 8), prolate grains along the disk major axis is viewed more pole-on, while those along the disk minor axis is viewed more edge-on. However, the largest polarized intensity and $p_{l}$ are no longer along the disk minor axis, but slightly offset clockwise due to $\delta=15^{\circ}$.

If the optical depth increases, such as $\tau_{0}=3$ of the third column of Fig. 16, scattering polarization dominates and provides a positive Stokes $Q^{\prime}$ (Fig. 16g). The polarization angles are mostly parallel to the disk minor axis near the center and the effects of spiral alignment is only visible in the outer optically thin regions as seen in Fig. 16p. The regions where Stokes $U^{\prime}=0$ remains offset from the disk major and minor axes (compare Fig. 16j and k). In contrast, if $\delta=0^{\circ}$ (rightmost column of Fig. 16), the lines where Stokes $U^{\prime}=0$ is completely along the disk major and minor axes (Fig. 16l).

With $\delta=15^{\circ}$ (third column of Fig. 16), all the polarization properties show modest differences to the perfect azimuthal alignment case (rightmost column of Fig. 16). However, showing $x^{\prime}=0$ and $y^{\prime}=0$ as guidelines on Stokes $U^{\prime}$ makes it easier to identify the deviation. Since scattering does not affect Stokes $U^{\prime}$ when the projected prolate grain is parallel or perpendicular to $\hat{\theta}$, Stokes $U^{\prime}$ will be 0 even in the optically thick regions where scattering dominates.

At face value, Stokes $U^{\prime}$ of Band 3 and 6 both exhibit a slight tilt that appears to suggest $\delta$ much less than the chosen $15^{\circ}$ example in Fig. 16 (we find that HL Tau is broadly consistent with $\delta\sim 5^{\circ}$). However, Band 7 does not resemble Bands 3 and 6 and also does not match any of the spiral alignment cases. In particular, the negative Stokes $U^{\prime}$ in the lower right quadrant almost vanished making it difficult to explain with just spirally aligned grains. Whether or not HL Tau does have spiral alignment will still need numerical confirmation, especially considering the elongated beam, which is beyond the scope of this paper. Higher angular resolution will help the simple diagnostic, since the lines where Stokes $U^{\prime}=0$ can be better resolved and differentiated from the beam convolution effects. The use of plane-parallel slab models as demonstrated in this section can constrain the alignment direction of the grains under the influence of scattering and may help identify the alignment mechanism of disks.

Note that folding across the minor axis and major axis described in Section 4.3 do not hold for spirally aligned grains. Instead, the image has rotational symmetry of order 2 if the disk is both geometrically thin and axisymmetric.

There are other sources that are similar to HL Tau which exhibit the scattering polarization pattern at the shorter wavelength, typically Band 7, and an azimuthal pattern at the longer wavelength, e.g., Band 3. The DG Tau disk observed at Band 3 (Harrison et al., 2019) shows a clear azimuthal pattern of polarized intensity with polarization angles that look elliptical. There is a bar of polarized intensity along the major axis with polarization angles parallel to the minor axis of the disk. At Band 7, the polarization angles are largely parallel to the disk minor axis (Bacciotti et al., 2018). The wavelength behavior matches what we expect from an increase of optical depth from Band 3 to 7. The Band 3 polarized intensity is roughly similar to the $\tau_{0}=0.05$ case of the reference convolved model image in Fig. 8 with a slightly stronger bar which would suggest a $\tau_{0}$ in between $0.05$ and $0.5$. In contrast, HL Tau does not have a bar at Band 3 which would suggest DG Tau is more optically thick than HL Tau at Band 3. Indeed at Band 7, the polarized intensity peak is offset from the Stokes $I$ peak which is expected if the dust is not settled and the disk is optically thick (Yang et al., 2017). Directly using the major and minor axis polarization fraction at Band 3, we have at $r=0.4\arcsec$, $q_{\text{maj}}\sim 4.5\%$ and $q_{\text{min}}\sim-6\%$. Applying Eq. (17) and (18), we get $S\sim 0.1\%$ and $T_{0}\sim 6.1\%$. Band 7 is notably non-settled which makes the technique inapplicable.

Haro 6-13 is a Class II disk with features similar to DG Tau in its Band 3 polarization image (Harrison et al., 2019). The outer region of the disk exhibits an azimuthal pattern with a bar of polarized intensity along the major axis. In the outer region, the polarization fraction along the minor axis is larger than that along the major axis which is expected from azimuthally aligned prolate grains modulated by little scattering contribution. The Band 7 image only has an obvious polarized intensity near the center with polarization parallel to the disk minor axis (Harrison et al. in prep.). The wavelength behavior is similar to HL Tau and can be explained by scattering of azimuthally aligned prolate grains.

AS 209 is another Class II source observed at two wavelengths. The Band 7 image has an outer azimuthal pattern with a polarized intensity bar across the major axis (Mori et al., 2019). In the region with an azimuthal pattern, the polarization fraction along the minor axis is just slightly larger than the polarization fraction along the major axis (see Fig. 4 in Mori et al. 2019). This would also suggest the presence of azimuthally aligned prolate grains with scattering. However, at the longer wavelength of Band 6, the central polarized intensity bar appears larger (Harrison et al., 2021) when we would expect the central bar to diminish and the thermal polarization to grow at the longer wavelength. The cause may be due to different beam sizes. The Band 7 image has a beam of $\sim 0.8\arcsec$ while the Band 6 images has a beam of $\sim 1.3\arcsec$. The larger beam of the longer wavelength image may have smeared the stronger polarized intensity due to scattering at the center to the outskirts of the disk where the thermal polarized intensity is low. Higher angular resolution observations of Band 6 will clarify if the contradiction from expectation is due to beam size or a scenario that is not explained by scattering of azimuthally aligned prolate grains like HL Tau.

As pointed out in Sadavoy et al. (2019), the Band 6 image of the Class I IRS 63 resembles Band 6 of HL Tau in that the polarization in the central region is largely parallel to the disk minor axis suggesting polarization due to scattering, while the outer regions are largely elliptical. Direct comparison with a perfect elliptical pattern shows clear deviations in which the polarization angles have a tendency to follow the direction of the disk minor axis. The deviation is expected if scattering polarization is comparable to the thermal polarization and the extra positive Stokes $Q^{\prime}$ pushes the polarization angle to follow the disk minor axis (see Fig. 5) as we found with Band 6 of HL Tau.

The accretion disk of the massive protostar, GGD27 MM1, that powers the HH 80-81 jet also has resolved polarized dust continuum at 1.14 mm (Girart et al., 2018). Similar to the previous sources, the polarization fraction shows two distinct regions: the inner region where most of the polarization is parallel to the disk minor axis with $\sim 0.6\%$ polarization and the outer region where the polarization is elliptical with $\sim 6\%$ polarization. The inner region exhibits a near/far-side asymmetry in the polarized intensity (and polarization fraction) which is expected from scattering if the disk is optically thick and geometrically thick (Yang et al., 2017). In the outer region, the minor axis polarization is larger than the major axis polarization which is expected from azimuthally aligned prolate grains. Given that the disk is likely geometrically thick, our plane-parallel model cannot capture the near/far-side asymmetry nor the radiation anisotropy. Nevertheless, the transition from the inner region to the outer region can also be explained based on a change in optical depth as our Fig. 8 demonstrates.

As a counter example, TMC-1A is a Class I source that cannot be explained by scattering azimuthally aligned prolate grains. The 1.3 mm polarization at the disk center is $\sim 0.7\%$ and parallel to the disk minor axis, while the polarization in the outer region is $\sim 10\%$ and mostly radial (Aso et al., 2021). The radial polarization pattern is better explained by azimuthally aligned oblate grains (Cho & Lazarian, 2007), and the central region is most likely due to scattering. At first glance, the decrease to low polarization could be due to dichroic extinction from large optical depth (like the previous cases). However, if the oblate grains are completely aligned azimuthally throughout the disk, we would not expect a depolarized region in between the outer region and the center. Along the disk minor axis, polarization of the outer region from oblate grains share the same polarization angle as that from the central scattering region. Thus, there is no cancellation involved to create the depolarization. Instead, the polarization fraction should transition smoothly (analogous to the $\phi=0^{\circ}$ curve of Fig. 4). The reasons behind a depolarized region is likely related to a true change in the grain properties or alignment efficiency. For example, the grains could be less aligned within the depolarized region.

The above discussions lead us to believe that azimuthally aligned scattering grains that appear to explain the multiwavelength polarization of the HL Tau disk may also exist in disks around other low-mass stars and possibly even massive stars. However, not all disks with evidence of aligned grains and scattering can be explained by the same way. Why grains in the HL Tau and other disks are aligned with their long axes along the azimuthal direction is an interesting question that deserves further investigation.