employs the linear model for the real part; thus it is pressure independent)
for: uncoated bubbles in water (a and b), coated bubbles in water (c $\&$ d)
and uncoated bubbles in tissue ($\rho$=1060 kg/m3, $C_{l}$=1540 m/s,
$\mu_{s}$= 0.00287 Pa.s, $G$=0.5 MPa, $\sigma$=0.056 N/m 39 ) (e and f).

###  D.2 Validation of the model at higher pressures against the semi-linear
Louisnard model

(a) (b)

(c) (d)

(e) (f)

Figure 8:  Comparison between the predictions the Louisnard $\&$ the nonlinear
model for sound speed and attenuation. Case of a bubbly medium with uncoated
MBs with $R_{0}$= 2 $\mu$ m and $\beta_{0}$=10-5. a) attenuation at $P_{a}$=40
kPa, b) sound speed at $P_{a}$=40 kPa, c) attenuation at $P_{a}$=100 kPa d)
sound speed at $P_{a}$=100 kPa, e) attenuation at $P_{a}$= 150 kPa and f)
sound speed at $P_{a}$=150 kPa.

As the pressure increases, assumptions (e.g. small amplitude MB oscillations)
on which the linear model is based on are no longer valid. To investigate the
effect of pressure, the radial oscillations of the MBs were simulated for
exposures of various acoustic pressure amplitudes. For the uncoated bubble
$P_{a}$= 40, 70, 100, 150 kPa, for the coated bubble $P_{a}$= 40, 70, 100,
150, 200 kPa and for the bubble in tissue $P_{a}$= 100, 500, 700, 1000 kPa
were chosen. The power dissipation expressions for nonlinear damping effects
which are given by Eq.50 (for uncoated bubble), Eq.51 for coated bubble and
Eq.52 for the bubble in tissue were used to calculate the total dissipated
power. The imaginary and real part of the wave number were then calculated
using Eqs.11 and 10 in case of the nonlinear model and Eq.53 and Eq.54 in case
of the Louisnard model 18 . The predictions of the two models are illustrated
in Fig. 7.
The left column of Fig. 7 shows that the
$\left\langle{}{\Im}\left(k^{2}\right)\right\rangle{}$ calculated by Eq. 11 is
in excellent agreement with the Louisnard model (Eq. 53) for all the acoustic
pressures and the bubble models that investigated. The simple approach
introduced here, only needs the radial oscillations of the bubble as input and
reduces the complexity of the Louisnard model where the equations for
different dissipation mechanisms must be derived for each bubble case.
Variations of $\left\langle{}{\Im}\left(k^{2}\right)\right\rangle{}$ with
pressure shows the importance of the considerations of the pressure effects as
the linear model fails to predict phenomena like the resonance shift (e.g. 14
), changes in the amplitude of the
$\left\langle{}{\Im}\left(k^{2}\right)\right\rangle{}$ with pressure and the
generation of SuH (e.g. sojahrood3 ) and subharmonic (SH) resonances (e.g. 54
; 55 ).
As an instance in case of the uncoated bubble in Fig. 7a
$\left\langle{}{\Im}\left(k^{2}\right)\right\rangle{}\approxeq$ 8.5$\times$108
m-2 at pressure dependent resonance( $\frac{f}{f_{r}}\approxeq 0.98$) when
$P_{a}$=40 kPa. However, as pressure increases to $P_{a}$=150 kPa resonance
shifts to $f/f_{r}\approxeq$ 0.64 and
$\left\langle{}{\Im}\left(k^{2}\right)\right\rangle{}\approxeq$ 7.3$\times$108
m-2. Moreover, a SuH occurs at ${f}/{f_{r}}\approxeq$ 0.34 with
$\left\langle{}{\Im}\left(k^{2}\right)\right\rangle{}\approxeq$ 2.4$\times$108
m-2. When $P_{a}=40kPa$ and at $\frac{f}{f_{r}}\approxeq$ 0.34,
$\left\langle{}{\Im}\left(k^{2}\right)\right\rangle{}\approxeq$ 3.9$\times$108
m-2. Thus, the pressure increase has a significant influence on the resonances
of the system and the magnitude of the
$\left\langle{}{\Im}\left(k^{2}\right)\right\rangle{}$.
The Louisnard model uses the linear assumptions (Eq.54) to calculate the
$\left\langle{}{\Re}\left(k^{2}\right)\right\rangle{}$. The predictions of the
nonlinear model Eq.10 for
$\left\langle{}{\Re}\left(k^{2}\right)\right\rangle{}-({\omega}/{C_{l}})^{2}$,
are compared with the predictions of Eq.54 in the right hand column of Fig. 7
and for 3 bubble cases (uncoated, coated and bubble in tissue). We have
subtracted the constant $({\omega}/{C_{l}})^{2}$ from
$\left\langle{}{\Re}\left(k^{2}\right)\right\rangle{}$ to better highlight the
pressure dependent changes. In each case, pressure increase leads to
significant changes in $\left\langle{}{\Re}\left(k^{2}\right)\right\rangle{}$,
and predictions of Eq. 3 significantly deviate from the linear values (Eq.
54). As an instance for the uncoated bubble (Fig. 2b) the linear model
predicts a maximum for
$\left\langle{}{\Re}\left(k^{2}\right)\right\rangle{}-({\omega}/{C_{l}})^{2}\approxeq$
4.1$\times$107 m-2 at ${f/f_{r}}\approxeq$ 0.9 and a minimum for
$\left\langle{}{\Re}\left(k^{2}\right)\right\rangle{}-({\omega}/{C_{l}})^{2}\approxeq$
-5$\times$107 m-2 at $\frac{f}{f_{r}}\approxeq 1.12$. However, when
$P_{a}=100kPa$ the maximum of
$\left\langle{}{\Re}\left(k^{2}\right)\right\rangle{}-({\omega}/{C_{l}})^{2}\approxeq$
1.9$\times$107 m-2 at ${f}/{f_{r}}\approxeq$ 0.761 and the minimum is
$\left\langle{}{\Re}\left(k^{2}\right)\right\rangle{}-({\omega}/{C_{l}})^{2}\approxeq$
-8.5$\times$107 m-2 at ${f}/{f_{r}}\approxeq$ 0.773.
The nonlinear model incorporates the pressure-dependent changes in
$\left\langle{}{\Re}\left(k^{2}\right)\right\rangle{}$ and thus can be used to
predict the changes of the
$\left\langle{}{\Re}\left(k^{2}\right)\right\rangle{}$ with pressure. To our
best knowledge this is the first time that the frequency-pressure dependence
of the $\left\langle{}{\Re}\left(k^{2}\right)\right\rangle{}$ in a bubbly
medium has been calculated. The ability of the nonlinear model to calculate
both the $\left\langle{}{\Im}\left(k^{2}\right)\right\rangle{}$ and
$\left\langle{}{\Re}\left(k^{2}\right)\right\rangle{}$ with pressure changes
increase the accuracy of the predictions of the medium attenuation and sound
speed changes.

(a) (b)

(c) (d)

(e) (f)

Figure 9: Influence of bubble-bubble interaction on the pressure dependent
sound speed and attenuation at 100kPa for a coated bubble with $R_{0}=2\mu m$
in Eq. A2: a-b) $\beta_{0}$=10-7, c-d) $\beta_{0}=5.1\times 10^{-6}$ and e-f)
$\beta_{0}=10^{-4}$. In each case 20 bubbles are considered and randomly
distributed in a cube. The side lengths of the cube were chosen to replicate
the $\beta_{0}$ in each case. The side length can be calculated as
$d=(20\times 4\pi R_{0}^{3}/3{}{\beta_{0}})^{1/3}$. The minimum distance
between neighboring MBs was chosen to be 10$\mu$m.

###  D.3 Importance of the accurate calculation of
$\left\langle{}{\Re}\left(k^{2}\right)\right\rangle{}$ in pressure dependent
attenuation and sound speed estimation

Fig. 8 compares the attenuation and sound speed that are calculated using the
nonlinear model and the Louisnard model. The values are calculated for the
uncoated bubble in Figs. 8a-b and at $P_{a}$=40 kPa, $P_{a}$= 100 kPa and
$P_{a}$= 150 kPa. At 40 kPa (Figs. 8a-b), the Louisnard model fails to capture
the sound speed fluctuation around ${f}/{f_{r}}\approxeq$ 0.5 due to the
occurrence of 2nd order superharmonic (SuH) regime. Moreover, the Louisnard
model over-estimates the attenuation at the resonance frequency by about 10
$\%$. The deviation in the predicted values between the two models increases
with increasing pressure. At $P_{a}$= 100 kPa (Figs. 8c-d), Louisnard model
overestimates the attenuation by about 40 $\%$. Moreover, Louisnard model can
not capture the the shift in the maximum sound speed to lower frequencies as
well as the $\approx$ 15$\%$ increase in its magnitude. At 150 kPa (Figs.
8e-f) the Louisnard model overestimates the attenuation peak by 77 $\%$ and
underestimates the sound speed peak by about 52 $\%$. The nonlinear model
predicts a shift in the frequency of the sound speed peak by about 42 $\%$.
Once again, the frequency at which the attenuation peaks (${f}/{f_{r}}$=0.65)
corresponds to the frequency at which ${C}/{C_{l}}$=1.
This, shows that pressure dependent effects of
$\left\langle{}{\Re}\left(k^{2}\right)\right\rangle{}-({\omega}/{C_{l}})^{2}$
can not be neglected and must be included in the calculation of sound speed
and attenuation. The proposed nonlinear model has the advantage of calculating
both of the pressure dependent
$\left\langle{}{\Re}\left(k^{2}\right)\right\rangle{}$ and
$\left\langle{}{\Im}\left(k^{2}\right)\right\rangle{}$.
As the pressure increases, the resonance frequency of the bubbles decreases 14
, which is observed as the peak of $\Im{(k^{2})}$ in Fig. 7 and attenuation
curve in Fig. 8 shift towards lower frequencies; this corresponds to the
frequencies at which the sound speed in the bubbly medium is equal to the
sound speed in the absence of the bubbles. This is seen in Fig. 8 where the
frequency in which attenuation peaks corresponds to the frequency in which
${C}/{C_{l}}$=1 in the blue curves that can only be captured by the nonlinear
proposed model. At pressure dependent resonances, the oscillations are in
phase with the driving acoustic pressure similar to the case of linear
resonance(when $f=f_{r}$ and at $P_{a}\approxeq<$ 1 kPa ${C}/{C_{l}}$=1 page
290 56 ). As the pressure increases, the maximum sound speed of the bubbly
medium increases and occurs at a lower frequency, which depends on the driving
acoustic pressure amplitude. The abrupt increases in the sound speed and
attenuation at particular frequencies in Figs. 8c-d and in
$\left\langle{}{\Im}\left(k^{2}\right)\right\rangle{}$ in Figs. 7a, 7c and 7d
are due to the pressure dependent resonance frequency which is described in
detail in 14 . We have previously shown that when MBs are sonicated with their
pressure dependent resonance frequency, the radial oscillation amplitude of
the MBs undergo a saddle node bifurcation (rapid increase in amplitude) as
soon as the pressure increases above a threshold 14  and the maximum stable
scattered pressure increases considerably.

