11institutetext: Laboratoire AIM, Paris-Saclay, CEA/IRFU/SAp - CNRS -
Université Paris Diderot, 91191 Gif-sur-Yvette Cedex, France 22institutetext:
Department of Physics, Nagoya University, Furo-cho, Chikusa-ku, Nagoya, Aichi
464-8602, Japan
22email<EMAIL_ADDRESS>
# Is the mm/submm dust polarization a robust tracer of the magnetic field
topology in protostellar envelopes ? A model exploration

Valeska Valdivia Is the mm/submm dust polarization a robust tracer of the
magnetic field topology in protostellar envelopes ? A model explorationIs the
mm/submm dust polarization a robust tracer of the magnetic field topology in
protostellar envelopes ? A model explorationIs the mm/submm dust polarization
a robust tracer of the magnetic field topology in protostellar envelopes ? A
model explorationIs the mm/submm dust polarization a robust tracer of the
magnetic field topology in protostellar envelopes ? A model exploration
Anaëlle Maury Is the mm/submm dust polarization a robust tracer of the
magnetic field topology in protostellar envelopes ? A model explorationIs the
mm/submm dust polarization a robust tracer of the magnetic field topology in
protostellar envelopes ? A model exploration Patrick Hennebelle Is the
mm/submm dust polarization a robust tracer of the magnetic field topology in
protostellar envelopes ? A model explorationIs the mm/submm dust polarization
a robust tracer of the magnetic field topology in protostellar envelopes ? A
model exploration

(Received Month dd, yyyy; accepted Month dd, yyyy)

###### Abstract

Context. High-resolution millimeter and submillimeter (mm and submm)
polarization observations have opened a new era in the understanding of how
magnetic fields are organized in star forming regions, unveiling an intricate
interplay between the magnetic fields and the gas in protostellar cores.
However, to assess the role of the magnetic field in the process of solar-type
star formation, it is important to understand to what extent these polarized
dust emission are good tracers of the magnetic field in the youngest
protostellar objects.

Aims. In this paper, we present a thorough investigation of the fidelity and
limitations of using dust polarized emission to map the magnetic field
topologies in low-mass protostars.

Methods. To assess the importance of these effects, we performed an analysis
of magnetic field properties in 27 realizations of magnetohydrodynamics (MHD)
models following the evolution of physical properties in star-forming cores.
Assuming a uniform population of dust grains the sizes of which follow the
standard MRN size distribution, we analyzed the synthetic polarized dust
emission maps produced when these grains align with the local $B$-field
because of radiative torques (B-RATs).

Results. We find that mm and submm polarized dust emission is a robust tracer
of the magnetic field topologies in inner protostellar envelopes and is
successful at capturing the details of the magnetic field spatial distribution
down to radii $\sim 100\leavevmode\nobreak\ \mathrm{au}$. Measurements of the
line-of-sight-averaged magnetic field line orientation using the polarized
dust emission are precise to $<15^{\circ}$ (typical of the error on
polarization angles obtained with observations from large mm polarimetric
facilities such as ALMA) in about $75\%-95\%$ of the independent lines of
sight that pass through protostellar envelopes. Large discrepancies between
the integrated $B$-field mean orientation and the orientation reconstructed
from the polarized dust emission are mostly observed in (i) lines of sight
where the magnetic field is highly disorganized and (ii) those that probe
large column densities. Our analysis shows that the high opacity of the
thermal dust emission and low polarization fractions could be used to avoid
using the small fraction of measurements affected by large errors.

###### Key Words.:

polarization – magnetic fields – protostars

##  1 Introduction

The fate of the evolution of protostellar cores is a result of the interplay
between their gravity, turbulence, magnetic fields, and chemical composition
through coupled processes that control the dynamics of the gas as well as its
thermo-dynamical state. In order to better understand the star formation
process and to untangle the different contributions, we rely on observations
on one hand and theoretical models on the other. The quality of the modeling
and the quality of the interpretation of actual observations depend on how
well we can connect these two approaches, making the study of how observed
data actually traces the physico-dynamical state of the gas a key topic in
modern astrophysics.

Particularly, magnetic fields ($B$-fields) are a fundamental component of the
star formation process at all scales and evolutionary stages. They are key
agents that regulate the star formation at scales of molecular clouds, by
providing a support against gravitational collapse (Mouschovias & Ciolek 1999;
Li et al. 2014), and are thought to channel the mass accretion to form
interstellar filaments (the preferred sites of star formation). At core scales
($\lower 2.15277pt\hbox{$\;\buildrel<\over{\sim}\;$}3000$ au), if $B$-fields
are efficiently coupled with the protostellar envelope made of dust and gas,
they may also play a key role in regulating the final outcome of the star
formation process, by transporting a significant fraction of the core angular
momentum outward, hence allowing the formation of stars, setting the pristine
disk size and shaping the properties of the upcoming star and planets.

While observing the $B$-fields in star-forming environments cannot be done
directly, various techniques have been developed to probe the magnetic
properties in astrophysical structures. Among them, the most widely used
method for the dense molecular conditions consists in interpreting the
polarized emission of ISM dust grains. Indeed, several theoretical works
suggest that nonspherical dust grains spinning around their axis of maximal
inertia (the shortest axis of the dust grain) tend to align with the magnetic
fields through dissipative processes producing a polarized dust emission that
is perpendicular to the magnetic field lines (Davis & Greenstein 1951b;
Purcell 1979; Roberge 2004). The mechanisms able to align the dust grains in
such a manner include the radiative alignment torques (B-RATs, Lazarian &
Hoang 2007a; Andersson et al. 2015), a new variant of RATs including super-
paramagnetic dust grains (MRAT, Hoang & Lazarian 2016), the supersonic
mechanical grain alignment (or Gold effect, Gold 1952), the mechanical
alignment torques (MATs, Davis & Greenstein 1951a; Lazarian & Hoang 2007b;
Reissl et al. 2022), and other mechanisms (not discussed in this paper)
operating at specific conditions that can produce a polarized signal but do
not necessarily trace the magnetic field. For instance, the dust self-
scattering (Kataoka et al. 2015) caused by the scattered light of anisotropic
continuum emission by large dust grains (with sizes comparable to the observed
wavelengths) can explain the polarization levels observed in the innermost
regions of protostellar disks independently of the magnetic field, as well as
another variant of the RAT theory (k-RATs, Lazarian & Hoang 2007b; Tazaki et
al. 2017), where dust grains in regions near a strong radiation source precess
around the radiation anisotropy vector instead of precessing around the
magnetic field lines. Current observations with single-dish instruments (such
as Planck, BLASTPOL, and the JCMT POL2), as well as interferometric facilities
(such as ALMA, CARMA and SMA), are providing a wealth of polarimetric data
that, if properly interpreted, can deepen our insight into the different
processes occurring in different environments. For instance, recent
interferometric observations with the SMA allowed Galametz et al. (2020) to
highlight the potentially crucial role of the relative orientation between the
magnetic field orientation inferred from the polarized dust emission and the
rotation axis inferred from the direction of the outflows and jets in the
fragmentation and multiplicity properties of protostars. Recently, a few works
used numerical models and radiative transfer methods to characterize the dust
polarization expected from physical models of clouds, cores, and disks
(Kuffmeier et al. 2020; Brunngräber & Wolf 2021; Lam et al. 2021) and
qualitatively compare them to observations. Most conclude that, considering
the known mechanisms that generate dust polarization at millimeter
wavelengths, such polarization likely originates from magnetically aligned
grains in inner envelopes and dust scattering in disks. However, while the
ability to reconstruct large-scale magnetic field structures from the
polarized signal arising from magnetically aligned dust grains is well
established for diffuse interstellar gas, hardly any work (see Kuffmeier et
al. 2020, focusing on bridge structures) has yet been dedicated to assessing
the robustness of polarization as a tracer of the magnetic field topology in
the case of objects displaying a higher degree of complexity such as young
protostars.

Specifically, polarimetric observations aiming at probing the magnetic fields
of protostellar cores face several challenges. As the collapse proceeds, the
magnetic field lines are dragged, bent, and distorted, making the definition
of a true mean orientation difficult to obtain. The increasing gas density
subsequently rises the opacity to levels where even at millimeter (mm)
wavelengths it is hard to trace the full information along the line of sight.
On the other hand, the properties of polarized dust emission heavily rely on
the ability of dust grains to align with the magnetic fields: local
irradiation and its properties may affect how accurately the polarized dust
emission represents the mean integrated orientation of the magnetic field
lines. Understanding how well and under which conditions the magnetic field
morphologies are constrained by polarimetric observations can help us to
understand the relative importance of the physical processes shaping the star
formation process.

Table 1: Parameters of the protostellar models

ID | $\mu$ | $\beta_{\mathrm{rot}}$ | $\theta$ | $\mathcal{M}$ |  $t_{\mathrm{form}}$ [yr] |  $t_{\mathrm{0.1M_{\odot}}}$ [yr] |  $t_{\mathrm{0.2M_{\odot}}}$ [yr] |  $t_{\mathrm{0.3M_{\odot}}}$ [yr]
---|---|---|---|---|---|---|---|---
R1 | 3.33 | 0.01 | $30^{\circ}$ | 0 | 55428 | 2493 | 5991 | 13484
R2 | 3.33 | 0.04 | $30^{\circ}$ | 0 | 60328 | 2472 | 8760 | 19509
R3 | 10.0 | 0.04 | $30^{\circ}$ | 0 | 50987 | 5179 | 14816 | 30282
R4 | 6.67 | 0.04 | $30^{\circ}$ | 0 | 51206 | 3693 | 8927 | 17943
R6 | 3.33 | 0.01 | $30^{\circ}$ | 1 | 66942 | 3618 | 12321 | 30486
R7 | 3.33 | 0.04 | $90^{\circ}$ | 0 | 56006 | 2461 | 6001 | 13532
R8 | 6.67 | 0.01 | $30^{\circ}$ | 1 | 57132 | 2460 | 7331 | 16413
R9 | 10.0 | 0.04 | $90^{\circ}$ | 0 | 49636 | 5112 | 15843 | 42005
R10 | 3.33 | 0.0025 | $30^{\circ}$ | 0 | 52577 | 2828 | 7015 | 12982
111 Selected numerical models from Hennebelle et al. (2020). $\mu$ is the
mass-to-flux ratio normalized by the critical value, $\beta_{\mathrm{rot}}$ is
the ratio between rotational energy and gravitational energy, $\theta$ is the
inclination angle of the initial rotation axis with respect to the magnetic
field axis, and $\mathcal{M}$ is the initial Mach number. $t_{\mathrm{form}}$
is the time at which the sink forms in the simulation, while
$t_{\mathrm{0.1M_{\odot}}}$ , $t_{\mathrm{0.2M_{\odot}}}$ and
$t_{\mathrm{0.3M_{\odot}}}$ are the ages of the sink at which it reaches
$0.1$, $0.2$, and $0.3\leavevmode\nobreak\ \mathrm{M_{\odot}}$, respectively.

##  2 Methods

###  2.1 Non-ideal MHD models of protostellar cores

Protostellar cores are expected to display a wide variety of physical
conditions, especially in terms of masses, level of turbulence, magnetization
degree, and level of rotation and relative orientation between the magnetic
field and the rotation axis. While an exploration of the complete parameter
space is out of the scope of this study, we attempt to use realistic models of
Class 0 protostars to investigate whether polarized dust emission is a
reliable tracer of the magnetic field topology inside protostellar envelopes.
We therefore use a subset of the physically motivated numerical simulations by
Hennebelle et al. (2020), which were initially run to understand how the
interplay between the turbulence, magnetization degree, initial rotation, and
the initial angle between rotation axis and the magnetic field orientation
define the outcome of the evolution of protostellar cores. Even though most
stars are formed in multiple stellar systems (see Offner et al. 2022, and
references therein) whose dynamics is still influenced by their environment,
these isolated cores are believed to grasp the complexity of the physical
process at work at envelope scales, and hence can be trusted as realistic 3D
models of the density, temperature, gas kinematics, and magnetic field
structures in prototypical solar-type Class 0 protostars.

The numerical simulations were performed using the Ramses code (Teyssier 2002;
Fromang et al. 2006), which is able to solve the MHD equations, including the
ambipolar diffusion, to treat the gas dynamics. The thermal evolution of the
gas is computed using an equation of state (EoS) that mimics the opacity
effects on the radiative transfer in dense gas (see Eq. 1 of Hennebelle et al.
2020). A sink particle approach is used to follow the formation and evolution
of the protostar (Bate et al. 1995; Bleuler & Teyssier 2014). The code
benefits from the adaptive mesh refinement (AMR) technique to adapt the
resolution. In these simulations, the refinement criterion is the local Jeans
mass, for which we impose that the Jeans length is resolved by at least 20
cells. The structure of the data is a fully threaded ocTree.
The setup of the simulations is the same as that of Hennebelle et al. (2016)
and Masson et al. (2016), with the only difference being that it includes sink
particles and covers a wider range of initial conditions.

All the model cores have an initial mass of $1\leavevmode\nobreak\
\mathrm{M_{\odot}}$, and the magnetic field is initially uniform and aligned
with the $z$ direction. The ratio of thermal to gravitational energy is
$\alpha=0.4$ for all the simulations. The selected models cover three levels
of rotation (set by the ratio of rotational to gravitational energy
$\beta_{\mathrm{rot}}=E_{\mathrm{rot}}/E_{\mathrm{grav}}=0.0025,0.01,0.04$),
three levels of magnetization (given by the mass-to-flux ratio normalized to
the critical mass-to-flux ratio
$\mu=(M/\Phi)/(M/\Phi)_{\mathrm{crit}}=3.33,6.67,10.0$). Most of the models
have a Mach number $\mathcal{M}=0$ and an inclination between the initial
rotation axis and the initial magnetic field orientation $\theta=30^{\circ}$,
but special cases including turbulence ($\mathcal{M}=1$) and cases using an
extreme tilt of $\theta=90^{\circ}$ are also included. The run identifiers
(ID), as well as their parameters, are summarized in the first part of Table
1.

Figure 1: Sink-mass evolution as a function of the time since the formation of
the sink $t_{\mathrm{age}}$.The color identifies the normalized mass-to-flux
ratio, the line-style identifies the rotational to gravitational energy ratio,
and the symbols indicate the presence of turbulence (circles) or an extreme
inclination of the rotation with respect to the magnetic field (squares). We
keep the same identifiers in the following figures.

