Evidence of Grain Alignment by Magnetically Enhanced Radiative Torques from Multiwavelength Dust Polarization Modeling of HL Tau

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 9\times10^2$. Maximum grain sizes probed at wavelengths $\lambda$ = 0.87, 1.3, and 3.1$\,$mm are constrained at $\sim$ 60, 80, and 90$\,\mu$m, respectively.


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 momen-tum 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 thermal emission from aligned grains has been known as the leading tool 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 ∼ 1 Myr) (White & Hillenbrand 2004) located in the Taurus molecular cloud, which is ∼ 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).Selfscattering 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 (grain size a ≳ 10 µ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 ( â1 ) with the angular momentum (J ), so-called internal alignment (IA); and the alignment of J with a preferred direction, so-called external alignment (EA).Internal alignment can be achieved by internal relaxation (INR), e.g., Barnett relaxation (Roberge et al. 1993) and inelastic relaxation (Molina et al. 2003).External alignment, on the other hand, can be accomplished by grain precession (e.g., Larmor precession, radiative precession, and mechanical precession).Grains achieving external alignment with Larmor precession, radiative precession, or mechanical precession have J aligned along the direction of the magnetic field (B), radiation beam (k), or gas flow (ν).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.
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 wellstudied ALMA Bands 3, 6, and 7 corresponding to the wavelengths λ = 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 modeling.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.

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 λ = 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-scattering1 , 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 λ = 3.1, 1.3, and 0.87 mm.The details of the information regarding the observational data used in this paper are summarized in Table 1.The observational data is taken directly from the ALMA archive, which has been reduced and calibrated by ALMA staff.We only further perform debiasing the polarized intensity and polarization degree via the method described in Hull et al. (2014) and Hull & Plambeck (2015).The ALMA minimum detectable polarization degree is expected at 0.1%, and defined as three times the systematic calibration uncertainty.Hence, we use the error of 0.033% when the polarization degree uncertainty is less than this value.
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 mm and (2) the polarization degree is generally at the level of ∼1% for the whole three wavelengths.

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 mix-ture 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 max , and a fixed minimum grain size a min = 0.05 µm.We assume the grain population to be the same throughout the disk.As discussed in Lin et al. (2022), the value of a 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 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πa/λ << 1.Indeed, later on, this is shown to be the case with the maximum grain sizes found to be ∼60, 80, and 90 µm for the three probing wavelengths 0.87, 1.3, and 3.1 mm, respectively.
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 and RAT efficiencies are then precalculated 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 max .Selfscattering polarization is most effective when a max ∼ λ/2π (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 (socalled low-J attractors).Grains with iron inclusion are dubbed paramagnetic material (PM) or superparamagnetic material (SPM).Considering paramagnetic grains, we assign the value of the fraction of iron atoms diffusely distributed within grains as f p = 0.1.For superparamagnetic grains, we fix the volume filling factor of iron clusters to ϕ sp = 0.1, which corresponds to an iron abundance of ∼30% presented in the form of iron clusters2 and vary the number of iron atoms per cluster N cl to describe different degree of grain alignment enhancement due to superparamagnetic relaxation.The possible value of N cl spans from ∼ 20 to 10 5 (Jones & Spitzer 1967).
Note that for internal relaxation of grains, although the case of RATs quickly bringing J to align with â1 (fast internal relaxation, Lazarian & Hoang 2021) has been well-studied, the case of internal relaxation timescale being much longer than gas damping timescale (slow internal relaxation) is not as well-defined (Hoang et al. 2022).Hence, in this study, for the case of slow internal relaxation, we fix the values of internal alignment efficiency Q X as follows.When grains are aligned at high-J attractors, the internal alignment efficiency is Q X,high−J = 0.4.On the other hand, when grains are at low-J attractors, they can be aligned at either right or wrong internal alignment (Hoang & Lazarian 2009).In the case of right internal alignment, J is aligned parallel to â1 , which coincides with the minor axis of the ellipsoidal grains.In the case of wrong internal alignment, their alignment direction is perpendicular instead.When Q X,low−J = −0.5, grains are defined to achieve perfect wrong internal alignment (Giang et al. 2023).We set Q X,low−J = 0.4 for the case of grains aligning at right internal alignment, and Q X,low−J = −0.4 at wrong internal alignment.The choice of the values of Q X,high−J and Q X,low−J at wrong internal alignment is taken as the median value of the constraint which will be given in Section 4.3.The value of Q X,low−J at right internal alignment is set to be similar to Q X,high−J based on the assumption that J would have a negligible effect on the grains' internal alignment efficiency when they are aligned with slow internal relaxation.Throughout the paper, we will call wrong internal alignment/right internal alignment for short in referring to the case of grains aligned with wrong or right internal alignment when they are at low-J attractors and aligned with slow internal relaxation.The procedure to calculate the alignment degree of dust grains 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 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, 1.3 mm, and 0.87 mm, which corresponds to ALMA observation band 3, 6, and 7, respectively.
First, we assume a single radiation source to be the protostar at the disk center.The protostar has the effective temperature T ⋆ = 4000 K and radius R ⋆ = 7R ⊙ , which corresponds to a luminosity L ⋆ ∼11L ⊙ (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 max = 100 µ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) where R is the distance to the disk center, R c = 78.9AU and γ = −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, where z is the height from the disk midplane and h d is the disk dust scale height whose profile is taken as with h 0 = 10 AU being the scale height at 100 AU from the center of the disk, and β = 1.15 the exponent factor characterizing the flaring effect of the disk (Pinte et al. 2016;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 selfscattering (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 d by a factor f settle = 10, which is defined as the dust-settling factor from the disk scale height, following Kataoka et al. (2016).In this calculation, the temperature distribution is taken identical to that from the radiative heating calculation with the original scale height.We assume that 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 max = 100 µm, we take the total dust mass of the system M d as a calibrating parameter to reproduce the three Stokes I emis-sions at the three wavelengths mentioned earlier.M d = 1.3 × 10 −2 M ⊙ is found to be the best-fitted value, which produces total emissions comparable to the observations as shown in Figure 2.This value is similar to the level of M d ∼ 1 × 10 −2 M ⊙ as estimated in Mori & Kataoka (2021) and Zhang et al. (2023), in the case where grain size at ∼ 100 µm is employed.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, under the common assumption of gas-to-dust ratio equal to 100 (i.e., the typical value in the ISM), the total disk mass (gas+dust) would be M disk ∼ 1.3 M ⊙ .This total disk mass is similar to that of the central star (M ⋆ = 1.7 M ⊙ , Pinte et al. 2016), which would imply that the disk is extremely unstable and therefore, unrealistic.On the other hand, Kwon et al. (2011), Liu et al. (2017), and Carrasco-González et al. (2019) give the value of total dust mass is M d = 1 × 10 −3 M ⊙ .Assuming gas-todust ratio equal to 100, the value of total disk mass is M disk = 0.1 M ⊙ , which is an order-of-magnitude lower than our estimated value and much more reasonable.Thus, to obtain the gas density, we normalize the density profile with M g = 0.1 M ⊙ , which is equivalent to a gas-to-dust ratio equal to ∼ 10.
The radiation field and temperature profile resulting from radiative transfer calculation with this disk model are shown in Appendix A. We assume gas temperature to be the same as the dust temperature for simplicity.
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.

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 λ = 3.1 mm with fiducial grain size a max = 100 µm.
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 internal alignment.The different grain alignment scenarios considered in this section are summarised in Table 3.
Initially, we examine the ideal scenario (top left panel) in which all grains have perfect internal alignment (i.e., the fraction of grains aligned at high-J attractors is f 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 ∼10%.
In the subsequent scenarios, we take into account the effect of slow internal relaxation and enhanced (super)paramagnetic alignment by MRAT.We consider both cases of right and wrong internal alignment (middle and right panels).The effect of paramagnetic grains (top panels) and superparamagnetic grains with N cl = 10 3 and 10 4 (middle and bottom panels) are also included for comparison.
In the case of right internal alignment (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 paramagnetic material, the polarization degree is ∼0.01%.In this case, the alignment of very large grains 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 superparamagnetic grains, the alignment is considerably more efficient (i.e., the polarization degree is two orders of magnitude higher than in the paramagnetic grains scenario).With increasing N cl , higher values of f 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 internal relaxation in the same direction, their internal alignment efficiencies are the same whether at high-J or low-J attractors.Thus, the polarization degrees observed at N cl = 10 3 and 10 4 are seemingly identical.
In the case of wrong internal alignment (right panels), the polarization morphology is more diverse.Grain alignment is weak for paramagnetic grains with polarization degrees as low as ∼0.01%.The polarization vectors show an azimuthal pattern owning to the inefficient internal alignment of grains' â1 perpendicular to J at low-J attractors with slow internal relaxation and external alignment of J along B. superparamagnetic grains, conversely, induce a higher polarization degree with increasing N cl .However, the polarization pattern alters when we keep increasing N cl .Indeed, the increased level of iron inclusion also causes the increased abundance of grains at high-J attractors (i.e., higher f high-J ) and more efficient Barnett relaxation, which would drive them to have more efficient internal alignment 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 cl = 10 4 .
Note that the features of disk polarization only induced by grain alignment observed above also hold for observations at both wavelengths λ = 1.3 and 0.87 mm.

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 internal alignment with superparamagnetic grains having N cl at the level of ∼10 3 Notes.χ 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 internal alignment/wrong internal alignment when they are at low-J attractors and slow internal relaxation.
Table 4. Summary of the correlation between the Models and Observations.
Grain Magnetic Property (SPM) Notes.We only consider superparamagnetic (SPM) grains in the model.Each column represents the model with a different value of the number of iron atoms per cluster (N cl ).Each row represents the model with different grain sizes at each wavelength.Each model with a couple value of N cl and amax has two validating factors AM and cP.AM is the average value of alignment measure across the image.cP is the correlation factor of pixel-by-pixel polarization degree between the models and observations.The closer the values of AM and cP to unity, the better the model compared to the observations.can produce the azimuthal polarization morphology and the polarization degree at ∼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 P ) while computing the Alignment Measure (AM ) as introduced in González-Casanova & Lazarian (2017) to compare the polarization orientation.
where θ r = |θ mod − θ 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 AM to represent the overall correlation between polarization vectors of the model and observations.Finally, the values of c P and AM are compared between each model to find the best model for the observations.Note that the closer the values of c P and AM to unity, the better the agree-ment between the model and the observations.Table 4 summarizes the values of c P and AM of the typical models that we considered.As shown in Table 4, the model that best reproduces the disk polarization observed by ALMA is found at N cl = 9 × 10 2 ; a max = 90, 80, and 60 µm for λ = 3.1, 1.3, and 0.87 mm, respectively.The disk polarization corresponding to this model is displayed in Figure 4.As expected, very large grains with wrong internal alignment 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

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 internal alignment for very large grains with slow internal relaxation, but their wrong internal alignment 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 very large grains with slow internal relaxation must have wrong internal alignment.Therefore, dust polarization can still trace the magnetic field within protoplanetary disks.Moreover, due to the wrong internal alignment, the polarization vectors are along the magnetic field direction, so that one does not need to rotate the polarization vectors by 90 • to infer the magnetic field, as in the case of right internal alignment.
However, Yang et al. (2019) and Mori & Kataoka (2021) modeled the HL Tau polarization observed at ALMA Band 3 (λ = 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 internal alignment for grains with slow internal relaxation, 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.Guillet et al. (2020) proposed that in the Mie regime (x = 2πa/λ ∼ 1), the polarized emission can turn into a negative value, which results in the polarization vectors parallel to the toroidal magnetic field for elongated spheroidal grains aligned with right internal alignment.This idea can provide an alternative to our model in explaining the azimuthal polarization pattern observed in the disk.However, the persistence of the azimuthal polarization pattern across a large range of wavelength, from λ = 0.87 mm up to 7 mm (Lin et al. 2023), means that they must all be well within the Mie regime, or the effective grain size must be up to a max ∼ 1 mm.This is in contradiction with our constrained grain size where they are only placed at ∼ 100 µm-scaled.Hence, the grain size constraint would have implications in determining whether the ALMA observations are in the Rayleigh or Mie regime.Alternatively, it is also worth noticing that as soon as λ > 2πa, the polarized emission of aligned grains is predicted to become positive.Therefore, observing whether there is a polarization flip to a radial pattern at a shorter wavelength than 0.87 mm or longer than 7 mm can help further clarify which is the correct explanation for the observed disk polarization.
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 degree and spectral index of thermal emission (Mori & Kataoka 2021;Ueda et al. 2021;Lin et al. 2023;Zhang et al. 2023).Earlier in the current study, a max was constrained to be different when probing at different wavelengths (90 µm, 80 µm, and 60 µm at λ = 3.1 mm, 1.3 mm, and 0.87 mm, respectively).However, they are all at ∼ 100 µm, which is still greatly differs from the 1 mm grain size constrained from spectral index SED fitting.Further discussions are presented in Section 4.5.

Constraints on Iron Abundance within Grains in the form of Clusters
Our results in Section 3 show that paramagnetic grains produce a very low polarization degree, ∼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 ϕ sp = 0.1 and varied N cl to obtain the different magnetic susceptibility of grains (see e.g.Equation 3in Giang et al. 2023).In this section, we treat ϕ sp as a free parameter to investigate its effect on disk polarization.For simplicity, we test the effect using q ′ 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 ′ 100 would have a positive value when the polarization vector is radial and negative when it is azimuthal.We demonstrate the effect at λ = 3.1 mm, where alignment is expected to be the dominant component over self-scattering.The value of −q ′ 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 ϕ sp , the curve describing the dependence of −q ′ 100 on N cl follows a parabolic shape.The reason behind this is as described in Section 3.1 in which increased N cl not only increases polarization degree through enhancing external alignment but also increases f high-J which drives polarization vector from azimuthal to radial direction.The value of −q ′ 100 is at its peak when the effect of superparamagnetic alignment from N cl is not negligible so that grains can efficiently align with the magnetic field, yet not too significant so that f high-J is small enough to keep the polarization vectors at the azimuthal direction.With the higher value of ϕ sp , this peak shifts to the position corresponding to the lower value of N cl .
This test puts a constraint on ϕ sp with its minimum value at ∼0.05, which corresponds to ∼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 cl ∼ 9 × 10 2 is consistent with the −q ′ 100 reference when considering ϕ sp = 0.1 as in our dust model, but the range of N cl is also found to be at ∼2 × 10 2 -2 × 10 3 within the valid values of ϕ sp .

Constraints on Internal Alignment Efficiency for Slow Internal Relaxation
Slow internal relaxation is generally responsible for internal alignment of large dust grains within protoplanetary disks since large grains above ∼10 µm and high density are ubiquitous in these environments (Hoang 2022).However, the efficiency of internal alignment for the case of slow internal relaxation is not yet available from theory (Hoang & Lazarian 2009;Hoang et al. 2022).Usually, the efficiency of internal alignment by slow internal relaxation is represented by two quantities Q X,low−J and Q X,high−J .In Section 3, our modeling results showed that Q X,low−J must have a negative value (i.e., wrong internal alignment) to reproduce the azimuthal polarization pattern observed toward the HL Tau disk.On the contrary, grains at high-J attractors always have the right internal alignment or Q X,high−J is always positive (Hoang & Lazarian 2009).It indicates that the directions of internal alignment of grains aligned at low-J and high-J attractors are opposite.
Here, we will attempt to empirically constrain the magnitude of Q X,low−J and Q X,high−J for slow internal relaxation using the observed polarization data.To do so, we consider the extreme scenarios where parameters involved in the alignment process are set to their respective extremes.We initially fix ϕ sp = ϕ sp,max = 0.3.To obtain the maximum value of |Q X,high−J |, we set |Q X,low−J | = |Q X,low−J,max | = 0.5.On the other hand, |Q X,high−J | = |Q X,high−J,min | ∼ 0.01 is set to constrain the minimum value of |Q X,low−J |.Following Section 4.2, we use −q ′ 100 at λ = 3.1 mm as the referenced polarization degree for the disk polarization.Its value is calculated for different values of Q X,high−J or Q X,low−J and over the range N 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 X,high−J | is ∼0.8, while the minimum value of |Q X,low−J | is at ∼0.15.

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 (socalled 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 internal alignment, 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 max ∼ 100 µm, and assuming a gas surface density of Σ g ∼ 5 g cm −2 , the Stokes number is estimated to be St ∼ 5.4 × 10 −3 .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 internal alignment, a circular polarization pattern is produced.On the other hand, if grains are aligned with the wrong internal alignment, 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 1 mm as constrained in Carrasco-González et al. (2019), which corresponds to St ∼ 5.4 × 10 −2 ≪ 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 current understanding on grain alignment theory that can explain the polarization signatures observed in the HL Tau disk by ALMA.

Grain Size Population
Grain size was constrained primarily from the dust polarization of the self-scattering component.Therefore, in order to address the detailed view of the grain size constraint given in Section 3, we provide the polarization efficiency from dust scattering of each dust model, as a function of wavelength in Figure 6.Polarization efficiency is represented by the polarization degree at 90 • scattering times the albedo Pω following Kataoka et al. (2015).
We could reproduce the observed dust polarization assuming 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 valwhen tracing the dust by longer wavelengths, namely a max = 60, 80, and 90 µm at λ = 0.87, 1.3, and 3.1 mm, respectively.
However, this argument is only valid 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.Irregular grain structures have been studied to alleviate the grain size conflict between polarization and spectral index SED fitting (Lin et al. 2023;Zhang et al. 2023).Moreover, it offers a promising explanation for the spectral dependence of the polarization degree across (sub)millimeter wavelengths as in Figure 6.According to the literature, grain size can be larger than 1 mm depending on the grain porosity.Furthermore, grains with large porosity have less steeper-slope spectral dependence of scattering efficiency, which generates a similar trend for polarization degree resulted from self-scattering.Therefore, it is crucial to incorporate the effect of grain porosity into detailed modeling of grain alignment to have a more complete picture of HL Tau polarization.

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 λ = 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 ∼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 ∼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 internal relaxation 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).One of which is the high polarization degree at the gaps, where their further modeling suggests that the intrinsic polarization could be up to 10-15 % without beam smearing effect.This is potentially 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 constant gas-to-dust ratio throughout the disk is no longer acceptable.The accumulation of magnetic field inside the gaps also needs to be considered if they were formed due to magnetohydrodynamic (MHD) effects (Flock et al. 2015).Moreover, the effect of mechanical torque alignment (Hoang et al. 2022) has not been taken into account.Its incorporation into modeling still faces difficulty due to its significant sensitivity to grain shape (Reissl et al. 2023).Hence, additional studies are needed to fully reproduce the high-resolution HL Tau dust polarization observed by ALMA.

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, especially its scattering property.
The assumed geometrically thick disk in intrinsic dust thermal emission calculation may cause unwanted nearfar 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.

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.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 internal relaxation, grains can either align with wrong or right internal alignment, which manifests in azimuthal or radial polarization patterns, respectively.Additionally, paramagnetic grains cannot reproduce the observed polarization data and thus cannot be the primary dust material in HL Tau.
2. A mixture of self-scattering and MRAT grain alignment mechanism where grains are aligned with wrong internal alignment at low-J attractors and slow internal relaxation can reproduce the dust polarization observed in HL Tau by ALMA at (sub)millimeter wavelengths.
3. Grains are shown to be made of superparamagnetic materials with the number of iron atoms in-side each cluster is N cl ∼ 9 × 10 2 .The abundance of embedded iron inside grains identified from our model is at ∼16% minimum.
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.
µ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 (A1)  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 λ = 0.87 mm as in Stephens et al. (2023), which yields the surface density profile where Σ a = 2.3, R 0 = 10 AU, p = 0.5, and the parameters for the rings are given in Table 5.In this investigation, instead of using temperature calculated from radiative transfer calculation, for simplicity, we adopt directly the temperature profile as in the literature: T (R) = 110 × [R/(10 AU )] −0.5 .We then normalize the dust density by the total dust mass of M d ∼ 5 × 10 −3 M ⊙ to produce a comparable Stokes I emission at λ = 0.87 mm as in Stephens et al. (2023).The gas density, on the other hand, is obtained by normalizing the density profile with the total gas mass of M g = 0.1 M ⊙ similar to that done in Section 2.2.This is equivalent to a gas-to-dust ratio of ∼ 20.
For the polarization calculation, we adopt the same dust models as constrained in Section 3.2 with N cl = 9 × 10 2 , and a max = 60 µm.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 3 .
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 high−J = 1).The central and right panels include the effects of slow internal relaxation and enhanced (super)paramagnetic relaxation.Central panels show the case of right internal alignment (rIA), while the right ones correspond to wrong internal alignment (wIA).Three different dust models are taken into consideration with paramagnetic grains (top) and superparamagnetic grains which have N cl = 10 3 (middle) and 10 4 (bottom), respectively.

Figure 4 .
Figure4.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 internal alignment, N cl ∼ 9 × 10 2 ; and amax = 90, 80, and 60 µm at λ = 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 pixelto-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 (AM , cP) of the model are ∼ (0.93, 0.68), (0.93, 0.63), and (0.93, 0.09) at λ = 3.1, 1.3, and 0.87 mm, respectively.

Figure 5 .
Figure 5. Dependence of polarization degree −q ′ of Stokes Q due to grain-alignment-only in the principal frame at 100 AU along the disk's minor axis on N 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 ϕsp with the original dust model set-up.The middle panel considers the varying values of Q X,high−J with the extreme scenario when ϕsp = 0.3 and Q X,low−J = −0.5.The bottom panel is plotted with different values of Q X,low−J when ϕsp = 0.3 and Q X,high−J = 0.01.Only the value of ϕsp ≳ 0.05; |Q X,high−J | ≲ 0.8, and |Q X,low−J | ≳ 0.15 can reach the threshold.

Figure 6 .
Figure 6.Polarization degree at 90 • scattering times the albedo Pω against the observing wavelength at the threegrain size models amax = 60 µm, 80 µm, and 90 µm.The three vertical lines correspond to the three ALMA wavelengths 0.87 mm, 1.3 mm, and 3.1 mm, where observations were used as references in our work.The three red dots correspond to the maximum grain size constrained at each observing wavelength.

Figure 7 .
Figure 7. From left to right: Radiation flux J rad , spectrum weighted wavelength λ, and spectrum weighted anisotropy parameter γ.The calculation is made assuming a single radiation source to be the central protostar with luminosity L⋆ ∼ 11L⊙.

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

Table 1 .
Information of the observational data used in this paper

Table 2 .
Key Parameters of the Dust Model

Table 5 .
Parameters for the identified rings in HL Tau.Ring number Σi Ci [AU] Wi [AU]