# Evidence of grain alignment by magnetically enhanced radiative torques from
multi-wavelength dust polarization modeling of HL Tau

Thang Nguyen Tat Department of Astrophysics, Vietnam National Space Center,
Vietnam Academy of Science and Technology, 18 Hoang Quoc Viet, Hanoi, Vietnam
Simons Astrophysics Group at IFIRSE (SAGI)/Institute For Interdisciplinary
Research in Science and Education (IFIRSE), ICISE, 07 Science Avenue, Ghenh
Rang Ward, 55121 Quy Nhon City, Binh Dinh Province, Vietnam  Pham Ngoc Diep
Department of Astrophysics, Vietnam National Space Center, Vietnam Academy of
Science and Technology, 18 Hoang Quoc Viet, Hanoi, Vietnam  Graduate
University of Science and Technology, Vietnam Academy of Science and
Technology, 18 Hoang Quoc Viet, Hanoi, Vietnam  Thiem Hoang Korea Astronomy
and Space Science Institute, Daejeon 34055, Republic of Korea  Korea
University of Science and Technology, 217 Gajeong-ro, Yuseong-gu, Daejeon,
34113, Republic of Korea  Le Ngoc Tram Max-Planck-Institut für
Radioastronomie, Auf dem Hügel 69, 53121, Bonn, Germany  Nguyen Bich Ngoc
Department of Astrophysics, Vietnam National Space Center, Vietnam Academy of
Science and Technology, 18 Hoang Quoc Viet, Hanoi, Vietnam  Graduate
University of Science and Technology, Vietnam Academy of Science and
Technology, 18 Hoang Quoc Viet, Hanoi, Vietnam  Nguyen Thi Phuong Department
of Astrophysics, Vietnam National Space Center, Vietnam Academy of Science and
Technology, 18 Hoang Quoc Viet, Hanoi, Vietnam  Simons Astrophysics Group at
IFIRSE (SAGI)/Institute For Interdisciplinary Research in Science and
Education (IFIRSE), ICISE, 07 Science Avenue, Ghenh Rang Ward, 55121 Quy Nhon
City, Binh Dinh Province, Vietnam  Bao Truong Korea Astronomy and Space
Science Institute, Daejeon 34055, Republic of Korea  Korea University of
Science and Technology, 217 Gajeong-ro, Yuseong-gu, Daejeon, 34113, Republic
of Korea

###### Abstract

Atacama Large Millimeter/Submillimeter Array (ALMA) has revolutionized the
field of dust polarization in protoplanetary disks across multiple
wavelengths. Previous observations and empirical modeling suggested multiple
mechanisms of dust polarization toward HL Tau, including grain alignment and
dust scattering. However, a detailed modeling of dust polarization based on
grain alignment physics is not yet available. Here, using our updated POLARIS
code, we perform numerical modeling of dust polarization arising from both
grain alignment by Magnetically Enhanced Radiative Torque (MRAT) mechanism and
self-scattering to reproduce the HL Tau polarization observed at three
wavelengths 0.87, 1.3, and 3.1 mm. Our modeling results show that the observed
multi-wavelength polarization could be reproduced only when large grains
contain embedded iron inclusions and those with slow internal relaxation must
have wrong internal alignment (i.e., the grain’s major axis parallel to its
angular momentum). The abundance of iron embedded inside grains in the form of
clusters is constrained to be $\gtrsim 16$%, and the number of iron atoms per
cluster is $N_{\rm cl}\sim 2\times 10^{3}$. Maximum grain sizes probed at
wavelengths $\lambda$ = 0.87, 1.3, and 3.1 mm are constrained at $\sim$ 60,
90, and 130$\,\mu$m, respectively. Assuming a dust differential settling
effect with grain sizes from the constraint gives the value of maximum grain
size at the disk mid-plane to be millimeter-scaled.

Protoplanetary disks; Polarimetry; Radio astronomy; Circumstellar dust;
Magnetic field

††journal: ApJ††software: POLARIS

##  1 Introduction

Magnetic fields are thought to play a crucial role in protoplanetary disks in
driving magnetorotational instability (Velikhov, 1959; Balbus & Hawley, 1991),
disk formation (e.g., magnetic braking and angular momentum transport) (Balbus
& Hawley, 1991, 1998; Hawley et al., 1994), protostellar outflow launching
(Blandford & Payne, 1982; Sharkawi et al., 2021), and planet formation
(Quillen & Trilling, 1998; Fendt, 2003; Johansen, 2009). Polarized emission
from aligned grains has been known as the leading channel to investigate the
magnetic fields in the interstellar medium (ISM) and molecular clouds (MCs).
Therefore, dust polarization is expected to be a powerful tool in studying
magnetic fields (Stephens et al., 2014) and dust properties (Kataoka et al.,
2015) in disks. Since the emergence of the Atacama Large
Millimeter/submillimeter Array (ALMA), the field of (sub)millimeter-wavelength
disk polarization has progressed significantly with the polarization
signatures being observed to incredible detail and accuracy (Stephens et al.,
2017; Kataoka et al., 2017; Cox et al., 2018; Harrison et al., 2019; Sadavoy
et al., 2019; Tang et al., 2023).

HL Tau is a perfect example to illustrate the evolutionary process of the
field. It is a young protoplanetary disk (SED class I/II, age $\sim\,$1 Myr)
(White & Hillenbrand, 2004) located in the Taurus molecular cloud, which is
$\sim\,140\,$pc from Earth (Rebull et al., 2004). It is one of the earliest
disks observed with polarimetry by Combined Array for Research in Millimeter-
wave Astronomy (CARMA) and Submillimeter Array (SMA). Its polarized emission
was used to study the disk’s magnetic fields, supposing grains are
magnetically aligned (Stephens et al., 2014). However, the study concluded
that a simple magnetic field structure could not fully explain the detected
polarization. Coming to the age of ALMA observations, the disk has been
revealed to have a transition of polarization morphology from azimuthal to
parallel to the disk’s minor axis when probing at longer to shorter
wavelengths (Kataoka et al., 2017; Stephens et al., 2017). These new
observations essentially raise questions about the validity of tracing
magnetic fields by dust polarization within the disk.

The complex polarization patterns indicate the possibility of multiple
mechanisms involved in producing the dust-polarized emission of HL Tau at
(sub)millimeter wavelengths. Indeed, previous studies have suggested that a
superposition of self-scattering and thermal emission from azimuthally-aligned
prolate grains can explain such observed polarized emission (Yang et al.,
2019; Mori & Kataoka, 2021; Lin et al., 2022). Self-scattering is caused by
dust grains with sizes comparable to the observed wavelengths scattering and
polarizing incident lights from the thermal emission of other grains (Kataoka
et al., 2015). Grain alignment, on the other hand, has elongated grains being
systematically aligned with a preferred direction and emitting polarized
thermal electromagnetic waves. Nevertheless, which mechanisms exactly induce
grain alignment in such a way is still a mystery since most of the proposed
mechanisms (i.e., magnetic alignment, radiative alignment, aerodynamic
alignment) fail to interpret the observed disk polarization (Tazaki et al.,
2017; Yang et al., 2019; Mori & Kataoka, 2021).

Grain alignment in the ISM and MCs has been studied with great success and
made dust polarization an invaluable tool to trace magnetic fields within the
environment (Andersson et al., 2015; Pattle & Fissel, 2019). However, its
applicability to the protostellar environment is complex due to the presence
of very large grains (VLGs, grain size $a\,\gtrsim\,\text{10}\,\mu$m) and
significant gas randomization effect. Recently, Hoang et al. (2022) considered
a more realistic grain alignment model within such an environment in contrast
to previous modelings where only perfect alignment is included. The process of
grain alignment therein consists of the alignment of the axis of maximum
moment of inertia ($\boldsymbol{\hat{a_{1}}}$) with the angular momentum
($\boldsymbol{J}$), so-called internal alignment (IA); and the alignment of
$\boldsymbol{J}$ with a preferred direction, so-called external alignment
(EA). IA can be achieved by internal relaxation (INR), e.g., Barnett
relaxation (Roberge et al., 1993) and inelastic relaxation (Molina et al.,
2003). EA, on the other hand, can be accomplished by grain precession (e.g.,
Larmor precession, radiative precession, and mechanical precession). Grains
achieving EA with Larmor precession, radiative precession, or mechanical
precession have $\boldsymbol{J}$ aligned along the direction of the magnetic
field ($\boldsymbol{B}$), radiation beam ($\boldsymbol{k}$), or gas flow
($\boldsymbol{\nu}$). The mechanisms behind such alignments are called
magnetic alignment (Hoang, 2022), radiative alignment ($k$-RAT; Lazarian &
Hoang, 2007a), or mechanical alignment ($v$-MET; Lazarian & Hoang, 2007b),
respectively.

Recently, Giang et al. (2023) incorporated such grain alignment model (i.e.,
realistic IA and EA model) into POLARIS (Reissl et al., 2016) taking into
account only the grain alignment mechanism of Magnetically-enhanced Radiative
Torque (MRAT) (Hoang & Lazarian, 2016) which was specifically applied to study
thermal dust polarization towards protostellar cores. With this mechanism,
grains can be aligned with the disk magnetic fields by radiative torques
(RATs; Dolginov & Mitrofanov, 1976; Draine & Weingartner, 1997; Hoang &
Lazarian, 2008) with enhanced efficiencies by (super)paramagnetic relaxation
(Davis & Greenstein, 1951; Hoang & Lazarian, 2008; Hoang & Lazarian, 2016).