###  2.2 Protostar mass evolution and output selection

As the gravitational collapse proceeds, the density increases towards the
center of the simulation box. When the local density reaches a threshold of
$n_{\mathrm{thres}}=3\times 10^{13}\leavevmode\nobreak\ \mathrm{cm^{-3}}$, a
sink particle is created at the highest resolution level, at the peak of the
local clump formed by the neighboring cells whose density is higher than
$n_{\mathrm{thres}}/10$. The time at which the sink forms in the simulation,
shown in the sixth column of Table 1, indicates how fast these conditions are
reached. The sink particle mimics the formation and evolution of the
protostellar embryo, allowing us to follow its mass evolution and dynamics.
The sink formation happens faster for the less magnetized models (R9 and R3),
while the slowest formation corresponds to the R6, which in addition to being
more magnetized ($\mu=3.33$) and having a moderate rotation
($\beta_{\mathrm{rot}}=0.01$), is turbulent. The importance of the additional
support against collapse provided by the rotation itself is reflected on the
delayed formation time of the sink in models R1, R2, and R10 that, in spite of
having the same initial parameters except for the rotation, increases with
$\beta_{\mathrm{rot}}$.

The sink particle can further increase its mass by accreting gas from
neighboring cells within a radius of four times the size of the most refined
grid nearby the sink position. At each time-step and from each of its
neighboring cells, the sink accretes $10\%$ of the gas whose density is above
$n_{\mathrm{thres}}/3$. The sink mass evolution shown in Fig. 1 seems to
indicate that the most important parameter in terms of the accretion is the
mass-to-flux ratio. But contrary to the trend seen for the sink formation
time, the initial accretion proceeds faster in the more magnetized cases as a
consequence of the formation of more massive disks when the magnetic field is
weaker and the rotation stronger. All in all, an increasing $B$-field delays
the formation of the protostellar embryo at the center of the collapsing core,
but allows the early mass accretion on the protostar to proceed faster.
Although other processes may also be playing a role, we note that the two
least magnetized simulations (R9 and R3) have the earliest sink formation, and
show the lowest mass-accretion rates. The strength and geometry of the
magnetic field can also play a more episodic role during the protostellar
formation process: for example simulation R2 exhibits a reduced mass-accretion
rate when its sink particle has reached roughly $0.15\leavevmode\nobreak\
\mathrm{M_{\odot}}$. This can be explained by the development of a magnetic
cavity (the loop-like structure seen in the $x-y$ map of Fig. 2), which is a
consequence of the interchange instability (Krasnopolsky et al. 2012) often
seen in simulated disks (Matsumoto et al. 2017; Machida et al. 2014; Joos et
al. 2012; Seifried et al. 2011).

In order to analyze and compare simulations with such different evolution, we
use the fraction of the initial-mass core accreted onto the stellar embryo
(traced by the sink mass) as an indicator of the evolutionary stage of the
protostar–envelope system. For each simulation, we selected three snapshots
that correspond to the moment when $10\%$, $20\%,$ and $30\%$ of the initial
core mass has been transferred to the embryo, or equivalently, the sink mass,
$m_{\mathrm{sink}}$, has reached $0.1$, $0.2$, and $0.3\leavevmode\nobreak\
\mathrm{M_{\odot}}$, respectively, thus sampling the Class 0 stage
($M_{\mathrm{env}}>M_{\mathrm{star}}$). The time at which the sink forms in
each simulation ($t_{\mathrm{form}}$), as well as the age at which the sink
reaches the selected masses relative to its formation time are shown in the
right-hand side of Table 1. The general evolution of the sink mass for all the
simulations is shown in Fig. 1. We decided not to use the very weakly
magnetized run R5 since it fragments early in the evolution and none of the
fragments reach $0.3\leavevmode\nobreak\ \mathrm{M_{\odot}}$.

###  2.3 Magnetic field mean orientation using bidirectional statistics

As dust grains with opposite spinning directions produce an identical
polarized pattern, the polarized dust emission is insensitive to the direction
of the magnetic field, and, as a consequence, the magnetic field orientation
can be considered as a bidirectional quantity. When the data are diametrically
bimodal, the mean angle obtained by computing the mean projected sinus and
cosinus is orthogonal to the actual mean orientation due to canceling effects.
For example, toroidal magnetic field lines will produce a poloidal field
geometry when they are computed edge-on because lines in the foreground and
background have opposite directions that will cancel out.

Even though the usual integrated magnetic field in simulations of interstellar
gas seems to satisfactorily trace the mean orientation, it is not well adapted
to describing the mean orientation in regions where the magnetic field lines
are particularly tangled, as in the case of protostellar cores because of the
aforementioned canceling effects. As this effect extends into the envelope of
protostellar cores, it is necessary to find a suitable description of the
magnetic field mean orientation that is compatible with the polarized dust
emission.

Here, and in the following, we compute the orientation angle, and any other
angle, as the position angle with respect to the northern (vertical)
direction, in the sense that increases anti-clockwise, as recommended by the
IAU222https://www.iau.org/static/archives/announcements/pdf/ann16004a.pdf, in
the range $-90^{\circ}$ to $+90^{\circ}$. To compute the mean orientation of
the magnetic field, we borrow a statistical technique widely used in geography
to treat bimodal orientation data (Jones & James 1969; Mardia 1975). In order
to calculate the mean orientation of the magnetic field lines along the line
of sight regardless of the direction of the $B$-vectors, we proceed to double
the individual angles and to compute the density-weighted mean sinus and
cosinus of $2\phi_{B}$, where $\phi_{B}$ is the local magnetic field angle
projected onto the plane of the sky, as follows:

$\displaystyle\langle\sin(2\phi_{B})\rangle$ $\displaystyle=$
$\displaystyle\frac{\int_{los}\rho\sin(2\phi_{B})\,\mathop{}\\!\mathrm{d}l}{\int_{los}\rho\,\mathop{}\\!\mathrm{d}l},$
(1) $\displaystyle\langle\cos(2\phi_{B})\rangle$ $\displaystyle=$
$\displaystyle\frac{\int_{los}\rho\cos(2\phi_{B})\,\mathop{}\\!\mathrm{d}l}{\int_{los}\rho\,\mathop{}\\!\mathrm{d}l,},$
(2)

where $\rho$ is the local gas number density, and the averages are computed
along the line of sight los. The mean orientation of the magnetic field
$\langle\phi_{B}\rangle$ is then calculated as:

$\langle\phi_{B}\rangle=0.5\arctan(\langle\sin(2\phi_{B})\rangle,\langle\cos(2\phi_{B})\rangle).$
(3)

The angle-doubling procedure is consistent with analytical formulae to compute
the Stokes parameters $Q$ and $U$, such as described in Landi Degl’Innocenti &
Landolfi (2004) or Martínez González et al. (2011). To estimate the intrinsic
level of organization of the magnetic field lines along the line of sight we
define the circular variance $V_{\phi}$ as follows:

$V_{\phi}=1-\sqrt{\langle\sin(2\phi_{B})\rangle^{2}+\langle\cos(2\phi_{B})\rangle^{2}}.$
(4)

Clustered values of $\phi_{B}$ will produce circular variance values close to
0, while an increasing dispersion of angles (or in other words an increasing
level of disorganization) will produce values closer to $1$. Figure 29 shows
the relation between the circular variance $V_{\phi}$ and the typical standard
deviation $\sigma_{\phi}$ in degrees.

Table 2: Summary of the parameters used in Polaris for all the models
described in Table 1.

Quantity | Symbol | Value
---|---|---
Radiation sources
Central source luminosity | $L_{\mathrm{star}}$ | $1\leavevmode\nobreak\ \mathrm{L_{\odot}}$
Central source radius | $R_{\mathrm{star}}$ | $1\leavevmode\nobreak\ \mathrm{R_{\odot}}$
Photon packages per $\lambda$ | $N_{\mathrm{ph,star}}$ | $500000$
External radiation field | $G_{0}$ |  $1$ (Mathis et al. 1983)
Photon packages per $\lambda$ | $N_{\mathrm{ph,ISRF}}$ | $10000$
Synthetic observation parameters
Distance to observer | $D$ | $250\leavevmode\nobreak\ \mathrm{pc}$
Wavelengths | $\lambda$ | $0.8,1.3,3.0\leavevmode\nobreak\ \mathrm{mm}$
Map size | $L_{\mathrm{map}}$ | $8000\leavevmode\nobreak\ \mathrm{au}$
Resolution | $\mathop{}\\!\mathrm{d}x$ | $5\leavevmode\nobreak\ \mathrm{au}$
Dust parameters
Gas to dust ratio | $\rho_{\mathrm{g}}/\rho_{\mathrm{d}}$ | 100
Minimum grain size | $a_{\mathrm{min}}$ | $5\leavevmode\nobreak\ \mathrm{nm}$
Maximum grain size | $a_{\mathrm{max}}$ | $20\leavevmode\nobreak\ \mathrm{\mu m}$
Aspect ratio | $\epsilon$ |  $0.5$ (oblate)
Power-law index (a) | $\alpha$ | -3.5
Astronomical silicates | $f_{\mathrm{sil}}$ | $62.5\leavevmode\nobreak\ \%$
Graphite | $f_{\mathrm{graph}}$ | $37.5\leavevmode\nobreak\ \%$
Suprathermal dust (b) | $f_{\mathrm{high}-J}$ | $0.25$
333(a) Power-law index of the size distribution
$\mathop{}\\!\mathrm{d}n(a)\propto a^{\alpha}$ (Mathis et al. 1977), (b)
Fraction of dust grains spinning suprathermally.

###  2.4 Polarized dust emission with the dust radiative transfer code POLARIS

####  2.4.1 Ramses datacubes

We use the subset of the simulations presented in Hennebelle et al. (2020)
summarized in Table 1. Each Ramses output provides a datacube containing the
gas density, the gas pressure (or equivalently the gas temperature), the three
components of the velocity field, as well as the three components of the
magnetic field. These gas properties are extracted and put in a suitable
format to be post-processed with Polaris. To facilitate the propagation of
photons and avoid missed photon packages from the central source due to long
computational times, we artificially emptied a small central sphere of radius
$4\leavevmode\nobreak\ \mathrm{au}$ around the sink particle as in Valdivia et
al. (2019). This procedure (only performed for the radiative transfer in post
processing) is a pragmatic choice to allow a faster radiative transfer despite
the very optically thick layers surrounding the sink, minimizing computational
artifacts and thus improving the quality of the radiative transfer at the
scales studied in this paper ($50-4000\leavevmode\nobreak\ \mathrm{au}$). As
most of the reprocessing of the protostellar photons into infrared (IR)
photons happens at radii $4-10\leavevmode\nobreak\ \mathrm{au}$, this
procedure does not noticeably affect the temperature and irradiation spectrum
at the scales we are investigating in this study.

To describe the evolution of the simulated collapse, we provide integrated
maps of the total column density ($N$), total mean magnetic pressure
($P_{\mathrm{mag}}$), and mean magnetic field orientation ($\phi_{B}$), as
well as the associated circular variance ($V_{\phi}$) at key time-steps in
Sect. 3. We point out that the advantage of providing the mean magnetic
pressure ($P_{\mathrm{mag}}$) rather than the mean magnetic field is that it
reflects the strength of the magnetic field without canceling out the
contributions of vectors with opposite directions. The mean magnetic pressure
is computed as the density-weighted average of the local magnetic pressure,
$P_{\mathrm{mag}}$, which is given by:

$P_{\mathrm{mag}}=\frac{B^{2}}{8\pi}=\left(\frac{B}{5\leavevmode\nobreak\
\mathrm{\mu G}}\right)^{2}\times 10^{-12}\leavevmode\nobreak\
\mathrm{erg\leavevmode\nobreak\ cm^{-3}},$ (5)

where $B=\sqrt{B_{x}^{2}+B_{y}^{2}+B_{z}^{2}}$ is the local magnetic field
strength in $\mathrm{\mathrm{\SIUnitSymbolMicro}G}$.

####  2.4.2 Radiation sources and dust properties

As the protostellar luminosity as well as the efficiency of the conversion
from gravitational energy into accretion luminosity are not well constrained
and can vary significantly during the Class 0 stage (Offner & McKee 2011;
Cheng et al. 2022; Kuffmeier et al. 2018; Jensen & Haugbølle 2018), we adopt a
simple description. For all the simulations, and regardless of the protostar
mass, we model the central radiation source as a blackbody of
$1\leavevmode\nobreak\ \mathrm{L}_{\odot}$ of radius of $1\leavevmode\nobreak\
\mathrm{R}_{\odot}$ located at the sink position. For a more detailed study
investigating the dependency of dust polarized emission with varying accretion
luminosity, and associated radiation field penetrating the envelopes of Class
0 protostars, see for example Le Gouellec et al. (2020) and Le Gouellec et al.
in prep. We include an external isotropic radiation field of strength
$G_{0}=1$ following the Mathis et al. (1983) description.
The dust properties are assumed uniform throughout the simulation box and are
the same as those in Valdivia et al. (2019), namely a gas-to-dust ratio of
$100$, with a composition of $62.5$ per cent astronomical silicates and $37.5$
per cent graphite grains, which matches the Galactic extinction curve (Mathis
et al. 1977). The dust grains are assumed to be oblate with an aspect ratio of
$0.5$ (Hildebrand & Dragovan 1995) and following a power-law size distribution
(Li & Draine 2001) of

$\mathop{}\\!\mathrm{d}n(a)\propto a^{-3.5}\mathop{}\\!\mathrm{d}a,$ (6)

where $a$ is the effective dust radius of a spherical grain of equivalent
volume, and $n(a)$ is the number of dust grains of effective radius $a$. The
lower and upper cutoff of the size distribution are set as
$a_{\mathrm{min}}=5\leavevmode\nobreak\ \mathrm{nm}$ and
$a_{\mathrm{max}}=20\leavevmode\nobreak\
\mathrm{\mathrm{\SIUnitSymbolMicro}m}$, respectively. The choice of
$a_{\mathrm{max}}$ stems from a compromise between our previous study
(Valdivia et al. 2019), suggesting that the polarization fractions observed
towards protostellar cores can only be approached when including a population
of large grains, and the need to avoid the regime of very large grains where
the DDSCAT (Draine & Flatau 2000) introduces artifacts on the absorption
coefficients (see Reissl et al. 2016).All the parameters are summarized in
Table 2 and are the same used for all the selected snapshots.

