Probing Magnetic Fields in Protoplanetary Disk Atmospheres through Polarized Optical/IR Light Scattered by Aligned Grains

Magnetic fields play essential roles in protoplanetary disks. Magnetic fields in the disk atmosphere are of particular interest, as they are connected to the wind-launching mechanism. In this work, we study the polarization of the light scattered off of magnetically aligned grains in the disk atmosphere, focusing on the deviation of the polarization orientation from the canonical azimuthal direction, which may be detectable in near-IR polarimetry with instruments such as VLT/SPHERE. We show with a simple disk model that the polarization can even be oriented along the radial (rather than azimuthal) direction, especially in highly inclined disks with toroidally dominated magnetic fields. This polarization reversal is caused by the anisotropy in the polarizibility of aligned grains and is thus a telltale sign of such grains. We show that the near-IR light is scattered mostly by $\mu$m-sized grains or smaller at the $\tau=1$ surface and such grains can be magnetically aligned if they contain superparamagnetic inclusions. For comparison with observations, we generate synthetic maps of the ratios of $U_\phi/I$ and $Q_\phi/I$, which can be used to infer the existence of (magnetically) aligned grains through a negative $Q_\phi$ (polarization reversal) and/or a significant level of $U_\phi/I$. We show that two features observed in the existing data, an asymmetric distribution of $U_\phi$ with respect to the disk minor axis and a spatial distribution of $U_\phi$ that is predominantly positive or negative, are incompatible with scattering by spherical grains in an axisymmetric disk. They provide indirect evidences for scattering by aligned non-spherical grains.