In this paper, our main goal is to test the validity of grain alignment with
the MRAT mechanism and put a strong constraint on grain magnetic properties as
well as grain growth in HL Tau. To do so, we perform multi-wavelength modeling
of disk polarization using the updated POLARIS (Giang et al., 2023) for three
well-studied ALMA Bands 3, 6, and 7 corresponding to the wavelengths $\lambda$
= 3.1, 1.3, and 0.87 mm.

The structure of this paper is as follows. In Section 2, we describe the disk
and dust model used in our modelling. Subsequently, our main results including
synthetic polarization from thermal aligned dust emission with MRAT mechanism;
and from the mixture of grain alignment and self-scattering are shown in
Section 3.1 and 3.2, respectively. We further discuss the implications, grain
alignment properties, and grain growth in Section 4. Finally, the main
findings are summarized in Section 5.

Figure 1: Left to right: Stokes $I$ intensity overlaid with ALMA polarimetric
observations of HL Tau at Band 3, Band 6, and Band 7. The line segments
represent the polarization vectors and are set to follow the scale of
polarization degree. The ones with polarized intensity below signal-to-noise
of 3 are not plotted. A reference vector of 1$\%$ is shown in the lower-left
corner of each panel. The contours represent the intensity levels of (10, 25,
50, 100, 200) $\times$ the noise level, which is $\sigma_{\rm I}$ = 58, 170,
and 409$\,\mu$Jy/beam at $\lambda$ = 3.1, 1.3, and 0.87 mm, respectively. The
beam sizes and their position angles are ($0.45\arcsec\times 0.29\arcsec$,
$-12.8^{\circ}$), ($0.33\arcsec\times 0.24\arcsec$, $-53.4^{\circ}$), and
($0.2\arcsec\times 0.13\arcsec$, $-80.9^{\circ}$), respectively.

##  2 Dust and Disk Models

The outline of our model is briefly described as follows. Firstly, we
construct a dust model and a best-fitted disk model to the Stokes $I$
emissions obtained by ALMA at the three wavelengths $\lambda$ = 3.1, 1.3, and
0.87 mm. Then, by separately measuring direct thermal emission from grains
aligned by MRAT and scattered emission from self-scattering, we combine their
Stokes parameters to obtain the combined linear polarization from the dust
emission.

Figure 1 shows ALMA Stokes $I$ intensity overlaid with polarization vectors of
HL Tau disk at $\lambda$ = 3.1, 1.3, and 0.87 mm. The notable features of the
observations are (1) the transition of the polarization patterns from
azimuthal to parallel to the disk’s minor axis from 3.1 to $0.87\,\rm mm$ and
(2) the polarization degree is generally at the level of $\sim$1% for the
whole three wavelengths.

###  2.1 Dust Model

We adopt a dust model with compositions from Birnstiel et al. (2018), where
dust grains are a mixture of astronomical silicate, troilite, refractory
organics, and water ice. The refractive index used in the calculation is from
Draine (2003) for astronomical silicate, Henning & Stognienko (1996) for
troilite, and Warren & Brandt (2008) for refractory organics and water ice.
Assuming compact grains, the mean refractive index of the mixture is then
calculated with the effective medium theory using Bruggeman rule (Bohren &
Huffman, 1998).

The grain sizes are assumed to follow a power-law distribution with a power
index $q=-3.5$ (Mathis et al., 1977), a maximum grain size $a_{\rm max}$, and
a fixed minimum grain size $a_{\rm min}$ = 0.05 $\mu$m. We assume the grain
population to be the same throughout the disk. As discussed in Lin et al.
(2022), the value of $a_{\rm max}$ does not converge for the three ALMA
wavelengths (i.e., there is no single value could satisfy the polarization
constraint at the three wavelengths at once). Therefore, we leave $a_{\rm
max}$ as a wavelength-dependent parameter. We assume the Rayleigh regime (Lee
& Draine, 1985) for the polarization calculation. This approximation holds
when the size parameter $x\,=\,2\pi a/\lambda\,<\,1$. Indeed, later on, this
is shown to be the case with the maximum grain sizes found to be $\sim$60, 90,
and 130$\,\mu$m for the three probing wavelengths 0.87, 1.3, and 3.1 mm,
respectively.

Table 1: Key Parameters of the Dust Model Parameter | Denotation | Value
---|---|---
Dust mixture | - | silicate, organics, water ice, troilite
Grain size distribution | dn/da | $Ca^{-3.5}$
Minimum grain size | $a_{\rm min}$ | 0.05$\,\mu$m
Maximum grain size | $a_{\rm max}$ |  $\lambda$ \- dependent
Grain aspect ratio | $s$ | 1/2 (ray-tracing) or 1 (MCRT)
Iron fraction of PM grains | $f_{\rm p}$ | 0.1
Iron cluster volume filling factor of SPM grains | $\phi_{\rm sp}$ | 0.1
Number of iron atoms per cluster of SPM grains | $N_{\rm cl}$ | [ 20, $10^{5}$ ]
IA efficiency at low-$J$ and slow INR | $Q_{\rm X,low-J}$ |  $-0.25$ (wrong IA) or 0.2 (right IA)
IA efficiency at high-$J$ and slow INR | $Q_{\rm X,high-J}$ | 0.2

We assume two different grain shapes for direct thermal dust emission ray
tracing and Monte-Carlo radiative transfer (MCRT) calculation (i.e., radiation
field, radiative heating, and dust scattering calculations). For dust emission
ray tracing, which is used to trace aligned grain emission, oblate grains with
aspect ratio $s$ ($0<s<1$) are assumed. Note that grain shape has a strong
effect on the degree of polarization. Grains with smaller $s$ are expected to
produce higher polarization degrees. However, since the exact grain geometry
in the disk environment is not yet constrained, we choose an average value of
$s=1/2$ for simplicity. The grains’ optical properties including extinction
cross-sections and RAT efficiency are then pre-calculated with DDSCAT (Draine
& Flatau, 1994). For the MCRT calculation, we adopt spherical grains (i.e.,
$s\,=\,$1), whose dust optical properties can be analytically calculated with
Mie theory (Wriedt, 2012). Though the two-grain shapes are not consistent, we
neglect the effect of scattering with non-spherical grains for simplicity. As
pointed out by Kirchschlager & Bertrang (2020), assuming perfect compact
spherical grains yields quantitative errors compared to non-spherical grains.
However, the errors are negligible in the Rayleigh regime (i.e., spherical
grains are still an acceptable representation).

Self-scattering-induced polarization is achieved by treating dust as photon
emitters in dust scattering radiative transfer simulation. The controlling
parameter of scattering is the maximum grain size, $a_{\rm max}$. Self-
scattering polarization is most effective when $a_{\rm
max}\,\sim\,\lambda/2\pi$ (Kataoka et al., 2015).

In the present work, grain alignment is considered to be induced only by MRAT.
Under the framework of MRAT, grains can efficiently be aligned with magnetic
fields owing to the enhanced efficiency by RATs and (super)paramagnetic
relaxation due to iron inclusion within the grains (Lazarian & Hoang, 2008;
Hoang & Lazarian, 2016). Recent observations show more evidence for iron
inclusion in the form of clusters in different environments such as in
protostellar core (Giang et al., 2023), filament (Ngoc et al., 2023), and the
envelope of evolved stars (Truong et al., 2023). RATs can spin grains up to
suprathermal rotation (so-called high-$J$ attractors) or de-spin them down to
thermal rotation (so-called low-$J$ attractors). Grains with iron inclusion
are dubbed paramagnetic material (PM) or superparamagnetic material (SPM).
Considering PM grains, we assign the value of the fraction of iron atoms
diffusely distributed within grains as $f_{\rm p}=0.1$. For SPM grains, we fix
the volume filling factor of iron clusters to $\phi_{\rm sp}$ = 0.1, which
corresponds to an iron abundance of $\sim$30% presented in the form of iron
clusters111In reality, $\phi_{\rm sp}$ can be as large as 0.3, which
corresponds to an iron abundance of $\sim\,100\,\%$ presented in the form of
iron clusters (Hoang & Lazarian, 2016) and vary the number of iron atoms per
cluster $N_{\rm cl}$ to describe different degree of grain alignment
enhancement due to superparamagnetic relaxation. The possible value of $N_{\rm
cl}$ spans from $\sim$ 20 to $10^{5}$ (Jones & Spitzer, 1967).

Note that for INR of grains, although the case of RATs quickly bringing
$\boldsymbol{J}$ to align with $\boldsymbol{\hat{a_{1}}}$ (fast INR, Lazarian
& Hoang, 2021) has been well-studied, the case of INR timescale being much
longer than gas damping timescale (slow INR) is not as well-defined (Hoang et
al., 2022). Hence, in this study, for the case of slow INR, we fix the values
of IA efficiency $Q_{\rm X}$ as follows. When grains are aligned at high-$J$
attractors, the IA efficiency is $Q_{\rm X,high-J}$ = 0.2. On the other hand,
when grains are at low-$J$ attractors, they can be aligned at either right or
wrong IA (Hoang & Lazarian, 2009). In the case of right IA, $\boldsymbol{J}$
is aligned parallel to $\boldsymbol{\hat{a_{1}}}$, which coincides with the
minor axis of the ellipsoidal grains. In the case of wrong IA, their alignment
direction is perpendicular instead. When $Q_{\rm X,low-J}=-0.5$, grains are
defined to achieve perfect wrong IA (Giang et al., 2023). We set $Q_{\rm
X,low-J}=0.2$ for the case of grains aligning at right IA, and $Q_{\rm
X,low-J}=-0.25$ at wrong IA. The choice of the values of $Q_{\rm X,high-J}$
and $Q_{\rm X,low-J}$ at wrong IA is taken as the average value of the
constraint which will be given in Section 4.3. The value of $Q_{\rm X,low-J}$
at right IA is set to be similar to $Q_{\rm X,high-J}$ based on the assumption
that $\boldsymbol{J}$ would have a negligible effect on the grains’ IA
efficiency when they are aligned with slow INR. Throughout the paper, we will
call wrong IA/right IA for short in referring to the case of grains aligned
with wrong or right IA when they are at low-$J$ attractors and aligned with
slow INR. The procedure to calculate the alignment degree of with magnetic
fields in POLARIS is nicely summarized in Figure 2 of Giang et al. (2023).