####  2.4.3 Temperatures and dust alignment with Polaris

We post-processed each one of the snapshots described in Sect. 2.1 with the
Polaris code (Reissl et al. 2016). The photon propagation is based on a Monte
Carlo scheme and, along with the local values of the density and the magnetic
field, it is used to compute the dust temperature and the grain alignment
efficiency on a cell-by-cell basis. As the gas temperature computed by Ramses
using an EoS is the result of an approximation and not the result of an actual
radiative transfer, we recompute the gas temperature as well.
In addition to the gas and dust heating, the radiation field induces a
radiation torque on dust grains, spinning them up to suprathermal rotation
velocities, at which they can align with the magnetic field (Barnett 1915;
Dolginov & Mitrofanov 1976; Lazarian & Hoang 2007a; Hoang & Lazarian 2014).
The grain alignment efficiency is computed according to the B-RATs theory as
implemented in Polaris (Reissl et al. 2016). In principle, the fraction of
dust grains spinning at suprathermal velocities ($f_{\mathrm{high}-J}$, which
is the stable population of aligned grains) depends on the dust grain
properties such as the specific shape, composition, and magnetic properties.
As this quantity cannot be calculated exactly, it is set to
$f_{\mathrm{high}-J}=0.25$. This rather conservative choice is able to produce
a detectable amount of polarization. Recently, Herranen et al. (2021) found
that the alignment efficiency might be higher than predicted, which in the
case of aggregates can reach $f_{\mathrm{high}-J}\sim 0.35-1$. Also, the
presence of iron inclusions can increase this quantity to nearly unity (Hoang
& Lazarian 2016). Additionally, we take into account the fact that, due to
internal thermal fluctuations, nonspherical dust grains are not expected to be
perfectly aligned, but precessing around their rotation axes; we do this by
enabling the imperfect internal alignment mechanism (II, Barnett 1915; Purcell
1979), as implemented in Polaris. To include the effect of random gas
collisions, we include the imperfect Davis-Greenstein alignment (IDG, or
paramagnetic relaxation Davis & Greenstein 1951b; Jones & Spitzer 1967;
Purcell 1979; Spitzer & McGlynn 1979), which despite being rather inefficient
(Reissl et al. 2017) can contribute to aligning a part of the small grains.
Silicate grains and graphite grains have very different paramagnetic
properties, which differ by roughly six orders of magnitude (Draine 1996;
Hoang & Lazarian 2014). This results in a very randomized distribution of
graphite grain orientations, preventing them from contributing to the
polarized emission. Finally, we do not include the effects of the self-
scattering (Kataoka et al. 2015), because it operates in the innermost regions
of the disk (typically within a radius of $\sim 10\leavevmode\nobreak\
\mathrm{au}$), which are excluded from our analysis, or the k-RAT mechanism
(alignment with the radiation anisotropy), which requires radiation conditions
much stronger than those considered in this work to be efficient (Lazarian &
Hoang 2007b; Tazaki et al. 2017).

####  2.4.4 Synthetic observations of the Stokes maps with Polaris

For each one of the selected snapshots, we set our model at a distance of
$250\leavevmode\nobreak\ \mathrm{pc}$ and observe a squared region of
$8000\leavevmode\nobreak\ \mathrm{au}$ on a side centered at the position of
the protostellar embryo (materialized by the sink particle) at a resolution of
$5\leavevmode\nobreak\ \mathrm{au}$ per pixel. Each snapshot is observed
synthetically along each one of the main axes $x$, $y$, and $z$ ($y$–$z$,
$x$–$z$, and $x$–$y$ maps respectively) at three wavelengths: $0.8$, $1.3,$
and $3.0\leavevmode\nobreak\ \mathrm{mm}$.

For each synthetic observation, Polaris provides perfectly observed maps of
the four components of the Stokes vector $S=(I,Q,U,V)^{T}$, where $I$ is the
total dust emission, $Q$ and $U$ are the linearly polarized dust emission
components, and $V$ is the circularly polarized emission.

Finally, we compute the total linearly polarized intensity maps ($PI$) and the
polarization fraction $p_{\mathrm{frac}}$ from the Stokes maps as follows:

$\displaystyle PI$ $\displaystyle=$ $\displaystyle\sqrt{Q^{2}+U^{2}}$ (7)
$\displaystyle p_{\mathrm{frac}}$ $\displaystyle=$ $\displaystyle PI/I.$ (8)

At far-infrared (FIR) to millimiter wavelengths, the magnetic field
orientation inferred from the polarized emission ($\phi_{\mathrm{pol}}$) is
obtained by rotating the polarization vector $\psi=\frac{1}{2}\arctan(U,Q)$ by
$90^{\circ}$.

Figure 2: Column density maps of all the selected outputs described in Table
1. Each block corresponds to a given simulation (the identifier is indicated
in the lower left corner of the block). Each block displays column density
maps integrated along the $x$, $y,$ and $z$ axes (from top to bottom: $y$-$z$,
$x$-$z$, and $x$-$y$ maps, respectively), and for a snapshot at
$m_{\mathrm{sink}}=0.1,0.2,$ and $0.3\leavevmode\nobreak\
\mathrm{M_{\odot}}$(from left to right). The projection and evolutionary step
are shown at the top of each individual plot. Higher resolution figures for
individual simulations are available online. Figure 3: Mean magnetic pressure
$P_{\mathrm{mag}}$ maps and the orientation of the magnetic field lines along
each projection for all the simulations. Each block corresponds to a given
simulation and the plots are organized using the same structure as in Fig. 2.
Higher resolution figures for individual simulations are available online.
Figure 4: Circular variance $V_{\phi}$ of the $B$-field lines in the Ramses
data cubes (See Sect. 2.3). The figure is organized using the same structure
as in Figs. 2 and 3. The simulation identifier is given in the bottom left
corner of the block, while the projection and evolutionary step are shown at
the top of each individual plot. $0.1\leavevmode\nobreak\
\mathrm{M_{\odot}}$$0.2\leavevmode\nobreak\
\mathrm{M_{\odot}}$$0.3\leavevmode\nobreak\ \mathrm{M_{\odot}}$ Figure 5:
Radial mean column density profiles. From left to right: profiles at
$m_{\mathrm{sink}}=0.1,0.2,$ and $0.3\leavevmode\nobreak\ \mathrm{M_{\odot}}$.
$0.1\leavevmode\nobreak\ \mathrm{M_{\odot}}$$0.2\leavevmode\nobreak\
\mathrm{M_{\odot}}$$0.3\leavevmode\nobreak\ \mathrm{M_{\odot}}$ Figure 6:
Radial mean profiles of the circular variance $V_{\phi}$. From left to right:
profiles at $m_{\mathrm{sink}}=0.1,0.2,$ and $0.3\leavevmode\nobreak\
\mathrm{M_{\odot}}$.

##  3 From 3D models to 2D maps: evolution of the integrated properties in the
models

The general evolution of the numerical simulations and the dependence on the
initial conditions are detailed in Hennebelle et al. (2020). In the present
paper, we focus on the description of the properties relevant to the polarized
dust emission, particularly the protostar evolution (represented by the sink
particle), the column densities, and the integrated mean orientation of the
magnetic field lines as well as its intrinsic degree of disorganization
(quantified by the variance of the angle orientation along the line of sight).

In this section we present integrated maps along the $x$, $y,$ and $z$ axes
for the selected simulations and snapshots described in Table 1. We chose to
present these maps at the same resolution and scale as those of the
synthetically observed ones in order to allow a direct comparison between
both444Higher resolution figures for individual panels of Figs. 2, 3, 8, 16,
and 17 are available online.. In Fig. 2, we present the total column density
maps, which correspond to the total number of hydrogen atoms integrated along
the line of sight. Figure 3 shows the density-weighted mean magnetic pressure
($P_{\mathrm{mag}}$, as described in Sect. 2.4.1) with overlaid mean
orientation ($\phi_{B}$, computed following Eq. 3), aiming at providing a
description of the magnetic field strength and mean orientation. Figure 4
shows the circular variance ($V_{\phi}$, as defined in Eq. 4), which
quantifies the level of disorganization of the magnetic field lines along the
line of sight. Finally, we summarize these results by computing the average
radial profiles of the three axes of $N$ and $V_{\phi}$ in Figs. 5 and 6,
respectively.

###  3.1 Integrated 2D maps

The column densities shown in Fig. 2 display a wide variety of gas
distribution features including outflows, spiral structures, loops, and more
complex density distributions depending on the initial physical conditions.
More precisely, the simulations including turbulence display very asymmetrical
features, while those of the remaining simulations appear more uniform. Figure
3 shows how the magnetic field lines are dragged after the initial
gravitational collapse, forming the characteristic hourglass shape, but the
further evolution of the gas distribution and the final morphology of the
magnetic field are controlled by the interplay of the specific physical
conditions. As in the case of the evolution of the sink particle, the most
important parameter is the mass-to-flux ratio ($\mu$). Figure 3, and more
particularly simulations R2, R3, and R4, shows that as the magnetic field
lines are more prone to bending and distortion in the less magnetized cases
($\mu=10$), the magnetic field lines evolve into a helicoidal morphology in
the outflow, which develops more easily. As the magnetization level increases,
the outflow remains more confined and the helicoidal feature is less evident.
The effect of the magnetization level can also be seen in simulations R7 and
R9. Both simulations have the same rotation parameters ($\beta=0.04$ and
$\theta=90^{\circ}$), but the outcome of the evolution is very different.
Simulation R9 ($\mu=10$) displays a well-developed spiral structure seen in
the $y-z$ projection (nearly face-on projection) of the column density map
shown in Fig. 2 as well as in the projected magnetic field lines (streamplot)
of Fig. 3. On the other hand, in simulation R7 the stronger magnetic tension
inhibits the development of high-contrast density waves seen as spiral arms,
and the associated warped $B$-field lines. Simulations R6 and R8 also share
the same rotation ($\beta=0.04$ and $\theta=30^{\circ}$) and turbulent
($\mathcal{M}=1$) parameters, but in this case it seems that the turbulence
dominates the evolution of the gas dynamics as well as the evolution of the
magnetic field lines. Nevertheless, the circular variance (Fig. 4) seems to be
higher and more extended for R6.

The circular variance maps (Fig. 4) describe the level of organization of the
projected orientation of the magnetic field lines onto the plane of the sky,
while traveling along the line of sight. As all the magnetic field lines are
initially set parallel to the $z$ direction, all the $V_{\phi}$ maps are
initially uniform and equal to zero. As the gas collapses inwards, dragging
the magnetic field lines, the projected orientations of the magnetic field are
no longer parallel and the $V_{\phi}$ starts to increase. At larger scales
($r>4000\leavevmode\nobreak\ \mathrm{au}$, not shown in the figure), the
magnetic field orientation remains extremely well organized and the hourglass
pattern is barely affected by the rotation or the outflow.

At smaller scales, as those shown in the figure, the magnetic field remains
relatively ordered along the lines of sight exhibiting yellow to light blue
colors, which corresponds to a typical angular dispersion of the orientation
of the magnetic field of $\sigma_{\phi}=20^{\circ}-30^{\circ}$ (see Appendix
C), while regions with purple to red have larger angular dispersions
($\sigma_{\phi}>30^{\circ}$). These regions of larger $V_{\phi}$ are
ubiquitously seen towards the midplanes of the nearly edge-on projections, as
well as in regions where the early hourglass structure of the magnetic field
is perturbed either by the development of the outflow (very prominent in
simulations R2 and R3, for example) or by the turbulent gas motions
(simulations R6 and R7).

We note that the projection that displays the most organized $B$-field is $z$
($x-y$ maps on the bottom row of each block,) which is the direction of the
initial magnetic field lines ($V_{\phi}$ remains small even for the cases R7
and R9 where the rotation axis is tilted by $90^{\circ}$ with respect to the
$z$ direction).

###  3.2 Radial profiles

To summarize the global behavior of the simulations and to compare the time
evolution of both the collapse and the stress it produces on the magnetic
field topology, we compare the mean radial profiles of the total column
density and of the circular variance of the $B$-field orientation (Figs. 5 and
6, respectively). For each evolutionary step and simulation, we stack the
integrated maps of the three projections, and compute the mean value of the
variable in radial logarithmic bins. From the evolution of the column density
radial distribution (Fig. 5), we clearly distinguish the usual development of
a high-density inner region peaking around the protostar (column densities
higher than $\sim 2\times 10^{24}\leavevmode\nobreak\ \mathrm{cm^{-2}}$ within
a radius of $r\lesssim 60\leavevmode\nobreak\ \mathrm{au}$). We also notice a
flattening of the density profile at intermediate radii with time, which is
due to the development of a flattened envelope sometimes also referred to as
the ”pseudo-disk” (Galli & Shu 1993; Li & Shu 1996; Hennebelle & Ciardi 2009),
while the outer layers of the protostellar envelopes
($r>1000\leavevmode\nobreak\ \mathrm{au}$) remain mostly unaffected. All the
simulations, except for the simulation with a lower degree of magnetization,
present a comparable global evolution. Simulation R9, and to a lesser extent
simulation R3 (low magnetization cases and $\theta=90^{\circ}$ and
$\theta=30^{\circ}$, respectively), display a wider overdensity region due to
the development of the spiral structures and a thicker equatorial region with
respect to the rotation axis.

The behavior of the mean circular variance radial profiles of Fig. 6 seems
quite different from that of the mean column density radial profile. The
$V_{\phi}$ profiles have slightly lower values in the inner region, remaining
more uniform with typical values of $V_{\phi}\sim 0.2-0.6$ at intermediate
radii, while in the outer region they decrease to unperturbed values. The
regions displaying the largest variance are located at intermediate distances,
following the development of the outflow. In these regions, the development of
a helicoidal magnetic field, along with the projection effects, induces a
higher circular variance in the cavity walls of the outflow, as remarkably
shown in the case of R3. However, even at advanced Class 0 evolution
($0.3\leavevmode\nobreak\ \mathrm{M_{\odot}}$ in the central embryo), the
stress applied by the gravitational collapse and outflowing motions to the
$B$-field only result in rather low $V_{\phi}$ with values $<0.6$ (equivalent
to typical dispersion of the $B$-field angle along the line of sight $\lower
2.15277pt\hbox{$\;\buildrel<\over{\sim}\;$}30^{\circ}$). A particular case is
the simulation R9, which displays a decreasing $V_{\phi}$ with time in the
inner $\sim 500\leavevmode\nobreak\ \mathrm{au}$ with time due to the
development of a more organized spiral pattern guided by the rotation due to a
lower magnetic tension, and that remains unperturbed due to the lack of
outflow. There is not a clear dependence between the mean $V_{\phi}$ and the
degree of magnetization, but the presence of turbulence seems to produce a
more uniform and more stable level of disorganization of the magnetic field
lines.