(a) (b)

Figure 10: Case of a bubbly medium with $\beta_{0}=5.3\times 10^{-3}$ and
$R_{0}$=2.07 mm silberman  (a- attenuation and b-sound speed curves). Blue
circles are constructed by solving the nonlinear model (NM) without bubble-
bubble interaction. Blue solid line is constructed by the linear Commander and
Prosperetti model 17 . The red line-circle is constructed by solving the
nonlinear model and incorporating bubble-bubble interaction. Green diamonds
are experimentally measured values by Silberman silberman . For the
simulations, similar to trujio  pressure amplitude of 10Pa is used.

##  Appendix E Attenuation and sound speed changes at higher void fractions

Some applications of MBs (e.g. pre-clinical drug delivery applications) employ
high concentration of MBs. At higher void fractions, bubble-bubble
interactions become significant. The resonance frequency and maximum
oscillation amplitude has been shown to decrease with MB-MB interaction guedra
; 55 ; Qin1 ; Qin2 . In order to shed light on the pressure dependent changes
of the attenuation and sound speed at higher void fractions we have performed
numerical simulations on a monodisperse population of MBs with and without
considering MB-MB interactions. Figure 9 shows the results of the numerical
simulations for a case of MBs with $R_{0}$=2$\mu$m at 100kPa of driving
acoustic pressure at void fractions of 10-7 (very dilute suspension),
5.1$\times 10^{-6}$ (similar to the void fractions in experiments in Figs.
4-5) and 10-4 (highly concentrated suspension). Radial oscillations in the
absence of interaction at each frequency were computed using Eq. A2. In the
presence of MB-MB interactions radial oscillations were calculated by
considering the pressure radiated by each MB in the location of other MBs by
adding the term $-\rho\sum_{j\neq
i}^{20}{R_{j}}/{d_{ij}}(R_{j}\ddot{R_{j}}+2\dot{R_{j}}^{2})$ to the right side
of Eq. A2. An approach similar to Eq. 14 and 55  is used for the solution of
the large number of the coupled ordinary differential equations. In case of
interacting bubbles, for simplicity we considered a case of randomly
distributed 20 MBs in a cube with a side length of $d$. The side length can
simply be calculated from the void fraction ($d=(20\times 4\pi
R_{0}^{3}/3{}{\beta_{0}})^{1/3}$)). For a sinusoidal acoustic excitation, the
attenuation and sound speed can be calculated using Eqs. 55 and 56 as follows:

$\langle\Re(k^{2})\rangle=\frac{{\omega{}}^{2}}{C_{l}^{2}}+\frac{2{\rho_{l}{}}}{T{P_{a}}}\frac{N}{20}\sum_{i=1}^{20}\int_{0}^{T}{\sin(2\pi
ft)}\frac{{\partial{}}^{2}{\beta{}}_{i}}{\partial{}t^{2}}dt$ (55)
$\langle\Im(k^{2})\rangle=-\frac{2{\rho_{l}{}}}{T{P_{a}}}\frac{N}{20}\sum_{i=1}^{20}\int_{0}^{T}{\cos(2\pi
ft)}\frac{{\partial{}}^{2}{\beta{}}_{i}}{\partial{}t^{2}}dt$ (56)

where $N$ can be calculated from $N$=$\beta$/(4/3$\pi$$R_{0}^{3}$).
At the lower void fraction of $\beta_{0}$=1$\times$10-7 (Fig. 9a-b), MB-MB
interactions are negligible. As concentration increases, the attenuation and
sound speed peak do not linearly scale. At the void fraction of
$\beta_{0}$=5.1$\times$10-6, MB-MB interactions lead to a $\approx$5.6$\%$ and
$\approx$8.5$\%$ decrease in the peak attenuation and the frequency of the
peak respectively (Fig. 9c). Sound speed peak increases by 3$\%$, while the
frequency of the sound speed peak decreased by 9.3$\%$ (Fig. 9d). When in a
solution, each MB receives a sum of the pressures of the multiple neighboring
bubbles, plus the incident acoustic pressure field. Since the scattered
pressure by each MB is not negligible, the interaction effects can not be
neglected. Moreover, these effects are stronger near the MB resonance
frequency. The effective pressure amplitude (pressure of the sound filed plus
the pressure radiated by the other bubbles) at the location of one of the
bubbles in Fig. 9c-d was calculated. Compared to the case in the absence of
interaction, the maximum pressure amplitude felt by the bubble increased by
14$\%$ (114kPa) suggesting a strong MB-MB interaction which can not be
neglected. At the higher void fraction of $\beta_{0}$=1$\times$10-4, MB-MB
interactions lead to a 13.6$\%$ and $\approx$16.2$\%$ decrease in the peak
attenuation and the frequency of the peak respectively (Fig. 9e). Sound speed
peak increases by 10.6$\%$, while the frequency of the sound speed peak
decreases by 16.3$\%$ (Fig. 9f). The results of the numerical simulations
emphasize the influence of the interactions within MBs. Analysis of the
experimental results at large void fractions must include these effects.
Bubble-bubble interactions are even important in the linear regime of the
oscillations. This is emphasized in Commander and Prosperetti 17  as one of
the reasons behind the disagreement between the linear model and the
experiments. For a large enough void fraction, the distance between bubbles
decreases such that the average pressure field exciting a bubble is smaller
than, or comparable with, the pressure wave scattered by a neighboring bubble,
and thus the linear model fails 17 . Since bubble-bubble interactions were
neglected, further analysis by Trujillo trujio ; trujio2 , also did not
achieve good agreement with experiential results near resonance (attenuation
was over-estimated by at least an order of magnitude).
Here, we briefly investigate this important effect, and extend the
experimental verification of our nonlinear model to higher void fractions (in
case of bigger uncoated bubbles). We considered the experimental results by
Silberman silberman  (green diamonds) in Figure 10 for air bubbles with
$R_{0}$=2.07mm and $\beta_{0}=5.3\times 10^{-3}$. Similar to Trujillo trujio
we used an excitation pressure of 10Pa to solve equation A1 at each given
frequency in Fig. 10. Equation A1, was solved with and without bubble
interactions effects. In case of interaction, we considered 20 randomly
distributed bubbles in a cube with a side length $d$. The minimum distance
between bubbles was set to be 10mm. Attenuation and sound speed were then
calculated from Eq. 55 and 56. Figure 10 shows that the linear model and the
nonlinear model in the absence of bubble-bubble interactions overestimate the
near resonance attenuation by about and order of magnitude (blue line in Fig.
10a). However, when bubble-bubble interactions are considered, there is a
significant improvement in the the prediction of the attenuation curve (red
line in Fig. 10a). Similar to the case of the coated MB in Fig. 9d, bubble-
bubble interaction increases the sound speed peak and reduces the frequency of
the peak (Fig. 10b).
Comparing Fig. 9 and Fig. 10 we see that for comparable void fractions
attenuation of MBs are orders of magnitude higher than mm sized bubbles.
Attenuation of a population of MBs with ${}{\beta_{0}}=1\times 10^{-4}$ is
about 800 dB/cm (Fig. 9c) while this is only about 10dB/cm for bubbles with
$R_{0}=2.07$mm at ${}{\beta_{0}}=5.4\times 10^{-3}$ (Fig. 10a). For the sound
speed changes, one can see the exact opposite relationship (Fig. 9d and Fig.
10b).

##  Appendix F Attenuation and sound speed changes at higher pressures

(a) (b)

(c) (d)

(e) (f)

Figure 11: Influence of increasing the pressure amplitude on the sound speed
and attenuation at $\beta_{0}=5.1\times 10^{-6}$ for a coated bubble with
$R_{0}=2\mu m$ in Eq. A2 when $P_{a}$ is: a-b) 200kPa c-d)500kPa and
e-f)1000kPa . In each case 20 bubbles are considered and randomly distributed
in a cube. The dimensions of the cube were chosen to replicate the $\beta_{0}$
in each case. The dimension can be calculated as $d=(20\times 4\pi
R_{0}^{3}/3{\beta_{0}})^{1/3}$. The minimum distance between neighboring MBs
was chosen to be 50$\mu$m to eliminate the possibility of MBs collisions at
higher pressures .

Many applications employ pressures above the Blake threshold of the MBs. To
investigate the frequency dependent attenuation at higher pressures, in this
section we considered a coated MB with $R_{0}$=2$\mu$m with
$\beta_{0}$=5.1$\times$10-6. It is assumed that MB integrity is maintained for
all exposures. Figure 11a-b shows the attenuation and sound speed when $P_{a}$
is 200kPa. The fundamental frequency of the attenuation peak further decreases
(30$\%$ compared to when $P_{a}$=100kPa in Fig. 9c) when bubble-bubble
interaction is considered. The attenuation of the 2nd order superharmonic
(SuH) resonance frequency exceeds that of the main resonance ($\approx$27$\%$)
in Fig. 11a. In the vicinity of the 2nd order SuH frequency (550kHz), the
sound speed peak which was below the medium sound speed at 100kPa
(0.95$C_{water}$) becomes larger ($\approx$1.04$C_{water}$). Bubble-bubble
interactions reduce the attenuation maxima and frequencies of the attenuation
peaks while increasing the maximum sound speed (Figs.11a-f).
The attenuation peak increase with pressure in the studied frequency range
(0.5-3.5MHz) and the influence of bubble-bubble interactions becomes stronger
with increasing pressure. The maximum sound speed increases with increasing
pressure and at 500kPa(625kHz) and 1MPa(575kHz) becomes approximately 2.5 and
4.5 times the $C_{water}$. It is interesting to note that at a $P_{a}$ of
100kPa, and at the same frequencies, the sound speed was below that of
$C_{water}$. The results indicate strong nonlinear changes in the attenuation
and sound speed even at the simulated low void fraction. We will classify
these changes during major non-linear regimes of the oscillations (e.g. pof1 )
in future studies.

##  Appendix G Importance of pressure dependent measurements in shell
characterization of lipid coated MBs undergoing buckling and rupture

(a) (b)

(c) (d)

(e) (f)

(g) (h)

Figure 12: Influence of increasing the pressure amplitude on the sound speed
and attenuation of a lipid coated MB with $R_{0}$=2.7$\mu$m,
$\beta_{0}=5.1\times 10^{-6}$ and different sets of shell parameters A,B,C and
D: a-b)12.5kPa c-d)25kPa and e-f)50kPa, g-h)100kPa. For each group, the shell
parameters are given in table 2. Radial oscillations are calculated using Eq.
12 and bubble-bubble interactions and transducer response are neglected for
simplicity. The duration of the sonicating pulse is 3$\mu$s.