A summary of the critical parameters used in our dust model is given in Table
1.

Figure 2: From left to right: Stokes $I$ emissions of model, observation, and
residual (observation$-$model) at 3.1 mm (top), 1.3 mm (middle), and 0.87 mm
(bottom), respectively. The left and middle panels show contours identical to
that in Figure 1. The right panels show residual contours at
$(-35,-5,5,35)\times\sigma_{\rm I}$

###  2.2 Disk Model

In this section, we construct an axisymmetric disk model that best reproduces
Stokes $I$ images of the polarization observations at 3.1 mm (Kataoka et al.,
2017; Stephens et al., 2017), 1.3 mm (Stephens et al., 2017), and 0.87 mm (Lin
et al., 2023a).

First, we assume a single radiation source to be the protostar at the disk
center. The protostar has the effective temperature $T_{\star}$ = 4000 K and
radius $R_{\star}=7R_{\odot}$, which corresponds to a luminosity
$L_{\star}\sim$$11L_{\odot}$ (Men’shchikov et al., 1999; Pinte et al., 2016;
Liu et al., 2017). The radiation field and dust temperature within the disk
are then calculated using radiative transfer calculation. Since temperature is
a strong function of grain size, we fix $a_{\rm max}$ = 100 ${\mu}$m for
temperature calculation to fix the disk’s temperature profile.

Next, we build the surface density profile for the disk. It is commonly
estimated by using the spectral energy distribution (SED) fitting. However, it
is only simple when the temperature profile is fixed or can be analytically
approximated. In our study, the main source of the disk temperature is
radiative heating requiring radiative transfer calculation which is
computationally costly. Thus, for simplicity, we assume a surface density
profile taken from Kwon et al. (2011),

$\displaystyle\Sigma\propto\left(\frac{R}{R_{c}}\right)^{-\gamma}\text{exp}\left[\left(\frac{R}{R_{\text{c}}}\right)^{2-\gamma}\right]\text{,}$
(1)

where $R$ is the distance to the disk center, $R_{\rm c}$ = 78.9 AU and
$\gamma=-0.2$. Note that the dust surface density profile is assumed to have a
smooth radial distribution and ignore all the gaps and rings observed with
higher spatial resolution in ALMA Partnership et al. (2015); this is the case
for the lower resolution of ALMA polarimetric observations that we are
attempting to reproduce.

The vertical dust density is assumed to follow Gaussian distribution,

$\displaystyle\rho=\frac{\Sigma}{\sqrt{2\pi}h_{\rm
d}}\text{exp}\left(\frac{-z^{2}}{2h_{\rm d}^{2}}\right),$ (2)

where $z$ is the height from the disk midplane and $h_{\rm d}$ is the disk
dust scale height whose profile is taken as

$\displaystyle h_{\rm
d}(R)=h_{0}\left(\frac{R}{100\;\text{AU}}\right)^{\beta},$ (3)

with $h_{0}$ = 15 AU being the scale height at 100 AU from the center of the
disk, and $\beta$ = 1.15 the exponent factor characterizing the flaring effect
of the disk (see Model A in Liu et al., 2017). However, HL Tau was studied to
have a high degree of dust settling based on the appearance of rings and gaps
on millimeter observations (Pinte et al., 2016). Dust settling is expected to
have a significant effect on polarization due to self-scattering (Yang et al.,
2017). Hence, in the calculation for self-scattering only, we take into
account the effect of dust-settling by dividing the scale height $h_{\rm d}$
by a factor $f_{\rm settle}=10$, which is defined as the dust-settling factor
from the disk scale height, following Kataoka et al. (2016). In this
calculation, the radial temperature profile is taken identical to that at the
disk midplane of radiative heating calculation with the original scale height.
The disk is also assumed to be vertically isothermal in this case. We assume
the disk geometrical thickness would have a minor effect on the intrinsic dust
thermal emission at (sub)millimeter wavelength for simplicity.

Assuming the fiducial grain size of $a_{\rm max}$ = 100 ${\mu}$m, we take the
total dust mass of the system $M_{\rm d}$ as a calibrating parameter to
reproduce the three Stokes $I$ emissions at the three wavelengths mentioned
earlier. $M_{\rm d}$ = $3.8\times 10^{-3}M_{\odot}$ is found to be the best-
fitted value, which produces total emissions comparable to the observations as
shown in Figure 2. Although the residuals are still significant, we stress
that our aim is to reproduce axisymmetric continuum emissions for the whole
three wavelengths with a simple model and ignore the asymmetries of the disk
in the later modeling. Therefore, we will not further search for more optimal
parameter sets. For the gas part, for simplicity, we take the gas-to-dust
ratio equal to 100 to derive the gas density profile.

The radiation field and temperature profile resulting from radiative transfer
calculation with this disk model is shown in Appendix A.

Concerning the magnetic field, we assume that it is dominated by toroidal
component as indicated in Flock et al. (2015) with the strength’s profile as
estimated in Bai (2011) assuming the accretion process within the disk is
solely driven by magnetorotational instability in the active layer.

$\displaystyle
B=1.0\;\text{G}\times\left(\frac{\dot{M}}{10^{-8}\;M_{\odot}\;\text{yr}^{-1}}\right)^{-1/2}\left(\frac{R}{1\;\text{AU}}\right)^{-11/8}$
(4)

where $\dot{M}=1.6\times 10^{-7}M_{\odot}\;\text{yr}^{-1}$ is the accretion
rate measured by White & Hillenbrand (2004).

Figure 3: Top left image represents the polarization map of the ideal scenario
where all of the grains achieve perfect internal and external alignment (i.e.,
$f_{\rm high-J}\,=\,1$). The central and right panels include the effects of
slow INR and enhanced (super)paramagnetic relaxation. Central panels show the
case of right IA (rIA), while the right ones correspond to wrong IA (wIA).
Three different dust models are taken into consideration with PM grains (top)
and SPM grains which have $N_{\rm cl}=10^{3}$ (middle) and $10^{4}$ (bottom),
respectively. Table 2: Summary of grain models considered in Section 3.1.

|  |  Grain Alignment Model |  |
---|---|---|---|---
|  |  perfect IA | right IA |  wrong IA
Grain Magnetic Property | $\chi\,\sim\,\infty$ |  Ideal |  |
| PM grains |  | rIA, PM |  wIA, PM
| SPM grains ($N_{\rm cl}\,=\,10^{3}$) |  | rIA, SPM, $N_{\rm cl}\,=\,10^{3}$ |  wIA, SPM, $N_{\rm cl}\,=\,10^{3}$
| SPM grains ($N_{\rm cl}\,=\,10^{4}$) |  | rIA, SPM, $N_{\rm cl}\,=\,10^{4}$ |  wIA, SPM, $N_{\rm cl}\,=\,10^{4}$

Notes. $\chi$ is the grain magnetic susceptibility, which takes a very large
value for the ideal case. Grain alignment models right IA/wrong IA denotes the
models where grains are aligned with right IA/wrong IA when they are at
low-$J$ attractors and slow INR.

Table 3: Summary of the correlation between the Models and Observations.

|  |  | Grain Magnetic Property (SPM) |  |  |  |  |
---|---|---|---|---|---|---|---|---
|  |  | $N_{\rm cl}\,=\,10^{2}$ |  | $N_{\rm cl}\,=\,2\,\times\,10^{3}$ |  | $N_{\rm cl}\,=\,10^{4}$ |
Grain Size |  $\lambda\,=\,3.1\,$mm |  $a_{\rm max}\,=\,100\,\mu$m | 0.62 |  $-0.06$ | 0.89 | 0.64 | $-0.73$ | 1.12
|  |  $a_{\rm max}\,=\,130\,\mu$m | 0.42 |  $-0.09$ | 0.94 | 0.62 | $-0.68$ | 1.11
|  |  $a_{\rm max}\,=\,200\,\mu$m | 0.11 |  $-0.15$ | 0.75 | 0.17 | $-0.44$ | 0.96
|  $\lambda\,=\,1.3\,$mm |  $a_{\rm max}\,=\,50\,\mu$m | 0.90 |  0.10 | 0.79 | 0.19 | $-0.19$ | $-0.49$
|  |  $a_{\rm max}\,=\,90\,\mu$m | 0.81 |  0.27 | 0.94 | 0.63 | 0.31 | $-0.80$
|  |  $a_{\rm max}\,=\,100\,\mu$m | 0.79 |  0.25 | 0.91 | 0.57 | 0.39 | $-0.82$
|  $\lambda\,=\,0.87\,$mm |  $a_{\rm max}\,=\,50\,\mu$m | 0.97 |  0.21 | 0.98 | 0.26 | 0.70 | $-1.22$
|  |  $a_{\rm max}\,=\,60\,\mu$m | 0.97 |  0.31 | 0.97 | 0.44 | 0.82 | $-1.09$
|  |  $a_{\rm max}\,=\,100\,\mu$m | 0.96 |  0.01 | 0.90 | $-0.28$ | 0.94 | $-0.38$
|  |  | $\overline{AM}$ |  $c_{\rm P}$ | $\overline{AM}$ | $c_{\rm P}$ | $\overline{AM}$ | $c_{\rm P}$