$0.1\leavevmode\nobreak\ \mathrm{M_{\odot}}$$0.2\leavevmode\nobreak\
\mathrm{M_{\odot}}$$0.3\leavevmode\nobreak\ \mathrm{M_{\odot}}$ Figure 7: Mean
circular variance $V_{\phi}$ per logarithmic bin of column density. From left
to right: results at $m_{\mathrm{sink}}=0.1,0.2,$ and $0.3\leavevmode\nobreak\
\mathrm{M_{\odot}}$.

To better understand how the circular variance is related to the column
density, we compute the mean value of $V_{\phi}$ in logarithmic bins of total
column density. Figure 7 shows these results for the stacked data. We notice
that initially only the intermediate column densities present high levels of
disorganization of the magnetic field lines. These regions are related to the
equatorial region, where the column densities increase due to the collapse and
the further accretion is guided along the magnetic field lines, at the same
time as the lines are distorted by the gas rotation. At later times, as the
gas becomes denser, regions with high $V_{\phi}$ develop at higher column
densities for the simulations with a higher degree of magnetization (purple
lines indicating $\mu=3.33$), for all the cases (with turbulence, or rotation
orthogonal to the magnetic field). It is interesting to notice that
simulations R3 and R9 (yellow lines), despite having the same level of
magnetization and rotation, differing only in the initial rotation angle,
display a very different behavior at intermediate column densities,
particularly at around $N\sim 10^{24}\leavevmode\nobreak\ \mathrm{cm^{-2}}$,
where the simulation R3 have values of $V_{\phi}>0.6$, while simulation R9
barely surpasses $0.2$.

Figure 8: Results from the synthetic observation at
$\lambda=1.3\leavevmode\nobreak\ \mathrm{mm}$ showing the polarization
fraction (background image), the polarized intensity contours, and the
inferred magnetic field orientation vectors. The contours correspond to the
linearly polarized intensity normalized to the peak value at levels $10^{-2}$
(black), $10^{-3}$ (purple), $10^{-4}$ (red), and $10^{-5}$ (yellow). Higher
resolution figures for individual simulations are available online.
---
$0.1\leavevmode\nobreak\ \mathrm{M_{\odot}}$$0.2\leavevmode\nobreak\
\mathrm{M_{\odot}}$$0.3\leavevmode\nobreak\ \mathrm{M_{\odot}}$ Figure 9:
Radial profiles of the synthetic observation at
$\lambda=1.3\leavevmode\nobreak\ \mathrm{mm}$. From top to bottom: Total dust
emission (Stokes $I$), total linearly polarized dust emission ($PI$), and
polarization fraction ($p_{\mathrm{frac}}$). From left to right: Results at
$m_{\mathrm{sink}}=0.1,0.2,$ and $0.3\leavevmode\nobreak\ \mathrm{M_{\odot}}$.

##  4 Synthetic maps of the polarized dust emission

We produced perfectly integrated synthetic observations of the polarized dust
emission for the central $8000\leavevmode\nobreak\ \mathrm{au}$, observed at a
distance of $250\leavevmode\nobreak\ \mathrm{pc}$, and at a resolution of
$5\leavevmode\nobreak\ \mathrm{au}$ per pixel (see Table 2). Fig. 8 presents
synthetic maps of the polarization fraction observed at
$1.3\leavevmode\nobreak\ \mathrm{mm}$ (background color maps), with overlaid
polarized intensity contours (normalized to the peak value) and the magnetic
field vectors (inferred from the polarized emission as described in Sect.
2.4). To summarize and simplify the comparison between all the simulations, we
computed the radial profiles, stacking the data for the integrated quantities
along the three projections. We provide the radial profiles of the Stokes $I$,
the intensity of the total linearly polarized dust thermal emission $PI,$ and
the polarization fraction $p_{\mathrm{frac}}$ in Fig. 9. We present similar
maps, but for $\lambda=0.8$ and $3.0\leavevmode\nobreak\ \mathrm{mm}$ (Figs.
16 and 17), as well as similar radial profiles (Figs. 18 and 19) in the
Appendix.
The total dust emission maps (Stokes $I$ maps, not shown here) peak towards
the central region neighboring the protostar, closely following the column
density as expected. This is also seen by comparing the radial profiles of the
Stokes $I$ maps (top rows of Fig. 9 and 18 and 19) to the total column density
maps (Fig. 5), which show good agreement between these two quantities, but
show slight departures at high column densities due to the large associated
opacities.
The total linearly polarized intensity (contour lines of Figs. 8, 16 and 17,
and middle row panels of Figs. 9, 18 and 19) shows higher values toward the
central region (due to higher densities and temperatures), but the associated
polarization fractions are the lowest in these regions (see the bottom panels
of Figs. 9, 18 and 19). The highest polarization fractions are found in the
external regions, typically reaching $5\%$ to $10\%$, and gradually decrease
in a nonmonotonical manner toward the central regions. More precisely, Fig. 8
shows that $p_{\mathrm{frac}}$ drops more strongly in the midplane and the
outflow regions, particularly near the cavity walls. Even though the outflow
regions show lower polarization levels, the simulation R3 shows a rather high
$p_{\mathrm{frac}}$ in the central part of the outflow region. This is likely
due to the fact that the outflow cavity is well developed, allowing the
transversal component of the helicoidal magnetic field to dominate the
polarization pattern. A qualitative comparison of Fig. 6 and Fig. 8 suggests
that the polarization fraction drops in the lines of sight where the $B$-field
happens to have a large variance ($V_{\phi}$). In the following section, we
quantify and discuss this behavior, and its consequences as to our ability to
trace the $B$-field in low-mass protostellar envelopes with dust polarization
observations.

It is interesting to note that, even though the simulations analyzed in this
paper correspond to an initial core mass of only $1\leavevmode\nobreak\
\mathrm{M_{\mathrm{\odot}}}$, the less magnetized cases, and particularly
simulation R3, produce a polarization pattern remarkably similar to the one
observed by Hull et al. (2017) for Ser-emb 8, particularly in the $x-z$ maps
(see Fig. 16) at the time when the protostellar mass reached
$0.3\leavevmode\nobreak\ \mathrm{M_{\mathrm{\odot}}}$, consistent with the
scenario proposed by Hull et al. (2017), but without invoking the presence of
initial turbulence.

##  5 Do polarized dust emission maps capture the underlying topology of
protostellar magnetic fields?

###  5.1 Fidelity of the inferred magnetic field orientation

In this work, we aim to assess whether or not (sub)mm polarized dust emission
maps of protostellar envelopes —which are routinely used to discuss the
affects of the magnetic field during the collapse of protostellar cores into a
star–disk system— are indeed robust tracers of the magnetic field lines
threading those cores. To do so, we compared the two quantities introduced in
Sect. 2: the magnetic field orientation reconstructed from the synthetic
observation of polarized dust emission ($\phi_{\mathrm{pol}}$, Sect.2.4.4) and
the mean orientation of the magnetic field in the model integrated along the
line of sight and computed using the circular statistics ($\phi_{B}$, Eq. 3).
We compute the difference between model and observations as the following
difference: $\Delta\phi=\phi_{\mathrm{pol}}-\phi_{B}$, and the discrepancy as
the absolute value of the difference, $|\Delta\phi|$.

####  5.1.1 Dust polarized emission: an overall good tracer of $B$-fields in
protostellar environments

Figure 10: $B$-field recovery rate for all the synthetic observations at
$\lambda=0.8\leavevmode\nobreak\ \mathrm{mm}$. Individual points correspond to
those of Table 5, while the cross symbols and vertical bars correspond to the
the averaged values and their standard variation for the two cases defined in
Table 3, respectively: moderate ($\mu=3.33$) and weak ($\mu>6$) magnetic
fields. The values are separated in three spatial scales (indicated at the
bottom of the figure), and slightly shifted for the sake of clarity for the
three evolutionary stages given by the protostellar mass indicated at the top
of the figure. Table 3: Fidelity of $B$-field mapping from polarized dust
emission at $\lambda=0.8\leavevmode\nobreak\ \mathrm{mm}$ in protostellar
envelopes at different scales.

$m_{\mathrm{sink}}$ |  $B$-field recovery rates and $\sigma$ $[\%]$
---|---
$[\mathrm{M_{\odot}}]$ | All scales |  $30-100$ au |  $100-500$ au |  $500-1500$ au
$\mu=3.33$
0.1 | 96.9 3.1 | 70.1 14.8 | 82.3 18.8 | 92.9 9.5
0.2 | 95.6 4.2 | 74.3 11.0 | 75.8 18.5 | 89.2 12.2
0.3 | 93.3 5.3 | 73.0 10.3 | 78.2 12.7 | 85.9 11.7
$\mu=6.67,10.$
0.1 | 97.3 2.1 | 83.7 7.6 | 86.4 9.8 | 95.4 4.5
0.2 | 95.7 2.7 | 79.3 9.0 | 90.4 6.4 | 89.9 7.3
0.3 | 90.7 6.3 | 75.2 13.2 | 88.9 8.6 | 89.5 6.3
555 The recovery rate corresponds to the fraction of lines of sight where
measurement of the $B$-field orientation from polarized dust emission recovers
the integrated $B$-field orientation within an error of $<15^{\circ}$ in the
modeled Class 0 envelopes. The percentages shown in this table correspond to
the averaged recovery rate, while the standard deviation of each measurement
is given in a smaller font next to each mean value. The results are separated
in two cases: moderate magnetic field ($\mu=3.33$, simulations R1, R2, R6, R7,
and R10) and weak magnetic field ($\mu=6.67,10.0$, simulations R3, R4, R8, and
R9). Percentages have been averaged over all the corresponding simulations and
projections at a given evolutionary step (given by the protostar mass
$m_{\mathrm{sink}}$). The scales correspond to the whole map including all the
scales (second column), and in concentric rings containing all the lines of
sight with radial distances to the protostar in the interval indicated in the
header of the table. All the individual recovery rates are given in Table 5.

Table 4: Fidelity of $B$-field mapping from polarized dust emission at
$\lambda=0.8\leavevmode\nobreak\ \mathrm{mm}$ in protostellar envelopes as a
function of the polarized intensity.

$m_{\mathrm{sink}}$ |  $B$-field recovery rates and $\sigma$ $[\%]$
---|---
$[\mathrm{M_{\odot}}]$ | $PI_{\mathrm{norm}}\geq 10^{-2}$ | $PI_{\mathrm{norm}}\geq 10^{-3}$ | $PI_{\mathrm{norm}}\geq 10^{-3}$
$\mu=3.33$
0.1 | 84.4 17.7 | 93.3 8.8 | 98.4 2.0
0.2 | 83.7 14.2 | 92.7 8.4 | 97.7 2.5
0.3 | 86.6 8.9 | 92.7 6.3 | 96.5 3.3
$\mu=6.67,10.$
0.1 | 92.9 4.3 | 95.6 3.7 | 99.0 1.0
0.2 | 93.6 4.4 | 94.0 3.8 | 98.0 1.6
0.3 | 92.2 6.0 | 93.8 4.7 | 95.4 3.6
666 The recovery rates are defined as in Table 3 for an error of
$<15^{\circ}$. The recovery rates shown in this table correspond to the mean
values found for regions with a total linear polarized intensity normalized to
the peak value ($PI_{\mathrm{norm}}=PI/PI_{\mathrm{peak}}$) higher than
$10^{-2}$, $10^{-3}$, and $10^{-4}$, respectively. These values correspond to
the first three contours of Fig. 16. As in Table 3, the corresponding standard
deviations are given in a smaller font.

We built difference maps between the model magnetic field line orientations
and the reconstructed magnetic field orientations from polarized dust emission
observed at $\lambda=0.8,1.3$, and $3.0\leavevmode\nobreak\ \mathrm{mm}$ (see
Figs. 20, 21 and 22, respectively). We note that there is generally good
agreement between the inferred ($\phi_{\mathrm{pol}}$) and the actual
($\phi_{B}$) mean orientations, with localized regions displaying large
differences. From these maps, we can quantify —in a statistical way— the
ability of the polarized dust emission to recover the true magnetic field
topology. We define the recovery rate as the fraction of sight lines in a map
for which the magnetic field orientation is retrieved within a given error, or
in other words, have a discrepancy $|\Delta\phi|$ lower than a given
threshold. This recovery rate is found to be similar at the three wavelengths
for each map, slightly increasing with the wavelength. At
$\lambda=0.8\leavevmode\nobreak\ \mathrm{mm}$, which is the frequency that
displays the lowest recovery rates, the magnetic field is typically recovered
within an error of less than $15^{\circ}$ (typical of most observations using
dust polarized emission as a probe for magnetic topologies) for more than
$95\%$ of the lines of sight over the full maps, with only some cases with no
clear dependence on the magnetization level, rotation, or inclination,
recovering less than $90\%$. These cases correspond to late evolutionary steps
and projections along the $y$-axis ($x-z$ maps). The worst recovery rate
($79.7\%$) is obtained for simulation R3,
$m_{\mathrm{sink}}=0.3\leavevmode\nobreak\ \mathrm{M_{\odot}}$ and projected
along the $y$-axis. For this map, most of the lines of sight for which the
orientations were not recovered within an error of $15^{\circ}$ are associated
to higher circular variances (the recovery rate for the lines of sight with
$V_{\phi}<0.5$ is $99.7\%$). Table 3 summarizes the mean recovery rates for
the synthetic maps at $\lambda=0.8\leavevmode\nobreak\ \mathrm{mm}$ at
different evolutionary stages and at different scales, averaged for all the
projections and simulations. The detailed results are given in Table 5.
Additionally, we provide a synthesis of these tables in Fig. 10. These tables
show that the magnetic field orientation reconstructed from polarized dust
emission observations is remarkably good at large scales and that even at
intermediate scales the recovery rate is typically higher than $70\%$ (see
also Fig. 10). Table 7 and Table 8 show even better results at longer
wavelengths ($\lambda=1.3\leavevmode\nobreak\ \mathrm{mm}$ and
$\lambda=3.0\leavevmode\nobreak\ \mathrm{mm}$, respectively), consistent with
the findings of Kuffmeier et al. (2020). To check whether or not the polarized
emission that is susceptible to being observed recovers the magnetic field
topology, we computed the recovery rates for different levels of total
linearly polarized dust emission normalized to the peak
($PI_{\mathrm{norm}}=PI/PI_{\mathrm{peak}}$) at
$\lambda=0.8\leavevmode\nobreak\ \mathrm{mm}$. We use the same values of the
contour levels of Fig. 16. We computed the mean recovery rates over all the
simulations and projections separated in the same two cases as in Table 3. We
find that even for the first contour, which corresponds to the most intense
emission and which is therefore the most likely to be detected, the recovery
rates are found to be higher than $83\%$ for the moderate magnetic field
simulations, and higher than $92\%$ for the less magnetized simulations. As
the contour level decreases, the recovery rate increases, consistently with
the results of Table 3 seen for larger scales.