Here we demonstrate one of the applications of the introduced nonlinear model
which is in the accurate characterization of the shell parameters of lipid
coated MBs. The nonlinear behavior of the lipid coating including bucking and
rupture 16  intensifies the nonlinear changes of the resonance frequency.
The linear resonance frequency of the Marmottant modelxia  is given by:

$f_{r}=\frac{1}{2\pi
R_{0}}\sqrt{\frac{1}{\rho}\left(3kP_{0}+(3k-1)\frac{2\sigma_{0}}{R_{0}}+\frac{4\chi}{R_{0}}\right)}$
(57)

This equation has limited sensitivity to the $\sigma_{0}$ value. For a MB with
$R_{0}$=2.7$\mu$m and $\chi$ of 1N/m, eq. 57 predicts only 80kHz changes in
the resonance frequency for $\sigma_{0}$ between 0 and 0.072N/m. More
importantly, the attenuation measurements are often performed using pressures
of $\approx>$20kPa 43 ; 44 ; 45 ; 46 ; faez . Thus, these MBs are already in
their nonlinear regime, and since the resonance frequency of the lipid coated
MBs shifts with increasing pressure (as small as 5kPa) thus, these methods may
underestimate the shell elasticity. Moreover, using the linear model to fit
the pressure dependent shell parameters 17  casts doubts on the accuracy of
the claims such as stiffness softening and shear thinning with increasing
pressure 17 .
An accurate fit requires interrogation of the frequency dependent attenuation
at multiple pressures in increasing steps such as what was done in this study
or in 26 . At the lower pressures (e.g. $\approx$$<$50kPa) and in the absence
of the shell rupture, the changes in the resonance frequency and attenuation
peak are majorly affected by the $\sigma_{0}$, $\chi$ and $k_{s}$ (Fig. 12).
At pressures where the rupture occurs, the resonance frequency undergoes a
sudden decrease. The magnitude of the shift in the resonance frequency and the
attenuation are largely affected by the $R_{r}$, $\chi$ and $k_{s}$. At lower
pressures (e.g. 12.5kPa) several parameters of $\sigma_{0}$, $\chi$, $R_{r}$
and $k_{s}$ can provide a good fit to the attenuation and sound speed data
(Fig. 12a-b). However, as the pressure increases (Figs. 12c-h), the
predictions of each group diverge and only one group provides the best fit to
the experimental curves at all excitation pressures. These behaviors can not
be captured by the linear model.

Shell properties
---
Shell |  $\chi$(N/m) |  $\sigma_{0}$(N/m) |  $k_{s}$ ($\times$10-9kg/s) |  $R_{r}$ ($\times$$R_{0}$)
A |  0.5 |  0.03 |  10 |  1.1
B |  1 |  0.016 |  7 |  1.1
C |  1.8 |  0.01 |  6.5 |  1.1
D |  0.4 |  0.063 |  10 |  1.1

Table 2: Shell properties of the lipid coated bubble in Fig. 12

.


$N_{\rm cl}\phi_{\rm sp}\sim 3$ to produce $p\sim 1\%$ at $d_{\rm proj}=200$
au, but if they grow to $a_{\rm max}=100\,{\mu\rm{m}}$, the level of iron
inclusions must increase to $N_{\rm cl}\phi_{\rm sp}\sim 30$ to produce $p\sim
1\%$ in the central region. The value of $N_{\rm cl}\phi_{\rm sp}$ must
increase more to produce $p\sim 1\%$ if grains with slow internal relaxation
have wrong IA.

Besides the fraction of iron locked inside grains, magnetic field strength
also affects on the alignment of grains with B . The stronger magnetic field
strength enhances the Larmor precession and the magnetic relaxation (Section
2), allowing more large grains to be efficiently aligned with B in the densest
central region. In our modeling, we have adopted the uniform
$B=134\,{\mu\rm{m}}$ in the entire protostellar core to emphasize the effects
of grain alignment physics on synthetic dust polarization. However, the
magnetic field strength is amplified with increasing gas density inward of the
protostar as a result of the freezing of B-fields. Observations toward Class
0/I YSO L1448 IRS 2 by Kwon et al. (2019) reveal the highly order hourglass-
shape magnetic field in the protostellar system, and they estimated the
magnetic field strength in the inner 100 au region of $B\sim 600\mu\rm G$.
Myers & Basu (2021) also found the positive correlation between the magnetic
field strength measured by the David Chandrasekhar Fermi (DCF) method
(Chandrasekhar & Fermi 1953) and the gas density of $B_{\rm pos}\sim n_{\rm
H}^{0.66\pm 0.05}$ in the sample of 17 low-mass prestellar and protostellar
cores. To understand how the magnetic field strength affects the magnetic
alignment in the densest central region, we increase the magnetic field
strength in the central 1000 au region into $B=1$ mG while keeping the same
$B=134\mu G$ at $r>1000$ au, denoted as model Realistic$-$rIA$-$amax$-$B and
Realistic$-$wIA$-$amax$-$B (Table LABEL:tab:model_table). The variation of
$p(\%)$ with $N_{\rm cl}\phi_{\rm sp}$ in the central region at $d_{\rm
proj}=200$ au is shown in Figure 19.

The left panel of Figure 19 shows the results for model
Realistic$-$rIA$-$amax$-$B and the right panel is for model
Realistic$-$wIA$-$amax$-$B. Comparing to the upper panels of Figure 17, one
gets higher $p(\%)$ for all values of $N_{\rm cl}\phi_{\rm sp}$. For low
$N_{\rm cl}\phi_{\rm sp}<20$, the difference is clearer as a result of the
extend of the range of grain alignment to larger sizes, i.e., higher $a_{\rm
max,JB}^{\rm Lar}$, due to the enhanced Larmor precession by stronger magnetic
field strength. For high $N_{\rm cl}\phi_{\rm sp}>20$ in which most large
micron-sized grains are aligned with B , the increase of $p(\%)$ is not clear
because the stronger B -field strength cannot enhance the magnetic relaxation
over the strong gas randomization here. It means that grains within few
hundreds of au around the protostar are mainly aligned with B by RATs (instead
of MRAT alignment as grains with a high level of iron inclusions in the
envelope).

Figure 19: Dependence of $p$ with $N_{\rm cl}\phi_{\rm sp}$ for different
values of $a_{\rm max}$ for models with stronger B -fields in the central
region Realistic$-$rIA$-$amax$-$B (left column) and model
Realistic$-$wIA$-$amax$-$B (right column) calculated in the central region at
$d_{\rm proj}=200$ au. Similar as the left panel of Figure 17, only SPM grains
can reproduce the observed degree of polarization of $p>1\%$, and higher iron
fractions are required to align grains with B if grains grow to VLGs. The
polarization degree increases with increasing B -field strength for low
$N_{\rm cl}\phi_{\rm sp}\leq 20$ as a result of the enhanced Larmor
precession. However, stronger B -fields does not help to increase $p(\%)$ for
grains with high $N_{\rm cl}\phi_{\rm sp}$ because the MRAT alignment is not
the major grain alignment mechanism in the dense central region.

###  10.4 Implications for the high polarization degree observed by ALMA
toward protostellar envelopes

In dense environments such as protostellar cores, one usually expects to
observe a much lower dust polarization degree than in the diffuse ISM due to
the inefficient alignment of grains with magnetic fields, according to the RAT
theory (Lazarian & Hoang 2007). However, SMA and ALMA observations toward
Class 0/I YSOs reveal a high degree of dust polarization up to $p\sim 10-40\%$
in the envelope (Hull et al., 2014; Cox et al., 2018; Galametz et al., 2019),
much higher than the maximum values of $p\sim 20\%$ observed in the diffuse
ISM (Planck Collaboration et al. 2015). The polarization degree observed in
the central region is low but still detectable of $p>1\%$. Although the RAT
alignment theory expects the rise of $p$ toward the protostar due to
increasing radiation field (Hoang et al. 2021), the reduction of polarization
degree toward the central region is still possible if the randomization effect
by gas-grain collisions is more effective than the grain alignment by RATs.
However, the numerical modeling with POLARIS by Le Gouellec et al. (2020)
showed that the RAT alignment mechanism implemented in the current version of
POLARIS produces a much lower polarization degree than observations toward
both low and high-mass protostars by ALMA. They found that ALMA data only can
be successfully reproduced by the model in which all grains are perfectly
aligned with B .

Theoretical studies in Hoang (2022) and Hoang et al. (2022) and our detailed
modelling of grain alignment by the MRAT mechanism in Section 5.2 clearly show
that large micron-sized grains can achieve the efficient magnetic alignment in
protostellar environments if they contain a high level of iron inclusions.
This effect is more prominent in the envelope, where VLGs can have perfect
alignment and produce the high polarization degree up to $p\sim 40\%$. Thus,
the detection of high $p\sim 40\%$ toward Class 0/I YSOs by ALMA and the
implication of perfect grain alignment by Le Gouellec et al. (2020) can be
successfully reproduced by the MRAT mechanism. Moreover, this finding
indicates that dust grains in protostellar environments must contain iron
inclusions.

###  10.5 Toward constraining iron abundance embedded in dust using dust
polarization

Iron is among the most abundant element in the Universe. Observations reveal
that more than $90\%$ of iron abundance are missing from the gas (Jenkins,
2009), which implies that iron must be dominantly present inside dust grains
(Dwek, 2016). Yet, observational study of iron in dust is rather difficult.
Hoang & Lazarian (2016b) first showed that one can constrain the level of iron
in dust using sub(mm) polarization data. Lazarian & Hoang (2019) obtained the
lower and upper limit for magnetic susceptibility using the ALMA polarization
patterns toward protostellar disks.

A synthetic modeling of dust polarization by Lam et al. (2021), which assumes
the perfect alignment for all grains with Larmor precession faster than the
gas randomization, found that the magnetic susceptibility of grains of $a\geq
1\,{\mu\rm{m}}$ must be increased by a factor of 20 compared to PM to
reproduce $p\geq 1\%$ in the inner 100 au of a dozen of Class 0/I YSOs
observed by ALMA (Cox et al., 2018). Moreover, a synthetic modelling with
POLARIS by Valdivia et al. (2019), which does not consider the grain magnetic
properites, concluded that grains must grow to $a\sim 10-50\,{\mu\rm{m}}$ to
reproduce $p\geq 1\%$ observed in the samples of Class 0 YSOs by Hull et al.
(2014). Recently, numerical studies of grain alignment in protostellar cores
and protoplanetery disks by Hoang (2022) and Hoang et al. (2022) showed that
large micron-sized grains must contain the high level of iron inclusions to be
efficiently aligned with B in such dense environments and produce polarized
thermal emission.