Notes. We only consider SPM grains in the model. Each column represents the
model with a different value of the number of iron atoms per cluster ($N_{\rm
cl}$). Each row represents the model with different grain sizes at each
wavelength. Each model with a couple value of $N_{\rm cl}$ and $a_{\rm max}$
has two validating factors $\overline{AM}$ and $c_{\rm P}$. $\overline{AM}$ is
the average value of alignment measure across the image. $c_{\rm P}$ is the
correlation factor of pixel-by-pixel polarization degree between the models
and observations. The closer the values of $\overline{AM}$ and $c_{\rm P}$ to
unity, the better the model compared to the observations.

##  3 Numerical Results

###  3.1 Effect of Grain Magnetic Properties on Disk Polarization

Here, we investigate the effect of grain magnetic properties on disk
polarization due to the MRAT mechanism at $\lambda$ = 3.1 mm with fiducial
grain size $a_{\rm max}=130\;\mu$m as measured in Mori & Kataoka (2021). Note
that although the same disk and dust models as described in Section 2 are
applied, we temporarily ignore the effect of self-scattering for this study.

In Figure 3, we show the resulting polarization degree maps overlaid with
polarization vectors considering three different scenarios of IA. The
different grain alignment scenarios considered in this section are summarised
in Table 2.

Initially, we examine the ideal scenario (top left panel) in which all grains
have perfect IA (i.e., the fraction of grains aligned at high-$J$ attractors
is $f_{\rm high-J}\,=\,1$). In this case, grains are expected to have their
longer axis aligned perpendicular to the toroidal magnetic field with optimal
efficiency, which results in the radial polarization pattern with a very high
polarization degree of up to the level of $\sim$10%.

In the subsequent scenarios, we take into account the effect of slow INR and
enhanced (super)paramagnetic alignment by MRAT. We consider both cases of
right and wrong IA (middle and right panels). The effect of PM grains (top
panels) and SPM grains with $N_{\rm cl}\;=\;10^{3}$ and $10^{4}$ (middle and
bottom panels) are also included for comparison.

In the case of right IA (middle panels), the polarization vectors are the same
as in the ideal case, yet the polarization degree is significantly lower. When
grains are made of PM, the polarization degree is $\sim$0.01%. In this case,
the alignment of VLGs is extremely weak and easily dwarfed by other mechanisms
such as $k$-RAT (Tazaki et al., 2017). Hence, it would be impossible to
observe such a disk polarization in reality. For SPM grains, the alignment is
considerably more efficient (i.e., the polarization degree is two orders of
magnitude higher than in the PM grains scenario). With increasing $N_{\rm
cl}$, higher values of $f_{\rm high-J}$ (i.e., more efficient grain alignment)
are expected due to increasing magnetic relaxation (Hoang & Lazarian, 2016).
However, when grains are aligned with slow INR in the same direction, their IA
efficiencies are the same whether at high-$J$ or low-$J$ attractors. Thus, the
polarization degrees observed at $N_{\rm cl}\,=\,10^{3}\;\text{and}\;10^{4}$
are seemingly identical.

In the case of wrong IA (right panels), the polarization morphology is more
diverse. Grain alignment is weak for PM grains with polarization degrees as
low as $\sim$0.01%. The polarization vectors show an azimuthal pattern owning
to the inefficient IA of grains’ $\boldsymbol{\hat{a_{1}}}$ perpendicular to
$\boldsymbol{J}$ at low-$J$ attractors with slow INR and EA of
$\boldsymbol{J}$ along $\boldsymbol{B}$. SPM grains, conversely, induce a
higher polarization degree with increasing $N_{\text{cl}}$. However, the
polarization pattern alters when we keep increasing $N_{\text{cl}}$. Indeed,
the increased level of iron inclusion also causes the increased abundance of
grains at high-$J$ attractors (i.e., higher $f_{\text{high-J}}$), which would
drive them to have more efficient IA with the grains’ longer axis being
aligned perpendicular to the magnetic field. As a result, the polarization
vectors would eventually have the radial direction when $N_{\rm
cl}\;=\;10^{4}$.

Figure 4: The top panels show the total dust polarization induced by both
self-scattering and aligned grains via the MRAT mechanism. These results are
obtained for grains aligned with wrong IA, $N_{\rm cl}\sim 2\times 10^{3}$;
and $a_{\rm max}$ = 130, 90, and 60 $\mu$m at $\lambda$ = 3.1, 1.3, and 0.87
mm, respectively. Second-row panels show the ALMA polarization observations at
the corresponding wavelengths. The third-row panels show the maps of alignment
measure ($AM$), which is an indicator of the similarities between the
polarization pattern from the model and that from the observations. The bottom
panels are the pixel-to-pixel correlations between the polarization degree
from the model and that from the observations. The blue dots correspond to the
polarization degree data points of each pixel. The red lines are linear fits
to the correlations. The values of ($\overline{AM},\,c_{\rm P}$) of the model
are $\sim$ (0.94, 0.62), (0.94, 0.63), and (0.97, 0.44) at $\lambda$ = 3.1,
1.3, and 0.87 mm, respectively.

Note that the features of disk polarization only induced by grain alignment
observed above also hold for observations at both wavelengths $\lambda=1.3$
and 0.87 mm.

###  3.2 Mixture of MRAT Grain Alignment and Self-scattering: Final Results

Mori & Kataoka (2021) concluded that a superposition of azimuthally aligned
grains and self-scattering can produce the HL Tau disk polarization. In the
previous Section 3.1, we demonstrate that only the grain alignment model of
wrong IA with SPM grains having $N_{\rm cl}$ at the level of $\sim$$10^{3}$
can produce the azimuthal polarization morphology and the polarization degree
at $\sim$1%, which are ones of the notable features of the observations and in
line with the description from the literature. Hence, in this section, we
limit our grain alignment model to this scenario and add the effect of self-
scattering to reproduce the observational data.

To validate the agreement between models and observations, we compare
separately their polarization degree and polarization patterns. We check the
pixel-by-pixel correlation for the polarization degree ($c_{\rm P}$) while
computing the Alignment Measure ($AM$) as introduced in González-Casanova &
Lazarian (2017) to compare the polarization orientation.

$\displaystyle AM=2\left(\text{cos}^{2}\theta_{\text{r}}-\frac{1}{2}\right),$
(5)

where $\theta_{\text{r}}$ = $|\theta_{\text{mod}}\;-\;\theta_{\text{obs}}|$ is
the angle difference between the model and observation. The scale of $AM$ is
from $-1$ to 1, where $AM=1$ corresponds to a perfect match of polarization
vectors and $AM=-1$ indicates that the two vectors are perpendicular to each
other. We then take its average value across the image $\overline{AM}$ to
represent the overall correlation between polarization vectors of the model
and observations. Finally, the values of $c_{\rm P}$ and $\overline{AM}$ are
compared between each model to find the best model for the observations. Note
that the closer the values of $c_{\rm P}$ and $\overline{AM}$ to unity, the
better the agreement between the model and the observations. Table 3
summarizes the values of $c_{\rm P}$ and $\overline{AM}$ of the typical models
that we considered.

As shown in Table 3, the model that best reproduces the disk polarization
observed by ALMA is found at $N_{\rm cl}=2\times 10^{3}$; $a_{\rm max}=$ 130,
90, and 60 $\mu$m for $\lambda$ = 3.1, 1.3, and 0.87 mm, respectively. The
disk polarization corresponding to this model is displayed in Figure 4. As
expected, VLGs with wrong IA produce azimuthal polarization vectors with
strong polarization degrees along the major axis while self-scattering
produces both polarization vectors and strong polarization degrees along the
disk minor axis. Furthermore, as the probing wavelength decreases, the optical
depth of the disk increases, which induces a weaker polarization signature
from grain alignment due to dichroic extinction (Lin et al., 2022). As a
result of their superposition at different wavelengths, the polarization shows
a transition from azimuthal to parallel to the disk minor axis patterns as in
the observations.

##  4 Discussions

###  4.1 Comparison to Previous Studies and Implications

Since ALMA’s first polarimetric detection of HL Tau (Stephens et al., 2014),
there have been many observational and theoretical studies aiming to
understand the diversity of dust polarization observed toward the disk. Using
the RAT alignment theory, Tazaki et al. (2017) first performed a theoretical
study of grain alignment in protoplanetary disks and concluded that magnetic
alignment is impossible for mm-sized grains in this environment, even with a
high degree of iron inclusion, consistent with predictions in (Hoang &
Lazarian, 2016). The authors suggested that such large grains can be aligned
via radiation direction (i.e., $k$-RAT, Lazarian & Hoang 2007a), which can
produce the observed azimuthal polarization pattern. However, Tazaki et al.
(2017) assumed the right IA for VLGs with slow INR, but their wrong IA cannot
be ruled out (Hoang et al., 2022). In this paper, our detailed modeling for HL
Tau using the updated POLARIS shows that the disk polarization can be
reproduced with magnetic alignment via the MRAT mechanism, but VLGs with slow
INR must have wrong IA. Therefore, dust polarization can still trace the
magnetic field within protoplanetary disks. Moreover, due to the wrong IA, the
polarization vectors are along the magnetic field direction, so that one does
not need to rotate the polarization vectors by $90^{\circ}$ to infer the
magnetic field, as in the case of right IA.