Finally, note that the models analyzed here implement nonideal MHD processes
under ionization conditions that may not be prototypical at the smallest
scales around protostars: for example, Cabedo et al. (2022) estimated cosmic-
ray ionization rates in B335 from observations probing the $\sim 100-500$ au
scales which are up to two orders of magnitude higher than the predicted
ionization generated by external cosmic rays alone. If the ionization was
indeed to be much stronger at small envelope radii due to locally accelerated
cosmic rays, that would make the coupling of B fields to the local gas much
more efficient and models that implement both ideal and nonideal MHD with
realistic local ionization conditions would need to be explored in the future.

####  5.1.2 Physical and observational effects responsible for poor
measurements in disks and outflows

We now present our investigation of the lines of sight associated to large
errors in the reconstructed magnetic field maps, where our goal is to identify
the local physical properties responsible for the mismatch. As the local gas
density and its temperature both scale with radii in protostellar envelopes,
we first looked for the distribution of inaccurate measurements as a function
of distance to the central protostellar embryo. Figure 12 shows the percentage
of lines of sight that show a discrepancy (expressed in terms of the absolute
value of the difference $|\Delta\phi|$) of higher or equal to $30^{\circ}$ as
a function of radius (using $100$ logarithmic bins). This figure shows
distinctive radial distances at which the points with large discrepancies
represent an important contribution. These regions develop outwards in time,
diluting this fraction over a larger volume. The more magnetized simulations
($\mu=3.33$) display percentages of higher than $10\%$ over larger regions.
Simulations with $\theta=90^{\circ}$ start to show regions of high discrepancy
at a much larger radius. Once the sinks in the simulations have reached
$0.3\leavevmode\nobreak\ \mathrm{M_{\odot}}$, most of the simulations show a
region of high discrepancy within a distance of $30-100\leavevmode\nobreak\
\mathrm{au}$ with respect to the sink position. Even though Figs. 12 and 6 do
not seem to display the same behavior, these regions of higher discrepancy
seem to be connected to larger levels of disorganization of the field line
orientations along the line of sight.

To better understand the role of the intrinsic degree of organization of the
magnetic field lines and the reliability of the orientation inferred from dust
polarization observations, we present 2D histograms showing the absolute value
of the difference $|\Delta\phi|$ as a function of the circular variance
$V_{\phi}$ at $\lambda=1.3\leavevmode\nobreak\ \mathrm{mm}$ in Figs. 13 and
14. Analogous results for $\lambda=0.8\leavevmode\nobreak\ \mathrm{mm}$ and
$\lambda=3.0\leavevmode\nobreak\ \mathrm{mm}$ are shown in the Appendix in
Figs. 23 and 24, and Figs. 25 and 26, respectively. All the results correspond
to the stacked values for the three projections ($y-z$, $x-z$ and $x-y$) of
Fig. 21 (Fig. 20 and Fig. 22 for $\lambda=0.8\leavevmode\nobreak\ \mathrm{mm}$
and $\lambda=3.0\leavevmode\nobreak\ \mathrm{mm}$). The contour lines
corresponds to the 2D histogram contour levels of the distribution of the
lines of sight at $10^{5}$ (yellow), $10^{4}$ (red), $10^{3}$ (purple), and
$10^{2}$ (black) counts, while the color-coded background displays the mean
column density for a given bin of $V_{\phi}$ and $|\Delta\phi|$. We present in
Fig. 13 the distribution for the more standard conditions (no initial
turbulence and $\theta=30^{\circ}$, namely simulations R1, R2, R3, R4 and
R10), while Fig. 14 shows the same results but for the simulations including
turbulence (R6, R8), or with a rotation axis perpendicular to the initial
magnetic field orientation (R7, R9). In both cases, the results have been
organized by blocks of decreasing level of magnetization (increasing $\mu$),
and increasing level of rotation (increasing values of
$\beta_{\mathrm{rot}}$). In all the simulations, most of the points in the map
(contours higher or equal to $10^{4}$ counts) are characterized by low levels
of $|\Delta\phi|$ (typically $\leq 10^{\circ}$) and low to moderate values of
$V_{\phi}$, with low values of the mean column density. These points
corresponds to regions that, other than the hourglass, have not been
significantly disturbed by the infall, the turbulence, or the outflows and
cover a great surface of the maps (yellowish regions of Fig. 4 for example).
Figures 13 and 14 show a general behavior for points at low to moderate column
densities ($\langle N\rangle\leq 10^{25}\leavevmode\nobreak\
\mathrm{cm^{-3}}$), where the discrepancy between the inferred magnetic field
orientation and the actual mean orientation increases with the intrinsic
degree of disorganization of the field lines and is consistent with the
expected angular dispersion (see Fig. 29).

On the other hand, at higher column densities, the behavior differs from what
is expected, displaying $|\Delta\phi|$ values higher than what can be
attributed to the dispersion only. Particularly, points with high
$|\Delta\phi|$ ($\geq 30^{\circ}$) and low values of $V_{\phi}$ ($\leq 0.2$)
are puzzling, but they remain extremely rare. Even the abundance of points
with $|\Delta\phi|>30^{\circ}$ and $V_{\phi}\leq 0.5$ remain extremely low at
the three wavelengths (typically less than $0.01\%$, with only the very
particular case of the simulation R7, with projection $y-z$ at
$m_{\mathrm{sink}}=0.3\leavevmode\nobreak\ \mathrm{M_{\odot},}$ reaching
$0.1\%$). These points can be only partially explained by a lower efficiency
of the grain alignment (see Appendix B for a comparison with the case of
perfectly aligned dust grains), a higher opacity of the layers surrounding the
protostar —which only allows us to recover the foreground information—, and a
dust temperature gradient, which are not included in the computation of the
mean orientation of the magnetic field. To assess whether or not the opacity
is responsible of the lower recovery rates, we stack the maps at
$\lambda=0.8\leavevmode\nobreak\ \mathrm{mm}$ for all the simulations, time-
steps, and projections and compute the recovery rate in three opacity ranges.
For all the lines of sight with $\tau_{\mathrm{0.8\leavevmode\nobreak\
mm}}\geq 1,$ we obtain a recovery rate of only $42.5\%$. For
$\tau_{\mathrm{0.8\leavevmode\nobreak\ mm}}$ in the range $0.1-1$, we obtain a
recovery rate of $79.4\%$, while for $\tau_{\mathrm{0.8\leavevmode\nobreak\
mm}}$ in the range $0.01-0.1$ the recovery rate reaches $81.7\%$. This clearly
shows that the reconstructed $B$-fields of regions with high column densities
most likely suffer from opacity effects, and that these regions might benefit
from multi-wavelength observations.
Figure 13 suggests that, in otherwise similar initial conditions
($\beta=0.04$, $\mathcal{M}=0$, $\theta=30^{\circ}$), the more magnetized case
(R2) has initially (at $m_{\mathrm{sink}}=0.1\leavevmode\nobreak\
\mathrm{M_{\odot}}$) higher levels of $|\Delta\phi|$, but as the simulations
evolve, the overall levels of discrepancy become comparable (similar contour
levels). Nevertheless, the most puzzling measurements being the sight lines
where the magnetic field is relatively well organized (low variance), but
displaying large discrepancies (high $|\Delta\phi|$), are more important for
the least magnetized simulation (R3) at all evolutionary stages (shown in Fig.
13), and they are mostly associated with high-column-density lines of sight
which correspond to the very small radii (disk scales) where radiative
transfer effects cannot easily be neglected, even at mm wavelengths.
For the nonstandard configurations shown in Fig. 14, those including
turbulence show relatively similar levels of discrepancy (simulations R6 and
R8), indicating that the relative importance of the magnetization level is
reduced in the presence of turbulence. Nevertheless, the presence of
turbulence does not seem to be the main driver of the low recovery rates. For
instance, Table 6 shows that simulations R1 and R6 display similar recovery
rates, within the error bars, at different scales.
On the other hand, simulations with $\theta=90^{\circ}$ (R7 and R9) show an
increased importance of the magnetization level. Simulation R9 (the least
magnetized one) exhibits a relatively peculiar behavior, where although the
overall $|\Delta\phi|$ values remain low, it develops a region where the
discrepancy is maximal (even for values of $V_{\phi}$ lower than $0.2$),
reaching values closer to $90^{\circ}$.
Figure 13 also shows that the accuracy of polarized dust emission maps in
tracing the magnetic field topology in protostars has weak and rather unclear
dependence on the initial values of the rotational energy
$\beta_{\mathrm{rot}}$ in the more magnetized case ($\mu=3.33$, first three
rows), displaying overall a larger level of discrepancy between the intrinsic
magnetic field mean orientation and the reconstructed orientation from the
polarized dust emission for the simulation with larger rotational energy (R2),
and the lowest one for the simulation with intermediate rotational energy
level (R1).

Figure 11: Two-dimensional histograms of the discrepancy between the mean
$B$-field orientation and the orientation inferred from the polarized dust
emission, and the intensity of the polarized dust emission at three
wavelengths. From left to right $\lambda=0.8,1.3$, and
$3.0\leavevmode\nobreak\ \mathrm{mm}$. Each histogram has been constructed
using the data stacked for all simulations, time-steps, and projections. The
color histograms show the results using only the lines of sight that have a
polarization fraction above $1\%$, while the gray-scale part corresponds to
the rest of the points. $0.1\leavevmode\nobreak\
\mathrm{M_{\odot}}$$0.2\leavevmode\nobreak\
\mathrm{M_{\odot}}$$0.3\leavevmode\nobreak\ \mathrm{M_{\odot}}$ Figure 12:
Percentage of pixels with $|\Delta\phi|>30^{\circ}$ at
$\lambda=1.3\leavevmode\nobreak\ \mathrm{mm}$ in logarithmic radial bins.
---
$\mu=3.33$$\mu=6.67$$\mu=10$$\beta=0.0025$$\beta=0.01$$\beta=0.04$$\beta=0.04$$\beta=0.04$
Figure 13: Distribution of the discrepancy between the $B$ angle inferred from
the synthetic polarized dust emission at $\lambda=1.3\leavevmode\nobreak\
\mathrm{mm}$ and the mean orientation of the $B$ lines in the simulation as a
function of the circular variance $V_{\phi}$ for the simulations with standard
conditions. All these simulations do not include any initial turbulence
($\mathcal{M}=0$) and have an initial inclination angle between the initial
rotation axis and the magnetic field $\theta$ of $30^{\circ}$. The contour
lines show the smoothed 2D histogram contour levels at $10^{5}$ (yellow),
$10^{4}$ (red), $10^{3}$ (purple), and $10^{2}$ (black) counts. The color-
coded background corresponds to the mean column density.
---
$\mu=3.33$$\mu=6.67$$\mu=10$$\beta=0.01,\mathcal{M}=1$$\beta=0.04,\theta=90^{\circ}$$\beta=0.01,\mathcal{M}=1$$\beta=0.04,\theta=90$
Figure 14: Distribution of $|\Delta\phi|$ at $\lambda=1.3\leavevmode\nobreak\
\mathrm{mm}$ as a function of the circular variance $V_{\phi}$ for the
simulations with nonstandard conditions. The contours and the background image
are the same as in Fig. 13. Simulations R6 and R8 include turbulence
($\mathcal{M}=1$), while the simulations R7 and R9 have an initial angle
between the rotation axis and the magnetic field axis of $\theta=90^{\circ}$.
---
B-RAT+II+IDGPA Figure 15: Two-dimensional histograms of the polarization
fraction as a function of the circular variance for all the stacked maps at
three wavelengths (from left to right $\lambda=0.8,1.3$ and
$3.0\leavevmode\nobreak\ \mathrm{mm}$). The top row shows the results for the
case where the alignment is computed according to the B-RAT theory along with
the II and IDG mechanisms (B-RAT+II+IDG), as described in Section 2.4.3, while
the bottom row shows the results for the case where dust grains are assumed to
be perfectly aligned (PA) with the magnetic field.

###  5.2 Depolarization and magnetic field organization

Finally, to better understand the influence of the magnetic field topology on
the polarization fraction, we compared the distribution of the polarization
fraction ($p_{\mathrm{frac}}$) as a function of the circular variance
($V_{\phi}$) obtained for our synthetic maps to a test case. In the case of
our synthetic maps, the alignment was computed according to the B-RAT theory,
including the imperfect internal alignment and the imperfect Davis-Greenstein
mechanisms, as described in Sect. 2.4.3 (top row of Fig. 15, labeled as
B-RAT+II+IDG), while for the test case the alignment is imposed in such a
manner that all the dust grains susceptible to align are perfectly aligned
with the $B$-field lines (bottom row of Fig. 15, labeled as PA). Fig. 15 shows
2D histograms of $p_{\mathrm{frac}}$ versus $V_{\phi}$ by stacking the results
of all the simulations, time-steps, and projections at three wavelengths. This
figure shows a clear anti-correlation between these two quantities, indicating
that the degree of polarization decreases with an increasing level of
disorganization of the magnetic field along the line of sight, which is
consistent with the anticorrelation between the polarization fraction and the
dispersion of the position angles found in ALMA observations of Class 0
objects (Le Gouellec et al. 2020).
The fact that both cases display the same behavior supports the idea that a
geometrically depolarized signal due to the nonorganized component of the
magnetic field could be greatly responsible for the lower levels of
polarization observed in denser regions down to the scales of protostellar
envelopes, which is in agreement with observations at larger scales (Doi et
al. 2020; Planck Collaboration et al. 2015b, a). However, while at large
scales the turbulence of the gas may be the main physical cause of the
magnetic field disorganization, inside protostellar cores it seems the
gravitational field and the launching of protostellar outflows are more likely
responsible, as models with turbulence do not show significantly more
disorganized $B$-field, nor poorer fidelity, than models without, as shown in
Table 6.