In Section 10.3, we show the dependence of $p$ with the level of iron
inclusions $N_{\rm cl}\phi_{\rm sp}$ in Figure 17. The green shaded area in
the left and right panel of the figure shows the polarization degree observed
by ALMA, JVLA, SMA, etc, toward Class 0/I YSOs (Hull et al. 2014, Ko et al.
2020, Le Gouellec et al. 2020, Cox et al. 2018). It is very clear that grains
must be SPM with the magnetic susceptibility enhanced by a factor of $\geq 5$
compared with PM grains to reproduce $p>1\%$ in the central region and $p\sim
5-40\%$ in the envelope. Using our modeling results of polarization degree and
observed polarization, we can infer again the level of iron inclusions locked
inside dust grains. For example, to reproduce $p\sim 1\%$ at $d_{\rm
proj}=200$ au with $a_{\rm max}=10\,{\mu\rm{m}}$, one needs $N_{\rm
cl}\phi_{\rm sp}\geq 1$. It corresponds to $N_{\rm cl}=200$ if $\phi_{\rm
sp}=0.005$ (or $1.67\%$ of iron abundance locked inside grains under the form
of cluster), or $N_{\rm cl}=10$ if $\phi_{\rm sp}=0.1$ (or $30\%$ of iron
abundance locked inside grains). If grains grow to VLGs with $a_{\rm
max}=100\,{\mu\rm{m}}$, they need more embedded iron of $N_{\rm c}\phi_{\rm
sp}\geq 20$ to be aligned with B in the central region. It corresponds to
$N_{\rm cl}=4000$ and $N_{\rm cl}=20$ if $\phi_{\rm sp}=0.005$ and $\phi_{\rm
sp}=0.1$, respectively.

Note that the main goal of our present study is to demonstrate the key effects
of grain alignment physics on synthetic dust polarization in protostellar
environments, thus we adopt the simple spherical symmetry for the protostellar
core and uniform B-fields along the $z$ direction. We also perform the
synthetic modeling assuming that the magnetic field only lies in the plane of
the sky to understand how the grain alignment efficiency affects properties of
polarized dust emission. In realistic situations, magnetic fields in
protostellar cores usually have the hourglass shape driven by gravitational
collapse toward protostars, which are theoretically predicted by Fiedler &
Mouschovias (1993), Allen et al. (2003), and widely confirmed by observations,
e.g., Girart et al. (2006), Rao et al. (2009), Hull et al. (2014). The complex
B -field structure along the line of sight combined with the tangling of B
-fields by turbulence results in a smaller polarization degree than obtained
by our modeling with uniform B -fields. In addition, the gas density in the
central region of our model is $n_{\rm H}=10^{7}\,{\rm{cm}}^{-3}$, which is
much lower than the measured gas density in protostellar disks of Class 0/I
YSOs of $n_{\rm H}\sim 10^{8}-10^{9}\,{\rm{cm}}^{-3}$ (e.g., Sandell et al.
1991, Takahashi et al. 2012) and the gas density from MHD simulations of
collapsing protostellar cores (Lam et al. 2019). The magnetic field strength
is also correlated with the gas density instead of being uniform as our
assumption. Thus, the higher gas density will decrease the efficiency of grain
alignment, but higher B-field strength would help to enhance MRAT alignment
efficiency in the central region. Detailed synthetic observations of non-ideal
MHD simulations of collapsing protostellar cores are required to accurately
understand how grains around protostars get aligned with B and how much iron
inclusions are needed to reproduce the observed polarization level of $p>1\%$
from ALMA, JVLA, SMA, etc, observations. Such a study also allows us to better
interpret observational data and clarify the role of B-fields the star
formation process. We will present the results from this study in our followup
paper.

###  10.6 Effect of iron inclusions on polarization by dichroic extinction

In protostellar environments, the polarization by dichroic extinction is
expected to be activated at sub-mm (Brauer et al. 2016) due to the presence of
VLGs (Kwon et al. 2009, Testi et al. 2014, Miotello et al. 2014). VLGs can
strongly attenuate polarized thermal emission at sub-mm, and this feature is
suggested to be the origin of the polarization hole (Brauer et al. 2016, Liu
et al. 2016, Ko et al. 2020, Liu 2021). Polarization by dichroic extinction is
suggested to become the main source of polarization at optically thick
wavelengths and cause the $90^{\circ}$ flipping of the polarization pattern
between mm and sub-mm wavelengths observed in IRAS4A (Ko et al. 2020) and
OMC-3/MMS 6 (Liu 2021).

However, the study by Brauer et al. (2016) and the discussion in Ko et al.
(2020) and Liu (2021) assumed that VLGs are efficiently aligned with B . In
Section 5.2, we show that this situation only can happen if grains are SPM
with a high level of iron inclusions (fourth column of Figures 13 and 14). If
it is not the case, VLGs are not efficiently aligned with B , then, they less
attenuate polarized thermal emission and cannot become the main source of
polarization at optically thick wavelengths. Thus, one will not obtain any
change of the polarization pattern between optically thin and optically thick
wavelengths (Figures 13 and 14).

Furthermore, even if VLGs are efficiently aligned B in protostellar
environments, the extinction by aligned VLGs is not efficient in reducing
polarized thermal emission in the central region compared to the effect of the
inefficient grain alignment on the polarization degree (Section 8.2). Brauer
et al. (2016) showed that the efficiency of polarization by dichroic
extinction increases with increasing gas density. But the high gas density
makes the magnetic alignment of VLGs more difficult, which mimics the effect
of polarization by dichroic extinction in dense environments. Thus, simply
attributing the extinction by aligned dust grains to any anomalous features
observed in dense environments without examining in detail the level of iron
inclusions may cause the wrong interpretation of dust polarization.

###  10.7 Implications for the origin of polarization holes and the role of
iron inclusions

Polarization observations toward Class 0 YSOs usually report the reduction of
the polarization degree at sub-mm/mm wavelengths toward the central region of
protostellar cores. For example, the existence of polarization hole is widely
reported on 30 star-forming cores observed at 1.3mm by CARMA (Hull et al.
2014), 10 Class 0 YSOs observed at $870\,{\mu\rm{m}}$ by ALMA (Cox et al.
2018), 6 objects in Bok globules observed at $870\,{\mu\rm{m}}$ by SCUBA
(Henning et al. 2001, Wolf et al. 2003). Multi-wavelength observations toward
NGC 1333 IRAS4A, OMC-3/MMS 6, and NGC 2071 also show the similar feature from
sub-mm to mm wavelengths (Ko et al. 2020, Liu 2021, Fanciullo et al. 2022).

Many scenarios are proposed to explain the polarization hole toward
protostars, including the inefficient grain alignment by RATs (Hoang &
Lazarian 2016b; Hoang et al. 2021), the dichroic extinction by aligned VLGs
(Brauer et al. 2016), or the inclination angle and the geometrical effect of
magnetic fields due to the self-collapse protostellar cores (Tazaki et al.
2017).

Polarization by dichroic extinction suggested by Brauer et al. (2016) is the
promising explanation of the polarization hole. But as we discuss in Section
10.6, dichroic extinction by aligned grains only plays a minor role on
reducing polarized thermal emission around protostars. We show that the
reduced internal alignment and external alignment by the MRAT mechanism in the
dense central region is the more effective mechanism producing the
polarization hole in protostellar cores (Section 8.2). Furthermore, this
mechanism does not require the presence of aligned VLGs in very dense
environments. Thus, we suggest that this scenario as a main mechanism for
explaining the origin of polarization holes in different class 0/I YSOs.

We note that the inefficient grain alignment in dense environments depends on
the magnetic properties of grains. In particular, the polarization hole will
be due to the loss of grain alignment in dense central region if grains are
PM. In contrast, for SPM grains, it is because of the increase of grains with
inefficient IA by slow internal relaxation and the reduced efficiency of MRAT
alignment with increasing gas randomization. Since the magnetic properties of
grains can be inferred from the observed polarization degree (Section 10.5),
one then can know in detail the behind mechanism which is responsible for the
reduction of $p$ toward the central region of protostellar cores.

###  10.8 On the $90^{\circ}$ rotation of polarization vectors

Another interesting feature observed toward Class 0 YSOs is the change of the
polarization pattern with wavelengths. Particularly, the polarization patterns
in the inner 100 au region of Class 0 YSO NGC 1333 IRAS4A (Ko et al. 2020) and
OMC-3/MMS (Liu 2021) change $90^{\circ}$ from mm observed by JVLA to sub-mm
observed by ALMA. Another case is the MC NGC2071 which shows the $90^{\circ}$
difference of P in the central region between $870\,{\mu\rm{m}}$ observed by
JCMT and $158\,{\mu\rm{m}}$ and $214\,{\mu\rm{m}}$ observed by SOFIA.

Fanciullo et al. (2022) argued that the twist of the polarization pattern
observed in NGC 2071 may be due to the presence of grains with wrong IA
suggested by Hoang (2022). In Figure 10, we show that if SPM grains have a
moderate level of iron inclusions and have wrong IA by slow internal
relaxation, the polarization pattern at the outer boundary of the envelope can
change from P $\|$ B at $870\,{\mu\rm{m}}$ (see also the map at mm in Figure
24) to P $\perp$ B at $450\,{\mu\rm{m}}$ due to the change in IA of the
emission source. In our calculations, this feature happens at optically thin
wavelengths. The flipping of P between $214\,{\mu\rm{m}}$ and
$870\,{\mu\rm{m}}$ observed in the MC NGC 2071 also appears in optically thin
central region. Thus, we suppose that this anomalous feature be due to the
change in IA of the emission source, in which JCMT and SOFIA measure polarized
emission from large grains with wrong IA with P $\|$ B and from smaller grains
with right IA with P $\perp$ B , respectively. Indeed, the gas density and the
maximum grain size of the MC NGC 2071 are different from our model adopted for
the protostellar core. However, we expect the similar results for dense
regions inside MCs because the inefficient IA is naturally produced if the gas
randomization is stronger than the internal relaxation.