INTRODUCTION
The evolution of protoplanetary disks is generally thought to be determined by magnetic fields through either magnetorotational instability (Balbus & Hawley 1991) or magnetized disk wind (Blandford & Payne 1982). The magnetic field structure in the disk atmosphere 1 is of particular interest in understanding the wind launching mechanisms. Strongly magnetized disks tend to launch magneto-centrifugal wind (MCW; Bland-ford & Payne 1982) with rigid and mostly poloidal magnetic field lines in the disk atmosphere. Weakly magnetized disks tend to launch magneto-thermal disk winds , and rely on the vertical gradient of the magnetic pressure from the toroidal field to launch the wind. The magnetic field strength is also related to other interesting issues of protoplanetary disks, such as accretion rates.
Dust grains can trace the magnetic field if they are magnetically aligned (Andersson et al. 2015;Lazarian 2007). While polarized (sub)millimeter dust thermal emission has been proven to be a powerful tool to study magnetic fields on scales larger than the disks (e.g., Planck Collaboration et al. 2020; Hull & Zhang 2019 and references therein), its application on the disk scale has not been as successful because we see mostly scattering-induced polarization at shorter (e.g., 870 µm) wavelengths Hull et al. 2018;Bacciotti et al. 2018;Dent et al. 2019) and complicated non-magnetic-origin patterns at longer (e.g., 3 mm) wavelengths (Kataoka et al. 2017;Harrison et al. 2019). Yang (2021) showed that the Larmor precession in the disk midplane is likely too slow to ensure magnetic alignment, which is likely the reason behind the failure of polarized dust thermal emission to trace magnetic fields in disks. However, Yang (2021) proposed that the micronsized dust grains in the disk atmosphere can potentially be magnetically aligned. This idea was partially supported by Li et al. (2016), who found polarized radiation at 10.3 µm, which may be explained in part by thermal emission from grains aligned with the magnetic field in the disk atmosphere.
The Spectro-Polarimetric High-contrast Exoplanet REsearch (SPHERE) at the Very Large Telescope (VLT) has been used for high-resolution polarimetric observations of protoplanetary disks in scattered near-IR light (Benisty et al. 2015;Avenhaus et al. 2018;Garufi et al. 2020Garufi et al. , 2022. These studies usually focus on the azimuthal Stokes parameter Q φ component to produce high resolution images of protoplanetary disks. The U φ component is often observed to be small and its divergence from zero is a sign of deviation from the simplest single Rayleigh scattering. Theoretically, Canovas et al. (2015) used a generic transition disk model with large dust grains to show that the scattered light from moderately inclined disks can possess a significant U φ . Whitney & Wolff (2002) studied the scattering by aligned grains, but focused on the circular polarization in the protostellar envelope. Our focus in this paper is on the near-IR photons scattered by aligned dust grains in the atmosphere of protoplanetary disks, which can potentially have polarization patterns different from the commonly expected pure azimuthal ones.
The structure of the paper is as follows. In Section 2 and Appendix A, we discuss the polarization of scattered light in the grain's frame, focusing on the difference with spherical dust grains. In Section 3, we calculate the polarization pattern in a disk configuration. In Section 4, we discuss the grain size distribution and grain alignment at the optical depth τ = 1 surface in a generic protoplanetary disk model. In Section 5, we discuss our results, including detectability of the U φ produced by the scattering of magnetically aligned grains, the differences between the patterns produced by aligned grains and multiple scattering of spherical grains, and implications for observations. We summarize our results in Section 6.

BASIC PHYSICS
Before considering the general case in a disk environment. we first illustrate the basic physics of why the polarization orientation of the light scattered by aligned grains can be different from that by spherical grains in a simple setup in the grain's frame. We consider a Cartesian coordinate system xyz (see Figure 1). Letẑ be the propagation direction for incoming light. Without loss of generality, we fix the scattered light in the xOz plane, with O being the particle location. Let θ be the scattering angle, the scattering directional vector is then (sin θ, 0, cos θ). Since light is transverse wave, its E vector can be decomposed into two components that are perpendicular to the scattering direction. We call the component perpendicular to both incoming and scattering lights E 1 (in theê 1 direction), and the other component E 2 (in theê 2 direction). We can easily see that the E 1 direction is also theŷ direction of our coordinate system. If the dust particle is spherical, the scattered light can either be polarized along E 1 or along E 2 direction, due to the symmetry of this scattering geometry. If we define the Stokes parameters such that fully polarized light with polarization along E 1 has Q = I, the Stokes U is always zero for spherical dust grains. If the dust particle is not spherical, and ifx andŷ are not the principle axes of the dust particle, the light will no longer be polarized along either E 1 or E 2 direction, which leads to non-zero Stokes U component. The break of the symmetry is the reason why the polarization can deviate from the spherical case.
It turns out that the deviation is always maximized in the forward and backward scattering directions. For the forward scattering, we have the scattering direction being the same as the incoming light direction (ẑ). If the dust particle is spherical, the scattered light will always be non-polarized, because of the symmetry. If the dust particle is non-spherical, say being elongated alongx direction, then the scattered light would be polarized along thex direction, which is qualitatively different from the spherical case. If we change the scattered light slightly away from the forward scattering direction with a small scattering angle θ in the xOz plane, then the light scattered by a spherical dust grain would be polarized alongŷ direction, in the often assumed Rayleigh scattering regime. For dust grains elongated alongx direction, on the contrary, this small deviation in scattering angle is not enough to change the polarization state of the scattered light, and the scattered light is still polarized alongx direction. Hence the angle difference between polarization orientation of light scattered by spherical dust grains and the polarization orientation of light scattered by elongated dust grains near the forward scattering direction can always be as large as 90 • , as they are perpendicular to each other in the set-up we discussed above.  Figure 1. The geometry of a simple setting in the grain's frame. The photon propagating in theẑ direction is scattered by the dust grain sitting at the origin O towards thê ns direction in the xOz plane. The scattering angle is θ. The scattered light is decomposed into E1 and E2 directions. Note that it is often assumed that the scattered light is polarized along E1 direction.
In Appendix A, we discuss the angle difference in the grain's frame in more detail. We show that the angle difference between the spherical dust grain and the elongated dust grain can easily reach 10 • in the grain's frame. We also show that the angle difference increases with the angle between the symmetry axis of the dust particle and the incoming light direction and with the aspect ratio of the dust particle. It also depend on the compositions of the dust grains. In what follows, we focus on the angle difference in a disk environment where most of the relevant observations are carried out.

Model prescription
Calculations in Section 2 and Appendix A focus on the polarization of the scattered light in the grain's frame. While more physically intuitive, it is not directly connected to the observed polarization. Here we study the scattered light in a protoplanetary disk, focusing on the deviation of polarization orientation from the direction perpendicular to the stellar light, the expected polarization orientation in the small spherical particle regime.
To calculate the polarization orientation in the scattered light from the surface of a protoplanetary disk, we consider a set-up shown in Figure 2. The black horizontal arrow represents the disk midplane. The red dot The black dot represents the central star. The red dot represents the scattering dust grain. Hencen1 connecting these two dots is the incoming light direction. The scattered light propagates alongn2, making an angle of i with the z axis. The local magnetic field direction is prescribed by the angles θB and φB. See text for more details.
represents the dust particle that scatters light fromn 1 direction towards then 2 direction, which makes an angle i with the z direction, the direction perpendicular to the disk midplane; i is simply the disk inclination angle (i = 0 is face-on). Note that only one dust grain is plotted in the figure, but it represents a ring of dust grains all of which have the same cylindrical radius R from the star and height H above the disk midplane. We assume that the local magnetic field B makes an angle θ B with the z direction, and an azimuthal angle φ B from the x direction, such that φ B = 0 implies a pure poloidal magnetic field. As usual, we assume that the grains are aligned with their shortest axis along the magnetic field direction. If the grains are spinning around the B field, they would be effectively oblate after assemble-averaging independent of their intrinsic shapes.
We adopt the same dust composition as used in Appendix A, which is the same as the one adopted by Birnstiel et al. (2018), and assume the dipole approximation for simplicity. We consider 5 parameters: H/R, i, θ B , φ B , and the dust aspect ratio s.

A limiting case: purely toroidal magnetic field
Before studying the polarization pattern in a generic model, we will first consider a limiting case that will help our understanding: a purely toroidal magnteic field alongŷ direction, with θ B = 90 • and φ B = 90 • . The incoming light propagating alongn 1 direction can be decomposed along two directions:ŷ andn 1 ×ŷ. We will denote the dipoles excited by these two components n 1 n 2 B P 1 P 2 Figure 3. The geometry of an oblate grain aligned by a pure toroidal magnetic field. By definition, the magnetic field B is perpendicular to incoming light directionn1. The incoming light excites two dipoles in the grain, denoted as P1 and P2. The scattered light directionn2 makes an angle of ι with P2. Note that we have |P2| > |P1| under this configuration. P 1 and P 2 , respectively. See Figure 3 for a schematic illustration of this setting.
The angle between the scattered light directionn 2 and the P 2 , defined as ι, will be important to our analysis that follows. Note that in the disk frame, we have ι = i + (−)tan −1 (H/R) for the near (far) side of the ring considered in Figure 2.
Because the dust grains are aligned with the magnetic field in P 1 direction, we have P 1 = α 3 E 1 . Since P 1 is always perpendicular ton 2 , we have E s1 ∝ P 1 ∝ α 3 , where E s1,2 is the E vector of scattered light induced by P 1,2 . Similarly, we have P 2 = α 1 E 2 . Since P 2 is making an angle ι withn 2 , we have E s2 ∝ P 2 sin ι ∝ α 1 sin ι. If E s1 > E s2 , the scattered light is polarized along E s1 direction, i.e. the direction perpendicular to bothn 1 andn 2 , which is the canonical (azimuthal) polarization direction for spherical particles. If E s1 < E s2 , the scattered light is polarized along E s2 direction, i.e. the P 2 projected towards the sky plane, the plane perpendicular ton 2 . This is also the direction of the projected n 1 (radial) direction and is perpendicular to canonical polarization direction for spherical particles. The polarization direction is changed by 90 • and the Q φ , the azimuthal Q parameter (de Boer et al. 2020), changes from positive to negative. We call this "polarization reversal". Define the critical angle ι c as ι c ≡ sin −1 (|α 3 |/|α 1 |). For ι > ι c (ι < ι c ), we have E s2 > E s1 (E s2 < E s1 ), and the polarization is (not) reversed. The dependence of the critical angle with dust aspect ratio is shown in Figure 4. In Section 3.3, we discuss a fiducial case of a moderately inclined disk with ι > ι c but not a purely toroidal  in the x-axis is the aspect ratio of dust grains. The ι is the angle between the scattered light propagating direction n2 and the second dipole direction P2 (c.f. Figure 3). In the gray parameter space, the polarization is potentially reversed, i.e. the light is polarized along that radial, rather than the canonical azimuthal direction. magnetic field. We will demonstrate that the polarization reversal still exists (i.e., it is not limited to the pure toroidal magnetic field configuration) and is likely once we reach the critical angle. In Section 3.4, we present a less inclined disk model without polarization reversal. We will show that there is still an appreciable angle deviation from the azimuthal pattern that is potentially detectable with VLT/SPHERE.

Fiducial case
For the fiducial case, we consider H/R = 0.2, i = 60 • , s = 1.5, θ B = 45 • , and φ B = 90 • (Model 1 hereafter). The last two angles imply that the toroidal to poloidal magnetic field ratio is B tor /B pol = 1. We can easily calculate the ι at the near side as ι = i + tan −1 (H/R) = 71 • > ι c = 59 • . The results are shown in Figure 5. In the left panel, we show the polarization orientation at each location on a dust ring of constant cylindrical radius R and height H in the sky plane, with black unilength line segments. To guide the eyes, we also use dotted lines to connect the central star and the scattering grains. In the often assumed case, the polarization is perpendicular to the dotted line, shown as red lines, with only the Q φ component being non-zero. We can see that the polarization is completely reversed near azimuthal angle of 330 • . To show the deviation more quantitatively, we plot the angle difference as a function of the azimuthal angle (of dust grains in the disk frame) in the upper right panel. The polarization fraction, the ratio p ≡ PI/I between the polarized intensity and the total intensity, as a function of the azimuthal angle is shown in the lower right panel, as a solid line. The polarization fraction for scattering by small spherical grains is also shown as a dashed line. We can see that the polarization is still maximized at 100%, while the phase function deviates from the spherical curve slightly. We can see that the polarization reversal location coincides with low polarization points, because it relies on the usually subdominant component P 2 to overwhelm the P 1 to have polarization reversal (see Section 3.2).
The same calculation is repeated on a finer azimuthal grid (360 points to sample the full azimuthal extent) with different combinations of the angles θ B and φ B that specify the magnetic field configurations. The maximum angle difference for each combination is shown as a colormap in Figure 6, with θ B and φ B as the x and y axis, respectively. We can see that the φ B has a very strong effect on the angle difference. At the upper right part of the Figure 6, the maximum angle difference is essentially 90 • (i.e., polarization reversal). The fluctuations are due to the finite resolution of the azimuthal grid, and we have tested that the fluctuations become smaller (with the angle deviation closer to 90 • ) with an increasing number of grid points. Away from the upper right polarization reversal region, the angle difference is still appreciable, and can easily be above 20 • for our fiducial disk inclination of 60 • .
To study the dependence of the polarization reversal region on the observing inclination angle i, we first find for the fiducial i = 60 • all pairs of θ B and φ B that would make the angle difference ∆η = 70 • 2 . The resulting constant angle difference ∆η = 70 • contour is labeled in Figure 6. We then repeat the same calculation for several inclination angles. The resulting contours are labeled in the figure as well. We can see that the parameter space with polarization reversal as marked by the ∆η = 70 • contour increases with increasing inclination angle. In the most inclined case, the polarization reversal is almost inevitable, with only a small region in the lower right corner being not completely reversed (the maximum angle difference is still large). We would like to note that highly inclined disks are also subject to strong forward scattering and potentially multiple scattering (Canovas et al. 2015). These effects, not considered in this simple first study, may change the results substantially.

Less inclined case
We have shown that for highly inclined disks, the polarization can be reversed, with an orientation along the radial, rather than the canonical azimuthal, direction. We now focus on less inclined cases, particularly the dependence of the maximum angle difference with the observing inclination angle i in the disk frame.
As an example of smaller inclination angle, we assume H/R = 0.2, s = 1.5 and i = 45 • . We also assume a less extreme configuration for the magnetic field with θ B = 30 • and φ B = 45 • . The results for this model (Model 2 hereafter) are shown in Figure 7. We can see that for this less inclined disk model, the maximum angle difference is 11 • . Note that the toroidal to poloidal magnetic field ratio is B tor /B pol = 1/ √ 7 ≈ 0.38. In Section 5.1, we derive a rough error estimate formula for the angle deviation as δη ≈ 1/(2SNR), with SNR being the signal to noise ratio. The angle difference of 11 • can be detected at a signal to noise level of δη/∆η = 0.38 SNR. If we ask for 3σ detection for the angle difference, we need only an SNR of 7.8, easily achievable with VLT/SPHERE.
If we allow the magnetic field configuration to change while fixing the other parameters, we get the maximum angle difference in the (θ B , φ B ) map shown in Figure 8. We can see that the trend is similar to our fiducial model, except that there is no polarization reversal in this map. The maximum angle difference increases as we increase the φ B . Even for this moderate inclination angle of i = 45 • , as φ B approaches 90 • , the angle difference can easily reach 20 • or even 30 • . The azimuthal angle of the magnetic field, φ B , also determines the ratio of the toroidal component to the poloidal component of the magnetic field, which is very important in determining the wind launching mechanisms. Disks with large magnetization have rigid magnetic field lines, and tend to launch magnetocentrifugal winds with a small toroidal component. Weakly magnetized disk, on the other hand, will have magnetic field lines winded up into mostly toroidal configuration first. The disk wind is then launched due to a vertical gradient of magnetic pressure .
The maximum angle differences in the (θ B , φ B ) map can be calculated for different inclination angles i while fixing H/R = 0.2. The results for three different dust aspect ratios s = 0.1, 1.5, and 2.0 are shown in Figure 9. We can see that the behaviors are similar among different s: the maximum angle difference gradually increases before reaching about 30 • , then it suddenly jumps to 90 • and enters the polarization reversal regime, the gray region in Figure 4. So for disks with small inclination angles, say i < 20 • , the deviation from the azimuthal polarization pattern due to grain alignment is likely negligible.

Disk model adopted
In this section, we perform a simple analysis of scattering in the disk atmosphere problem. To do so, we adopt a Minimum Mass Solar Nebular model (Weidenschilling 1977) with a column density Σ(R) = Σ 0 R −1.5 au , where Σ 0 = 10 3 g/cm 2 and R au is the cylindrical radius in units of au. We assume a vertically isothermal temperature profile with T (R) = T 0 R −0.5 au , with T 0 = 300 K. This results in a mildly flared disk with (H g /R) = 0.045 R 1/4 au , where H g is the gas scale height. The dust scale height is different from H g and depends on the grain size, as (Youdin & Lithwick 2007): where St = ρ s a/Σ is the Stokes number that determines how well the dust grains are coupled with the gas, and α is the turbulence parameter, which we take to be α = 10 −4 . We assume dust grains have a powerlaw distribution (Mathis et al. 1977): N (a) ∝ a −3.5 after vertical integration, between a min = 0.01µm and a max = 1 mm. The total column density of the dust grains is 0.01Σ, where a gas-to-dust ratio of 100 is assumed. In practice, we use 100 bins of dust grains distributed evenly in logarithmic space, each represented by the center of the bin. For each bin, the dust grains follows vertical Gaussian distributions according to a scale height from Equation (1).

Grain size at τ = 1 surface
The τ = 1 surface is the surface where radial optical depth τ reaches 1. This is where the stellar light is scattered by the dust grains, and the properties of dust grains at this surface is very important.
In order to calculate the τ = 1 surface, we calculate the extinction cross section at 1.5 µm for grains of dif- Figure 6. The maximum angle difference for different configurations of magnetic fields, assuming an inclination angle of 60 • , H/R of 0.2, and a dust aspect ratio of s = 1.5. The color map represent the maximum angle difference ∆η. Also plotted are contours of constant angle difference ∆η = 70 • for different inclination angles. The inclination angle i is labeled next to each line. ferent sizes using Mie theory (Bohren & Huffman 1983) through the miepython module 3 . We then integrate the optical depth at λ = 1.5 µm radially outward from the center of the disk. The τ = 1 contour is shown in the upper panel of Figure 10. Also shown in the upper panel as a dashed line is the H/R of the surface. We can see that the τ = 1 surface is largely flat, despite the fact that our disk model is mildly flared. There are two reasons for this flat surface. On the one hand, the dust settles towards the midplane more at large radii due to the less turbulent stirring from the more diffuse gas. This effect can also be viewed in the lower panel of Figure 10, where the z τ =1 /H g is plotted against R as a solid line. We can see that the τ = 1 surface in terms of the gas scale height H g gradually decreases. On the other hand, and probably more importantly, the outer regions are blocked or "shadowed" by the inner regions, so that the τ = 1 surface can never bend towards midplane.
Note that Avenhaus et al. (2018) finds that the protoplanetary disks are moderately flared, as opposed to being flat in our model. In our turbulent stirring model, if we fix the grain size, the Stokes number goes as St ∼ R γ , where −γ = −1.5 is the power-law index for the adopted column density profile. In the limit of α Available at https://github.com/scottprahl/miepython H g /R ∼ R (1−q)/2 , with −q = −0.5 being the power-law index for the temperature profile. So we have, for fixing grain size a, H d /R ∼ R (1/2)(1−q−γ) . Since we have adopted q = 0.5, γ = 1.5, the H d /R decreases with R, and the dust at larger radii is shadowed by the dust at inner radii. If we require the near-IR scattering surface to be flared as well, we need q + γ < 1. This is very hard to achieve. The requirement may be alleviated with the introduction of radial variation of turbulent α. If α goes as R ω , i.e. the disk is more turbulent at a larger radius, the above constraint becomes q+γ < 1+ω. Another potential and more likely way to make a flaring scattering surface is to abandon the turbulent stirring model and introduce a disk wind, which may entrain small grains and make the disk appear flaring. We will not discuss these alternatives in more detail in what follows. The turbulent stirring model we introduce here is only for illustrative proposes and to lay the foundations for the following discussions on synthetic maps and grain alignments, both of which are not sensitive to whether the disk is flared or not.
In the lower panel of Figure 10, we also plot the grain size with St = α as a function of R. It characterizes the maximum size of the grains that can be stirred up to a height comparable to the gas scale height (c.f. Equation 1).
To view the dust distribution at τ = 1 surface more clearly, we plot the mass, number density, and extinction cross section at R = 30 au for each grain size bin in Figure 11. The mass ∆M is the mass density of dust grains between a and a + ∆a. The number density is ∆N/∆a = ∆M/(m a ∆a), with m a = ρ s (4π/3)a 3 being the mass of a grain with radius a. The extinction cross section is defined as ∆σ ext = ∆M κ ext (a), with κ ext (a) being the extinction opacity for the grains of radius a.
At R = 30 au, we have z τ =1 = 6.4 au, and a(St = α) = 2.3µm. When the grain size increases beyond 1µm, the mass and the number density drop very quickly. The extinction is also dominated by grains with radius of ∼ 0.6µm, or size parameter of 2πa/λ ∼ 2.5. This justifies our discussion on grain size in Section A.6 and the adoption of dipole approximation in most of this work.

Synthetic maps
With the τ = 1 surface and the corresponding H/R obtained, we can calculate the Stokes parameters at each location and generate synthetic maps. In order to focus on the deviation from the azimuthal pattern, we choose to show maps of the ratios of azimuthal Stokes parameters, U φ /I and Q φ /I. The results for Model 1 and Model 2 are shown in Figure 12. In order to make small values more visible while using the same colormap for all pan-  also smaller compared to the more inclined Model 1 and stays below ∼ 0.05 over most of the disk.
Last but not the least, the sign of U φ depends on the direction of the toroidal magnetic field. If we change the φ B to −φ B , the Q φ /I maps are unaffected but the U φ /I maps will change signs everywhere while keeping the magnitude the same. This behaviour, alongside information on disk rotation directions, can be used to distinguish it from other mechanisms that produce U φ , such as through multiple scattering (Canovas et al. 2015). We will discuss how to distinguish our U φ -producing mechanism from others in more detail in Section 5.2.

Grain alignment at τ = 1 surface
In this section, we briefly discuss the grain alignment at the τ = 1 surface based on timescales of several most relavent processes. The discussion is similar to Tazaki et al. (2017) and Yang (2021).
1. The gaseous damping timescale: it determines how fast random collisions with gas particles disalign dust grains and is given by: where ρ s is the solid density of dust grains, n g is the number density of gas molecules (assuming mean molecular weight of 2.3). For our adopted disk model at τ = 1 surface and 1 µm dust grains, the gaseous damping timescale is plotted as a blue curve in Figure 13.
2. The Larmor precession timescale: it is the timescale that rotating dust grains with magnetic moment due to the Barnett effect (Barnett 1915) precess around an external magnetic field: where T d is the dust temperature which we take to be the same as the midplane gas temperature prescribed above (vertically isothermal). For the magnetic field, following (Yang 2021), we adopt the estimate from (Bai 2011): where we have assumed the mass accretion rateṀ = 10 −8 M yr −1 , typical for classical T Tauri stars (Hartmann et al. 2016).
The Larmor precession timescale is plotted as an orange curve in Figure 13. It is possible for dust grains to carry superparamagnetic inclusions (Jones & Spitzer 1967), which can shorten the Larmor precession timescale by up to a factor ofχ ∼ 10 3 (Yang 2021). We also plot the Larmor precession timescale with superparamagnetic inclusions withχ = 10 3 as a green curve.
3. RAT precession timescale: it is the timescale for the precession due to the Radiative Alignment Torque Figure 12. Synthetic maps. The left two panels are U φ /I and Q φ /I for Model 1 with i = 60 • . The right two panels are the same but for Model 2 with i = 45 • . Note that the colormap is the same for all panels and uses symmetric logarithmic normalization while being linear between −0.1 and 0.1.  (Lazarian & Hoang 2007). It can be estimated as (Tazaki et al. 2017): whereλ is the energy weighted averaged wavelength of the radiation, γ is the anisotropy of the radiation, and u ISRF = 8.64 × 10 −13 is the interstellar radiation energy density (Mathis et al. 1983). We consider only the stellar light, hence γ = 1. For radiation energy density u rad , we assume solar parameters, with effective temperature of ∼ 6000 K and bolometric luminosity L = L . This yields: andλ = 0.89µm. Sinceλ ≤ 1.8a, we have |Q Γ | ≈ 0.4 (Lazarian & Hoang 2007). The RAT precession timescale is plotted as a red curve in Figure 13. From Figure 13, we can see that t RAT,p is always smaller than t d , indicating efficient Radiative Alignment torque. If the dust grains are of regular paramagnetic materials (orange curve), we also have t L > t RAT,p . In this case, we have k-RAT, i.e. grains aligned with radiation flux. As a result, the dust grains are aligned with short axes along the stellar light direction. As a result, the dust grains would look round from the star and they scatter stellar light exactly the same as spherical dust grains in the dipole regime with small particles. Hence we expect no deviation from the azimuthal pattern if dust grains are aligned with k-RAT.
If the dust grains possess superparamagnetic inclusions (SPIs), the Larmor precession timescales can be reduced by a factor up to about 10 3 (Yang 2021). In this case (the green curve), we have t L > t RAT,p outside a radius of 0.7 au. So it is likely that superparamagnetic dust grains are aligned with the magnetic field rather than the radiation flux for the majority of the disk at tens of au scale. As SPI-candidates are seen in meteorites (Goodman & Whittet 1995), it is possible that dust grains in the disk atmosphere possess SPIs as well. The near-IR wavelength scattering polarimetry of protoplanetary disks can be an excellent probe for the existence of SPIs, which will help understanding magnetic alignment of large dust grains in other environments.
Lastly, it is worth mentioning that the internal relaxation can also be problematic in disk atmosphere. Tazaki et al. (2017) estimated the internal relaxation timescale as a function of grain radius in their Figure 2 and showed that the internal relaxation timescale of 1 µm dust grains can be tens of years, longer than any of the timescales considered in Figure 13. In this case, the degree of alignment may be reduced due to the lack of internal relaxation (Hoang & Lazarian 2009;Tazaki et al. 2017). More detailed and quantitative discussion and modeling on grain alignment is beyond the scope of this paper and will be deferred to future investigations.

Detectability
In previous sections, we have shown that the angle deviation from azimuthal polarization can easily reach 10 • or more for a moderately inclined disk. Here we discuss the detectability of such angle deviations.
The angle with azimuthal direction can be calculated through η = (1/2)arctan2(U φ , Q φ ). Without loss of generality, we limit our discussion here to the quadrant where Q φ > 0, U φ > 0, so that η = (1/2)arctan(U φ /Q φ ). The total differential is then: We can see that the above expression is on the order of 1/(2SNR), where SNR is the signal to noise ratio of polarized intensity defined as δPI/PI, with PI being the polarized intensity. This estimate can be made more accurate if we focus on the deviation from the azimuthal pattern where Q φ = PI, U φ = 0, so that δη = δU φ /2Q φ = 1/(2SNR). For the 10 • angle deviation we obtained before, we need SNR≤ 7.8. In the survey presented by Avenhaus et al. (2018), the SNR is better than 20 for most cases, which translates into an error in angle of δη ∼ 1.4 • . So if the dust is aligned with the magnetic field in the atmosphere of a moderately inclined disk, we should be able to detect the deviation from the azimuthal polarization pattern as predicted by this work.

Distinguishing different mechanisms
The most important feature of our polarization mechanism is that we relies on elongated dust grains that are aligned with magnetic fields. If a mechanism that produces near-IR scattering polarization accounts for only spherical dust grains or non-spherical dust grains but without grain alignment, the optical properties of the ensemble of dust grains will have spherical symmetry, i.e. the scattering matrix is solely a function of scattering angle. Under this assumption, for light last scattered at the near-side (the right point in Figure 5) or at the far-side (the left point in Figure 5), the geometry of the scattering problem is symmetric between up and down 4 . As a result, the end polarization orientation would either be along the radial direction or along the azimuthal direction with no other possible outcomes, which means U φ = 0. In contrast, the U φ is maximized near the near-side and the far-side points in our models (c.f. Figure 12).
The main alternative discussed in the literature so far is the multiple scattering at high optical depth and high inclination (Canovas et al. 2015). Here we give a heuristic argument on how this mechanism works. If we take single scattering of spherical dust grains as the zeroth order problem, the first order problem will be the photons scattered twice before they reach our telescope. As discussed above, the zeroth order problem considering only single scattering cannot produce U φ . The U φ is then produced primarily by the first order problem with double scattering. Since the disk atmosphere is optically thick at near-IR, the first scattering site cannot be too far from the second scattering site. We shall refer to the particle at first scattering site as particle A and the particle at the second scattering site as particle B. The light coming from particle B is then what we observe. In this first order double scattering problem, the local anisotropy of radiation field as viewed at the location of particle B is what determines the polarization state of the scattered light. If we further ignores the polarization of the light between particle A and B and treat the light scattered by the particle A as non-polarized, the problem reduces to a problem that is the same as the self-scattering problem at high optical depth at (sub)millemeter wavelengths. The particle A is comparable with the original source of the dust thermal emission at (sub)millemeter wavelengths, and the particle B is the scattering particle of the self-scattering problem. As discussed by Yang et al. (2017), the polarization orientation is along the "minor axis" of the local disk surface, the direction that is coplanar with both the final scattering direction and normal direction of the local disk surface (see Figure 1 Figure 14. The U φ /I profiles. The red curve shows the U φ /I for multiple scattering, using a simple doublescattering approximation, assuming an inclination angle of 60 • . See text for discussions. The blue solid curve comes from the fiducial model, Model 1, with i = 60 • . The blue dashed curve is similar to Model 1, but with φB = −90 • . This changes the helicity of the assumed magnetic field configurations, but nothing else, resulting in U φ flips its sign in the whole disk. This is a unique feature for polarization produced by magnetic alignments. suming i = 60 • are shown as the red curves in Figure 14. Note that the absolute values of U φ /I from the above simple double scattering model is arbitrary in the sense that the contribution from single scattering is not taken into account in this simple double scattering model. The single scattering does not produce U φ , so it should not affect the overall profile, if azimuthal variation in single scattering is ignored. Higher order scattering events may have more significant contributions which are not taken into account here. We choose to multiply the whole curve by 0.2 to make the results comparable with those from the magnetic alignment. Despite the simplicity of this model, it still captures most of the physics and the result agrees with Canovas et al. (2015)'s moderate opacity model very well. The most important features of U φ generated by multiple scattering is that the U φ is opposite between symmetric points with respect to the disk minor axis (0 • − 180 • vs. 180 • − 360 • ; c.f. Figure  2 of Canovas et al. 2015).
For comparison, we also plot the the U φ for Model 1 with i = 60 • as blue curves. We can see that the U φ is maximized at the near and far sides, with near side having positive 0.15I and the far side having −0.15I. If we change the helicity of the magnetic field by changing φ B to −φ B , the U φ at the near and far sides will also flip to their opposite values. The curve for −φ B is plotted as blue dashed curve in Figure 14. We can see that the solid and dashed blue curves are symmetric with re-spect to the middle point. The U φ maps presented in Figure 12 can be changed to the magnetic field configuration with opposite helicity by multiplying the whole map by −1 and then flip upside down. This dependence on the helicity of magnetic fields can be very important in distinguishing our mechanisms from others. If we can infer the helicity of the magnetic fields through the rotation curves in the outflows or jets, we can then check against near-IR scattered polarimetry and see if there is any magnetic field signatures and see if the predicted U φ map for given helicity of magnetic fields agrees with observations.

Potential sources with aligned grains
Before going into specific systems, we summarize our discussions on the difference between our mechanism that relies on elongated dust grains aligned with magnetic fields and the other mechanisms that rely only on scattering by spherical dust grains. For axis-symmetric system, alternative mechanisms with only spherical dust grains will have U φ map being symmetric with respect to the minor axis of the disk. The symmetry is to be considered in the sense of opposite signs. That is to say, if there is a structure of positive U φ in one side of the disk, there has to exist the exact same structure of negative U φ in the other side of the disk, with these two structures being mirror symmetric with respect to the minor axis. For our mechanism that relies on elongated dust grains, there is no such requirement, and the U φ map can even be dominated by either positive or negative U φ (c.f. Figure 12). The asymmetry in the U φ map with respect to minor axis and/or predominant positive/negative U φ in axis-symmetric systems are both signs of our mechanism.
With this in mind, we find that CU Cha, HD 169142, MWC 614, Hen 3-365, and HD 142527 are some good candidates in the Gemini-LIGHTS survey (Rich et al. 2022). They all show clear deviations from mirror symmetry expected for scattering by only spherical dust grains. In addition, the outer disk of HD 142527 are predominantly positive. Similarly HD 169142 is also predominantly positive. HD 34700 A has U φ maximized along the minor axis, which cannot be explained by spherical dust grains. MWC 614 has clear asymmetry in U φ and slightly more positive than negative U φ . We include Hen 3-365 (HD 87643) as a good candidate, despite its non-axisymmetric structures which complicate the interpretation. In addition to the substantial asymmetry in the U φ image, Hen 3-365 has a wedge of negative Q φ (see also Laws et al. 2020), similar to the Model 1 (c.f. Figure 12). In the DARTTS-S survey (Avenhaus et al. 2018), we find V4046 Sgr and DoAr 44 as candidates based on their asymmetric U φ images.
We would like to note that the calibration of near-IR polarimetry data is a very complicated process. Part of the calibration involves correction for instrumental polarization, stellar polarization, and/or foreground interstellar contamination. Since there is no priori knowledge on what the stellar polarization and interstellar polarization are, a parameterized approach is usually adopted (Avenhaus et al. 2018). How these effects affect the detection of the signals predicted in this work remains to be determined.

SUMMARY
In this paper, we have studied the scattering of the near-IR stellar light by small dust grains, that are aligned with respect to the magnetic fields, in the atmosphere of a protoplanetary disk. We focused on the polarization orientation of the scattered light and showed that the deviation from the often assumed azimuthal pattern can be significant. The main findings are as follows.
1. We calculated the polarization pattern in a disk frame (DF). We focused on two models: Model 1 with a relatively large inclination angle (i = 60 • ) and a large toroidal magnetic field component and Model 2 with moderate parameters (i = 45 • ). The Model 1 has polarization reversal, i.e. the scattered light is polarized in the radial rather than azimuthal direction, at certain locations. The Model 2 doesn't have polarization reversal, but still has a maximum angle difference of 11 • , detectable if we have SNR > 7.8 for Stokes parameters.
2. We gave a geometric explanation of the polarization reversal in Section 3.2. We showed that the polarization reversal is almost inevitable for disks with large inclination angles i > 75 • , regardless of the magnetic field configuration.
3. The angle difference strongly depends on the field configuration. In particular, it increases with an increasing toroidal component of the magnetic field. Hence it can be used to probe the launching mechanism of magnetized disk wind.
4. With a simple minimum mass solar nebular model, we studied the τ = 1 surface for scattering near-IR stellar light and found that the maximum grain size there is on the order of 1 µm, assuming a turbulent viscosity of α = 10 −4 . This justifies the focus of this initial study on relatively small grains. 5. We calculated synthetic U φ /I and Q φ /I maps for the two disk models on the τ = 1 surface. The peak U φ /I is on the order of ∼ 0.15 and ∼ 0.05 for Model 1 and 2, respectively. Interestingly, the U φ /I is reversed with magnitude unaffected if we change the azimuthal direction of the magnetic field (through φ B → −φ B ). This effect, together with disk rotation information, can be used to distinguish our mechanism from other U φ -producing mechanisms.
6. We conducted a grain alignment analysis at the τ = 1 surface. We found that Radiative Alignment Torque should be operating. For regular paramagnetic dust grains, our model favors k-RAT, i.e. grains aligned with radial stellar light. If grains possess substantial superparamagnetic inclusions, B-RAT becomes likely.
7. We compared the azimuthal profiles of U φ between our model and an alternative model that relies on multiple scattering of spherical dust grains. We argue that a spatial distribution of U φ that is predominantly positive or negative and/or asymmetric with respective to the minor axis of an intrisically axisymmetric disk are are signals of aligned elongated dust grains. We identified a handful of systems in the existing literature that are potential targets to look for magnetically aligned grains in future studies. In this appendix, we expand on the discussion of the scattering-induced polarization in the grain's frame presented in § 2.
The geometry of the set-up is shown in Figure A.1. The scattering particle is placed at the center of the frame, with z direction along the symmetry axis of the grain. The blue arrow denotes the incoming lightn i , which is placed in the xz plane, without loss of generality. It makes an angle i with the z axis. The red arrow denotes the scattered lightn s , defined by two position angles θ and φ. For the incoming light, its polarization is defined withê 1 andê 2 . A positive Q corresponds to polarization alongê 1 and a positive U corresponds to polarization along a direction bisectingê 1 and e 2 . For the scattered light, its polarization is defined withθ andφ. Here, a positive Q corresponds to polarization alongθ and a positive U corresponds to polarization along a direction bisectingθ andφ. The polarization orientation angle η is defined in theθ-φ plane as the angle starting fromθ and increasing towardsφ counterclockwise, with values between 0 • to 180 • .
In Sections A.1-A.5, we will limit our discussions to the dipole approximation (for small grains) to gain a better understanding. The calculation is done with the electrostatic approximation or dipole approximation (Bohren & Huffman 1983). See Appendix B for a brief description of this method. The impact of large dust grains is briefly discussed in Section A.6.

A.1. Scattering by Small Spherical Particles
Before calculating the polarization from scattering off elongated dust grains, let's first look at the simpler case for spherical particles when i = 45 • . In the left two panels of Figure A.2, we show the polarization fraction p (first left; defined as Q 2 + U 2 /I) and the polarization orientation angle η (second left) as we change the direction of the scattered light. Note that since the particle is perfectly spherical, the polarization is completely determined by the angle betweenn i andn s and the results contain no new information but the well known polarization profile p = sin 2 <n i ,n s > /(1 + cos 2 <n i ,n s >) and the fact that the polarization direction is perpendicular to bothn i and n s . Nonetheless, the spherical case in the left two panels of Figure A.2 will serve as a benchmark to help us understand the cases for aspherical dust grains later.
In the first left panel of Figure A.2, we overlay a curve which has p = 1, i.e. the scattered light is fully polarized. In this small spherical dust grain case, this is the direction that is perpendicular to the incoming light direction.
In the second left panel of Figure A.2, we can clearly see that there exist two singular points, one with (θ = 45 • , φ = 0 • ), and one with (θ = 135 • , φ = 180 • ). These two points correspond to the forward scattering and backward scattering, correspondingly. In these two directions, the polarization fraction in scattered light is 0, and the polarization direction is ill defined. If we walk around the singular point while fixing θ, the polarization will be alongθ direction, in order to be perpendicular to both incoming and scattered light. Similarly, if we walk around the singular point while fixing φ, the polarization will be alongφ direction. Note that η = 0 • and η = 180 • correspond to the same polarization orientation.
It is worth mentioning that the polarizability matrix (see Appendix B and Equation B7) in this spherical case is isotropic and diagonal:ᾱ = diag{α s , α s , α s }, where with a being the grain size, and is the complex dielectric function.

A.2. Fiducial case
Now let's move on the the more interesting case with aspherical dust grains. For our fiducial case, we consider a small dust grain in the dipole regime. In this work, we consider only oblate dust grains characterized by an aspect ratio s > 1, which we set to s = 1.5 in our fiducial case. In the grain's frame depicted in Figure A.1, the symmetry axis of the dust grain is placed along the z direction. The light makes an angle i = 45 • with the z direction, the same as the spherical case discussed above. For the fiducial model and most of the models in this paper, we assume the composition from Birnstiel et al. (2018). It is a mixture of 20% water ice (Warren & Brandt 2008), 33% astronomical silicates (Draine 2003), 7% troilite (Henning & Stognienko 1996), and 40% refractory organics (Henning & Stognienko 1996) by mass. Thoughout this paper, we assume an observing wavelength of 1.5 µm. The results are largely independent of the specific choice of the wavelength. The results are shown in the right panels of Figure A.2.
In the grain's frame, the polarizability matrix is always diagonal as P = diag{α 1 , α 1 , α 3 }, and |α 1 | > |α 3 | because we assume oblate dust grains. See Appendix B for more details. For an incoming radiation propagating along then i direction, we can decompose the light to two components: E i = E 1ê1 + E 2ê2 . The dipole excited in response to these two components are: and We can see that P 2 is always along y direction, i.e. theê 2 direction. At the same time P 1 is not along theê 1 direction any more, due to the difference between α 1 and α 3 . This is the very reason why scattering by aligned aspherical grains is different from scattering by spherical grains.
In the third left panel of Figure A.2, we can see that the maximum polarization is still p = 1, i.e. fully polarized. The location where p = 1 is achieved is plotted as a solid curve in the figure. The p = 1 curve for the spherical case is also plotted in the figure as a dashed line. We can see that the p = 1 locations are slightly different between these two cases. The difference is zero at θ = 90 • , φ = 90 • or 270 • directions. This is because these two directions correspond to the ±ŷ direction, which is along the dipole P 2 . As a result, P 2 does not contribute to the scattered light, and the scattered light is fully polarized.
Along the φ = 0 line in the right most panel of Figure A.2, the fully polarization (p = 1) is achieved at θ = 130.3 • . It differs from 135 • in the spherical case, by 4.7 • , due to the difference in direction between P 1 andê 1 . This is also the difference in the polarization orientation η (Right panel of Figure A.3) at θ = 90 • and φ = 90 • or 270 • , which is for the same reason.
In the right panel of Figure A.2, we first notice that the singular points for the spherical grains ((θ = 45 • , φ = 0 • ) and (θ = 135 • , φ = 180 • )) are no longer singular. There are two reasons for this behavior. Firstly, the forward scattering directionn 1 is no longer perpendicular to P 1 . As a result, the emission is no longer maximized for the dipole radiation from P 1 , making it inferior to the radiation from P 2 . Secondly and more importantly, |P 1 | < |P 2 | because α 1 > α 3 . These two reasons combine to make the dipole radiation from P 2 dominate over that from P 1 in the forward and backward scattering directions. At these two singular points in the spherical cases, the polarization is thus along P 2 direction, with η = 90 • . Because η near these points can be 0 • when varying along constant θ, the difference in polarization orientation between the spherical case and aligned aspherical case is always as large as 90 • near these points. There still exist singular points in the diagram for the aspherical case. They are located symmetrically around the previous singular points with the same θ but different values of φ between the spherical and fiducial cases.
In the left two panels of Figure A.3, we show the difference in polarization degree and the difference in η for different scattering directions. We can see that the difference in polarization degree can reach up to 13%. The difference in η strongly depends on the scattered light direction. Near the forward and backward scattering direction, i.e. the singular points in the spherical case, ∆η ≈ 90 • , which applies up to the new singular points and forms ribbon-like structures in the ∆η plot.
In the right panel of Figure A.3, we show the histogram of the angle difference ∆η. We use a vertical dashed line to show ∆η = 4.7 • , which is the angle difference betweenê 1 and P 1 . We can see that while most scattering directions have ∆η 4.7 • , a fraction of them have substantially larger ∆η of 10 • − 20 • .

A.3. Dependence on the dust aspect ratio
To compare the angle difference with different dust models and/or different inclination angles, we propose the following two metrics. The first one is the ∆η at the scattering angle of θ = 90 • and φ = 90 • , which we call ∆η y . This is always the direction along P 2 (andŷ) and hence the scattered light comes purely from P 1 , so that the ∆η y equals the angle difference between P 1 andê 1 .
In the upper panels of Figure A.4, we show the results for an extreme case with s = 100 and its comparison with the spherical case, again assuming i = 45 • . We can clearly see in the upper left panel that the p = 1 curve moves closer to θ = 90 • line. This is due to the fact that |α 1 | |α 3 |, so that P 1 is close tox. In the limit that P 1 x, it can be easily verified that only scattered lights in the xy plane, with θ = 90 • , are fully polarized.
In the lower panel of Figure A.4, we show the histogram for the s = 100 case. The vertical dashed line represents ∆η y = 22.3 • . We can see that ∆η y nicely characterizes the angle difference in this case as well, with most scattered light having an angle difference comparable to or less than ∆η y .
While ∆η y sets a scale for the most probable angle difference, it does not provide a good description for the spread in the histogram beyond ∆η y . To characterize the spread, we propose a second metric max(p∆η), the maximum value of p∆η for all scattering directions. It is motivated by the fact that large ∆η directions tend to have low polarization fractions (see Figure A.2 and Figure A.4). In the case of s = 100, we have max(p∆η) = 29.7 • . Note that we always have max(p∆η) ≥ ∆η y , because when scattered towardsŷ, we have p = 1 and ∆η = ∆η y . The difference between max(p∆η) and ∆η y is a measure of the spreading beyond ∆η y .
In the left panel of Figure A.5, we show the two metrics, ∆η y and max(p∆η), as a function of the aspect ratio s. We can clearly see that more flattened (oblate) grains have larger deviations in polarization orientation due to scattering compared to spherical grains. For s = 2, the angle difference is typically on order of 10 • , while for extremely elongated grains, the difference can be as large as 30 • for a large fraction of scattering angles.

A.4. Dependence on the inclination angle
In the right panel of Figure A.5, we show the two metrics, ∆η y and max(p∆η), as a function of the inclination angle i. We can see that ∆η y increases as we increase the inclination angle i initially, then falls back to 0 as we go towards i = 90 • . In the two limiting cases, i = 0 and i = 90 • , we have ∆η y = 0. This is becauseê 1 becomes aligned with one of the principle axes (ẑ for i = 0 andx for i = 90 • ), hence P 1 ê 1 . However, this doesn't imply the angle difference also goes to 0 as the inclination angle approaches 90 • . In fact, for larger inclination angles, the difference between |P 1 | and |P 2 | increases with i, and the polarization orientation deviates from the spherical case (with |P 1 | = |P 2 |) progressively. These differences are not captured by the first metric ∆η y .
At a larger inclination angle, there is a stronger need for the second metric. To see what exactly happens at a larger inclination angle, we show the polarization fraction and angle difference for i = 75 • in upper panels of Figure A.6. To compare, we show the same for i = 15 • in lower panels of Figure A.6. These two models have similar ∆η y (2.5 • and 2.1 • for i = 75 • and i = 15 • , respectively). We can see that even though the deviations of the p = 1 locations from corresponding spherical models are similar for these two inclination angles, the ribbon-like structures are much larger in i = 75 • than in i = 15 • . This results in a larger spread in the histogram of angle differences in i = 75 • than i = 15 • , shown in the right panels of Figure A.6. This spread is nicely captured by our second metric, max(p∆η), which are 13.2 • and 2.1 • at i = 75 • and i = 15 • , respectively. According to max(p∆η), there is a substantial fraction of scattering directions with angle differences as large as ∼ 13 • when i = 75 • , even though the ∆η y is only 2.5 • .
With a better understanding of the two metrics, ∆η y and max(p∆η), we now come back to the right panel of Figure A.5. We conclude that the angle difference with the spherical model increases monotonically as the inclination angle i increases. There is a large fraction of scattering directions with ∆η of ∼ 15 • or larger when i is close to 90 • .
Interestingly, max(p∆η) = ∆η y for small inclination angles. In these cases, the ribbon-like structures are very thin and the angle difference is dominated by the difference between P 1 andê 1 (see upper lower panels of Figure A.6).

A.5. Dependence on dust composition
The composition of dust has a strong impact on the optical properties of dust grains. To show the dependence of angle difference on dust composition, we calculated ∆η y and max(p∆η) assuming i = 45 • and s = 1.5 for five illustrative dust compositions. The results are tabulated in Table A.5. We also listed the real (n) and the imaginary part (k) of the refractive index for each composition at the wavelength λ = 1.5µm.
Among the materials considered, troilite produces the largest angle difference. This may be related to its absorptive nature, characterized by its large imaginary part of the refractive index k.

A.6. Results for moderately large dust grains
The size of grains has a strong impact on dust scattering. To relax the previous small grain size assumption, we use the PyTMatrix 5 module (Leinonen 2014), which is a wrapper for the TMatrix code (Mishchenko & Travis 1994). The results for MRN-distributed dust grains with s = 1.5 and a max = 1µm, corresponding to a size parameter of x max = 4.2, are shown in Figure A.7. We can see that the polarization fraction decreases significantly from 100% at the peak curve, but the distribution of η is still similar to the one in dipole regime in right panels of Figure A.2. The histogram of the angle difference is also larger for the 1µm grains. The dipole approximation is reasonable for 1µm grains or smaller, if we focus on the polarization orientation η.
We note that even larger grains (a max ≥ 2µm) can produce more complicated polarization patterns and distributions of η (results not presented in this paper), which are harder to use for interpreting observational results. Larger grains may account for the strong forward scattering observed in some systems (e.g. Avenhaus et al. 2018). We will postpone a full exploration of larger dust grains to a future investigation and focus on the dipole approximation in this work.