###  5.3 Significance of other polarization-producing processes in
protostellar envelopes

Our analysis relies on the hypothesis that dust grains imperfectly align with
the local magnetic field in a way that reflects the local radiation conditions
(B-RATs) and grain properties. Dedicated studies of other physical effects
producing polarization of the dust thermal emission, along with better dust
models, are necessary to shed light on the precise conditions in which
polarization seen at mm wavelengths could trace other processes. However, here
we reinforce our hypothesis by recalling a few of them, and give simple
qualitative arguments as to why they are not expected to be dominating in
typical conditions reigning in protostellar envelopes.

While at the protostar disk scales the alignment of dust grains with local
magnetic field lines can become more difficult to realize because of the
higher rate of collisions (Tazaki et al. 2017), the gaseous damping timescale
is inversely proportional to gas density and gas temperatures. Hence, such
loss of efficiency is not expected in protostellar envelopes where densities
are lower than $10^{6}\leavevmode\nobreak\ \mathrm{cm^{-3}}$ and temperatures
are quite low ($T<20\leavevmode\nobreak\ \mathrm{K}$). It has also been
suggested that in conditions of high opacity and anisotropic radiation field,
usually met in disks, the dust self-scattering of the radiation field from
large grains ($\sim 100\leavevmode\nobreak\ \mathrm{\mu m}$) could dominate
the polarized signal at mm wavelengths (Kataoka et al. 2015, 2017, see also
Kirchschlager & Bertrang 2020; Brunngräber & Wolf 2020). Despite the
indications of large grains at the protostellar scales probed in our analysis,
not only are the conditions not met to produce significant polarization from
self-scattering, but the levels of observed polarization in embedded
protostars (a few percent to a few tens of percent) cannot be produced with
self-scattering polarization.

A highly anisotropic radiation field could reduce the radiative precession
timescale to values shorter than the Larmor timescale, inducing polarization
of dust thermal emission due to the k-RAT mechanism, which is related to the
radiation anisotropy and not to the magnetic field. This mechanism is favored
for the large grains (for a discussion, see e.g. the work of Pattle et al.
2021 in the Orion Source I), and necessitate extremely irradiated conditions
at wavelengths similar to the dust grain sizes, which are not typical of
embedded low-mass protostars (Le Gouellec et al. in prep).

Finally, mechanical alignment of dust grains with the local gas flow has been
proposed as an extra mechanism to produce polarized dust emission from
irregular grains (Hoang et al. 2018), with a rotation axis parallel to the
velocity drift direction. While this mechanism could in principle be at work
for example around massive protostars where radiative pressure drives the
outflow (Mignon-Risse et al. 2021), its conditions seem very sensitive to the
flow properties (high velocities to enhance the alignment efficiency, but not
as high as to disrupt the dust grains) and quite unable to explain the large
range of observations showing polarized dust emission in protostars. We also
stress that no convincing observational signature of such alignment process
has yet been found in protostellar environments (Cortes et al. 2021; Aso et
al. 2021).

##  6 Conclusions

In this paper, we present a thorough investigation of the fidelity and
limitations of using dust polarized emission as a tracer of the magnetic field
topologies in the dense regions of protostellar cores forming solar type
stars. We focus here mostly on quantifying the errors induced by (i) widely
varying physical conditions along each individual line of sight used to probe
the dust grain emission of the protostellar objects, and (ii) its potentially
devastating combination with the averaging effects of complex $B$-field line
orientations along those same lines of sight.

To assess the importance of these effects, we performed an analysis of
magnetic field properties in 27 realizations of MHD models following the
evolution of physical properties in star-forming cores. Assuming a uniform
population of oblate dust grains whose sizes follow the standard MRN
distribution (Mathis et al. 1977), we produced synthetic polarized dust
emission maps for dust grains aligned with the local $B$-field thanks to
B-RATs mechanism. We provide a detailed comparison of these synthetic maps and
their relative models.

We find that, in most cases, (sub-)mm polarized dust emission is a robust
tracer of the magnetic field topologies in inner protostellar envelopes and is
successful at capturing the details of the magnetic field spatial distribution
down to radii $\sim 100$ au. Measurements of the line-of-sight-averaged
magnetic field line orientation using the polarized dust emission are precise
to $<15^{\circ}$ (typical of the error on polarization angles obtained with
observations from large mm polarimetric facilities such as ALMA) in about
$95\%$ of the independent lines of sight that pass through protostellar
envelopes at all radii. When focusing on the smaller envelope radii
$<500\leavevmode\nobreak\ \mathrm{au}$, where the magnetic field lines are
more likely perturbed from the initially smooth configuration by infall and
outflow, $75\%$ of the lines of sight still give robust results.

Large discrepancies between the integrated $B$-field mean orientation and the
orientation reconstructed from the polarized dust emission are mostly observed
in lines of sight probing large column densities. Our analysis shows that, at
disk scales, physical conditions producing thermal dust emission with high
opacity are mostly responsible for the small fraction of mediocre measurements
(inaccuracy of $50\%$ or more on the recovered $B$-field position angle at
opacities higher than $1\leavevmode\nobreak\ \mathrm{mag}$). While the use of
longer wavelengths allows a slightly better recovery of the true orientation
of the magnetic fields, our study suggest that, at those scales, most
measurements of magnetic field line orientation will suffer from large
($>30^{\circ}$) error bars. Discrepancy is also found associated to a small
fraction of lines of sight where the magnetic field is highly disorganized,
such as outflow cavity walls and accretion shocks: in such lines of sight the
concept of integrated $B$-field mean orientation does not make sense, and no
good measurement can therefore be done. We find that the polarization fraction
is correlated with the degree of organization of the magnetic field along the
line of sight, confirming that the behavior observed at larger scales in less
dense gas holds true down to the scales of the envelopes of protostars. Hence,
low levels of polarization fraction in the observations could be used as a
flag to avoid these measurements affected by large errors.

###### Acknowledgements.

We thank M. Küffmeier for constructive comments to improve the readability of
this paper. This research has received funding from the European Research
Council (ERC) under the European Union’s Horizon 2020 research and innovation
programme (MagneticYSOS project, Grant Agreement No 679937).

## References

  * Andersson et al. (2015) Andersson, B.-G., Lazarian, A., & Vaillancourt, J. E. 2015, ARA&A, 53, 501
  * Aso et al. (2021) Aso, Y., Kwon, W., Hirano, N., et al. 2021, ApJ, 920, 71
  * Barnett (1915) Barnett, S. J. 1915, Physical Review, 6, 239
  * Bate et al. (1995) Bate, M. R., Bonnell, I. A., & Price, N. M. 1995, MNRAS, 277, 362
  * Bleuler & Teyssier (2014) Bleuler, A. & Teyssier, R. 2014, MNRAS, 445, 4015
  * Brunngräber & Wolf (2020) Brunngräber, R. & Wolf, S. 2020, A&A, 640, A122
  * Brunngräber & Wolf (2021) Brunngräber, R. & Wolf, S. 2021, A&A, 648, A87
  * Cabedo et al. (2022) Cabedo, V., Maury, A., Girart, J. M., et al. 2022, arXiv e-prints, arXiv:2204.10043
  * Cheng et al. (2022) Cheng, Y., Gutermuth, R. A., Offner, S., et al. 2022, MNRAS, 512, 960
  * Cortes et al. (2021) Cortes, P. C., Le Gouellec, V. J. M., Hull, C. L. H., et al. 2021, ApJ, 907, 94
  * Davis & Greenstein (1951a) Davis, Leverett, J. & Greenstein, J. L. 1951a, ApJ, 114, 206
  * Davis & Greenstein (1951b) Davis, Jr., L. & Greenstein, J. L. 1951b, ApJ, 114, 206
  * Doi et al. (2020) Doi, Y., Hasegawa, T., Furuya, R. S., et al. 2020, ApJ, 899, 28
  * Dolginov & Mitrofanov (1976) Dolginov, A. Z. & Mitrofanov, I. G. 1976, Ap&SS, 43, 291
  * Draine (1996) Draine, B. T. 1996, in Astronomical Society of the Pacific Conference Series, Vol. 97, Polarimetry of the Interstellar Medium, ed. W. G. Roberge & D. C. B. Whittet, 16
  * Draine & Flatau (2000) Draine, B. T. & Flatau, P. J. 2000, DDSCAT: The discrete dipole approximation for scattering and absorption of light by irregular particles
  * Fromang et al. (2006) Fromang, S., Hennebelle, P., & Teyssier, R. 2006, A&A, 457, 371
  * Galametz et al. (2020) Galametz, M., Maury, A., Girart, J. M., et al. 2020, A&A, 644, A47
  * Galli & Shu (1993) Galli, D. & Shu, F. H. 1993, ApJ, 417, 243
  * Gold (1952) Gold, T. 1952, MNRAS, 112, 215
  * Hennebelle & Ciardi (2009) Hennebelle, P. & Ciardi, A. 2009, A&A, 506, L29
  * Hennebelle et al. (2016) Hennebelle, P., Commerçon, B., Chabrier, G., & Marchand, P. 2016, ApJ, 830, L8
  * Hennebelle et al. (2020) Hennebelle, P., Commerçon, B., Lee, Y.-N., & Charnoz, S. 2020, A&A, 635, A67
  * Herranen et al. (2021) Herranen, J., Lazarian, A., & Hoang, T. 2021, ApJ, 913, 63
  * Hildebrand & Dragovan (1995) Hildebrand, R. H. & Dragovan, M. 1995, ApJ, 450, 663
  * Hoang et al. (2018) Hoang, T., Cho, J., & Lazarian, A. 2018, ApJ, 852, 129
  * Hoang & Lazarian (2014) Hoang, T. & Lazarian, A. 2014, MNRAS, 438, 680
  * Hoang & Lazarian (2016) Hoang, T. & Lazarian, A. 2016, ApJ, 831, 159
  * Hull et al. (2017) Hull, C. L. H., Mocz, P., Burkhart, B., et al. 2017, ApJ, 842, L9
  * Jensen & Haugbølle (2018) Jensen, S. S. & Haugbølle, T. 2018, MNRAS, 474, 1176
  * Jones & Spitzer (1967) Jones, R. V. & Spitzer, Jr., L. 1967, ApJ, 147, 943
  * Jones & James (1969) Jones, T. A. & James, W. R. 1969, Journal of the International Association for Mathematical Geology, 1, 129
  * Joos et al. (2012) Joos, M., Hennebelle, P., & Ciardi, A. 2012, A&A, 543, A128
  * Kataoka et al. (2015) Kataoka, A., Muto, T., Momose, M., et al. 2015, ApJ, 809, 78
  * Kataoka et al. (2017) Kataoka, A., Tsukagoshi, T., Pohl, A., et al. 2017, ApJ, 844, L5
  * Kirchschlager & Bertrang (2020) Kirchschlager, F. & Bertrang, G. H. M. 2020, A&A, 638, A116
  * Krasnopolsky et al. (2012) Krasnopolsky, R., Li, Z.-Y., Shang, H., & Zhao, B. 2012, ApJ, 757, 77
  * Kuffmeier et al. (2018) Kuffmeier, M., Frimann, S., Jensen, S. S., & Haugbølle, T. 2018, MNRAS, 475, 2642
  * Kuffmeier et al. (2020) Kuffmeier, M., Reissl, S., Wolf, S., Stephens, I., & Calcutt, H. 2020, A&A, 639, A137
  * Lam et al. (2021) Lam, K. H., Chen, C.-Y., Li, Z.-Y., et al. 2021, MNRAS, 507, 608
  * Landi Degl’Innocenti & Landolfi (2004) Landi Degl’Innocenti, E. & Landolfi, M. 2004, Polarization in Spectral Lines, Vol. 307
  * Lazarian & Hoang (2007a) Lazarian, A. & Hoang, T. 2007a, MNRAS, 378, 910
  * Lazarian & Hoang (2007b) Lazarian, A. & Hoang, T. 2007b, ApJ, 669, L77
  * Le Gouellec et al. (2020) Le Gouellec, V. J. M., Maury, A. J., Guillet, V., et al. 2020, A&A, 644, A11
  * Li & Draine (2001) Li, A. & Draine, B. T. 2001, ApJ, 554, 778
  * Li et al. (2014) Li, Z.-Y., Banerjee, R., Pudritz, R. E., et al. 2014, Protostars and Planets VI, 173
  * Li & Shu (1996) Li, Z.-Y. & Shu, F. H. 1996, ApJ, 472, 211
  * Machida et al. (2014) Machida, M. N., Inutsuka, S.-i., & Matsumoto, T. 2014, MNRAS, 438, 2278
  * Mardia (1975) Mardia, K. V. 1975, Journal of the Royal Statistical Society: Series B (Methodological), 37, 349
  * Martínez González et al. (2011) Martínez González, M. J., Manso Sainz, R., Asensio Ramos, A., & Belluzzi, L. 2011, Monthly Notices of the Royal Astronomical Society, 419, 153
  * Masson et al. (2016) Masson, J., Chabrier, G., Hennebelle, P., Vaytet, N., & Commerçon, B. 2016, A&A, 587, A32
  * Mathis et al. (1983) Mathis, J. S., Mezger, P. G., & Panagia, N. 1983, A&A, 128, 212
  * Mathis et al. (1977) Mathis, J. S., Rumpl, W., & Nordsieck, K. H. 1977, ApJ, 217, 425
  * Matsumoto et al. (2017) Matsumoto, T., Machida, M. N., & Inutsuka, S.-i. 2017, ApJ, 839, 69
  * Mignon-Risse et al. (2021) Mignon-Risse, R., González, M., & Commerçon, B. 2021, A&A, 656, A85
  * Mouschovias & Ciolek (1999) Mouschovias, T. C. & Ciolek, G. E. 1999, in NATO Advanced Study Institute (ASI) Series C, Vol. 540, The Origin of Stars and Planetary Systems, ed. C. J. Lada & N. D. Kylafis, 305
  * Offner & McKee (2011) Offner, S. S. R. & McKee, C. F. 2011, ApJ, 736, 53
  * Offner et al. (2022) Offner, S. S. R., Moe, M., Kratter, K. M., et al. 2022, arXiv e-prints, arXiv:2203.10066
  * Pattle et al. (2021) Pattle, K., Lai, S.-P., Wright, M., et al. 2021, MNRAS, 503, 3414
  * Planck Collaboration et al. (2015a) Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2015a, A&A, 576, A105
  * Planck Collaboration et al. (2015b) Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2015b, A&A, 576, A104
  * Purcell (1979) Purcell, E. M. 1979, ApJ, 231, 404
  * Reissl et al. (2022) Reissl, S., Meehan, P., & Klessen, R. S. 2022, arXiv e-prints, arXiv:2201.03694
  * Reissl et al. (2017) Reissl, S., Seifried, D., Wolf, S., Banerjee, R., & Klessen, R. S. 2017, A&A, 603, A71
  * Reissl et al. (2016) Reissl, S., Wolf, S., & Brauer, R. 2016, A&A, 593, A87
  * Roberge (2004) Roberge, W. G. 2004, in ASP Conf. Series, Vol. 309, Astrophysics of Dust, ed. A. N. Witt, G. C. Clayton, & B. T. Draine, 467
  * Seifried et al. (2011) Seifried, D., Banerjee, R., Klessen, R. S., Duffin, D., & Pudritz, R. E. 2011, MNRAS, 417, 1054
  * Spitzer & McGlynn (1979) Spitzer, Jr., L. & McGlynn, T. A. 1979, ApJ, 231, 417
  * Tazaki et al. (2017) Tazaki, R., Lazarian, A., & Nomura, H. 2017, ApJ, 839, 56
  * Teyssier (2002) Teyssier, R. 2002, A&A, 385, 337
  * Valdivia et al. (2019) Valdivia, V., Maury, A., Brauer, R., et al. 2019, MNRAS, 488, 4897