For the case of IRAS4A and OMC-3/MMS 6, Ko et al. (2020) and Liu (2021) found
that the cores of two objects are optically thick at sub-mm measured by ALMA
and optically thin at mm measured by JVLA. Therefore, they assign the
polarization flipping in the inner 100 au region to the change of the
polarization mechanism from dichroic emission to dichroic extinction. However,
as we show in Section 8.1 and discuss in Section 10.6, the change of the
polarization pattern from optically thin to optically thick wavelengths in
protostellar environments only can happen if VLGs contain a very high level of
iron inclusions (Figures 13 and 14). In our calculations, the inner 1000 au
region becomes optically thick at $250\,{\mu\rm{m}}$, and VLGs must have
$N_{\rm cl}\phi_{\rm sp}=50$ ($\sim 5000$ larger than PM) to produce the
polarization flipping between $250\,{\mu\rm{m}}$ and $89\,{\mu\rm{m}}$. For
IRAS4A and OMC-3/MMS 6, Ko et al. (2020) and Liu (2021) found that inner 100
au region becomes optically thick at $\leq 1$ mm, indicating that these
objects are very dense. Thus, to activate the effect of dichroic extinction at
$\leq 1$ mm, grains in IRAS4A and OMC-3/MMS 6 must have much higher iron
inclusions than our studied case. However, the amount of iron locked inside
grains still be constrained by the crystalline structure of grains (Yang
2021). Thus, one cannot increase amount of iron too much to make the scenario
proposed by Ko et al. (2020) and Liu (2021) come true. In case that
polarization by dichroic extinction is not effective, we suppose that the
change in IA of the emission source from inefficient IA at mm to efficient IA
at sub-mm to be the reason of the $90^{\circ}$ flipping of P with wavelengths.
However, IRAS4A and OMC-3/MMS 6 have higher luminosity and higher gas density
(Sandell et al. 1991, Takahashi et al. 2012) than our studied parameters. The
strong radiation field can allow more large grains near the source to be
aligned with B with efficient IA (Figure 6, left panel) and activate
polarization by dicroic extinction at sub-mm. However, the strong radiation
field also can disrupt large grains by RAdiative Torque Disruption (RAT-D)
mechansim proposed by Hoang et al. (2019) (Hoang et al. 2021). If it is the
case, the polarization by dichroic extinction cannot cause the change of P
with wavelengths. The detail study of grain alignment by the MRAT mechanism
and grain disruption by RAT-D around high-luminosity protostar is required to
accuracy figure out the physics behind the twist of the polarization pattern
of IRAS4A and OMC-3/MMS 6.

One interesting feature is that if the change of the polarization pattern with
wavelengths is caused by the change in IA of the emission, it implies that
grains with slow internal relaxation in the studied objects have wrong IA.
Until now, the detailed study about the dynamics of grains with slow internal
relaxation is still missing. Thus, it is necessary to study basic properties
and effects of environments on establishing the alignment direction of grains
with slow internal relaxation. Knowing this information allows us to
accurately model grain alignment in dense environments, which is the key for
better analyzing dust polarization.

###  10.9 From polarization vectors to the magnetic field direction

Here we discuss the conditions for which the polarization vectors can be
rotated by $90^{\circ}$ to infer the magnetic field direction in dense
environments. As discussed in Section 10.2, the orientation of P with B
strongly depends on the level of iron locked inside grains and the internal
alignment direction. Furthermore, multi-wavelength observations from mm to
sub-mm should be done to accurately answer the question whether we should
rotate P by $90^{\circ}$ or not to infer B .

In particular, if the polarization pattern does not change from mm to sub-mm,
it can be classified to three cases: SPM grains with the high level of iron
inclusions, SPM grains with a moderate level of iron inclusions which have
right IA with slow internal relaxation, and PM grains. For the first and
second cases, grains always align with their longest axis perpendicular to B
and thus emit polarized thermal radiation with P $\perp$ B . This case is
accompanied by a high polarization degree of $p>5\%$ in the envelope and
$p>1\%$ in the central region. With this case, we can confidently rotate P by
$90^{\circ}$ to reconstruct B -fields morphology on the plane of sky. For the
last case, P can be perpendicular or parallel to B , but the expected
polarization degree is rather low of $p\sim 1\%$ in the envelope and $p\ll
1\%$ in the central region. If it is the case, we should keep P because we
don’t exactly know whether the grain longest axis is perpendicular to parallel
to B , or dust polarization does not trace the magnetic field direction. Until
now, all observations show the high polarization level of $p>5\%$ in
protostellar envelopes, which requires dust grains to be SPM rather than PM.
Therefore, one can infer B-fields in the envelope by rotating the polarization
orientation by $90^{\circ}$.

If the polarization pattern changes with wavelengths, it can originate from
the change in IA of the emission source (Figures 10 and 9), or due to the
change of the polarization mechanism. The latter case happens at optically
thick wavelengths, but it requires the additional information of iron
fractions locked inside dust grains before giving the conclusion as discussed
in Section 10.8. Given this case, the polarization signal at optically thin
wavelengths is due to polarization by dichroic emission while signal at
optically thick wavelengths is due to polarization by dichroic extinction.
Thus, we can rotate P at optically thin wavelengths by $90^{\circ}$ to infer B
, and keep P at optically thick wavelengths. If the origin is the former
scenario, large grains tend to have wrong IA due to slow internal relaxation
and produce P $\|$ B , while small grains can have efficient IA and produce P
$\perp$ B . For this case, we should keep P at long wavelengths and rotate P
at shorter wavelengths by $90^{\circ}$ to reconstruct magnetic fields
structure.

In conclusion, we suggest that multi-wavelength polarimetric observations
toward protostellar cores are required to accurately probe the dust, gas, and
magnetic properties around protostars. A combined analysis of multi-wavelength
polarization patterns and polarization degrees from synthetic observations
with real observations will help us to first constrain amount of iron locked
inside grains, then the grain alignment efficiency in studied environments,
and lastly magnetic field geometry.

##  11 Conclusions

In our study, we have incorporated the magnetic properties of dust grains and
the advanced grain alignment within the RAT paradigm into the POLARIS code.
Then, we used the updated version of POLARIS to study the effects of iron
inclusions on grain alignment and synthetic polarimetric observations of the
protostellar core. Our findings are summarized as follows:

  1. 1.

Paramagnetic (PM) grains are not aligned with B in the central region with
high gas density of $n_{\rm H}\geq 10^{7}\,{\rm{cm}}^{-3}$ due to the grain
randomization by gas collisions that occurs faster than the Larmor precession.
They can be aligned with B in the envelope but always have inefficient IA due
to slow internal relaxation. In contrast, superparamagnetic (SPM) grains can
have the magnetic alignment with efficient IA in the protostellar core due to
the enhanced Barnett relaxation and Larmor precession by iron inclusions. The
grain alignment efficiency increases with the increasing level of iron
inclusions locked inside dust grains. However, in the protostellar central
region, large grains of $a\geq 1\,{\mu\rm{m}}$ aligned with B at low-J
attractors always have slow internal relaxation and can align their longest
axis perpendicular or parallel with B .

  2. 2.

The polarization patterns produced by both PM and SPM grains are uniform with
P $\perp$ B at optically thin wavelengths when grains with slow internal
relaxation have right IA. In contrast, when grains with slow internal
relaxation have wrong IA, the polarization pattern produced by PM grains will
be uniform in space but with P $\|$ B . SPM grains with a moderate level of
iron inclusions produce the uniform pattern with P $\|$ B at mm, but they
produce complex ones at sub-mm with P $\perp$ B in the boundary of the
envelope and in the inner $\sim 300$ au region, and P $\|$ B in the middle
region of the core. The change in polarization pattern with the radial
distance and wavelengths at optically thin wavelengths is due to the change in
the internal alignment (IA) of the emission source. For SPM grains with a high
level of iron inclusions, the polarization pattern is uniform, with P $\perp$
B from mm to sub-mm wavelengths because of the efficient grain alignment by
the MRAT mechanism.

  3. 3.

The polarization degree increases with increasing the abundance of iron
present in dust, $N_{\rm cl}\phi_{\rm sp}$, but always decreases in the
central region near the protostar. PM grains produce a low polarization of
$p\sim 1\%$ in the envelope because of their inefficient grain alignment and a
negligible polarization level in the inner 1000 au because they are not
aligned with B there. In contrast, SPM grains with a high level of iron
inclusions can achieve perfect alignment by the MRAT mechanism and produce a
very high polarization degree up to $p\sim 40\%$ in the envelope. It
successfully explains the detection of high $p\sim 40\%$ by ALMA toward Class
0 YSOs. The polarization fraction produced by SPM grains slightly decreases to
$p<10\%$ in the central region due to the reduced MRAT alignment efficiency in
very dense environments.

  4. 4.

The polarization by dichroic extinction at optically thick sub-mm wavelengths
is less efficient than previously thought due to the weak alignment of VLGs in
protostellar environments. We found that the inefficient grain alignment in
dense regions induces much stronger reduction of $p(\%)$ with increasing
intensity than the extinction effect by VLGs. Furthermore, dichroic extinction
cannot become the main source of polarization at sub-mm and cannot induce the
polarization flipping from P $\perp$ B at optically thin wavelengths to P $\|$
B at optically thick wavelengths, excepts for SPM grains with very high levels
of iron inclusions. The estimated amount of iron inclusions locked inside dust
is required to accurately quantify the efficiency of grain alignment and
interpret polarimetric data observed at optically thick wavelengths.

  5. 5.

We found that if grains with wrong IA exist in protostellar cores, the change
in IA of the polarized emission source can cause the $90^{\circ}$ flipping of
the polarization pattern between the envelope and the central region, and
between mm and sub-mm. This feature happens at optically thin wavelengths, and
be a potential mechanism for the variation of P with distances and
wavelengths. In addition, it also helps us to study the dynamic of grains with
slow internal relaxation.

  6. 6.

The polarization degree increases with increasing iron abundance in the dust
but decreases with increasing the maximum grain size. VLGs must be SPM with
the magnetic susceptibility enhanced by a factor of $\geq 10^{3}$ compared to
PM to reproduce the observed level of polarization of $p\geq 1\%$ in the
central region and high $p\sim 5-40\%$ in the envelope. Thus, it opens the new
tool for tracing the level of iron inclusions locked inside dust grains via
dust polarimetry. Knowing the magnetic properties of grains in dense
environments is very important for probing the reability of studying magnetic
fields by using dust polarization.

  7. 7.

We find that, in addition to the increase of gas randomization toward
protostars, the decrease of grain magnetic susceptibility with increasing dust
temperature can also decrease the efficiency of external alignment (via the
Larmor precession) and internal alignment (via Barnett relaxation). Therefore,
although the RAT alignment increases toward the protostar, the decrease of
grain magnetic susceptibility and increase of gas randomization leads to the
rapid decrease of $p$ toward the protostar. This is a new effect that can
reproduce the polarization hole in the central region.

## Acknowledgements