However, Yang et al. (2019) and Mori & Kataoka (2021) modeled the HL Tau
polarization observed at ALMA Band 3 ($\lambda$ = 3.1 mm) and concluded that
$k$-RAT could not explain its polarization pattern, which turns out to be more
delicate being elliptical, namely the polarization radiation is required to be
emitted only in the disk plane. Therefore, the authors suggest that the HL Tau
polarization could be reproduced by effective prolate grains that are aligned
with their major axis along the azimuthal direction. However, the reason why
the grains could be aligned in such a way is not discussed. In our detailed
modeling, such a scenario is the outcome of grain alignment physics in which
the grains’ longer axes are aligned along their angular momenta due to wrong
IA for grains with slow INR, and the grains’ angular momenta are aligned along
the direction of the toroidal magnetic field as a result of enhanced
superparamagnetic relaxation and Larmor precession for grains with embedded
iron inclusions.

In addition to the polarization properties, grain growth in HL Tau has also
been in the spotlight due to the conflict of its constraint from polarization
and from spectral index of thermal emission (Mori & Kataoka, 2021; Ueda et
al., 2021; Lin et al., 2023b; Zhang et al., 2023). Earlier in the current
study, we found that $a_{\rm max}\,=\,130\,\mu$m at $\lambda\,=\,3.1\,$mm.
This result is consistent with the value measured by using polarization with
azimuthally-aligned prolate grain in Mori & Kataoka (2021). Additionally,
$a_{\rm max}$ is constrained to differ from this value when probing at the
other two wavelengths (90$\,\mu$m and 60$\,\mu$m at
$\lambda\,=\,\text{1.3}\,$mm and 0.87 mm, respectively). The variation of the
constrained grain size proves to be a potentially convenient tool to further
constrain the true scale of grain size at the disk midplane using the
differential dust settling effect assumption. Further details are presented in
Section 4.5.

###  4.2 Constraints on Iron Abundance within Grains in the form of Clusters

Our results in Section 3 show that PM grains produce very low polarization
degree, $\sim$0.01%, which is much lower than the observed polarization level.
On the contrary, grains with embedded iron clusters can produce the observed
polarization level. Therefore, iron atoms are unlikely distributed diffusely
inside dust grains, but must be present in the form of iron clusters.

In Section 3, for convenience, we fixed the value of $\phi_{\rm sp}$ = 0.1 and
varied $N_{\text{cl}}$ to obtain the different magnetic susceptibility of
grains (see e.g. Equation 3 in Giang et al. 2023). In this section, we treat
$\phi_{\rm sp}$ as a free parameter to investigate its effect on disk
polarization. For simplicity, we test the effect using $q^{\prime}_{100}$,
which is the polarization degree of Stokes $Q$ in the principal frame at the
radius of 100 AU along the disk minor axis. A principal frame is a frame whose
ordinate coincides with the disk minor axis. In this way, $q^{\prime}_{100}$
would have a positive value when the polarization vector is radial and
negative when it is azimuthal. We demonstrate the effect at $\lambda$ = 3.1
mm, where alignment is expected to be the dominant component over self-
scattering. The value of $-q^{\prime}_{100}=1.95$ due to the alignment
component as decomposed in Lin et al. (2022) will be chosen as the reference
value. An error of 0.5$\%$ is assigned to represent the uncertainty of the
method.

The demonstration is shown in the top panel of Figure 5. Generally, at each
value of $\phi_{\text{sp}}$, the curve describing the dependence of
$-q^{\prime}_{100}$ on $N_{\text{cl}}$ follows a parabolic shape. The reason
behind this is as described in Section 3.1 in which increased $N_{\text{cl}}$
not only increases polarization degree through enhancing EA but also increases
$f_{\text{high-J}}$ which drives polarization vector from azimuthal to radial
direction. The value of $-q^{\prime}_{100}$ is at its peak when the effect of
superparamagnetic alignment from $N_{\text{cl}}$ is not negligible so that
grains can efficiently align with the magnetic field, yet not too significant
so that $f_{\text{high-J}}$ is small enough to keep the polarization vectors
at the azimuthal direction. With the higher value of $\phi_{\text{sp}}$, this
peak shifts to the position corresponding to the lower value of
$N_{\text{cl}}$.

Figure 5: Dependence of polarization degree $-q^{\prime}$ of Stokes $Q$ due to
grain-alignment-only in the principal frame at 100 AU along the disk’s minor
axis on $N_{\rm cl}$. The dashed line corresponds to the value measured in Lin
et al. (2022) using a plane-slab model with an image of HL Tau at Band 3. The
filled area corresponds to an uncertainty of 0.5%. The top panel takes into
consideration different values of $\phi_{\rm sp}$ with the original dust model
set-up. The middle panel considers the varying values of $Q_{\rm X,high-J}$
with the extreme scenario when $\phi_{\rm sp}=0.3$ and $Q_{\rm X,low-J}=-0.5$.
The bottom panel is plotted with different values of $Q_{\rm X,low-J}$ when
$\phi_{\rm sp}$ = 0.3 and $Q_{\rm X,high-J}$ = 0.01. Only the value of
$\phi_{\rm sp}\gtrsim$ 0.05; $|Q_{\rm X,high-J}|\lesssim$ 0.4, and $|Q_{\rm
X,low-J}|\gtrsim$ 0.1 can reach the threshold.

This test puts a constraint on $\phi_{\rm sp}$ with its minimum value at
$\sim$0.05, which corresponds to $\sim$16% of iron locked inside grains in the
form of cluster. Finally, it is crucial to point out that not only our
constrained value of $N_{\rm cl}\sim 2\times 10^{3}$ is consistent with the
$-q^{\prime}_{100}$ reference when considering $\phi_{\rm sp}$ = 0.1 as in our
dust model, but the value range of $N_{\rm cl}$ is also found to be
$\sim$$2(\pm 1)\times 10^{3}$ for all of the valid values of $\phi_{\rm sp}$.

###  4.3 Constraints on Internal Alignment Efficiency for Slow Internal
Relaxation

Slow INR is generally responsible for IA of large dust grains within
protoplanetary disks since large grains of more than $\sim$10 $\mu$m and high
density are ubiquitous in these environments (Hoang, 2022). However, the
efficiency of IA for the case of slow INR is not yet available from theory
(Hoang & Lazarian, 2009; Hoang et al., 2022). Usually, the efficiency of IA by
slow INR is represented by two quantities $Q_{\rm X,low-J}$ and $Q_{\rm
X,high-J}$. In Section 3, our modeling results showed that $Q_{\rm X,low-J}$
must have a negative value (i.e., wrong IA) to reproduce the azimuthal
polarization pattern observed toward the HL Tau disk. On the contrary, grains
at high-$J$ attractors always have right IA or $Q_{\rm X,high-J}$ is always
positive (Hoang & Lazarian, 2009). It indicates that the directions of IA of
grains aligned at low-$J$ and high-$J$ attractors are opposite.

Here, we will attempt to empirically constrain the magnitude of $Q_{\rm
X,low-J}$ and $Q_{\rm X,high-J}$ for slow INR using the observed polarization
data. To do so, we consider the extreme scenarios where parameters involved in
the IA process are set to their respective extremes. We initially fix
$\phi_{\rm sp}=\phi_{\text{sp,max}}=0.3$. To obtain the maximum value of
$|Q_{\rm X,high-J}|$, we set $|Q_{\rm X,low-J}|=|Q_{\rm X,low-J,max}|=0.5$. On
the other hand, $|Q_{\rm X,high-J}|=|Q_{\rm X,high-J,min}|\sim$ 0.01 is set to
constrain the minimum value of $|Q_{\rm X,low-J}|$. Following Section 4.2, we
use $-q^{\prime}_{100}$ at $\lambda=3.1$ mm as the referenced polarization
degree for the disk polarization. Its value is calculated for different values
of $Q_{\rm X,high-J}$ or $Q_{\rm X,low-J}$ and over the range $N_{\rm cl}$ =
[$10^{2}$, $10^{4}$]. The plots are shown in the middle and bottom panels of
Figure 5. By considering the values satisfying the reference value, we can
deduce the maximum possible value of $|Q_{\rm X,high-J}|$ is $\sim$0.4, while
the minimum value of $|Q_{\rm X,low-J}|$ is at $\sim$0.1.

###  4.4 Other Grain Alignment Mechanisms

Besides MRAT, grain alignment due to interaction of the gas flow and grains,
including Gold alignment (Gold, 1952) and mechanical torque (MET) alignment
with the grain’s long axis along the flow velocity (so-called $v$-MET; Hoang
et al. 2022) can also produce an azimuthal polarization pattern.

Gold alignment was thought to be the most promising mechanism to explain the
polarization morphology caused by grain alignment in HL Tau (Yang et al.,
2019). However, the possibility of the Gold mechanism causing grain alignment
in disk environments is highly questionable (Mori & Kataoka, 2021). Since the
mechanism requires supersonic gas speed to have effectively aligned grains,
whereas it is generally subsonic in protoplanetary disks (Cho & Lazarian,
2007), it is unlikely to induce grain alignment within this environment.

$v$-MET, on the other hand, offers a rather promising alternative. If grains
with Keplerian motions align with the drift direction with wrong IA, the
polarization pattern is also expected to be azimuthal (Hoang et al., 2022).
However, it may not be the case if we take a more thorough inspection. The
polarization direction depends on the direction of gas flow in the frame of
dust grains, which is a strong function of Stokes number denoting how well
dust grains couple to ambient gas. With $a_{\rm max}\sim 100$ $\mu$m, and
assuming a gas surface density of $\sim$30 gcm-2, the Stokes number is
estimated to be $St\sim 8.8\times 10^{-4}$. This value is much smaller than
unity which indicates that the direction of gas velocity with respect to the
grains will be radial (see Figure 2 in Kataoka et al., 2019). In the case of
right IA, a circular polarization pattern is produced. On the other hand, if
grains are aligned with wrong IA, the pattern will be radial. Both cases are
not compatible with the elliptical patterns from observations (Yang et al.,
2019). Note that the argument still holds even when considering grain size up
to 1000 $\mu$m which corresponds to $St\sim 8.8\times 10^{-2}\ll 1$.