##  Appendix A Additional material

In this section, we provide complementary material on synthetic observations
and their comparisons to properties measured from the models. Figures 16 and
17 show the maps from the synthetic observations of the polarized dust
emission of the 27 models, at $\lambda=0.8$ and $3.0\leavevmode\nobreak\
\mathrm{mm}$. Figures 18 and 19 show the azimuthally averaged profiles of the
total dust emission flux density, linearly polarized flux density, and the
polarization fraction (averaged for the three projections) at $\lambda=0.8$
and $3.0\leavevmode\nobreak\ \mathrm{mm}$, respectively.
Figures 20, 21, and 22 are maps of the difference in measured mean $B$-field
orientation between the model value and the value inferred from the polarized
dust emission at $\lambda=0.8$, $1.3,$ and $3.0\leavevmode\nobreak\
\mathrm{mm}$. We find mostly small differences (light colors or white), and
highlight in dark colors the regions where the $B$-field orientation inferred
from the polarized dust emission is very different from the mean magnetic
field orientation.
Figures 23 and 24 (Figures 25 and 26) show the distribution of the discrepancy
between the inferred orientation of the magnetic field and the true mean
orientation at $\lambda=0.8\leavevmode\nobreak\ \mathrm{mm}$
($\lambda=3.0\leavevmode\nobreak\ \mathrm{mm}$) for the simulations with
standard conditions and the ones starting from less prototypical conditions
(either $\mathcal{M}=1$ or $\theta=90^{\circ}$), respectively.
Table 5 gives the detailed recovery rates at $\lambda=0.8$ for all the
simulations (ID), projections ($los$), and evolutionary stages
($m_{\mathrm{sink}}$) at different scales. The fourth column shows the
recovery rate over the whole map, while columns 5, 6, and 7 correspond to the
recovery rate in concentric rings of radial distances comprised in the
interval given in the header of the table.
We verified the influence of the presence of turbulence in Table 6. This table
shows the mean recovery rates for simulations R1 and R6. These two otherwise
identical simulations ($\mu=3.33$, $\beta_{\mathrm{rot}}=0.01$ and
$\theta=30^{\circ}$), differing only in the degree of turbulence (R6 has
$\mathcal{M}=1$), show that the degree of turbulence does not have a strong
impact on the recovery rate.
Finally, Tables 7 and 8 show the mean recovery rates separated into two cases
(moderate and low magnetization degree) similar to Table 3, but at
$\lambda=1.3$, and $3.0\leavevmode\nobreak\ \mathrm{mm}$, respectively.

Figure 16: Results from the synthetic observation at
$\lambda=0.8\leavevmode\nobreak\ \mathrm{mm}$ showing the polarization
fraction (background image), the polarized intensity contours, and the
inferred magnetic field orientation vectors. The contours correspond to the
linearly polarized intensity normalized to the peak value at levels $10^{-2}$
(black), $10^{-3}$ (purple), $10^{-4}$ (red), and $10^{-5}$ (yellow). Higher
resolution figures for individual simulations are available online. Figure 17:
Results from the synthetic observation at $\lambda=3.0\leavevmode\nobreak\
\mathrm{mm}$ showing the polarization fraction (background image), the
polarized intensity contours, and the inferred magnetic field orientation
vectors. The contours correspond to the linearly polarized intensity
normalized to the peak value at levels $10^{-2}$ (black), $10^{-3}$ (purple),
$10^{-4}$ (red), and $10^{-5}$ (yellow). Higher resolution figures for
individual simulations are available online.
---
$0.1\leavevmode\nobreak\ \mathrm{M_{\odot}}$$0.2\leavevmode\nobreak\
\mathrm{M_{\odot}}$$0.3\leavevmode\nobreak\ \mathrm{M_{\odot}}$ Figure 18:
Radial profiles of the synthetic observation at
$\lambda=0.8\leavevmode\nobreak\ \mathrm{mm}$. From top to bottom: Total dust
emission (Stokes $I$), total linearly polarized dust emission ($PI$), and
polarization fraction ($p_{\mathrm{frac}}$).
---
$0.1\leavevmode\nobreak\ \mathrm{M_{\odot}}$$0.2\leavevmode\nobreak\
\mathrm{M_{\odot}}$$0.3\leavevmode\nobreak\ \mathrm{M_{\odot}}$ Figure 19:
Radial profiles of the synthetic observation at
$\lambda=3.0\leavevmode\nobreak\ \mathrm{mm}$. From top to bottom: Total dust
emission (Stokes $I$), total linearly polarized dust emission ($PI$), and
polarization fraction ($p_{\mathrm{frac}}$). Figure 20: Difference maps
$\Delta\phi$ at $\lambda=0.8\leavevmode\nobreak\ \mathrm{mm}$. The simulation
identifier is shown in the bottom left corner of each block, while the
projection and evolutionary step (given by the sink mass) are indicated at the
top of each individual map. Figure 21: Difference maps $\Delta\phi$ at
$\lambda=1.3\leavevmode\nobreak\ \mathrm{mm}$. The structure of the figure is
the same as for Fig. 20. Figure 22: $\Delta\phi$ at
$\lambda=3.0\leavevmode\nobreak\ \mathrm{mm}$. The structure of the figure is
the same as for Figs. 20 and 21.
---
$\mu=3.33$$\mu=6.67$$\mu=10$$\beta=0.0025$$\beta=0.01$$\beta=0.04$$\beta=0.04$$\beta=0.04$
Figure 23: Distribution of the discrepancy $|\Delta\phi|$ between the $B$
angle inferred from the synthetic polarized dust emission at
$\lambda=0.8\leavevmode\nobreak\ \mathrm{mm}$ and the mean orientation of the
$B$ lines in the simulation as a function of the circular variance $V_{\phi}$
for the simulations with standard conditions. None of these simulations
includes any initial turbulence ($\mathcal{M}=0$) and all have an initial
inclination angle between the initial rotation axis and the magnetic field
$\theta$ of $30^{\circ}$. The contour lines show the smoothed 2D histogram
contour levels at $10^{5}$ (yellow), $10^{4}$ (red), $10^{3}$ (purple), and
$10^{2}$ (black) counts. The color coded background corresponds to the mean
column density.
---
$\mu=3.33$$\mu=6.67$$\mu=10$$\beta=0.01,\mathcal{M}=1$$\beta=0.04,\theta=90^{\circ}$$\beta=0.01,\mathcal{M}=1$$\beta=0.04,\theta=90$
Figure 24: Distribution of $|\Delta\phi|$ at $\lambda=0.8\leavevmode\nobreak\
\mathrm{mm}$ as a function of the circular variance $V_{\phi}$ for the
simulations with nonstandard conditions. The contours and the background image
are the same as in Fig. 23. Simulations R6 and R8 include turbulence
($\mathcal{M}=1$), while the simulations R7 and R9 have an initial angle
between the rotation axis and the magnetic field axis of $\theta=90^{\circ}$.
---
$\mu=3.33$$\mu=6.67$$\mu=10$$\beta=0.0025$$\beta=0.01$$\beta=0.04$$\beta=0.04$$\beta=0.04$
Figure 25: Distribution of the discrepancy between the $B$ angle inferred from
the synthetic polarized dust emission at $\lambda=3.0\leavevmode\nobreak\
\mathrm{mm}$ and the mean orientation of the $B$ lines in the simulation as a
function of the circular variance $V_{\phi}$ for the simulations with standard
conditions. All these simulations do not include any initial turbulence
($\mathcal{M}=0$) and have an initial inclination angle between the initial
rotation axis and the magnetic field $\theta$ of $30^{\circ}$. The contour
lines show the smoothed 2D histogram contour levels at $10^{5}$ (yellow),
$10^{4}$ (red), $10^{3}$ (purple), and $10^{2}$ (black) counts. The color
coded background corresponds to the mean column density.
---
$\mu=3.33$$\mu=6.67$$\mu=10$$\beta=0.01,\mathcal{M}=1$$\beta=0.04,\theta=90^{\circ}$$\beta=0.01,\mathcal{M}=1$$\beta=0.04,\theta=90$
Figure 26: Distribution of $|\Delta\phi|$ at $\lambda=3.0\leavevmode\nobreak\
\mathrm{mm}$ as a function of the circular variance $V_{\phi}$ for the
simulations with non standard conditions. The contours and the background
image are the same as in Fig. 25. Simulations R6 and R8 include turbulence
($\mathcal{M}=1$), while the simulations R7 and R9 have an initial angle
between the rotation axis and the magnetic field axis of $\theta=90^{\circ}$.
Table 5: Fidelity of the $B$-field mapping from polarized dust emission at
$\lambda=0.8\leavevmode\nobreak\ \mathrm{mm}$ at different scales.