We thank Stefan Reissl for help with the Polaris code during the initial stage
of this work. We thank the members of Vietnam Astrophysics Research Network
(VARNET) for various discussion and comments. T.H. is supported by the
National Research Foundation of Korea (NRF) grant funded by the Korea
government (MSIT (No. 2019R1A2C1087045). J.-G.K acknowledges support from the
EACOA Fellowship awarded by the East Asia Core Observatories Association.

## Data Availability

The data underlying this article will be shared on reasonable request to the
corresponding author.

## References

  * Allen et al. (2003) Allen A., Li Z.-Y., Shu F. H., 2003, ApJ, 599, 363
  * Andersson et al. (2015) Andersson B. G., Lazarian A., Vaillancourt J. E., 2015, ARA& A, 53, 501
  * Barnett (1915) Barnett S. J., 1915, Phys. Rev., 6, 239
  * Bethell et al. (2007) Bethell T. J., Chepurnov A., Lazarian A., Kim J., 2007, ApJ, 663, 1055
  * Bok (1977) Bok B. J., 1977, PASP, 89, 597
  * Brauer et al. (2016) Brauer R., Wolf S., Reissl S., 2016, A&A, 588, A129
  * Chandrasekhar & Fermi (1953) Chandrasekhar S., Fermi E., 1953, ApJ, 118, 113
  * Cho & Lazarian (2007) Cho J., Lazarian A., 2007, The Astrophysical Journal, 669, 1085
  * Clemens et al. (1991) Clemens D. P., Yun J. L., Heyer M. H., 1991, ApJS, 75, 877
  * Clemens et al. (2020) Clemens D. P., et al., 2020, ApJS, 249, 23
  * Collaboration et al. (2020) Collaboration P., et al., 2020, Astronomy and Astrophysics, 641, A12
  * Cox et al. (2018) Cox E. G., Harris R. J., Looney L. W., Li Z.-Y., Yang H., Tobin J. J., Stephens I., 2018, The Astrophysical Journal, 855, 92
  * Crutcher (2012) Crutcher R. M., 2012, ARA& A, 50, 29
  * Davis & Greenstein (1951) Davis Leverett J., Greenstein J. L., 1951, ApJ, 114, 206
  * Dent et al. (2019) Dent W. R. F., Pinte C., Cortes P. C., Ménard F., Hales A., Fomalont E., de Gregorio-Monsalvo I., 2019, MNRAS, 482, L29
  * Dolginov & Mitrofanov (1976) Dolginov A. Z., Mitrofanov I. G., 1976, Ap&SS, 43, 291
  * Draine (1996) Draine B. T., 1996, in Roberge W. G., Whittet D. C. B., eds, Astronomical Society of the Pacific Conference Series Vol. 97, Polarimetry of the Interstellar Medium. p. 16 (arXiv:astro-ph/9603053)
  * Draine & Hensley (2021) Draine B. T., Hensley B. S., 2021, The Astrophysical Journal, 909, 94
  * Draine & Lazarian (1998) Draine B. T., Lazarian A., 1998, ApJ, 508, 157
  * Draine & Weingartner (1996) Draine B. T., Weingartner J. C., 1996, ApJ, 470, 551
  * Dwek (2016) Dwek E., 2016, The Astrophysical Journal, 825, 136
  * Fanciullo et al. (2022) Fanciullo L., et al., 2022, Monthly Notices of the Royal Astronomical Society, 512, 1985
  * Fiedler & Mouschovias (1993) Fiedler R. A., Mouschovias T. C., 1993, ApJ, 415, 680
  * Frank et al. (2014) Frank A., et al., 2014, in Beuther H., Klessen R. S., Dullemond C. P., Henning T., eds, Protostars and Planets VI. p. 451 (arXiv:1402.3553), doi:10.2458/azu_uapress_9780816531240-ch020
  * Galametz et al. (2018) Galametz M., et al., 2018, A&A, 616, A139
  * Galametz et al. (2019) Galametz M., Maury A. J., Valdivia V., Testi L., Belloche A., André P., 2019, A&A, 632, A5
  * Girart et al. (2006) Girart J. M., Rao R., Marrone D. P., 2006, Science, 313, 812
  * Greenberg (1968) Greenberg J. M., 1968, in Middlehurst B. M., Aller L. H., eds, , Nebulae and Interstellar Matter. p. 221
  * Guillet et al. (2020) Guillet V., Hennebelle P., Forêts G. P. d., Marcowith A., Commerçon B., Marchand P., 2020, Astronomy & Astrophysics, 643, A17
  * Henning et al. (2001) Henning T., Wolf S., Launhardt R., Waters R., 2001, ApJ, 561, 871
  * Hensley & Draine (2021) Hensley B. S., Draine B. T., 2021, The Astrophysical Journal, 906, 0
  * Herranen et al. (2021) Herranen J., Lazarian A., Hoang T., 2021, ApJ, 913, 63
  * Hirashita & Li (2013) Hirashita H., Li Z. Y., 2013, MNRAS, 434, L70
  * Hoang (2022) Hoang T., 2022, The Astrophysical Journal, 928, 102
  * Hoang & Lazarian (2008) Hoang T., Lazarian A., 2008, Monthly Notices of the Royal Astronomical Society, 388, 117
  * Hoang & Lazarian (2009) Hoang T., Lazarian A., 2009, The Astrophysical Journal, 697, 1316
  * Hoang & Lazarian (2014) Hoang T., Lazarian A., 2014, Monthly Notices of the Royal Astronomical Society, 438, 680
  * Hoang & Lazarian (2016a) Hoang T., Lazarian A., 2016a, The Astrophysical Journal, 821, 91
  * Hoang & Lazarian (2016b) Hoang T., Lazarian A., 2016b, ApJ, 831, 159
  * Hoang et al. (2019) Hoang T., Tram L. N., Lee H., Ahn S.-H., 2019, Nature, 3, 766–775
  * Hoang et al. (2021) Hoang T., Tram L. N., Lee H., Diep P. N., Ngoc N. B., 2021, ApJ, 908, 218
  * Hoang et al. (2022) Hoang T., Tram L. N., Minh Phan V. H., Giang N. C., Phuong N. T., Dieu N. D., 2022, AJ, 164, 248
  * Hull & Zhang (2019) Hull C. L. H., Zhang Q., 2019, Frontiers in Astronomy and Space Sciences, 6, 254
  * Hull et al. (2014) Hull C. L. H., et al., 2014, The Astrophysical Journal Supplement Series, 213, 13
  * Jenkins (2009) Jenkins E. B., 2009, The Astrophysical Journal, 700, 1299
  * Jones & Spitzer (1967) Jones R. V., Spitzer Lyman J., 1967, ApJ, 147, 943
  * Kataoka et al. (2012) Kataoka A., Machida M. N., Tomisaka K., 2012, ApJ, 761, 40
  * Kataoka et al. (2013) Kataoka A., Tanaka H., Okuzumi S., Wada K., 2013, A&A, 557, L4
  * Ko et al. (2020) Ko C.-L., Liu H. B., Lai S.-P., Ching T.-C., Rao R., Girart J. M., 2020, ApJ, 889, 172
  * Krumholz & Federrath (2019) Krumholz M. R., Federrath C., 2019, Frontiers in Astronomy and Space Sciences, 6, 7
  * Kwon et al. (2009) Kwon W., Looney L. W., Mundy L. G., Chiang H.-F., Kemball A. J., 2009, ApJ, 696, 841
  * Kwon et al. (2019) Kwon W., Stephens I. W., Tobin J. J., Looney L. W., Li Z.-Y., van der Tak F. F. S., Crutcher R. M., 2019, The Astrophysical Journal, 879, 25
  * Lam et al. (2019) Lam K. H., Li Z.-Y., Chen C.-Y., Tomida K., Zhao B., 2019, Monthly Notices of the Royal Astronomical Society, 489, 5326
  * Lam et al. (2021) Lam K. H., Chen C.-Y., Li Z.-Y., Yang H., Cox E. G., Looney L. W., Stephens I., 2021, MNRAS, 507, 608
  * Lazarian (2007) Lazarian A., 2007, J. Quant. Spectrosc. Radiative Transfer, 106, 225
  * Lazarian & Hoang (2007) Lazarian A., Hoang T., 2007, MNRAS, 378, 910
  * Lazarian & Hoang (2008) Lazarian A., Hoang T., 2008, ApJL, 676, L25
  * Lazarian & Hoang (2019) Lazarian A., Hoang T., 2019, The Astrophysical Journal, 883, 122
  * Lazarian & Hoang (2021) Lazarian A., Hoang T., 2021, The Astrophysical Journal, 908, 12
  * Lazarian & Roberge (1997) Lazarian A., Roberge W., 1997, The Astrophysical Journal, 484, 230
  * Lazarian et al. (2015) Lazarian A., Andersson B.-G., Hoang T., 2015, in Kolokolova L., Hough J., Levasseur-Regourd A.-C., eds, Polarimetry of stars and planetary systems. (New York: Cambridge Univ. Press), p. 81
  * Le Gouellec et al. (2020) Le Gouellec V. J. M., et al., 2020, A&A, 644, A11
  * Leung (1985) Leung C. M., 1985, in Black D. C., Matthews M. S., eds, Protostars and Planets II. pp 104–136
  * Liu (2021) Liu H. B., 2021, The Astrophysical Journal, 914, 25
  * Liu et al. (2016) Liu H. B., et al., 2016, ApJ, 821, 41
  * Lucy (1999) Lucy L. B., 1999, A&A, 344, 282
  * Martin (1974) Martin P. G., 1974, ApJ, 187, 461
  * Mathis et al. (1977) Mathis J. S., Rumpl W., Nordsieck K. H., 1977, ApJ, 217, 425
  * Mathis et al. (1983) Mathis J. S., Mezger P. G., Panagia N., 1983, A&A, 500, 259
  * McKee & Ostriker (2007) McKee C. F., Ostriker E. C., 2007, ARA& A, 45, 565
  * Miotello et al. (2014) Miotello A., Testi L., Lodato G., Ricci L., Rosotti G., Brooks K., Maury A., Natta A., 2014, A&A, 567, A32
  * Morrish (2001) Morrish A. H., 2001, The Physical Principles of Magnetism
  * Myers & Basu (2021) Myers P. C., Basu S., 2021, ApJ, 917, 35
  * Nakano & Nakamura (1978) Nakano T., Nakamura T., 1978, PASJ, 30, 671
  * Okuzumi & Tazaki (2019) Okuzumi S., Tazaki R., 2019, ApJ, 878, 132
  * Okuzumi et al. (2012) Okuzumi S., Tanaka H., Kobayashi H., Wada K., 2012, ApJ, 752, 106
  * Pattle & Fissel (2019) Pattle K., Fissel L., 2019, Frontiers in Astronomy and Space Sciences, 6, 15
  * Planck Collaboration et al. (2015) Planck Collaboration et al., 2015, A&A, 576, A104
  * Purcell (1979) Purcell E. M., 1979, ApJ, 231, 404
  * Rao et al. (2009) Rao R., Girart J. M., Marrone D. P., Lai S.-P., Schnee S., 2009, ApJ, 707, 921
  * Reissl et al. (2016) Reissl S., Wolf S., Brauer R., 2016, A&A, 593, A87
  * Reissl et al. (2017) Reissl S., Seifried D., Wolf S., Banerjee R., Klessen R. S., 2017, Astronomy and Astrophysics, 603, A71
  * Reissl et al. (2020) Reissl S., Guillet V., Brauer R., Levrier F., Boulanger F., Klessen R. S., 2020, A&A, 640, A118
  * Roberge et al. (1993) Roberge W. G., Degraff T. A., Flaherty J. E., 1993, ApJ, 418, 287
  * Sandell et al. (1991) Sandell G., Aspin C., Duncan W. D., Russell A. P. G., Robson E. I., 1991, ApJL, 376, L17
  * Seifried & Walch (2015) Seifried D., Walch S., 2015, MNRAS, 452, 2410
  * Seifried et al. (2018) Seifried D., Walch S., Reissl S., Ibáñez-Mejía J. C., 2018, arXiv.org, astro-ph.GA
  * Shu et al. (1987) Shu F. H., Adams F. C., Lizano S., 1987, ARA& A, 25, 23
  * Takahashi et al. (2012) Takahashi S., Saigo K., Ho P. T. P., Tomida K., 2012, The Astrophysical Journal, 752, 10
  * Tazaki et al. (2017) Tazaki R., Lazarian A., Nomura H., 2017, ApJ, 839, 56
  * Testi et al. (2014) Testi L., et al., 2014, in Beuther H., Klessen R. S., Dullemond C. P., Henning T., eds, Protostars and Planets VI. p. 339 (arXiv:1402.1354), doi:10.2458/azu_uapress_9780816531240-ch015
  * Tram & Hoang (2022) Tram L. N., Hoang T., 2022, arXiv e-prints, p. arXiv:2208.13195
  * Valdivia et al. (2019) Valdivia V., Maury A., Brauer R., Hennebelle P., Galametz M., Guillet V., Reissl S., 2019, MNRAS, 488, 4897
  * Wolf et al. (2003) Wolf S., Launhardt R., Henning T., 2003, The Astrophysical Journal, 592, 233
  * Wolf et al. (2004) Wolf S., Launhardt R., Henning T., 2004, Ap&SS, 292, 239
  * Yang (2021) Yang H., 2021, ApJ, 911, 125

##  Appendix A Polarized radiative transfer equation for the Stokes parameters

The polarization state of radiation can be described by the Stokes vector S
$=(I~{}Q~{}U~{}V)^{\rm T}$ where Stokes $I$ presents the total intensity,
Stokes $Q$ and $U$ present the linear polarization, and Stokes $V$ presents
the circular polarization. The radiative transfer equation of the Stokes
vector through a dusty environment is given by (Martin 1974, Reissl et al.
2016):

$\displaystyle\frac{d}{ds}\begin{pmatrix}I\\\ Q\\\ U\\\
V\end{pmatrix}=-\begin{pmatrix}\alpha_{\rm I}&&\alpha_{\rm Q}&&0&&0\\\
\alpha_{\rm Q}&&\alpha_{\rm I}&&0&&0\\\ 0&&0&&\alpha_{\rm I}&&\kappa_{\rm
Q}\\\ 0&&0&&-\kappa_{\rm Q}&&\alpha_{\rm I}\\\
\end{pmatrix}\begin{pmatrix}I\\\ Q\\\ U\\\
V\end{pmatrix}+\begin{pmatrix}j_{\rm I}\\\ j_{\rm Q}\\\ 0\\\ 0\end{pmatrix}$
(35)

where $\alpha_{\rm I}$ is the extinction coefficient over the grain size
distribution:

$\displaystyle\alpha_{\rm I}=\int_{a_{\rm min}}^{a_{\rm max}}C_{\rm
ext}(a,\lambda)(dn/da)da,$ (36)

$\alpha_{\rm Q}$ is the linear polarization coefficient due to extinction:

$\displaystyle\alpha_{\rm Q}=\int_{a_{\rm align}}^{a_{\rm max,JB}^{\rm
Lar}}C_{\rm pol}^{\rm ext}(a,\lambda)(dn/da)da,$ (37)

$\kappa_{\rm Q}$ is the circular polarization coefficient:

$\displaystyle\kappa_{\rm Q}=\int_{a_{\rm align}}^{a_{\rm max,JB}^{\rm
Lar}}C_{\rm circ}(a,\lambda)(dn/da)da,$ (38)

where $C_{\rm ext}(a,\lambda)$ is the extinction cross section, $C_{\rm
ext}^{\rm pol}(a,\lambda)$ is the linear polarization cross section, and
$C_{\rm circ}(a,\lambda)$ is the circular polarization cross section of the
grain of size $a$ and radiation wavelength $\lambda$.

Above, $j_{\rm I}$ and $j_{\rm Q}$ are the thermal dust emissivity from all
dust grains and aligned dust grains, respectively. For dust grains at thermal
equilibrium of temperature $T_{d}$, they are given by

$\displaystyle j_{\rm I}=\int_{a_{\rm min}}^{a_{\rm max}}C_{\rm
abs}(a,\lambda)B_{\lambda}(T_{\rm d}(a))(dn/da)da,$ (39)

and

$\displaystyle j_{\rm Q}=\int_{a_{\rm align}}^{a_{\rm max,JB}^{\rm Lar}}C_{\rm
pol}^{\rm abs}(a,\lambda)B_{\lambda}(T_{\rm d}(a))(dn/da)da,$ (40)

where $B_{\rm\lambda}(T_{\rm d}(a))$ is the Planck function at the grain
temperature $T_{\rm d}(a)$, $C_{\rm abs}(a,\lambda)$ is the absorption cross
section, and $C_{\rm pol}^{\rm abs}(a,\lambda)$ is the polarization cross
section by absorption of the grain of size $a$.

The extinction cross section $C_{\rm ext}$ in Equation (36) is calculated as :

$\displaystyle C_{\rm ext}=\frac{C_{\rm ext,x}+C_{\rm ext,y}}{2},$ (41)

and the linear polarization cross section $C_{\rm pol}^{\rm ext}$ in Equation
(37) is given by:

$\displaystyle C_{\rm pol}^{\rm ext}=\frac{C_{\rm ext,x}-C_{\rm ext,y}}{2}.$
(42)

$C_{\rm ext,x}$ and $C_{\rm ext,y}$ are the extinction cross section of grain
of size for the electric field E oscillating along x and y axis on the plane
of sky, which is:

$\displaystyle C_{\rm ext,x}=\langle C_{\rm
ext}\rangle+\frac{1}{3}~{}R~{}(C_{\rm ext,\parallel}-C_{\rm ext,\perp}),$ (43)

and

$\displaystyle C_{\rm ext,y}=\langle C_{\rm
ext}\rangle+\frac{1}{3}~{}R~{}(C_{\rm ext,\parallel}-C_{\rm
ext,\perp})(1-3\sin^{2}\psi)$ (44)

The term $\langle C_{\rm ext}\rangle$ is the extinction cross section over
grain orientation in the grain frame, which is:

$\displaystyle\langle C_{\rm ext}\rangle=\frac{2C_{\rm ext,\parallel}+C_{\rm
ext,\perp}}{3},$ (45)

where $C_{\rm ext,\parallel}$ and $C_{\rm ext,\perp}$ are the extinction cross
section produced when the electric field E oscillating along the major and
minor axis of the grain, respectively.

Putting all terms above into Equations (41) and (42), one obtains

$\displaystyle C_{\rm ext}$ $\displaystyle=$
$\displaystyle\frac{(4+3R\cos^{2}\psi-R)C_{\rm ext,\|}}{6}$ (46)
$\displaystyle+$ $\displaystyle\frac{(2-3\cos^{2}\psi+R)C_{\rm
ext,\perp}}{6},$

and $C_{\rm pol}^{\rm ext}$ is:

$\displaystyle C_{\rm pol}^{\rm ext}=\frac{(C_{\rm ext,\|}-C_{\rm
ext,\perp})R\sin^{2}\psi}{2},$ (47)

where $R$ is the Rayleigh reduction factor described in Section 4.1.3.

The values of $C_{\rm abs}$, $C_{\rm pol}^{\rm abs}$, and $C_{\rm circ}$ have
the same form with $C_{\rm ext}$ and $C_{\rm pol}^{\rm ext}$.

The term $dn/da$ in Equation (36) is the grain size distribution. Here we
assume the typical power law of $dn/da=Ca^{-3.5}$ as in the diffuse ISM (the
MRN distribution, Mathis et al. 1977) where $C$ the normalization constant is
determined by the gas-to-dust mass ratio:

$\displaystyle\eta=\frac{\int_{a_{\rm min}}^{a_{\rm max}}m_{\rm d}(a)n_{\rm
H}Ca^{-3.5}da}{n_{\rm H}m_{\rm H}},$ (48)

where $m_{\rm d}(a)=4\pi\rho sa^{3}/3$ is the mass of the grain of size $a$
and axial ratio $s$.

##  Appendix B Validity of the one-cell assumption in the azimuthal direction

Figure 20: Variation of the mean radiation field strength U (left panel) and
the anisotropic degree $\gamma_{\rm rad}$ with distances for the spherical
grid model with $N_{\rm r}\times N_{\vartheta}\times N_{\varphi}=250\times
101\times 1$ and $N_{\rm r}\times N_{\vartheta}\times N_{\varphi}=250\times
101\times 60$, respectively. The grid model with one cell along the azimuthal
direction gives a similar average $U$ and $\gamma_{\rm rad}$ as the model with
60 cells on the azimuthal direction due to the symmetry of the gas density
distribution inside the spherical core.

Figure 20 shows the comparison between the mean radiation field strength $U$
(left panel) and the anisotropic degree $\gamma_{\rm rad}$ (right panel) as a
function of distances to the central protostar using a spherical grid cell
with $N_{\rm r}\times N_{\vartheta}\times N_{\varphi}=250\times 101\times 1$
and $N_{\rm r}\times N_{\vartheta}\times N_{\varphi}=250\times 101\times 60$,
respectively. One can clearly see that $U$ and $\gamma_{\rm rad}$ calculated
from the model with $N_{\rm varphi}=1$ are nearly similar to the results from
the model with $N_{\rm varphi}=60$ due to the spherical symmetry of the gas
density model. Thus, using one cell along the azimuthal direction does not
affect results of our adopted model.

##  Appendix C Grain alignment along x and z direction

Figure 21: Variation along the $z$ direction of the range of aligned dust
grains $a_{\rm align}-a_{\rm max,JB}^{\rm Lar}$ (upper left panel), grains
with fast internal relaxation at high and low-J attractors $a_{\rm
min,aJ}-a_{\rm max,aJ}$ (upper right and lower left panels), and the maximum
size that grains are aligned with B by the MRAT mechanism, $a_{\rm
max,JB}^{\rm DG,0.5}$, and $a_{\rm max,JB}^{\rm DG,1}$. The results are for PM
and SPM grains with different $N_{\rm cl}$. Generally, large grains can be
aligned with B and have fast internal relaxation in the envelope due to the
reduced gas randomization. Iron inclusions enhance the internal and external
alignment and help more large grains in the envelope to have perfect alignment
with B by MRAT alignment.

Figures 21 and 22 shows in detail the grain alignment of PM and SPM grains
along the z and x direction, respectively (see the gas density map on x-z
plane in the left panel of Figure 4). The upper left panel shows the minimum
and maximum grain alignment size $a_{\rm align}$ and $a_{\rm max,JB}^{\rm
Lar}$. The upper right and lower left panels show the minimum and maximum size
for grains having fast internal relaxation at high-J and low-J attractors,
respectively. The lower right panel shows the maximum size for grains having
$f_{\rm high-J}=0.5$ and $f_{\rm high-J}=1$. In general, more large grains can
be aligned with B (i.e., larger range $a_{\rm align}-a_{\rm max,JB}^{\rm
Lar}$), and have fast internal relaxation at high-J attractors (i.e., larger
range $a_{\rm min,aJ}^{\rm high-J}-a_{\rm max,aJ}^{\rm high-J}$), in the
envelope due to the reduced gas randomization. Iron inclusions allow more
grains to have perfect magnetic alignment (i.e., $f_{\rm high-J}=1$) there due
to the enhanced MRAT alignment.

However, one can see that the alignment efficiency of grains on the $x$
direction is slightly smaller than grains on the $z$ direction. It is because
grains are spun up by RATs weaker in the area where the radiation field is
perpendicular with B-field (Equation 18). As a result, the range of aligned
dust grains having fast internal relaxation at high-J attraction reduces,
i.e., larger $a_{\rm align}$ and larger $a_{\rm max,aJ}^{\rm high-J}$,
inducing lower $p(\%)$ on the equatorial as shown in the polarization maps in
Section 6.

Figure 22: Similar as Figure 21 but for grains on the x direction. The
variation of grain alignment along $x$ is similar to along the $z$-direction,
but alignment efficiency is weaker, i.e., narrow alignment range of $a_{\rm
align}-a_{\rm max,JB}^{\rm Lar}$ (upper left panel) and narrow range of grains
with fast internal relaxation at high-J attractors $a_{\rm min,aJ}^{\rm
high-J}-a_{\rm max,aJ}^{\rm high-J}$ (upper right panel), due to the reduced
RATs efficiency in the area where k $\perp$ B .

##  Appendix D Polarization map at millimeter wavelengths for $a_{\rm
max}=100\,{\mu\rm{m}}$