Nevertheless, it is crucial to stress that we cannot rule out the possibility
of enhanced alignment with grain drift (Hoang et al., 2022), which is
neglected in our model. Since the direction of gas flow with respect to the
grains is radial, which coincides with the direction of the radiation beam,
the combined effect of grain drift and radiation can significantly improve the
efficiency of grain alignment within the disk.

In conclusion, it has been demonstrated that MRAT is the only grain alignment
mechanism in the constraint of the current grain alignment theory that can
explain the polarization signatures observed in the HL Tau disk by ALMA.

###  4.5 Grain Size Population

We could reproduce the dust polarization in accordance with the observations
with the initial assumption that the three wavelengths trace different layers
of the disk, which contain grains of discrete sizes. In the optically thick
region of HL Tau, this is the case because the differential dust-settling
effect can cause light at different wavelengths to trace different layers of
the disk with longer wavelengths leaning towards the disk midplane where
larger grains are expected to be present (Brunngräber & Wolf, 2020; Ueda et
al., 2021). Indeed, the grain size in our results is constrained at larger
values when tracing the dust by longer wavelengths, namely $a_{\rm max}$ = 70,
90, and 130$\,\mu$m at $\lambda$ = 0.87, 1.3, and 3.1 mm, respectively.
Motivated by this idea, it is possible to derive the value of maximum grain
size at the disk midplane by extrapolating from the values of $a_{\rm max}$.

To do so, it is essential to derive approximately the relationship between the
maximum grain size detected $a_{\rm max}$ and its tracing wavelength
$\lambda$. First of all, we assume a simple correlation between $\lambda$ and
penetration depth as such,

$\displaystyle\lambda=k[h_{\rm d}(a_{\rm min})-h_{\rm d}]$ (6)

where $k$ is a correlation factor ($k>0$), $h_{\rm d}$ is the dust scale
height as a function of dust size, $h_{\rm d}(a_{\rm min})$ is the scale
height of grains with the minimum grain size. Next, we assume a power-law
dependence of the dust scale height $h_{\rm d}$ on grain size a :

$\displaystyle h_{\rm d}(a)=h_{\rm d}(a_{\rm min})\left(\frac{a}{a_{\rm
min}}\right)^{-\alpha}$ (7)

where $\alpha$ is the exponent component. Combining the two equations above,
we can derive the effective grain size probed by a certain wavelength as
follows.

$\displaystyle a(\lambda)=a_{\rm min}(1-\beta\lambda)^{-1/\alpha}$ (8)

where $\beta=1/[kh_{\rm d}(a_{\rm min})]$ ($\beta>0$), $a_{\rm min}$ is the
minimum grain size. Note that in this approximation, we do not consider the
flaring shape of the disk, which means that the dust scale height is assumed
to be independent of the radius from the disk center. Hence, some familiar
parameters such as $a_{\rm min}$ may not be the same as specified earlier in
our model. With that in mind, fitting the above function using the three
values of $a_{\rm max}$ constrained at the three wavelengths in Section 3.2
gives us the parameter set ($a_{\rm min},\beta,\alpha$) = ($50\,\mu$m,
$0.01\,\text{mm}^{-1}$, 0.04). We choose $\lambda$ = 9 mm and $\lambda$ = 13
mm as lower and upper bounds for our constraint of maximum grain size as we
estimate the opacity index of the disk inner region to be at $\sim$ unity
(i.e., the disk becomes completely optically thin in the whole radius range)
following the model in Carrasco-González et al. (2019). This estimation is
based, in particular, on bands Ka and Q data of VLA (top panel of Figure 6 of
Carrasco-González et al. 2019). Assuming 9 mm $\lesssim\lambda_{a_{\rm
max}}\lesssim$ 13 mm, we can derive the value of maximum grain size to be 640
$\mu$m $\lesssim a_{\rm max}\lesssim 2000\;\mu$m. An illustration of the
extrapolation can be found in Figure 6. Interestingly, it indicates that dust
grains within the disk can be as large as millimeter-sized as measured by
Carrasco-González et al. (2019) using spectral index SED fitting.

Figure 6: The orange data points represent the maximum grain sizes constrained
at the three wavelengths (see Section 3.2). The blue line is the result of the
fit of the dependence of grain size on the tracing wavelength using the orange
data points. The two vertical dashed lines correspond to the lower and upper
bounds of the pink area where the disk is expected to go from optically thick
to completely optically thin (i.e., $\tau\sim 1$). This criterion is expected
to trace the grains with maximum size at the disk mid-plane.

However, the above argument is only correct at the optically thick inner
region of the disk. The optical depth would decrease with the radial distance
from the disk center. Thus, the dust differential settling effect would have a
weaker impact on the outer region, which requires an alternative explanation
for the grain size constraint inconsistency between polarization and spectral
index SED fitting. Irregular grain structures have also been studied to
alleviate this grain size conflict (Lin et al., 2023b; Zhang et al., 2023),
where grain size can be larger than 1 mm depending on the grain porosity.
Therefore, it is crucial that future studies implement the effect of grain
porosity into detailed modeling to have a more complete picture of HL Tau
polarization.

###  4.6 Disk Polarization with Substructures

Previously in our original disk model in Section 2.2, the disk is assumed to
have a smooth radial density distribution. However, ALMA Partnership et al.
(2015) has revealed that the disk contains an alternate distribution of rings
and gaps, which is a crucial sign for planet formation. Stephens et al. (2023)
are recently able to detect the dust polarization within the substructures at
$\lambda\,=\,0.87\,$mm and find the distinguishable features observed in gaps
compared to the rings. In the gaps, the authors find the polarization vectors
to be azimuthal with a higher polarization degree (up to $\sim$4%) than within
the rings. This is understandable considering within the gaps, the optical
depth is at unity, which induces a higher polarization degree due to self-
scattering (up to $\sim$3%, Kataoka et al., 2016). The lower dust opacity also
induces lower dichroic extinction, which means a higher polarization degree of
the intrinsic emission from aligned grains. The lower gas density in the gaps
is also expected to make the grain alignment process more efficient with less
effective gas damping (Hoang et al., 2022). However, too low gas density may
also induce more grains aligning with fast INR and lower the degree of
azimuthal polarization. This essentially signifies the importance of having
accurate gas density for polarization calculation. Furthermore, with lower gas
density, dust and gas are expected to couple more efficiently and the effect
of grain alignment due to grain drift would be non-negligible. We attempt to
model the disk with substructures in Appendix B. However, it is only a rough
estimation and cannot reproduce the entire detailed polarimetric observation
as in Stephens et al. (2023). This is because the dust vertical settling
effect has not been taken into account, which ignores the larger grain size
being probed in the gaps due to less opacity. The gas-to-dust ratio contrast
in these gaps and rings is also still uncertain (Pinte et al., 2016; Yen et
al., 2019) and the simple assumption of gas-to-dust ratio equals to 100 is no
longer acceptable. Moreover, the effect of mechanical torque alignment (Hoang
et al., 2022) has not been taken into account. Its incorporation into modeling
is still difficult due to its sensitivity to grain shape. Hence, additional
studies are needed to fully reproduce the high-resolution HL Tau dust
polarization observed by ALMA.

###  4.7 Caveats

Our model is subject to the unresolved grain shape inconsistency. While we
calculate the direct polarized thermal emission from aligned elongated dust
grains, spherical grains are used in the MCRT calculations. This is due to the
difficulties in computing the optical properties of non-spherical dust grains.

The assumed geometrically thick disk in intrinsic dust thermal emission
calculation may cause unwanted near-far side asymmetry in the disk emission at
optically thick region due to more dichroic extinction in the near side of the
disk. However, we ignore this effect due to the nature of our simple disk
model setup and the preferred need to use a precise approximation of the
radiation field and gas density for grain alignment calculation.

##  5 Summary

The origin of the polarized emission of HL Tau at different wavelengths has
been debated in previous studies. In this work, through our detailed multi-
wavelength modeling, the following results are found.

  1. 1.

The ideal scenario of magnetic alignment significantly overestimates the
polarization degree compared to the realistic cases. In the realistic
scenario, especially when grains are aligned at low-$J$ attractors with slow
INR, grains can either align with wrong or right IA, which manifests in
azimuthal or radial polarization patterns, respectively. Additionally, PM
grains cannot reproduce the observed polarization data and thus cannot be the
primary dust material in HL Tau.

  2. 2.

A mixture of self-scattering and MRAT grain alignment mechanism where grains
are aligned with wrong IA at low-$J$ attractors and slow INR can reproduce the
dust polarization observed in HL Tau by ALMA at (sub)millimeter wavelengths.

  3. 3.

Grains are shown to be made of superparamagnetic materials with the number of
iron atoms inside each cluster is $N_{\rm cl}\sim 2\times 10^{3}$. The
abundance of embedded iron inside grains identified from our model is at
$\sim$16% minimum.

  4. 4.

By assuming grain sizes probed to be $\lambda$-dependent, maximum grain sizes
are constrained at $a_{\rm max}\sim$ 60, 90, and 130$\,\mu$m when probing at
$\lambda$ = 0.87, 1.3, and 3.1 mm, respectively. Further extrapolation with
dust differential settling effect gives the value of grain size at the disk
midplane to be millimeter-scaled consistent with estimation from SED fitting.

  5. 5.

MRAT is found to be the only mechanism in the contemporary grain alignment
theory that can explain the polarization morphology caused by the grain
alignment component within the HL Tau disk.