ID | $m_{\mathrm{sink}}$ | $los$ | Recovery rate at different scales $[\%]$ | ID | $m_{\mathrm{sink}}$ | $los$ | Recovery rate at different scales $[\%]$
---|---|---|---|---|---|---|---
| $[\mathrm{M_{\odot}}]$ |  | All scales | $30-100$ | $100-500$ | $500-1500$ |  | $[\mathrm{M_{\odot}}]$ |  | All scales | $30-100$ | $100-500$ | $500-1500$
|  |  |  | au | au | au |  |  |  |  | au | au | au
R1 | 0.1 | $x$ | 99.1 | 82.2 | 96.8 | 99.0 | R7 | 0.1 | $x$ | 99.3 | 88.9 | 96.6 | 99.5
|  | $y$ | 96.8 | 60.9 | 76.3 | 90.5 |  |  | $y$ | 92.1 | 55.1 | 33.6 | 71.1
|  | $z$ | ¿99.9 | 68.7 | ¿99.9 | 100.0 |  |  | $z$ | ¿99.9 | 95.1 | ¿99.9 | 100.0
| 0.2 | $x$ | 97.7 | 80.5 | 72.2 | 95.5 |  | 0.2 | $x$ | 97.8 | 76.2 | 70.4 | 98.6
|  | $y$ | 94.4 | 66.3 | 70.3 | 82.3 |  |  | $y$ | 91.9 | 93.4 | 46.1 | 66.1
|  | $z$ | ¿99.9 | 73.4 | 99.0 | 100.0 |  |  | $z$ | ¿99.9 | 92.7 | 99.5 | 100.0
| 0.3 | $x$ | 97.2 | 72.1 | 90.4 | 87.0 |  | 0.3 | $x$ | 95.8 | 76.0 | 83.6 | 92.4
|  | $y$ | 93.3 | 61.5 | 76.0 | 84.4 |  |  | $y$ | 89.4 | 76.6 | 66.8 | 63.1
|  | $z$ | 99.5 | 72.5 | 90.3 | 99.5 |  |  | $z$ | 99.8 | 81.6 | 96.3 | ¿99.9
R2 | 0.1 | $x$ | 97.7 | 51.0 | 75.2 | 97.2 | R8 | 0.1 | $x$ | 94.8 | 80.4 | 70.9 | 88.0
|  | $y$ | 91.0 | 50.3 | 57.0 | 76.2 |  |  | $y$ | 97.4 | 70.8 | 96.6 | 96.9
|  | $z$ | 99.8 | 75.7 | 99.4 | ¿99.9 |  |  | $z$ | 93.2 | 75.3 | 90.8 | 96.6
| 0.2 | $x$ | 96.5 | 67.2 | 65.7 | 86.9 |  | 0.2 | $x$ | 93.5 | 72.7 | 79.9 | 86.5
|  | $y$ | 88.1 | 79.3 | 48.4 | 68.5 |  |  | $y$ | 96.9 | 78.3 | 88.5 | 97.8
|  | $z$ | 99.7 | 82.8 | 96.2 | 99.1 |  |  | $z$ | 93.9 | 92.0 | 83.7 | 95.6
| 0.3 | $x$ | 92.6 | 69.5 | 73.8 | 85.1 |  | 0.3 | $x$ | 92.5 | 86.8 | 95.6 | 87.3
|  | $y$ | 83.7 | 81.7 | 79.7 | 69.6 |  |  | $y$ | 94.7 | 85.7 | 87.7 | 88.7
|  | $z$ | 98.4 | 67.4 | 87.1 | 97.8 |  |  | $z$ | 95.5 | 85.1 | 92.2 | 92.7
R3 | 0.1 | $x$ | 98.7 | 89.4 | 87.8 | 94.3 | R9 | 0.1 | $x$ | 99.8 | 91.8 | 99.7 | 99.5
|  | $y$ | 96.4 | 93.0 | 89.7 | 86.2 |  |  | $y$ | 94.9 | 97.2 | 71.4 | 98.0
|  | $z$ | 98.2 | 80.2 | 95.2 | 99.2 |  |  | $z$ | 99.8 | 80.7 | 95.5 | 100.0
| 0.2 | $x$ | 94.5 | 86.5 | 92.8 | 87.4 |  | 0.2 | $x$ | 99.5 | 63.7 | 99.7 | 98.7
|  | $y$ | 90.5 | 92.9 | 97.6 | 78.7 |  |  | $y$ | 97.1 | 91.1 | 91.9 | 95.1
|  | $z$ | 95.4 | 77.1 | 91.1 | 86.9 |  |  | $z$ | 99.3 | 79.7 | 94.9 | 97.2
| 0.3 | $x$ | 88.7 | 68.9 | 96.0 | 89.1 |  | 0.3 | $x$ | 99.0 | 67.5 | 98.5 | ¿99.9
|  | $y$ | 79.7 | 84.3 | 97.9 | 94.6 |  |  | $y$ | 84.4 | 91.4 | 90.5 | 90.0
|  | $z$ | 81.1 | 67.3 | 87.6 | 91.4 |  |  | $z$ | 96.3 | 90.4 | 89.4 | 87.1
R4 | 0.1 | $x$ | 98.9 | 81.7 | 73.9 | 96.4 | R10 | 0.1 | $x$ | 99.3 | 54.7 | 84.5 | 99.2
|  | $y$ | 96.1 | 86.8 | 78.1 | 90.6 |  |  | $y$ | 93.1 | 82.1 | 62.1 | 86.3
|  | $z$ | 99.6 | 77.2 | 86.9 | 99.3 |  |  | $z$ | 99.9 | 78.0 | 98.5 | 99.9
| 0.2 | $x$ | 96.9 | 69.5 | 80.1 | 86.0 |  | 0.2 | $x$ | 98.9 | 69.6 | 80.6 | 97.3
|  | $y$ | 92.7 | 73.4 | 87.1 | 75.6 |  |  | $y$ | 91.6 | 79.0 | 46.2 | 83.7
|  | $z$ | 98.8 | 74.9 | 97.3 | 99.7 |  |  | $z$ | 99.7 | 59.4 | 94.5 | 100.0
| 0.3 | $x$ | 93.0 | 64.6 | 70.4 | 84.0 |  | 0.3 | $x$ | 97.1 | 61.4 | 65.5 | 90.3
|  | $y$ | 85.8 | 47.8 | 73.0 | 73.7 |  |  | $y$ | 89.1 | 90.4 | 44.4 | 65.2
|  | $z$ | 97.1 | 62.5 | 88.6 | 96.0 |  |  | $z$ | 99.0 | 49.8 | 84.7 | 98.6
R6 | 0.1 | $x$ | 95.5 | 81.9 | 91.5 | 96.0 |  |  |  |  |  |  |
|  | $y$ | 92.5 | 49.5 | 74.2 | 80.0 |  |  |  |  |  |  |
|  | $z$ | 96.8 | 76.7 | 89.3 | 98.4 |  |  |  |  |  |  |
| 0.2 | $x$ | 93.6 | 74.8 | 89.9 | 97.5 |  |  |  |  |  |  |
|  | $y$ | 87.1 | 50.5 | 66.3 | 68.3 |  |  |  |  |  |  |
|  | $z$ | 96.8 | 68.9 | 91.7 | 94.1 |  |  |  |  |  |  |
| 0.3 | $x$ | 89.5 | 89.5 | 82.0 | 93.2 |  |  |  |  |  |  |
|  | $y$ | 83.4 | 73.1 | 67.7 | 78.7 |  |  |  |  |  |  |
|  | $z$ | 91.9 | 71.5 | 85.1 | 84.0 |  |  |  |  |  |  |

777 Here we present the recovery rates obtained for each simulation (ID),
evolutionary stage ($m_{\mathrm{sink}}$), and projection ($los$) using the
synthetic observations at $\lambda=0.8\leavevmode\nobreak\ \mathrm{mm}$. The
fourth column shows the recovery rate over the whole map, while the rest of
the scales correspond to the lines of sight with radial distances to the
central source comprised in the interval indicated in the header. A synthesis
of this table is given in Table 3.

Table 6: Influence of the turbulence on the mean fidelity of $B$-field mapping
from polarized dust emission at $\lambda=0.8\leavevmode\nobreak\ \mathrm{mm}$
in protostellar envelopes at different scales.

$m_{\mathrm{sink}}$ |  $B$-field recovery rates and $\sigma$ $[\%]$
---|---
$[\mathrm{M_{\odot}}]$ | All scales |  $30-100$ au |  $100-500$ au |  $500-1500$ au
R1 | 97.5 2.4 | 70.9 7.5 | 85.7 12.0 | 93.1 7.2
R6 | 91.9 4.5 | 70.7 13.2 | 82.0 10.1 | 87.8 10.5
888 Mean recovery rates at different scales for an error bar of $15^{\circ}$.
The values have been averaged for all the evolutionary stages and projections
for simulations R1 (non-turbulent) and R6 (turbulent) given in Table 5. The
values given in a smaller font correspond to the standard deviation.

Table 7: Fidelity of $B$-field mapping from polarized dust emission in
protostellar envelopes at $\lambda=1.3\leavevmode\nobreak\ \mathrm{mm}$.

$m_{\mathrm{sink}}$ |  $B$-field recovery rates and $\sigma$ $[\%]$
---|---
$[\mathrm{M_{\odot}}]$ | All scales |  $30-100$ au |  $100-500$ au |  $500-1500$ au
$\mu=3.33$
0.1 | 97.1 3.0 | 71.9 14.9 | 83.7 18.4 | 93.4 9.1
0.2 | 95.9 4.0 | 76.2 11.3 | 77.3 17.9 | 90.2 11.5
0.3 | 93.7 5.1 | 74.6 10.4 | 79.4 12.5 | 86.9 11.5
$\mu=6.67,10.$
0.1 | 97.5 2.0 | 85.6 7.8 | 87.9 9.1 | 95.7 4.4
0.2 | 96.0 2.7 | 82.2 9.5 | 91.3 6.0 | 90.6 7.1
0.3 | 91.2 5.8 | 77.8 13.7 | 90.7 7.4 | 90.6 6.5
999 Mean recovery rates as in Table 3, but for maps at
$\lambda=1.3\leavevmode\nobreak\ \mathrm{mm}$.

Table 8: Fidelity of $B$-field mapping from polarized dust emission in
protostellar envelopes at $\lambda=3.0\leavevmode\nobreak\ \mathrm{mm}$.

$m_{\mathrm{sink}}$ |  $B$-field recovery rates and $\sigma$ $[\%]$
---|---
$[\mathrm{M_{\odot}}]$ | All scales |  $30-100$ au |  $100-500$ au |  $500-1500$ au
$\mu=3.33$
0.1 | 97.2 3.0 | 72.9 14.7 | 84.9 18.0 | 93.9 8.8
0.2 | 96.2 3.8 | 77.4 11.3 | 78.8 17.3 | 90.9 11.0
0.3 | 94.0 5.0 | 75.8 10.5 | 80.5 12.0 | 87.7 11.2
$\mu=6.67,10.$
0.1 | 97.6 2.0 | 86.4 7.9 | 89.1 8.4 | 95.9 4.3
0.2 | 96.2 2.7 | 83.4 10.0 | 92.1 5.6 | 91.3 7.1
0.3 | 91.7 5.6 | 79.3 13.6 | 92.2 6.4 | 91.3 6.4
101010 Mean recovery rates as in Table 3, but for maps at
$\lambda=3.0\leavevmode\nobreak\ \mathrm{mm}$.

##  Appendix B Additional tests for understanding the origin of points with
low dispersion and high discrepancy

Even though the great majority of discrepancies between the polarization
vectors and the mean orientation of the magnetic field can be understood as a
consequence of the intrinsic organization level, some points of high
discrepancy ($|\Delta\phi|>30^{\circ}$) and low disorganization levels
($V_{\phi}<0.5$) remain unexplained. To better understand this issue we
explore the effect of the opacity, which can cause a signal to be dominated by
the foreground external layers, and thus not tracing the full line of sight,
and the effect of the alignment efficiency, which if lowered can hinder the
ability of polarization to properly trace the magnetic field.
The opacity effect can be studied by analyzing the distribution of
$|\Delta\phi|$ as a function of the opacity and the projected radius at the
available wavelengths. In Fig. 27 we present the stacked (for all simulations,
projections and time-steps) distribution of the discrepancy $|\Delta\phi|$ as
a function of the intrinsic level of disorganization of the field lines
($V_{\phi}$), the opacity ($\tau$) and as a function of the radial distance
($r$). The colored region corresponds to the subset of points where
$|\Delta\phi|>30^{\circ}$ and $V_{\phi}\leq 0.5$, while the gray background
corresponds to all the points (given for reference). This figure shows that
the bulk of the points at $|\Delta\phi|\sim 30-40^{\circ}$ is still consistent
with the intrinsic disorganization level ($\sigma_{\phi}\sim 30^{\circ}$ when
$V_{\phi}\sim 0.5$), and that they are likely produced in regions disturbed by
the outflows at a distance of $200-1000\leavevmode\nobreak\ \mathrm{au}$.
Lines of sight with extremely high discrepancy levels
$|\Delta\phi|>70^{\circ}$ are associated to regions located within a radius of
roughly $100\leavevmode\nobreak\ \mathrm{au}$ with high opacities of about
$2-100$. This might indicate that in these regions the signal recovered in the
synthetic observations is just partial. This is coherent with the decrease of
the number of problematic lines of sight for longer wavelengths, that have a
lower optical depth.
To test the effect of the alignment efficiency we performed an extra set of
identical synthetic observations, except for the alignment properties. The
Polaris code allows us to simulate the polarized dust emission arising from
perfectly aligned dust grains. When this hypothesis is used, all the dust
grains susceptible to be aligned will contribute to the polarized dust
emission, and the resulting polarized intensity and polarization fraction will
correspond to the upper limits of these quantities. Figure 28 displays the
same plots as Fig.27, but for the perfectly aligned dust grains hypothesis.
From this figure we find that the number of lines of sight in the region
$|\Delta\phi|>30^{\circ}$ and $V_{\phi}<0.5$ is reduced by $25$, $22$, and
$25\%$ at $\lambda=0.8,1.3,$ and $3.0\leavevmode\nobreak\ \mathrm{mm}$,
respectively, when using the perfect alignment hypothesis, indicating that a
lower alignment efficiency is also contributing to a worse description of the
mean magnetic field orientation, especially in the inner
$200\leavevmode\nobreak\ \mathrm{au}$. Additionally, these results highlight
the role of the opacity effects, since even though dust grains are perfectly
aligned, there are remaining patches of problematic lines of sight. This is
likely due to the fact that the emission is dominated by the most external
layers in the foreground, which is consistent with the fact that the number of
conflicting lines of sight decreases with the observed wavelength (or
equivalently with the opacity).

|  |
---|---|---
|  |
|  |
Figure 27: Two-dimensional histograms for the RAT case of the stacked maps (all simulations, time-steps, and projections) showing the $\Delta\phi$ as a function of $V_{\phi}$, vs. optical depth and vs. the radial distance $r$ (from left to right), at three wavelengths $\lambda=0.8,1.3,3.0\leavevmode\nobreak\ \mathrm{mm}$ (top to bottom). |  |
---|---|---
|  |
|  |
Figure 28: Same as Fig. 27 but for the case where dust grains are perfectly
aligned.

##  Appendix C Relation between the circular variance and the angular
dispersion

To understand how the circular variance ($V_{\phi}$), defined in Sect. 2.3,
and the angular dispersion ($\sigma_{\phi}$) are related we simulated several
ensembles of random angles and computed the resulting $V_{\phi}$ and
$\sigma_{\phi}$. To avoid any bias related to the position angle, we set $180$
regularly spaced central angles ($\phi_{c}$) from $-90^{\circ}$ to
$89^{\circ}$. Secondly, around each central angle we defined $18$ opening
angles $\Delta\theta$, varying from $5^{\circ}$ to $90^{\circ}$ and we
randomly generated an ensemble of $N=500$ uniformly distributed angles in the
interval $[\phi_{c}-\Delta\theta,\phi_{c}+\Delta\theta)$.
As the mean angle does not represent the true orientation when dealing with
bidirectional data, we define the mean orientation $\langle\phi\rangle$ as in
Eq 3. The dispersion is thus defined as follows:

$\sigma_{\phi}=\sqrt{\frac{1}{N}\sum_{i=1}^{N}{(\phi_{i}-\langle\phi\rangle)^{2}}},$
(9)

where we use the smallest difference with respect to the mean mean orientation
($|(\phi_{i}-\langle\phi\rangle)|\leq 90^{\circ}$).
We computed $V_{\phi}$ and $\sigma_{\phi}$ for each ensemble, and repeated the
procedure $200$ times for each central angle and opening angle (totaling
$648000$ ensembles). We present the results for the relation between the
circular variance $V_{\phi}$ and the angular dispersion in Fig. 29. Each point
of the figure shows the circular variance as a function of the angular
dispersion of each individual simulated ensemble of $500$ data points. In this
figure we color coded the opening angle $\Delta\theta$ of each ensemble. We
note that the angular dispersion is necessarily smaller than $\Delta\theta$,
which only represents the maximum allowed variation with respect to the
central angle $\phi_{c}$. Figure 29 shows that $\sigma_{\phi}$ increases
monotonically with $V_{\phi}$ (and with $\Delta\theta$), but its maximum value
is limited to less than $55^{\circ}$.

Figure 29: Relation between the circular variance $V_{\phi}$ and the angular
dispersion.