Figure 23 shows the polarization patterns calculated at 2mm by model Ideal
(first column) and model Realistic$-$rIA. The upper panels show the results in
the envelope, and the lower panels show the results for the inner region of
1000 au around the protostar. The polarization patterns produced by both PM
and SPM are uniform with P $\perp$ B from the envelope to the central region,
assuming that grains with slow internal relaxation have right IA.

Figure 23: Polarization maps calculated at 2 mm for model Ideal (first column)
and model Realistic$-$rIA. The upper panels show the results in the envelope,
and the lower panels are for the inner 1000 au around the protostar. The
polarization patterns for both PM and SPM grains are uniform with P $\perp$ B
from the envelope to the center when grains with slow internal alignment have
right IA. The polarization degree in our model decreases inward due to the
inefficient alignment of grains with B . Figure 24: Similar as Figure 23 but
for model Realistic$-$wIA. At mm, PM and SPM grains with low $N_{\rm cl}=100$
produce the uniform polarization pattern with P $\|$ B in the entire
protostellar core due to the emission of grains with wrong IA. In contrast,
SPM grains with high $N_{\rm cl}=10^{4}$ produce polarization vectors P
$\perp$ B as model Ideal.

Figure 24 shows the similar results as Figure 23 but for model
Realistic$-$wIA. Similar as the results at $870\,{\mu\rm{m}}$ (Figures 10 and
11), PM and SPM grains with $N_{\rm cl}=100$ produce the uniform polarizaton
pattern with P $\|$ B from the envelope to the center due to the emission of
VLGs with wrong IA at mm wavelengths. In contrast, SPM grains with $N_{\rm
cl}=10^{4}$ produce the uniform pattern with P $\perp$ B in the entire
protostellar core because of the efficient alignment of VLGs by the MRAT
mechanism with B .

