© V. Valdivia et al. 2022

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 (≲ 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.

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 M_⊙, and the magnetic field is initially uniform and aligned with the z direction. The ratio of thermal to gravitational energy is α = 0.4 for all the simulations. The selected models cover three levels of rotation (set by the ratio of rotational to gravitational energy β_rot = E_rot/E_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 μ = (M/Φ)/(M/Φ)_crit = 3.33, 6.67, 10.0). Most of the models have a Mach number M = 0 and an inclination between the initial rotation axis and the initial magnetic field orientationθ= 30°, but special cases including turbulence (M = 1) and cases using an extreme tilt of θ = 90° are also included. The run identifiers (ID), as well as their parameters, are summarized in the first part of Table 1.

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_thres = 3 × 10¹³ cm⁻³, 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_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 (μ = 3.33) and having a moderate rotation (β_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 β_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_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 M_⊙. 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_sink, has reached 0.1, 0.2, and 0.3 M_⊙, respectively, thus sampling the Class 0 stage (M_env > M_star). The time at which the sink forms in each simulation (t_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 M_⊙.

Fig. 1 Sink-mass evolution as a function of the time since the formation of the sink tage.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.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 anticlockwise, as recommended by the IAU¹, in the range −90° to +90°. 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ϕ_B, where ϕ_B is the local magnetic field angle projected onto the plane of the sky, as follows: (1) (2)

where ρ is the local gas number density, and the averages are computed along the line of sight los. The mean orientation of the magnetic field 〈ϕ_B〉 is then calculated as: (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_ϕ as follows: (4)

Clustered values of ϕ_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 C.1 shows the relation between the circular variance V_ϕ and the typical standard deviation σ_ϕ in degrees.

Fig. 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 msink = 0.1, 0.2, and 0.3 M⊙ (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.

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 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 au). As most of the reprocessing of the protostellar photons into infrared (IR) photons happens at radii 4–10 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_mag), and mean magnetic field orientation (ϕ_B), as well as the associated circular variance (V_ϕ) at key time-steps in Sect. 3. We point out that the advantage of providing the mean magnetic pressure (P_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_mag, which is given by: (5)

where is the local magnetic field strength in μ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 L_⊙ of radius of 1 R_⊙ 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₀ = 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 (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_min = 5 nm and a_max = 20 um, respectively. The choice of a_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_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_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_high−J ~ 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 ~10 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 pc and observe a squared region of 8000 au on a side centered at the position of the protostellar embryo (materialized by the sink particle) at a resolution of 5 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 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_frac from the Stokes maps as follows: (7) (8)

At far-infrared (FIR) to millimiter wavelengths, the magnetic field orientation inferred from the polarized emission (ϕ_pol) is obtained by rotating the polarization vector by 90°.

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 both. 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_mag, as described in Sect. 2.4.1) with overlaid mean orientation (ϕ_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_ϕ 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_ϕ in Figs. 5 and 6, respectively.

Fig. 3 Mean magnetic pressure Pmag 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.

Fig. 4 Circular variance Vϕ of the B-fleld 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.

Fig. 5 Radial mean column density profiles. From left to right: profiles at msink = 0.1, 0.2, and 0.3 M⊙.

Fig. 6 Radial mean profiles of the circular variance Vϕ. From left to right: profiles at msink = 0.1, 0.2, and 0.3 M⊙.

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 (μ). 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 (μ = 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 (β = 0.04 andθ= 90°), but the outcome of the evolution is very different. Simulation R9 (μ = 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 β-field lines. Simulations R6 and R8 also share the same rotation (β = 0.04 andθ= 30°) and turbulent (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 μ 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_ϕ starts to increase. At larger scales (r > 4000 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 σ_ϕ = 20°–30° (see Appendix C), while regions with purple to red have larger angular dispersions (σ_ϕ > 30°). These regions of larger 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_ϕ remains small even for the cases R7 and R9 where the rotation axis is tilted by 90° 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-neld 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 ~ 2 × 10²⁴ cm⁻² within a radius of r ≲ 60 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 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θ= 90° andθ= 30°, 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_ϕ profiles have slightly lower values in the inner region, remaining more uniform with typical values of V_ϕ ~ 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 M_⊙ in the central embryo), the stress applied by the gravitational collapse and outflowing motions to the B-field only result in rather low with values <0.6 (equivalent to typical dispersion of the B-field angle along the line of sight ≲ 30°). A particular case is the simulation R9, which displays a decreasing V_ϕ with time in the inner ~500 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_ϕ 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.

To better understand how the circular variance is related to the column density, we compute the mean value of V_ϕ 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_ϕ develop at higher column densities for the simulations with a higher degree of magnetization (purple lines indicating μ = 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 ~ 10²⁴ cm⁻², where the simulation R3 have values of V_ϕ > 0.6, while simulation R9 barely surpasses 0.2.

Fig. 7 Mean circular variance Vϕ per logarithmic bin of column density. From left to right: results at msink = 0.1, 0.2, and 0.3 M⊙.

4 Synthetic maps of the polarized dust emission

We produced perfectly integrated synthetic observations of the polarized dust emission for the central 8000 au, observed at a distance of 250 pc, and at a resolution of 5 au per pixel (see Table 2). Figure 8 presents synthetic maps of the polarization fraction observed at 1.3 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_frac in Fig. 9. We present similar maps, but for λ = 0.8 and 3.0 mm (Figs. A.1 and A.2), as well as similar radial profiles (Figs. A.3 and A.4).

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 Figs. 9 and A.3 and A.4) 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, A.1 and A.2, and middle row panels of Figs. 9, A.3 and A.4) 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, A.3 and A.4). The highest polarization fractions are found in the external regions, typically reaching 5–10%, and gradually decrease in a nonmonotonical manner toward the central regions. More precisely, Fig. 8 shows that p_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_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 Figs. 6 and 8 suggests that the polarization fraction drops in the lines of sight where the B-field happens to have a large variance (V_ϕ). 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 M_⊙, 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. A.1) at the time when the protostellar mass reached 0.3 M_⊙, consistent with the scenario proposed by Hull et al. (2017), but without invoking the presence of initial turbulence.

Fig. 8 Results from the synthetic observation at λ = 1.3 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.

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 (ϕ_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 (ϕ_B, Eq. (3)). We compute the difference between model and observations as the following difference: Δϕ = ϕ_pol − ϕ_B, and the discrepancy as the absolute value of the difference, |Δϕ|.

Fig. 9 Radial profiles of the synthetic observation at λ = 1.3 mm. From top to bottom: total dust emission (Stokes I), total linearly polarized dust emission (PI), and polarization fraction (pfrac). From left to right: results at msink = 0.1, 0.2, and 0.3 M⊙.

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

We built difference maps between the model magnetic field line orientations and the reconstructed magnetic field orientations from polarized dust emission observed at λ = 0.8, 1.3, and 3.0 mm (see Figs. A.5, A.6 and A.7, respectively). We note that there is generally good agreement between the inferred (ϕ_pol) and the actual (ϕ_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 |Δϕ| 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 λ = 0.8 mm, which is the frequency that displays the lowest recovery rates, the magnetic field is typically recovered within an error of less than 15° (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_sink = 0.3 M_⊙ 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° are associated to higher circular variances (the recovery rate for the lines of sight with V_ϕ < 0.5 is 99.7%). Table 3 summarizes the mean recovery rates for the synthetic maps at λ = 0.8 mm at different evolutionary stages and at different scales, averaged for all the projections and simulations. The detailed results are given in Table A.1. 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). Tables A.3 and A.4 show even better results at longer wavelengths (λ = 1.3 mm and λ = 3.0 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_norm = PI/PI_peak) at λ = 0.8 mm. We use the same values of the contour levels of Fig. A.1. 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. This behavior is seen for the three studied wavelengths as shown in Fig. 11, which depicts the distribution of the discrepancy angle as a function of the polarized intensity.

Finally, note that the models analyzed here implement non-ideal 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 ~ 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.

Fig. 10 B-field recovery rate for all the synthetic observations at λ = 0.8 mm. Individual points correspond to those of Table A.1, while the cross symbols and vertical bars correspond to the averaged values and their standard variation for the two cases defined in Table 3, respectively: moderate (μ = 3.33) and weak (μ > 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.

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 |∆ϕ|) of higher or equal to 30° 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 (μ = 3.33) display percentages of higher than 10% over larger regions. Simulations withθ= 90° start to show regions of high discrepancy at a much larger radius. Once the sinks in the simulations have reached 0.3 M_⊙, most of the simulations show a region of high discrepancy within a distance of 30–100 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 |∆ϕ| as a function of the circular variance V_ϕ at λ = 1.3 mm in Figs. 13 and 14. Analogous results for λ = 0.8 mm and λ = 3.0 mm are shown in Figs. A.8 and A.9, and Figs. A.10 and A.11, respectively. All the results correspond to the stacked values for the three projections (y–z, x–z and x – y) of Fig. A.6 (Figs. A.5 and A.7 for λ = 0.8 mm and λ = 3.0 mm). The contour lines corresponds to the 2D histogram contour levels of the distribution of the lines of sight at 10⁵ (yellow), 10⁴ (red), 10³ (purple), and 10² (black) counts, while the color-coded background displays the mean column density for a given bin of V_ϕ and |∆ϕ|. We present in Fig. 13 the distribution for the more standard conditions (no initial turbulence andθ= 30°, 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 μ), and increasing level of rotation (increasing values of β_rot). In all the simulations, most of the points in the map (contours higher or equal to 10⁴ counts) are characterized by low levels of |∆ϕ| (typically ≤ 10°) and low to moderate values of V_ϕ, 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 (〈N〉 < 10²⁵ cm⁻³), 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. C.1).

On the other hand, at higher column densities, the behavior differs from what is expected, displaying |∆ϕ| values higher than what can be attributed to the dispersion only. Particularly, points with high |∆ϕ| (≥ 30°) and low values of V_ϕ (≤0.2) are puzzling, but they remain extremely rare. Even the abundance of points with |∆ϕ| > 30° and V_ϕ ≤ 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_sink = 0.3 M_⊙, 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 λ = 0.8 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 τ_{0.8 mm} ≥ 1, we obtain a recovery rate of only 42.5%. For τ_{0.8 mm} in the range 0.1–1, we obtain a recovery rate of 79.4%, while for τ_{0.8 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 (β = 0.04, M = 0, θ= 30°), the more magnetized case (R2) has initially (at m_sink = 0.1 M_⊙) higher levels of |∆ϕ|, 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 |∆ϕ|), 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 A.2 shows that simulations R1 and R6 display similar recovery rates, within the error bars, at different scales.

On the other hand, simulations with θ= 90° (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 |∆ϕ| values remain low, it develops a region where the discrepancy is maximal (even for values of V_ϕ lower than 0.2), reaching values closer to 90°.

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 β_rot in the more magnetized case (μ = 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).

Fig. 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 λ = 0.8, 1.3, and 3.0 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.

Fig. 12 Percentage of pixels with |∆ϕ| > 30° at λ = 1.3 mm in logarithmic radial bins.

Fig. 13 Distribution of the discrepancy between the B angle inferred from the synthetic polarized dust emission at λ = 1.3 mm and the mean orientation of the B lines in the simulation as a function of the circular variance Vϕ for the simulations with standard conditions. All these simulations do not include any initial turbulence (M = 0) and have an initial inclination angle between the initial rotation axis and the magnetic field 8 of 30°. The contour lines show the smoothed 2D histogram contour levels at 105 (yellow), 104 (red), 103 (purple), and 102 (black) counts. The color-coded background corresponds to the mean column density.

Fig. 14 Distribution of |∆ϕ| at λ = 1.3 mm as a function of the circular variance Vϕ 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 (M = 1), while the simulations R7 and R9 have an initial angle between the rotation axis and the magnetic field axis of θ= 90°.

Fig. 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 λ = 0.8, 1.3 and 3.0 mm). Top row: 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 Sect. 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_frac) as a function of the circular variance (V_ϕ) 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). Figure 15 shows 2D histograms of p_frac versus V_ϕ 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 Int. XIX 2015a; Planck Collaboration XX 2015b). 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 A.2.

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⁶ cm⁻³ and temperatures are quite low (T < 20 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 (~100 μ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 ~100 au. Measurements of the line-of-sight-averaged magnetic field line orientation using the polarized dust emission are precise to <15° (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 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 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^°) 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.

Supplementary Material

Higher resolution figures Access here

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).

Appendix A Additional material

In this section, we provide complementary material on synthetic observations and their comparisons to properties measured from the models. Figures A.1 and A.2 show the maps from the synthetic observations of the polarized dust emission of the 27 models, at λ = 0.8 and 3.0 mm. Figures A.3 and A.4 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 λ = 0.8 and 3.0 mm, respectively.

Figures A.5, A.6, and A.7 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 λ = 0.8, 1.3, and 3.0 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 A.8 and A.9 (Figures A.10 and A.11) show the distribution of the discrepancy between the inferred orientation of the magnetic field and the true mean orientation at λ = 0.8 mm (λ = 3.0 mm) for the simulations with standard conditions and the ones starting from less prototypical conditions (either M = 1 or θ = 90^°), respectively.

Table A.1 gives the detailed recovery rates at λ = 0.8 for all the simulations (ID), projections (los), and evolutionary stages (msink) 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 A.2. This table shows the mean recovery rates for simulations R1 and R6. These two otherwise identical simulations (μ = 3.33, β_rot = 0.01 and θ= 30°), differing only in the degree of turbulence (R6 has M = 1), show that the degree of turbulence does not have a strong impact on the recovery rate.

Finally, Tables A.3 and A.4 show the mean recovery rates separated into two cases (moderate and low magnetization degree) similar to Table 3, but at λ = 1.3, and 3.0 mm, respectively.

Fig. A.1 Results from the synthetic observation at λ = 0.8 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.

Fig. A.2 Results from the synthetic observation at λ = 3.0 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.

Fig. A.3 Radial profiles of the synthetic observation at λ = 0.8 mm. From top to bottom: total dust emission (Stokes I), total linearly polarized dust emission (PI), and polarization fraction (pfrac).

Fig. A.4 Radial profiles of the synthetic observation at λ = 3.0 mm. From top to bottom: total dust emission (Stokes I), total linearly polarized dust emission (PI), and polarization fraction (pfrac).

Fig. A.5 Difference maps ∆ϕ at λ = 0.8 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.

Fig. A.6 Difference maps ∆ϕ at λ = 1.3 mm. The structure of the figure is the same as for Fig. A.5.

Fig. A.7 ∆ϕ at λ = 3.0 mm. The structure of the figure is the same as for Figs. A.5 and A.6.

Fig. A.8 Distribution of the discrepancy |∆ϕ| between the B angle inferred from the synthetic polarized dust emission at λ = 0.8 mm and the mean orientation of the B lines in the simulation as a function of the circular variance Vϕ for the simulations with standard conditions. None of these simulations includes any initial turbulence (M = 0) and all have an initial inclination angle between the initial rotation axis and the magnetic fieldθof 30°. The contour lines show the smoothed 2D histogram contour levels at 105 (yellow), 104 (red), 103 (purple), and 102 (black) counts. The color coded background corresponds to the mean column density.

Fig. A.9 Distribution of |∆ϕ| at λ = 0.8 mm as a function of the circular variance Vϕ for the simulations with nonstandard conditions. The contours and the background image are the same as in Fig. A.8. Simulations R6 and R8 include turbulence (M = 1), while the simulations R7 and R9 have an initial angle between the rotation axis and the magnetic field axis of θ = 90°.

Fig. A.10 Distribution of the discrepancy between the B angle inferred from the synthetic polarized dust emission at λ = 3.0 mm and the mean orientation of the B lines in the simulation as a function of the circular variance Vϕ for the simulations with standard conditions. All these simulations do not include any initial turbulence (M = 0) and have an initial inclination angle between the initial rotation axis and the magnetic field θ of 30°. The contour lines show the smoothed 2D histogram contour levels at 105 (yellow), 104 (red), 103 (purple), and 102 (black) counts. The color coded background corresponds to the mean column density.

Fig. A.11 Distribution of |∆ϕ| at λ = 3.0 mm as a function of the circular variance Vϕ for the simulations with non standard conditions. The contours and the background image are the same as in Fig. A.10. Simulations R6 and R8 include turbulence (M = 1), while the simulations R7 and R9 have an initial angle between the rotation axis and the magnetic field axis of θ = 90°.

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 (|∆ϕ| > 30°) and low disorganization levels (V_ϕ < 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 |∆ϕ| as a function of the opacity and the projected radius at the available wavelengths. In Fig. B.1 we present the stacked (for all simulations, projections and time-steps) distribution of the discrepancy |∆ϕ| as a function of the intrinsic level of disorganization of the field lines (V_ϕ), the opacity (τ) and as a function of the radial distance (r). The colored region corresponds to the subset of points where |∆ϕ| > 30° and ≤ 0.5, while the gray background corresponds to all the points (given for reference). This figure shows that the bulk of the points at |∆ϕ| ~ 30 – 40° is still consistent with the intrinsic disorganization level (σ_ϕ ~ 30° when V_ϕ ~ 0.5), and that they are likely produced in regions disturbed by the outflows at a distance of 200 – 1000 au. Lines of sight with extremely high discrepancy levels |∆ϕ| > 70° are associated to regions located within a radius of roughly 100 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 B.2 displays the same plots as Fig. B.1, but for the perfectly aligned dust grains hypothesis. From this figure we find that the number of lines of sight in the region |∆ϕ| > 30° and V_ϕ < 0.5 is reduced by 25, 22, and 25% at λ = 0.8, 1.3, and 3.0 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 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).

Fig. B.1 Two-dimensional histograms for the RAT case of the stacked maps (all simulations, time-steps, and projections) showing the ∆ϕ as a function of Vϕ, vs. optical depth and vs. the radial distance r ( from left to right), at three wavelengths λ = 0.8, 1.3, 3.0 mm (top to bottom).

Fig. B.2 Same as Fig. B.1 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_ϕ), defined in Sect. 2.3, and the angular dispersion (σ_ϕ) are related we simulated several ensembles of random angles and computed the resulting V_ϕ and σ_ϕ. To avoid any bias related to the position angle, we set 180 regularly spaced central angles (ϕ_c) from −90° to 89°. Secondly, around each central angle we defined 18 opening angles Δθ, varying from 5° to 90° and we randomly generated an ensemble of N = 500 uniformly distributed angles in the interval [ϕ_c – Δθ, ϕ_c + Δθ).

As the mean angle does not represent the true orientation when dealing with bidirectional data, we define the mean orientation 〈ϕ〉 as in Eq 3. The dispersion is thus defined as follows: (C.1)

where we use the smallest difference with respect to the mean mean orientation (|(ϕ_i − 〈ϕ〉)| ≤ 90°).

We computed V_ϕ and σ_ϕ 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_ϕ and the angular dispersion in Fig. C.1. 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 ∆θ of each ensemble. We note that the angular dispersion is necessarily smaller than ∆θ, which only represents the maximum allowed variation with respect to the central angle ϕ_c. Figure C.1 shows that σ_ϕ increases monotonically with V_ϕ (and with ∆θ), but its maximum value is limited to less than 55°.

Fig. C.1 Relation between the circular variance and the angular dispersion.

References

All Tables

All Figures

Fig. 1 Sink-mass evolution as a function of the time since the formation of the sink tage.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.

In the text

Fig. 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 msink = 0.1, 0.2, and 0.3 M⊙ (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.

In the text

	Fig. 3 Mean magnetic pressure Pmag 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.
In the text

	Fig. 4 Circular variance Vϕ of the B-fleld 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.
In the text

	Fig. 5 Radial mean column density profiles. From left to right: profiles at msink = 0.1, 0.2, and 0.3 M⊙.
In the text

	Fig. 6 Radial mean profiles of the circular variance Vϕ. From left to right: profiles at msink = 0.1, 0.2, and 0.3 M⊙.
In the text

	Fig. 7 Mean circular variance Vϕ per logarithmic bin of column density. From left to right: results at msink = 0.1, 0.2, and 0.3 M⊙.
In the text

Fig. 8 Results from the synthetic observation at λ = 1.3 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.

In the text

	Fig. 9 Radial profiles of the synthetic observation at λ = 1.3 mm. From top to bottom: total dust emission (Stokes I), total linearly polarized dust emission (PI), and polarization fraction (pfrac). From left to right: results at msink = 0.1, 0.2, and 0.3 M⊙.
In the text

Fig. 10 B-field recovery rate for all the synthetic observations at λ = 0.8 mm. Individual points correspond to those of Table A.1, while the cross symbols and vertical bars correspond to the averaged values and their standard variation for the two cases defined in Table 3, respectively: moderate (μ = 3.33) and weak (μ > 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.

In the text

Fig. 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 λ = 0.8, 1.3, and 3.0 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.

In the text

	Fig. 12 Percentage of pixels with |∆ϕ| > 30° at λ = 1.3 mm in logarithmic radial bins.
In the text

Fig. 13 Distribution of the discrepancy between the B angle inferred from the synthetic polarized dust emission at λ = 1.3 mm and the mean orientation of the B lines in the simulation as a function of the circular variance Vϕ for the simulations with standard conditions. All these simulations do not include any initial turbulence (M = 0) and have an initial inclination angle between the initial rotation axis and the magnetic field 8 of 30°. The contour lines show the smoothed 2D histogram contour levels at 105 (yellow), 104 (red), 103 (purple), and 102 (black) counts. The color-coded background corresponds to the mean column density.

In the text

	Fig. 14 Distribution of |∆ϕ| at λ = 1.3 mm as a function of the circular variance Vϕ 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 (M = 1), while the simulations R7 and R9 have an initial angle between the rotation axis and the magnetic field axis of θ= 90°.
In the text

Fig. 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 λ = 0.8, 1.3 and 3.0 mm). Top row: 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 Sect. 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.

In the text

Fig. A.1 Results from the synthetic observation at λ = 0.8 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.

In the text

Fig. A.2 Results from the synthetic observation at λ = 3.0 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.

In the text

	Fig. A.3 Radial profiles of the synthetic observation at λ = 0.8 mm. From top to bottom: total dust emission (Stokes I), total linearly polarized dust emission (PI), and polarization fraction (pfrac).
In the text

	Fig. A.4 Radial profiles of the synthetic observation at λ = 3.0 mm. From top to bottom: total dust emission (Stokes I), total linearly polarized dust emission (PI), and polarization fraction (pfrac).
In the text

	Fig. A.5 Difference maps ∆ϕ at λ = 0.8 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.
In the text

	Fig. A.6 Difference maps ∆ϕ at λ = 1.3 mm. The structure of the figure is the same as for Fig. A.5.
In the text

	Fig. A.7 ∆ϕ at λ = 3.0 mm. The structure of the figure is the same as for Figs. A.5 and A.6.
In the text

Fig. A.8 Distribution of the discrepancy |∆ϕ| between the B angle inferred from the synthetic polarized dust emission at λ = 0.8 mm and the mean orientation of the B lines in the simulation as a function of the circular variance Vϕ for the simulations with standard conditions. None of these simulations includes any initial turbulence (M = 0) and all have an initial inclination angle between the initial rotation axis and the magnetic fieldθof 30°. The contour lines show the smoothed 2D histogram contour levels at 105 (yellow), 104 (red), 103 (purple), and 102 (black) counts. The color coded background corresponds to the mean column density.

In the text

	Fig. A.9 Distribution of |∆ϕ| at λ = 0.8 mm as a function of the circular variance Vϕ for the simulations with nonstandard conditions. The contours and the background image are the same as in Fig. A.8. Simulations R6 and R8 include turbulence (M = 1), while the simulations R7 and R9 have an initial angle between the rotation axis and the magnetic field axis of θ = 90°.
In the text

Fig. A.10 Distribution of the discrepancy between the B angle inferred from the synthetic polarized dust emission at λ = 3.0 mm and the mean orientation of the B lines in the simulation as a function of the circular variance Vϕ for the simulations with standard conditions. All these simulations do not include any initial turbulence (M = 0) and have an initial inclination angle between the initial rotation axis and the magnetic field θ of 30°. The contour lines show the smoothed 2D histogram contour levels at 105 (yellow), 104 (red), 103 (purple), and 102 (black) counts. The color coded background corresponds to the mean column density.

In the text

	Fig. A.11 Distribution of |∆ϕ| at λ = 3.0 mm as a function of the circular variance Vϕ for the simulations with non standard conditions. The contours and the background image are the same as in Fig. A.10. Simulations R6 and R8 include turbulence (M = 1), while the simulations R7 and R9 have an initial angle between the rotation axis and the magnetic field axis of θ = 90°.
In the text

	Fig. B.1 Two-dimensional histograms for the RAT case of the stacked maps (all simulations, time-steps, and projections) showing the ∆ϕ as a function of Vϕ, vs. optical depth and vs. the radial distance r ( from left to right), at three wavelengths λ = 0.8, 1.3, 3.0 mm (top to bottom).
In the text

	Fig. B.2 Same as Fig. B.1 but for the case where dust grains are perfectly aligned.
In the text

	Fig. C.1 Relation between the circular variance and the angular dispersion.
In the text