## Acknowledgements

We thank Ph.D. candidate Nguyen Chau Giang (KASI) for providing the updated
POLARIS and significant discussion about the code along with its embedded
grain alignment physics. N.T.T. is supported by a grant from the Simons
Foundation to IFIRSE, ICISE (916424, N.H.). N.B.N., N.T.P., and P.N.D. are
funded by Vingroup Innovation Foundation (VINIF) under project code
VINIF.2023.DA.057. ALMA is a partnership of ESO (representing its member
states), NSF (USA), and NINS (Japan), together with NRC (Canada), MOST and
ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic
of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO, and NAOJ.

##  Appendix A Radiation Field

With the radiation source at the disk center, the radiation flux $J_{\rm
rad}$, spectrum weighted wavelength $\overline{\lambda}$, and anisotropy
parameter $\overline{\gamma}$ (defined in Bethell et al., 2007) are calculated
and illustrated in Figure 7.

While stellar radiation dominates the upper layers of the disk, the inner disk
region is dominated by cold dust emission. The upper layers of the disk
correspond to the disk photosphere, where it can be observed with IR/NIR
observations. The inner disk region, on the other hand, can be traced by radio
observations, which is the case in our modeling. Calculations show that the
radiation at this region is dominated at $\overline{\lambda}\sim$ 140$\;\mu$m,
with anisotropy parameter $\overline{\gamma}\sim 0.1$. Therefore, only about
10$\%$ of the radiation flux primarily at $\overline{\lambda}$ $\sim$ 140
$\mu$m can contribute to RATs towards the grains within the disk because only
the anisotropic radiation is important for RAT alignment.

In our model, we also use radiative heating to calculate the temperature
profile of the disk. Our radiative transfer calculation yields a temperature
profile at the disk mid-plane of

$\displaystyle T(R)=311\left(\frac{R}{1\,\text{AU}}\right)^{-0.55}$ (A1)

Similar mid-plane temperature profiles were also derived from previous
radiative transfer models for HL Tau (Men’shchikov et al., 1999; Kwon et al.,
2015; Okuzumi & Tazaki, 2019).

Figure 7: From left to right: Radiation flux $J_{\text{rad}}$, spectrum
weighted wavelength $\overline{\lambda}$, and spectrum weighted anisotropy
parameter $\overline{\gamma}$. The calculation is made assuming a single
radiation source to be the central protostar with luminosity $L_{\star}\sim
11L_{\odot}$.

##  Appendix B Disk model with ring/gap structures

Throughout this work, we adopted the smoothed structure of the HL Tau disk.
However, Stephens et al. (2023) reported ring and gap structures with the new
ALMA observations. They found that the polarization degree is much higher in
the gaps than in the rings, and the polarization vectors are azimuthal.

This section aims to understand this newly observed phenomenon based on the
same physical model. We first mimic the disk structure at
$\lambda\,=\,0.87\,$mm as in Stephens et al. (2023), which yields the surface
density profile

$\displaystyle\Sigma\propto\Sigma_{a}\left(\frac{R}{R_{0}}\right)^{-p}\text{exp}\left[-\left(\frac{R}{R_{\text{c}}}\right)^{1.5}\right]+\sum_{i=1}^{N}\Sigma_{i}\,\text{exp}\left[-\frac{1}{2}\left(\frac{R-C_{i}}{W_{i}}\right)^{2}\right]\text{,}$
(B1)

where $\Sigma_{a}\,=\,2.3$, $R_{0}\,=\,10\,$AU, $p\,=\,0.5$, and the
parameters for the rings are given in Table 4. However, we still use the
temperature profile taken from radiative heating calculation rather than a
fixed profile as used in the literature. We then normalize the dust density by
the total dust mass of $M_{\rm d}\sim 5\times 10^{-3}M_{\odot}$ to produce a
comparable Stokes $I$ emission at $\lambda=0.87$ mm as in Stephens et al.
(2023).

Table 4: Parameters for the identified rings in HL Tau. Ring number | $\Sigma_{i}$ |  $C_{i}$ [AU] |  $W_{i}$ [AU]
---|---|---|---
1 | 8 | 24 | 4
2 | 5 | 39 | 3
3 | 5 | 49 | 3
4 | 8 | 59 | 2
5 | 3 | 73 | 4
6 | 8 | 88 | 3
7 | 8 | 102 | 2
8 | 8 | 116 | 2

For the polarization calculation, we adopt the same dust models as constrained
in Section 3.2 with $N_{\rm cl}=2\times 10^{3}$, and $a_{\rm max}=60\,\mu$m.
However, we set $Q_{\rm X,low-J}\,=\,-0.5$ to model the higher IA efficiency
in the gaps than in the rings due to less effective gas damping. The predicted
polarization degree and polarization angles with the ring/gap disk are shown
in Figure 8. Even with the simple dust model, we can spot differences between
the polarization observed in the rings and gaps. In the gaps, we found that
the polarization degree is higher, and the polarization caused by the self-
scattering component becomes less dominant, resulting in a more azimuthal
polarization pattern. These features are consistent with the new observations
in HL Tau.

Figure 8: Prediction of the toy model of the HL Tau disk with the ring-gap
structures at $\lambda\,=\,0.87\,$mm. The polarization degree is higher in the
gap, and the polarization orientation becomes azimuthal.