##  Appendix E Polarization results for $a_{\rm max}=10\,{\mu\rm{m}}$

###  E.1 Iron Inclusions and Polarization Pattern

Figure 25: Polarization maps calculated in the envelope (first to third rows)
and the inner region of 1000 au (fourth to sixth rows) for model Ideal$-$a10
(first column) and model Realistic$-$rIA$-$a10 (second to fourth columns)
(Table LABEL:tab:model_table). For model that all grains have right IA, the
polarization pattern is uniform in space with P $\perp$ B from 2mm to
$450\,{\mu\rm{m}}$.

Figure 26: Similar as Figure 25 but for model Realistic$-$wIA$-$a10. In the
presence of grains with wrong IA, at fixed wavelengths, the polarization
pattern produced by PM and SPM grains with low $N_{\rm cl}$ changes from P
$\perp$ B in the outer boundary of the envelope to P $\|$ B in the inner
region and flips again to P $\perp$ B around protostar. This complicated
polarization pattern is due to the complex distribution of grains with right
and wrong IA in the protostellar core (e.g., Figure 6). For SPM grains with
high $N_{\rm cl}=10^{4}$, P is uniform with P $\perp$ B . In addition, the
polarization patterns in the entire protostellar core do not change with
wavelengths for $a_{\rm max}=10\,{\mu\rm{m}}$.

Figure 27: Effect of iron inclusions on the variation of $p$ with $I/I_{\rm
max}$ at $870\,{\mu\rm{m}}$ (upper panels) and $450\,{\mu\rm{m}}$ (lower
panels) for model Realistic$-$rIA$-$a10 (left column) and
Realistic$-$wIA$-$a10 (right column). The polarization degree from our models
always decreases with increasing intensity due to the reduced grain alignment
efficiency in dense environments, and generally increases with increasing
$N_{\rm cl}$. If grains with wrong IA are present in the core, one will obtain
the "valley-polarization hole" on $p-I/I_{\rm max}$ curve caused by the
cancelling effect of polarized emission from grains with right and wrong IA .
Figure 28: Effect of iron inclusions on the polarization curve $p(\lambda)$
from 2 mm to $89\,{\mu\rm{m}}$ calculated at $d_{\rm proj}=200$ au (upper
panels) and $d_{\rm proj}=2000$ au (lower panels) for model
Realistic$-$rIA$-$a10 (left column) and Realistic$-$wIA$-$a10 (right column).
With $a_{\rm max}=10\,{\mu\rm{m}}$, the polarization degree is nearly
independent with mm/sub-mm wavelengths in both the central region and the
envelope (corresponds to the independence of the polarization patterns with
wavelengths in Figures 25 and 26).

Figure 25 shows the polarization maps obtained from the outer envelope (first
to third rows) to the central region (fourth to sixth rows) by model
Realistic$-$rIA$-$a10 (Table LABEL:tab:model_table). From the top to bottom,
the observed wavelengths decrease from 2mm to $870\,{\mu\rm{m}}$ and
$450\,{\mu\rm{m}}$, while from the left to right, the results are for model
Ideal$-$a10 (first column) and our models for PM and SPM grains with $N_{\rm
cl}=100$ and $N_{\rm cl}=10^{4}$, assuming $\phi_{\rm sp}=0.005$. When all
aligned dust grains have right IA, the polarization pattern (black segment) is
uniform with P $\perp$ B from the envelope to the central region and does not
change with wavelengths. The polarization degree (color code) decreases toward
the center due to the inefficient magnetic alignment of grains in dense core.

Figure 26 shows similar results as Figure 25 but for model
Realistic$-$wIA$-$a10. For PM grains (left column), the polarization pattern
is quite uniform with P $\|$ B within 10000 au where most grains have wrong IA
due to slow internal relaxation. However, beyond $\sim 10000$ au, P flips
$90^{\circ}$ to P $\perp$ B because a small amount of micron-sized grains can
have efficient IA by ISRF here (Figure 6, left panel). For SPM grains with
$N_{\rm cl}=100$ (center column), the polarization pattern changes from P
$\perp$ B beyond 5000 au (first three rows) to P $\|$ B in the inner region
and changes back to P $\perp$ B within 200 au from the protostar (fourth to
sixth rows) as a result of the complex distribution of grains with wrong and
right IA from the center to the envelope (Figure 6, left panel). In contrast,
SPM grains with $N_{\rm cl}=10^{4}$ (right column) produce much more simple
polarization pattern, with P $\perp$ B in the entire protostellar core due to
their efficient alignment with B -fields.

One can see that that in the absence of VLGs, the polarization vectors do not
rotate from P $\|$ B at mm to P $\perp$ B at sub-mm as in the case of model
Realistic$-$wIA for SPM grains with $N_{\rm cl}=100$ and $a_{\rm
max}=100\,{\mu\rm{m}}$ (Figures 10 and 11). It is because for the model with
$a_{\rm max}=10\,{\mu\rm{m}}$, most all SPM grains have efficient IA and thus
produce P $\perp$ B with wavelengths. In contrast, with the model of $a_{\rm
max}=100\,{\mu\rm{m}}$, the main source for therma dust emission at mm (VLGs)
arises from VLGs that have wrong IA due to slow internal relaxation, while the
source of dust emission at sub-mm arises from small (micron-size) grains that
have efficient IA. Consequently, the polarization pattern changes with
wavelengths due to the change in IA configuration of the emitting dust grains.
This is an interesting feature that (may) be classified as an imprint of the
grain growth process.

###  E.2 Iron Inclusions and Intensity-dependent Polarization Degree

Figure 27 shows the effect of iron inclusions on the variation of $p(\%)$ with
$I/I_{\rm max}$ at $870\,{\mu\rm{m}}$ (upper panels) and $50\,{\mu\rm{m}}$
(lower panels) for model Realistic$-$rIA$-$a10 (left column) and model
Realistic$-$wIA$-$a10 (right column). The result for the ideal model of grain
alignment is also plotted by the black dashed line for comparison. Similar as
the results with $a_{\rm max}=100\,{\mu\rm{m}}$, $p$ always decreases with
increasing $I/I_{\rm max}$ due to the decreased grain alignment efficiency in
dense regions. The polarization degree produced by SPM grains with higher
$N_{\rm cl}$ generally shows higher values thanks to the enhanced grain
alignment by the MRAT mechanism. In addition, if aligned dust grains have
wrong IA, one will obtain the "valley-polarization hole" caused by the
cancelling effect of polarized emission due to the co-existence of grains with
right and wrong IA (left column).

###  E.3 Iron Inclusions and Wavelength-dependent Polarization Degree

Figure 28 shows the effect of iron inclusions on the polarization curve
$p(\lambda)$ observed at $d_{\rm proj}=200$ au (upper panels) and $d_{\rm
proj}=12000$ au (lower panels) for model Realistic$-$rIA$-$a10 (left columns)
and Realistic$-$wIA$-$a10 (right columns). Without the presence of VLGs
(Section 9), the polarization degree is nearly independent with wavelengths
from $89\,{\mu\rm{m}}$ to 2mm in both the envelope and the central region. The
slight rise of $p(\%)$ from 2mm to $89\,{\mu\rm{m}}$ produced by SPM grains
with $N_{\rm cl}<1000$ for model Realistic$-$wIA$-$a10 is due to the increased
emission of micron-sized grains with efficient IA at short wavelengths.