## References

  * ALMA Partnership et al. (2015) ALMA Partnership, Brogan, C. L., Pérez, L. M., et al. 2015, ApJ, 808, L3, doi: 10.1088/2041-8205/808/1/L3
  * Andersson et al. (2015) Andersson, B. G., Lazarian, A., & Vaillancourt, J. E. 2015, ARA&A, 53, 501, doi: 10.1146/annurev-astro-082214-122414
  * Bai (2011) Bai, X.-N. 2011, ApJ, 739, 50, doi: 10.1088/0004-637X/739/1/50
  * Balbus & Hawley (1991) Balbus, S. A., & Hawley, J. F. 1991, ApJ, 376, 214, doi: 10.1086/170270
  * Balbus & Hawley (1998) —. 1998, Reviews of Modern Physics, 70, 1, doi: 10.1103/RevModPhys.70.1
  * Bethell et al. (2007) Bethell, T. J., Chepurnov, A., Lazarian, A., & Kim, J. 2007, ApJ, 663, 1055, doi: 10.1086/516622
  * Birnstiel et al. (2018) Birnstiel, T., Dullemond, C. P., Zhu, Z., et al. 2018, ApJ, 869, L45, doi: 10.3847/2041-8213/aaf743
  * Blandford & Payne (1982) Blandford, R. D., & Payne, D. G. 1982, MNRAS, 199, 883, doi: 10.1093/mnras/199.4.883
  * Bohren & Huffman (1998) Bohren, C. F., & Huffman, D. R. 1998, Absorption and Scattering of Light by Small Particles (Wiley-VCH)
  * Brunngräber & Wolf (2020) Brunngräber, R., & Wolf, S. 2020, A&A, 640, A122, doi: 10.1051/0004-6361/202037981
  * Carrasco-González et al. (2019) Carrasco-González, C., Sierra, A., Flock, M., et al. 2019, ApJ, 883, 71, doi: 10.3847/1538-4357/ab3d33
  * Cho & Lazarian (2007) Cho, J., & Lazarian, A. 2007, ApJ, 669, 1085, doi: 10.1086/521805
  * Cox et al. (2018) Cox, E. G., Harris, R. J., Looney, L. W., et al. 2018, ApJ, 855, 92, doi: 10.3847/1538-4357/aaacd2
  * Davis & Greenstein (1951) Davis, Leverett, J., & Greenstein, J. L. 1951, ApJ, 114, 206, doi: 10.1086/145464
  * Dolginov & Mitrofanov (1976) Dolginov, A. Z., & Mitrofanov, I. G. 1976, Ap&SS, 43, 291, doi: 10.1007/BF00640010
  * Draine (2003) Draine, B. T. 2003, ARA&A, 41, 241, doi: 10.1146/annurev.astro.41.011802.094840
  * Draine & Flatau (1994) Draine, B. T., & Flatau, P. J. 1994, Journal of the Optical Society of America A, 11, 1491, doi: 10.1364/JOSAA.11.001491
  * Draine & Weingartner (1997) Draine, B. T., & Weingartner, J. C. 1997, ApJ, 480, 633, doi: 10.1086/304008
  * Fendt (2003) Fendt, C. 2003, A&A, 411, 623, doi: 10.1051/0004-6361:20034154
  * Flock et al. (2015) Flock, M., Ruge, J. P., Dzyurkevich, N., et al. 2015, A&A, 574, A68, doi: 10.1051/0004-6361/201424693
  * Giang et al. (2023) Giang, N. C., Hoang, T., Kim, J.-G., & Tram, L. N. 2023, MNRAS, 520, 3788, doi: 10.1093/mnras/stad020
  * Gold (1952) Gold, T. 1952, MNRAS, 112, 215, doi: 10.1093/mnras/112.2.215
  * González-Casanova & Lazarian (2017) González-Casanova, D. F., & Lazarian, A. 2017, ApJ, 835, 41, doi: 10.3847/1538-4357/835/1/41
  * Harrison et al. (2019) Harrison, R. E., Looney, L. W., Stephens, I. W., et al. 2019, ApJ, 877, L2, doi: 10.3847/2041-8213/ab1e46
  * Hawley et al. (1994) Hawley, J. F., Gammie, C. F., & Balbus, S. A. 1994, in Astronomical Society of the Pacific Conference Series, Vol. 54, The Physics of Active Galaxies, ed. G. V. Bicknell, M. A. Dopita, & P. J. Quinn, 73
  * Henning & Stognienko (1996) Henning, T., & Stognienko, R. 1996, A&A, 311, 291
  * Hoang & Lazarian (2016) Hoang , T., & Lazarian, A. 2016, ApJ, 831, 159, doi: 10.3847/0004-637X/831/2/159
  * Hoang (2022) Hoang, T. 2022, ApJ, 928, 102, doi: 10.3847/1538-4357/ac5408
  * Hoang & Lazarian (2008) Hoang, T., & Lazarian, A. 2008, MNRAS, 388, 117, doi: 10.1111/j.1365-2966.2008.13249.x
  * Hoang & Lazarian (2009) —. 2009, ApJ, 697, 1316, doi: 10.1088/0004-637X/697/2/1316
  * Hoang et al. (2022) Hoang, T., Tram, L. N., Minh Phan, V. H., et al. 2022, AJ, 164, 248, doi: 10.3847/1538-3881/ac9af5
  * Johansen (2009) Johansen, A. 2009, in Cosmic Magnetic Fields: From Planets, to Stars and Galaxies, ed. K. G. Strassmeier, A. G. Kosovichev, & J. E. Beckman, Vol. 259, 249–258, doi: 10.1017/S1743921309030592
  * Jones & Spitzer (1967) Jones, R. V., & Spitzer, Lyman, J. 1967, ApJ, 147, 943, doi: 10.1086/149086
  * Kataoka et al. (2016) Kataoka, A., Muto, T., Momose, M., Tsukagoshi, T., & Dullemond, C. P. 2016, ApJ, 820, 54, doi: 10.3847/0004-637X/820/1/54
  * Kataoka et al. (2019) Kataoka, A., Okuzumi, S., & Tazaki, R. 2019, ApJ, 874, L6, doi: 10.3847/2041-8213/ab0c9a
  * Kataoka et al. (2017) Kataoka, A., Tsukagoshi, T., Pohl, A., et al. 2017, ApJ, 844, L5, doi: 10.3847/2041-8213/aa7e33
  * Kataoka et al. (2015) Kataoka, A., Muto, T., Momose, M., et al. 2015, ApJ, 809, 78, doi: 10.1088/0004-637X/809/1/78
  * Kirchschlager & Bertrang (2020) Kirchschlager, F., & Bertrang, G. H. M. 2020, A&A, 638, A116, doi: 10.1051/0004-6361/202037943
  * Kwon et al. (2011) Kwon, W., Looney, L. W., & Mundy, L. G. 2011, ApJ, 741, 3, doi: 10.1088/0004-637X/741/1/3
  * Kwon et al. (2015) Kwon, W., Looney, L. W., Mundy, L. G., & Welch, W. J. 2015, ApJ, 808, 102, doi: 10.1088/0004-637X/808/1/102
  * Lazarian & Hoang (2008) Lazarian , A., & Hoang, T. 2008, ApJ, 676, L25, doi: 10.1086/586706
  * Lazarian & Hoang (2007a) Lazarian, A., & Hoang, T. 2007a, MNRAS, 378, 910, doi: 10.1111/j.1365-2966.2007.11817.x
  * Lazarian & Hoang (2007b) —. 2007b, ApJ, 669, L77, doi: 10.1086/523849
  * Lazarian & Hoang (2021) —. 2021, ApJ, 908, 12, doi: 10.3847/1538-4357/abd02c
  * Lee & Draine (1985) Lee, H. M., & Draine, B. T. 1985, ApJ, 290, 211, doi: 10.1086/162974
  * Lin et al. (2022) Lin, Z.-Y. D., Li, Z.-Y., Yang, H., et al. 2022, MNRAS, 512, 3922, doi: 10.1093/mnras/stac753
  * Lin et al. (2023a) Lin, Z.-Y. D., Li, Z.-Y., Stephens, I. W., et al. 2023a, arXiv e-prints, arXiv:2309.10055, doi: 10.48550/arXiv.2309.10055
  * Lin et al. (2023b) Lin, Z.-Y. D., Li, Z.-Y., Yang, H., et al. 2023b, MNRAS, 520, 1210, doi: 10.1093/mnras/stad173
  * Liu et al. (2017) Liu, Y., Henning, T., Carrasco-González, C., et al. 2017, A&A, 607, A74, doi: 10.1051/0004-6361/201629786
  * Mathis et al. (1977) Mathis, J. S., Rumpl, W., & Nordsieck, K. H. 1977, ApJ, 217, 425, doi: 10.1086/155591
  * Men’shchikov et al. (1999) Men’shchikov, A. B., Henning, T., & Fischer, O. 1999, ApJ, 519, 257, doi: 10.1086/307333
  * Molina et al. (2003) Molina, A., Moreno, F., & Martínez-López, F. 2003, A&A, 398, 809, doi: 10.1051/0004-6361:20021742
  * Mori & Kataoka (2021) Mori, T., & Kataoka, A. 2021, ApJ, 908, 153, doi: 10.3847/1538-4357/abd08a
  * Ngoc et al. (2023) Ngoc, N. B., Diep, P. N., Hoang, T., et al. 2023, ApJ, 953, 66, doi: 10.3847/1538-4357/acdb6e
  * Okuzumi & Tazaki (2019) Okuzumi, S., & Tazaki, R. 2019, ApJ, 878, 132, doi: 10.3847/1538-4357/ab204d
  * Pattle & Fissel (2019) Pattle, K., & Fissel, L. 2019, Frontiers in Astronomy and Space Sciences, 6, 15, doi: 10.3389/fspas.2019.00015
  * Pinte et al. (2016) Pinte, C., Dent, W. R. F., Ménard, F., et al. 2016, ApJ, 816, 25, doi: 10.3847/0004-637X/816/1/25
  * Quillen & Trilling (1998) Quillen, A. C., & Trilling, D. E. 1998, ApJ, 508, 707, doi: 10.1086/306421
  * Rebull et al. (2004) Rebull, L. M., Wolff, S. C., & Strom, S. E. 2004, AJ, 127, 1029, doi: 10.1086/380931
  * Reissl et al. (2016) Reissl, S., Wolf, S., & Brauer, R. 2016, A&A, 593, A87, doi: 10.1051/0004-6361/201424930
  * Roberge et al. (1993) Roberge, W. G., Degraff, T. A., & Flaherty, J. E. 1993, ApJ, 418, 287, doi: 10.1086/173390
  * Sadavoy et al. (2019) Sadavoy, S. I., Stephens, I. W., Myers, P. C., et al. 2019, ApJS, 245, 2, doi: 10.3847/1538-4365/ab4257
  * Sharkawi et al. (2021) Sharkawi, M., Basu, S., & Machida, M. N. 2021, in American Astronomical Society Meeting Abstracts, Vol. 53, American Astronomical Society Meeting Abstracts, 329.01
  * Stephens et al. (2014) Stephens, I. W., Looney, L. W., Kwon, W., et al. 2014, Nature, 514, 597, doi: 10.1038/nature13850
  * Stephens et al. (2017) Stephens, I. W., Yang, H., Li, Z.-Y., et al. 2017, ApJ, 851, 55, doi: 10.3847/1538-4357/aa998b
  * Stephens et al. (2023) Stephens, I. W., Lin, Z.-Y. D., Fernandez-Lopez, M., et al. 2023, arXiv e-prints, arXiv:2311.08452, doi: 10.48550/arXiv.2311.08452
  * Tang et al. (2023) Tang, Y.-W., Dutrey, A., Koch, P. M., et al. 2023, ApJ, 947, L5, doi: 10.3847/2041-8213/acc45b
  * Tazaki et al. (2017) Tazaki, R., Lazarian, A., & Nomura, H. 2017, ApJ, 839, 56, doi: 10.3847/1538-4357/839/1/56
  * Truong et al. (2023) Truong, B., Hoang, T., Chau Giang, N., et al. 2023, arXiv e-prints, arXiv:2308.01215, doi: 10.48550/arXiv.2308.01215
  * Ueda et al. (2021) Ueda, T., Kataoka, A., Zhang, S., et al. 2021, ApJ, 913, 117, doi: 10.3847/1538-4357/abf7b8
  * Velikhov (1959) Velikhov, E. P. 1959, Soviet Journal of Experimental and Theoretical Physics, 9, 995
  * Warren & Brandt (2008) Warren, S. G., & Brandt, R. E. 2008, Journal of Geophysical Research (Atmospheres), 113, D14220, doi: 10.1029/2007JD009744
  * White & Hillenbrand (2004) White, R. J., & Hillenbrand, L. A. 2004, ApJ, 616, 998, doi: 10.1086/425115
  * Wriedt (2012) Wriedt, T. 2012, in The Mie Theory: Basics and Applications, ed. W. Hergert & T. Wriedt, Vol. 169 (Springer), 53, doi: 10.1007/978-3-642-28738-1_2
  * Yang et al. (2017) Yang, H., Li, Z.-Y., Looney, L. W., Girart, J. M., & Stephens, I. W. 2017, MNRAS, 472, 373, doi: 10.1093/mnras/stx1951
  * Yang et al. (2019) Yang, H., Li, Z.-Y., Stephens, I. W., Kataoka, A., & Looney, L. 2019, MNRAS, 483, 2371, doi: 10.1093/mnras/sty3263
  * Yen et al. (2019) Yen, H.-W., Gu, P.-G., Hirano, N., et al. 2019, ApJ, 880, 69, doi: 10.3847/1538-4357/ab29f8
  * Zhang et al. (2023) Zhang, S., Zhu, Z., Ueda, T., et al. 2023, ApJ, 953, 96, doi: 10.3847/1538-4357/acdb4e