A statistical analysis of dust polarization properties in ALMA observations of Class 0 protostellar cores

Recent observational progress has challenged the dust grain-alignment theories used to explain the polarized dust emission routinely observed in star-forming cores. In an effort to improve our understanding of the dust grain alignment mechanism(s), we have gathered a dozen ALMA maps of (sub)millimeter-wavelength polarized dust emission from Class 0 protostars, and carried out a comprehensive statistical analysis of dust polarization quantities. We analyze the statistical properties of the polarization fraction P_frac and dispersion of polarization position angles S. More specifically, we investigate the relationship between S and P_frac as well as the evolution of the product S*P_frac as a function of the column density of the gas in the protostellar envelopes. We find a significant correlation in the polarized dust emission from protostellar envelopes seen with ALMA; the power-law index differs significantly from the one observed by Planck in star-forming clouds. The product S*P_frac, which is sensitive to the dust grain alignment efficiency, is approximately constant across three orders of magnitude in envelope column density. This suggests that the grain alignment mechanism producing the bulk of the polarized dust emission in star-forming cores may not depend systematically on the local conditions such as local gas density. Ultimately, our results suggest dust alignment mechanism(s) are efficient at producing dust polarized emission in the various local conditions typical of Class 0 protostars. The grain alignment efficiency found in these objects seems to be higher than the efficiency produced by the standard RAT alignment of paramagnetic grains. Further study will be needed to understand how more efficient grain alignment via, e.g., different irradiation conditions, dust grain characteristics, or additional grain alignment mechanisms can reproduce the observations.


Introduction
Magnetic fields have been considered to play a key role in the formation of molecular clouds and in the regulation of star formation (Shu et al. 1987;McKee et al. 1993;McKee & Ostriker 2007). For example, fields are partially responsible for setting the star formation rate (Krumholz & Federrath 2019), as the gas motions tend to follow the orientations of magnetic fields, whose strengths can regulate the gravitational collapse of these structures (Mouschovias & Ciolek 1999). Past observations of molecular clouds have shown that the magnetic field seems to be a key player in the formation of parsec-scale density structures (Planck Collaboration et al. 2016;Soler 2019;Seifried et al. 2020), and appears regulate star formation inside these structures (Li et al. 2017). One of the main ways to characterize the spatial distribution of magnetic fields is to observe the polarized thermal emission from dust grains. Indeed, since dust grains are not perfectly spherical, they tend to align themselves with the ambient magnetic field under some conditions (Lazarian 2007;Andersson et al. 2015), resulting in polarized thermal emission that can be used to infer the magnetic field orientation integrated along the line of sight. This linear polarization emanating from this dust grain population is orthogonal to the magnetic field component projected on the plane of the sky.
Observations of magnetic fields via polarized dust emission are still subject to caveats due to the strong dependence of grain alignment on the local environmental conditions. Understanding the impact of the key factors enabling dust grain alignment via the Radiative Alignment Torques (RATs) theory-such as the degree of anisotropy of the radiation field, the dust grain size distribution, the gas temperature, and the density distribution-has been Article number, page 1 of 38 arXiv:2009.07186v1 [astro-ph.GA] 15 Sep 2020 A&A proofs: manuscript no. STAPs_arxiv the focus of numerical works where radiative transfer was performed on magneto-hydrodynamic (MHD) simulations (Padoan et al. 2001;Bethell et al. 2007;Pelkonen et al. 2009;Brauer et al. 2016). One of their main goals was to investigate the widespread phenomenon of depolarization, i.e., the drop in the ratio of linearly polarized dust emission to the total intensity emission toward high density zones in molecular clouds and cores; this is the so-called "polarization hole" phenomenon. Single-dish observations of molecular clouds (Poidevin et al. 2013;Fissel et al. 2016), single-dish observations of a starless core (Alves et al. 2014), and high-resolution interferometric observations of Class 0 protostellar cores (Hull et al. 2014;Galametz et al. 2018) found a significant decrease of the polarization fraction with increasing column density, and interpreted this drop as either depolarization caused by disorganized magnetic field lines smeared out in the synthesized beam (i.e., the resolution element of the observations), or as a possible loss of alignment efficiency of the dust particles caused by a lack of irradiation and/or changes in dust grain characteristics toward high column density regions. This depolarization phenomenon was analyzed in the scope of several possible physical explanations: the collisional de-alignment of dust grains due to high gas temperature and density (Reissl et al. 2020); the reddening of the radiation field when reprocessed during its propagation (Lazarian 2007); the change in grain size and shape due to coagulation and formation of icy mantles (Juárez et al. 2017); the lack of the necessary anisotropy in the radiation field as a result of high optical depth , studied the drop of polarization degree in dense regions of Bok globules); and the level of disorder of the magnetic field lines caused by turbulence (Falceta-Gonçalves et al. 2008). Most of the time, high angular resolution observations of dust polarization revealed that the drops of polarization fraction in observations of cores at coarse angular resolution were partly explained by beam smearing, i.e., the fact that the fluctuations of magnetic fields in the plane of the sky could not be resolved in high density zones. However, while these higher angular resolution observations did tend to detect polarized emission in the holes seen in the low-resolution data, these same high resolution observations saw their own polarization holes at smaller scales, as pointed out in Galametz et al. (2018).
In observations, this depolarization effect was also quantified thanks to the statistical analysis of both the polarization fraction and the local dispersion of polarization angles (Planck Collaboration et al. 2015a,b). Assuming perfect alignment efficiency of the dust grains with magnetic field lines, both quantities are correlated with the level of disorder in the magnetic field, and thus both vary as a function of the amount of local fluctuations of the magnetic field. The polarization fraction is sensitive to the cancellation of polarization along the line of sight and hence to the fluctuations of the apparent magnetic field along the line of sight. Conversely, the dispersion of polarization angles provides the fluctuation of the apparent magnetic field orientations in the plane of the sky. Assuming the fluctuations of the magnetic field lines are isotropic, the product of these two quantities gives access to the dust grain alignment efficiency: Planck Collaboration et al. (2018) analysed these statistical estimators as a function of column density from the diffuse interstellar medium (ISM) to molecular clouds, probing column density up to 10 22 cm −2 . They found no variation of dust grain alignment efficiency with varying conditions typi-cal of these environments. More recently, Reissl et al. (2020) applied this statistical tool to simulations of the diffuse ISM in order to quantify the relative influence that radiative torque intensity and gas pressure have on the dust grain alignment efficiency. They found no significant differences in dust grain alignment efficiency when analyzing polarization from perfectly aligned dust grains versus those aligned by RATs in these environments embedded in the interstellar radiation field.
Given the uncertain validity of RATs in high density environments where irradiation is much less homogeneous and is shifted to long wavelengths, the dense parts of protostellar cores represent regions of interest. While ALMA has recently produced a large number of high-sensitivity observations of young stellar objects, the thermal dust emission emitted by the youngest sources, known as prestellar cores, is heavily spatially filtered by ALMA and thus hardly detected (Dunham et al. 2016), rendering investigations of their polarized dust emission at high spatial resolution challenging. However, once these cores initiate their gravitational collapse, because their densest regions become warm enough to dissociate H 2 , a compact structure forms around the nascent embedded protostar, which enables interferometric observations. These youngest protostellar objects, known as Class 0 protostars (André et al. 2000, are engaged in a short but vigorous accretion phase during which the central protostar will gather most of its final mass, triggering also ejection of material in the form of bipolar outflows visible in molecular emission lines. These sources are ideal for our study because during this phase, most of the thermal dust emission is emitted by the envelope surrounding the central embryo; during the (later) Class I phase, the envelope has already been largely accreted/dissipated. Here, we focus on the ALMA dust polarization observations of Class 0 envelopes to optimize the number of detections. Most of these recent observations have shown that specific regions of the cores such as the walls of the bipolar outflow cavities, potential magnetized accretion streamers, and core equatorial plane are preferentially polarized (Hull et al. 2017b,a;Maury et al. 2018;Sadavoy et al. 2018a,b;Kwon et al. 2019;Le Gouellec et al. 2019;Takahashi et al. 2019;Hull et al. 2020). The specific locations of the recovered polarized emission raise questions regarding the local conditions required to align a significant fraction of dust grains along magnetic lines, and thus to produce the level of polarized emission observed. Our goal in this paper is to investigate the statistical behavior of polarized dust emission observed with ALMA within protostellar envelopes, by adapting and applying some of the tools previously developed to characterize polarized dust emission at cloud scales.
The paper is structured as follows. In Section 2, we introduce the statistical tools we use to analyse the polarized dust emission detected by ALMA in young protostars, and we present the sample of ALMA observations of Class 0 objects that are the focus of this statistical study. We present the methodology and results of the statistics in Section 3. Finally, we discuss the results we obtain regarding the dust grain alignment in young protostellar objects in Section 4, along with comparisons between the ALMA observations and synthetic observations of MHD simulations. We draw our conclusions in Section 5.

Statistical tools
Our objective is to characterize the polarized emission emanating from Class 0 protostellar cores, at envelope scales, targeting the emission from circumstellar material at radii of ∼ 10 − 2000 au. The properties of the linear polarization of thermal dust emission are expressed by the Stokes parameters Q and U . Stokes I represents the total intensity. We will denote the polarized intensity P (defined as P = Q 2 + U 2 , which we systematically debias, see Section 3.1), the polarization fraction P frac (defined as P frac = P/I), and the polarization position angle φ (defined as φ = 0.5 arctan U Q ). In the diffuse ISM and molecular clouds, Planck Collaboration et al. (2015aCollaboration et al. ( ,b, 2018 found a correlation between the local dispersion of the polarization position angle and the polarization fraction. Note that a similar correlation was found in Alves et al. (2008), using optical background-starlight polarization observations of the Pipe Nebula. The polarization angle dispersion function S, which quantifies the local (non)-uniformity of the polarization angle, is defined as follows: where the angle dispersion is calculated at a given position r and for a given neighborhood δ, which is also known as the lag. The lag describes the area over which the dispersion of polarization angles is derived, and thus corresponds to the characteristic length scale at which we quantify the disorganization of polarization position angles. The computation is performed on N neighboring pixels contained in an annulus centered on r, having inner and outer radii of δ/2 and 3δ/2, respectively; each of the N pixels is indexed by i, and located at r + δ i (Planck Collaboration et al. 2018). Planck Collaboration et al. (2018) developed an analytical model (briefly described in Appendix A) that relates the two quantities S and P frac . They found, among other results, that S ∝ P frac −1 in the diffuse ISM and molecular clouds: this correlation is shown as a red solid line in our plots. Exploring the evolution of the quantity S × P frac -which is a proxy for grain alignment efficiency-as a function of the column density and the dust temperature, Planck Collaboration et al. (2018) did not detect a significant drop of efficiency with increasing column density. We apply a modified version of this technique to ALMA observations assuming the dust grains are aligned with the ambient magnetic field at the typical scales of a protostellar core.
The dispersion of the polarization position angles S gives us information about the level of disorder in the magnetic field projected in the plane of the sky; the higher the value of S, the more disorganized the apparent magnetic field. S will saturate at the value of π/ √ N , as (φ = 0 • ) ≡ (φ = 180 • ) in a polarization map. Note that here we distinguish between the disorganization of the apparent magnetic field lines as seen by the observer and the actual (3-dimensional) turbulent component of the magnetic field. Indeed, if a given, moderately turbulent magnetic field is oriented closer to the line of sight, the observed dispersion S will be larger. In contrast, a uniform magnetic field will have dispersion S close to zero, regardless of the line of sight. These facts limit the capability of S to trace the turbulent component of the magnetic field. However, we should also note that the value of the polarized intensity P (and thus P frac ) directly depends on the orientation of the magnetic with respect to the line of sight. Thus, in the extreme line-of-sight cases where S reaches high values, P may drop below the detection limit.
The polarization fraction P frac is another tool linked with the disorganization of the magnetic field. A disorganized magnetic field along the line of sight will result in a low value of P frac as seen by the observer. Consequently, assuming an isotropic turbulent component of the magnetic field, and given the caveats listed above, S and P frac are directly linked to level of disorder in the magnetic field lines in a core.
The aim of the study presented here is to search for, compare, and interpret the possible physical causes for a correlation between S and P frac in Class 0 protostars, toward which polarized dust emission was observed with ALMA in recent years.

ALMA observations of polarized dust emission in
Class 0 protostellar cores In its polarized mode, ALMA produces visibility measurement sets of the three Stokes parameters I, Q, and U , which can be imaged and combined to produce maps of the polarized dust emission. We gathered publicly available ALMA dust polarization observations toward nearby, lowand intermediate-mass Class 0 protostars. Since the statistical tools we use require a large number of statistically independent measurements at the typical scale of the object studied, we selected observations with the most extended polarized dust emission. The regions of interest in these protostellar cores correspond to the inner envelope scales (∼10-2000 au). Therefore, we selected the ALMA datasets whose polarized dust emission was observed with combinations of sensitivity and angular resolution that allow us to detect low levels of polarized emission beyond the peak of continuum emission, at these inner envelopes typical scales.
We present the resulting sample in Table 1. We use the three Stokes maps provided by the authors of the corresponding publications (see Table 1) to create the polarized emission maps. In the case of NGC1333 IRAS4A, however, because these data were not yet published at the time we started our investigations, we calibrated and imaged these observations ourselves. We produced the polarized dust continuum images using the task tclean in version 5.4 of CASA (McMullin et al. 2007). We applied four rounds of consecutive phase-only self-calibration, using the total intensity (Stokes I) solutions as the model, with a Briggs weighting parameter of robust = 1. The three Stokes parameters I, Q, and U were cleaned separately after the last round of self-calibration using an appropriate residual threshold and number of iterations. In order to calculate appropriate thresholds for the data (see the method developed in Section 3.1), we require an homogeneous level of noise across the individual fields of view, and thus we do not perform any primary beam correction at this step of the analysis. However, the total intensity maps are primary beam corrected before deriving the column density whose ranges are reported in Table 1.
It is crucial, when building polarized dust emission maps from the combination of the Stokes maps, to have a robust assessment of the rms noise levels in each of these maps. This is particularly important to consider in our statistical measurements so that we do not introduce noise bias, since Article number, page 3 of 38 A&A proofs: manuscript no. STAPs_arxiv values of polarization fraction P frac can be affected significantly when dividing by Stokes I values that are uncertain. Here, we compute the rms noise values in each of the three Stokes maps I, Q, and U (σ I , σ Q , and σ U respectively) by measuring the root mean square in an area without strong emission. We notice that typically σ Q ≈ σ U , so we use a single value σ P ≡ σ Q ≈ σ U . We present in Figure 1 the distribution of P frac and polarization position angles φ in the region where the Stokes I is > 5 σ I , from all individual maps of all sources at each wavelength. In these histograms, the uncertainties in P frac , and φ (in radians) are showed as shaded areas, and are calculated as follows: Finally, assuming the dust emission recovered in the ALMA observations (at scales of 10-2000 au) is optically thin, we calculate the column density from the total intensity dust emission maps as follows: where S ν the flux density measured, d is the distance to the source (see Table 1), B ν (T d (r)) is the Planck function at the frequency ν of our observations for dust of a given temperature T d (r) (see below), κ ν is the opacity at a specific wavelength taken from Ossenkopf & Henning (1994), m H is the mass of a hydrogen atom, µ H2 is the mean molecular weight per hydrogen molecule (µ H2 = 2.8 for gas composed of 71% hydrogen, 27% helium, and 2% metals by mass; Kauffmann et al. 2008), and A is the area over which we calculate the flux density. We assume a gas-to-dust ratio of 100. The value of the dust temperature at a radius r from the position of the protostellar embryo (assumed to coincide with the peak position of the dust continuum emission in the ALMA Stokes I map), can be estimated assuming that only the central protostellar object heats the dust in the inner envelope, following Terebey et al. (1993) and Motte & André (2001): where L int is internal luminosity of the protostar, which is directly linked to the protostellar accretion luminosity. While for some of these sources-Serpens Emb 8, Serpens Emb 8(N), B335, L1448 IRS2, and NGC1333 IRAS4A-the  Table B.1). We find that the ALMA observations are sensitive to material in the inner envelope with typical column densities ∼ 10 22 − 10 25 cm −2 .
The individual ranges of column densities probed in each map are reported in Table 1.
Finally, while most of the polarized dust emission toward the sample of sources we present is caused by thermal emission of dust grains aligned with respect to the magnetic field, this may not be the case where the dust emission becomes optically thick and where the radiation from dust is highly anisotropic (such as protoplanetary disks); in these regions the polarized dust emission can be caused by the selfscattering of thermal dust emission (Kataoka et al. 2015). Within the sample of sources we present, two of them have been clearly identified as having polarized dust emission due to self-scattering in their inner region; these are the two Ophiuchus sources, IRAS 16293A/B and VLA 1623A/B (Sadavoy et al. 2018a(Sadavoy et al. ,b, 2019. We estimate that only 2% of the pixels could be contaminated in IRAS16293A/B, whereas up to 40% of the pixels could contain polarized emission mostly due to self-scattering (based on the pattern of polarization position angles) in VLA1623A/B. Thus, we exclude these pixels from our analysis. Moreover, in our sample the dust emission is also optically thick in the inner 100 au region of IRAS4A (Ko et al. 2020): this represents < 1% of the pixels in both our observations at 1.3 and 0.87 mm. We also exclude these pixels from our analysis. Note that although Kwon et al. (2019) and Takahashi et al. (2019) disfavoured self-scattering as the cause of the linear polarization at the very center of the L1448 IRS2 and OMC3 MMS6 cores, respectively, we cannot rule it out. However, if scattering were present in these sources, it would only affect the few pixels at the peak of the dust continuum emission.

Applying Planck statistical tools to interferometric observations
We aim to apply the statistical tools developed for the analysis of the Planck maps of the polarized ISM to interferometric ALMA observations. We compare the statistical properties of dust polarization in the dense regions of protostellar cores with the properties found in larger-scale star-forming clouds. However, using interferometric data requires us to adapt the Planck collaboration's methods for investigating large-scale maps. For example, unlike ALMA observations, Planck observations are not affected by spatial filtering. We treat the ALMA polarization products as follows. We regrid the maps of the three Stokes parameters I, Q, and U to a Nyquist sampling, with exactly 4 pixels per beam in terms of area. We then calculate the polarization angle dispersion function S in each pixel i of the Stokes maps, with respect to each of its n = 8 nearest neighbouring pixels j as follows: Considering the sampling described above, the equivalent δ parameter (see also Equation 1) is approximately 1/2 of a beam width (comparable to the value chosen in Planck Article number, page 5 of 38 A&A proofs: manuscript no. STAPs_arxiv Collaboration et al. 2018). Note that the measured value of S scales with the pixel gridding. Indeed, at a given angular resolution, changing the gridding pattern (i.e., how many pixels a beam contains) to a fewer number of pixels per beam leads to a measurement of S that covers a larger area, and thus the lag δ is larger. This in turn causes us to quantify the disorganization of the magnetic field across a larger physical area. As the angular resolution of the observations is fixed, increasing the lag will cause S to increase, as we lose spatial coherence in the apparent magnetic field, which in turn causes an increase in the calculated level of disorganization in the apparent magnetic field. We perform the same analysis with different gridding and choose the value of 4 pixels per beam area in order to strike a balance between statistical accuracy (i.e., using a large number of points) and independence of the individual points. Finally, as explained in Section 2.1, the way the dispersion S is derived causes the distribution to saturate, i.e., a completely random distribution of polarization angles will produce a maximum value of S of π/ √ n ∼ 63 • (as we chose n = 8).
While Planck Collaboration et al. (2018) produced covariance maps and were able to assess finely the noise properties at different spatial scales, interferometric maps are severely affected by imaging systematics such as the limited dynamic range of the images. Furthermore, Stokes I images tend to be much more dynamic range limited than Stokes Q and U images. In addition, the sources lie at different distances, which leads us to probe different angular extents. Finally, the data we analyze are heterogeneous in their uv-coverage and sensitivity. Therefore, the noise is neither spatially homogeneous nor correlated in the ALMA maps. Consequently, we compute the rms noise σ in each Stokes map, using regions close to the observation pointing center but devoid of emission. This is how we define the signal-to-noise ratio (S/N) of polarized intensity, i.e., P/σ P . Note when creating these polarized emission maps, one must to correct for the bias that occurs at low S/N levels: to do so and to construct fully sampled P maps, we follow the method from Wardle & Kronberg (1974) (see also Hull & Plambeck 2015 for an application of this method to interferometric data).
We follow the method introduced in Planck Collaboration et al. (2018) to compute a pixel-selection criterion in order to test appropriately the correlation between S and P frac in our objects. This pixel-selection is a cutoff based on Stokes I, which allows us to remove the noise-biased data. We obtain this cutoff by analyzing the average S/N of the polarized intensity P , which typically increases with increasing Stokes I. When this average S/N of P , plotted as a function of the total intensity Stokes I for each dataset, meets the value S/N = 5, we use the corresponding value of Stokes I as the pixel-selection cutoff for the given dataset. We show an example in Figure 2 (top panel ), where we plot the S/N of the polarized intensity map as a function of the total intensity for the 1.3 mm observations of the B335 core. The vertical dotted line, which denotes this cutoff value of Stokes I, thus corresponds to the value of Stokes I where P/σ P I ≥ 5. We then take all points lying above this cutoff in Stokes I and form the sample to which we will apply our statistical method. Note that if this method provides a threshold limit of Stokes I below 5 σ I , we chose 5 σ I as the cutoff for the dataset.
As shown in the lower panel of Figure 2, the values of S × P frac diverge at Stokes I values lower than the aforementioned threshold, as a consequence of the noise-bias of S and P frac at low values of Stokes I. It is important to note that we have not performed any selection based on the P values; we select only on the I values in order to keep the pixels exhibiting low polarized intensity that contain the information of depolarization, which is essential for our statistical analysis. As an example, in Figure 3 we show the maps of Stokes I, polarized intensity P , polarization fraction P frac , and dispersion of polarization angles S from the 1.3 mm observations of the B335 core. The contours indicate both the threshold in Stokes I found with the method introduced above, as well as the 5σ I level. Similar maps of all sources can be seen in Figure B

Results from the statistical analysis of the polarized dust emission in protostellar envelopes
We present here the outcome of the statistical analysis of the polarized dust emission from the sample of ALMA observations presented in Table 1. In Figure B.3 we show he distribution of the dispersion of the polarization angles S as a function of the polarization fraction P frac in the 15 maps probing the dust emission in protostellar cores. In these plots, the running mean of P frac (shown as black points and line) shows the average trend and evolution of S with P frac . In particular, one can clearly see the area of the distribution affected by the saturation effect of S described above, and that the distribution is linear in the logarithmic two-dimensional (2D) space outside of this saturated area.
In each distribution, the points are coloured based on their Stokes I value. The relationship between S and P frac observed by Planck at cloud scales is reported in each diagram as a red line, for reference. We find a global trend similar to the Planck findings, with high values of polarization fraction P frac associated with low dispersion in polarization angles S in regions of faint Stokes I values. Conversely, we see high S and low P frac in regions with bright Stokes I. We list in Table 2 the values of the power law indexes α derived from the fit to the S ∝ P frac −α relation in each individual core. These values range from α = 0.523 ± 0.094 to α = 0.866 ± 0.040.
In order to take full advantage of the statistical power of our methodology and to discuss global properties of the polarized dust emission in Class 0 protostars, we have merged all data from each of the 15 ALMA observations. Figure 4 shows the merged distribution of S as a function of P frac , along with the linear fit previously described, which is defined by the two parameters α and f such that S = f /P frac α . Note that at low values of P frac , the distribution of S flattens because of the saturation of S for high dispersion values. This is an artifact arising from the definition of S, and thus these points should be excluded from the linear fit. To do so, we establish a threshold in P frac of 1.3%, indicated by the vertical dot-dashed grey line in Figure 4, which denotes the P frac level beyond which the distribution is linear. We calculate this threshold by determining where the α value from the linear fit would no longer have changed if we had moved the threshold up in polarization fraction. We obtain a power law index α = 0.79 ± 0.03, which is flatter than the results and the analytical correlation found with Planck observations at larger scales, where S ∝ P frac −1 . Merging all the ALMA data does not enable us to investigate the relation between S and P frac with respect to Stokes I because of the heterogeneous properties of the sources and observations (e.g., the wavelength of the observations). Thus, we use the column density (calculated as described  Top: Evolution of the S/N of the polarized intensity, i.e., P/σP , as a function of the total intensity Stokes I in Jy str −1 . The dot-dashed horizontal black line is at the value of P/σP = 5. The dotted vertical line is the selected cutoff in Stokes I described in Section 3.1. The dot-dashed vertical line is the 5σI value. The solid line is the running mean, which is calculated along the Stokes I; the shaded area represents ± the standard deviation of the Gaussian fit performed on each bin. Bottom: S × P frac for the selected pixels as a function of the total intensity. To the left of the cutoff in Stokes I (the red dotted line, plotted as in the top panel), the points are no longer plotted and the running mean turns in a translucent dashed line. Each color corresponds to an angular resolution: red is the original resolution, whereas blue and green are 4 × and 9 × lower resolution (in terms of beam area), respectively. Note that, as expected, one see that decreasing the resolution, and thus increasing the spatial length of the lag, causes on the dispersion S to increase as well, on average (see Section 2.1). Results from the correlations presented in Figure B.3. We list the values of the power law indexes α and associated uncertainties obtained from the linear fits, as well as the wavelength of observations λ, the number of points selected in each case, and the cutoff in Stokes I applied .
in Section 2.2) of the individual lines of sight in order to collect all the data points in a single plot. Figure 5 presents the variation of S, P frac , and S × P frac as a function of the local column density N H2 in the envelopes. Figure 6 shows the same distribution of points, but where all column density values are normalized to the maximum column density in the map including all optically thin lines of sight (e.g., excluding highly extinct lines of sight where polarized dust emission could be severely contaminated by self-scattering) 1 . We notice a global trend in the merged data shown in Figure 6 of increasing S and decreasing P frac with increasing column density. Figure B.4 presents the results of linear regressions on the trends of S and P frac as a function of N H2 in individual cores: the resulting polar-law indices and R-squared values suggest a significant decrease of P frac with increasing N H2 across the sample, and hints of increasing S with column density, although these latter trends are noisier. This suggests that the behavior of the disorganization of the apparent magnetic field evolves in the same way from the outer to the inner core, despite the widely varying ranges of column density in each core (see Table 1). 00. Bottom-Left. Polarization fraction P frac in colorscale, shown where I > 3σI . Bottom-Right. Dispersion S of polarization position angles in color scale; the pixel size corresponds to the pixel size considered in the statistics. The dashed white contour represents the 5σI level. The solid white contour represents the threshold level of Stokes I calculated as described in Section 3.1, above which the mean S/N of P > 5. The beam size is 1 . 14 × 0 . 90, with a position angle of 89.1 • .

Comparing the statistical properties of the polarized dust emission in protostellar cores with those in star-forming clouds
The statistical analysis of the polarized dust emission from the sample of 15 datasets analyzed (see Table 1) reveals a significant correlation between the dispersion of polarization angles S observed in the plane of the sky and the polarization fraction P frac measured in each line of sight in these 11 protostellar envelopes. However, the S ∝ P frac −0.79 relationship we find at core scales is shallower than the S ∝ P frac −1 relationship found at larger scales in the Planck observations of star-forming clouds (Planck Collaboration et al. 2018). Moreover, we obtain on average higher values of S and P frac than those found in the lower density molecular clouds probed by Planck. Here we discuss possible origins of the different polarization properties at protostellar scales versus cloud scales. We start by investigating the nature of the disorganized component of the magnetic field (Section 4.1.1); we then address the different intrinsic spatial scales of total intensity versus polarized emission (Section 4.1.2), and how interferometric filtering may affect observed polarization properties (Section 4.1.3).

What physics governs the disorganized component of the magnetic field?
The correlation between S and P frac is governed by the level of disorganization of the apparent magnetic field lines projected on the plane of the sky. The magnetic field is also linked with the kinematics of the gas, assuming the gas is angles S as a function of the polarization fraction P frac from all of the datasets merged together. The points were selected according the method developed in Section 3.1. The color scale represents the number density of points in the plot. The solid black line and black points represent the running mean of P frac ; the associated black error bars are ± the standard deviation of each bin. We plot the linear fit in purple, which is a linear regression. We take into account the saturation of S in the derivation of the linear fit by applying a threshold in polarization fraction, indicated by the vertical dot-dashed grey line. This threshold denotes the P frac level beyond which the distribution is linear. The solid red line corresponds to the Planck correlation from Planck Collaboration et al. (2018), which we scaled down to the highest angular resolution of our ALMA observations. The red shaded area extends up to the same Planck correlation, scaled down this time to the largest field of view of our ALMA observations. As we gather all of the ALMA observations at their various angular resolutions, this red shaded area encompasses all of the corresponding scalings of the Planck correlation. The two parameters f and α are derived from the linear fit, where the analytical correlation is as follows: S = f /P α frac . The histograms in the two little subplots show the distributions of the values of α and f values derived from a large number of randomly chosen sub-samples of points. We calculate the uncertainties in f and α as standard deviations of Gaussian fits to those histograms.
well coupled to the field. The polarization is detected as long as the main orientation of the magnetic field is not along the line of sight (see Section 2.1), which is unlikely to be common considering the relatively high polarization fractions observed in the protostellar envelopes considered here.
The differences in the power law index α relating S and P frac between the Planck results at cloud scales and the correlation found at core scales with ALMA data (see our correlation and the red line in Figure 4) may be caused by different natures of the disorganized components of the magnetic field at these two spatial scales, where local physical conditions are very different. In the analytical model of Planck Collaboration et al. 2018, the function f m (δ), which depends on the lag, quantifies the disorganized component of the magnetic field relative to its uniform component. Using the dependence of f m (δ) as a function of δ, one cannot adequately perform the extrapolation of the correlation's intercept value we found (which we denote as f in Figure 4) between the Planck and ALMA scales, because the underlying analytical model used to express the dependence of this function on the scale relies on the hypothesis that the disorganized component of the magnetic field is isotropic, which in turn reflects the properties of the turbulent cascade at work in the diffuse ISM. 2 This model would predict values of f m (δ) (and thus levels of turbulence) that are too small at core scales; the vertical shift between the red line and our correlation confirms this point ( Figure 4). For typical lowmass cores, the contribution of the turbulence from the ISM is expected to be negligible. However, cores are observed to be turbulent at some level: typical linewidths are subsonic in the ∼ 1000 au-scale inner envelopes of Class 0 protostellar cores (Gaudel et al. 2020), to trans-sonic at low-mass starforming cores scales (Friesen et al. 2017;Keown et al. 2017). In addition, within these cores, it is expected that the turbulent component of an initially homogeneous magnetic field at core scales would originate from gravo-turbulence induced by collapse motions (Vazquez-Semadeni 2012; Mocz et al. 2017;Ballesteros-Paredes et al. 2018;Vázquez-Semadeni et al. 2019) and outflow phenomena (Zhang et al. 2005;Arce et al. 2007;Plunkett et al. 2013;Frank et al. 2014;Plunkett et al. 2015).
Moreover, note that an adaption of the Planck Collaboration et al. (2018) analytical model (the original version of which included multiple layers of turbulence along the line of sight; see Appendix A) using the specific case of a single-layer model of randomly oriented magnetic field predicts S ∝ P frac −0.5 . It is therefore possible that the flatter correlation between S and P frac observed in cores is due to a smaller number of contributing layers along the line of sight, resulting in an overall less turbulent component of the apparent magnetic field compared with that produced by the multi-scale turbulence at work in the ISM.
The two left panels of Figure 6 show the evolution of S and P frac as a function of the normalized envelope column density N H2 in the envelopes. It seems that S and P frac show opposite trends, which are the result of an increase in the fluctuations in the apparent magnetic field with increasing column density. Note that a similar trend in P frac was found with increasing column density in the diffuse ISM of the Vela C molecular cloud (Fissel et al. 2016). In spite of this intrinsic increase of complexity of the apparent magnetic field with increasing local column density, we still detect substantially organized magnetic fields, as shown by the relatively high values of polarization fraction observed in cores even at high column densities (> 3% at N H2 > 10 24 cm −2 ). In addition, despite of the fact that at the core scale, the main sources of the magnetic field disorder are the dynamical phenomena occurring in the core (e.g., gravitational collapse, outflows, A&A proofs: manuscript no. STAPs_arxiv  S × P frac (right), as a function of the column density NH 2 , where the data from all the cores are merged. The color scale represents number density of points in the plots. The solid black line and black points represent the running mean of S, P frac , and S × P frac ; the associated black error bars are ± the standard deviation of each bin. rotation), we tend to detect strongly polarized emission linked to organized magnetic fields in regions associated with infalling material.
The angular resolution remains an important factor in the statistical analysis of dust polarization observations, because depolarization effects can occur if the resolution of the observations is not high enough to resolve the characteristic length scales of the phenomena driving the small-scale magnetic field morphologies both along the line of sight as well as in the plane of the sky. Beyond the heterogeneity in the characteristics of the ALMA observations that we analyze (such as the angular resolutions and the dynamic range), at the scales we probe here, the magnetic field strength, ionization fraction, gravitational potential, and gas kinematics will affect how an initially uniform magnetic field at envelope scales will develop a complex topology. Given the simple assumption that the gravitational potential is isotropic, considering that the typical spatial resolution we have is on the order of or smaller than the typical Jeans length at the envelope densities we probe, the typical spatial scales at which gravity is expected to significantly distort the magnetic field lines are mostly resolved at the scales (a few beams) where we compute the dispersion S. Nevertheless, if the magnetic field is highly complex at smaller scales than the ones we probe, then indeed P frac drops and conversely S rises toward its highest values.

On the differences in the intrinsic scales of the total intensity (Stokes I) and polarized (Stokes Q and U ) emission
The spatial distributions of both the polarized and unpolarized emission in the plane of the sky show characteristics that are likely to affect the polarization fraction toward protostellar cores, and thus the statistical results we present in this paper. A qualitative view of typical ALMA maps of the dust emission from cores often reveals that the emission in Stokes Q and U looks sharper and more extended that in Stokes I. We have therefore examined the spatial power spectra, which quantify the power present at each spatial scale, of the Stokes I, Q, and U emission in each of our 15 datasets (see Appendix C). Each spectrum is normalized to its maximum value, which allows us to compare the relative power of the emission as a function of spatial scale. We find that, generally, once normalized to its maximum value, the power in the Stokes Q and U maps tends to be larger than the power in Stokes I sometimes by more than one order of magnitude. A larger fraction of the total polarized power resides at larger spatial scales, which explains why the polarized intensity maps appears less peaked than the total intensity maps. This effect could be due to severe dynamic-range limitations and image recovery artifacts affecting the Stokes I maps. However, we stress that the discrepancies in power between Q and U versus I, in the majority of the sources, are more and more significant as we probe larger spatial scales. Since all Stokes are from the same electromagnetic waves received by the same interferometer, if the polarized and unpolarized emission originally had similar spatial distributions, there would be no reason that interferometric filtering would create such differences in the power recovered at different angular scales for the different Stokes parameters. Thus, it is likely that these power spectra reflect the different intrinsic ("true") spatial distribution of polarized and total dust emission at typical core scales probed with ALMA observations. The differences in the spatial distribution of power between Stokes maps towards protostellar cores could be the underlying cause of the high values of polarization fractions observed at lowest observed column densities, which correspond to the largest radii in the envelopes (indeed, see the trend of Figure 6 and the core-by-core analysis of Figure B.4 where we see that the high polarization fraction values correspond to low column density values). While models do not predict such high levels of polarization (see Section 4.1.3 and the simulation and radiative transfer presented in Appendix D), this may contribute to the shallower correlation we find in the S versus P frac relation from the ALMA observations with respect to the trend obtained with large-scale cloud observations: indeed, a decrease in the the highest values of P frac would result in a relation closer to S ∝ P frac −1 , found in Planck and predicted by their analytical model. Confirmation of this result will require further investigation, as it is crucial to understand how much the scales of emission differ between polarized and total intensity emission. Quantifying this would allow us to remove the biases in the values of polarization fraction derived from interferometric observations.

On the effect of spatial filtering on statistical polarization properties
One major issue we face with the ALMA dust polarization observations is spatial filtering by the interferometer, which removes the scales of emission that are not included in the uvcoverage of the dataset. In contrast to the statistical analysis of dust polarization performed with single dish instruments such as Planck (Planck Collaboration et al. 2015a,b, 2018, BLASTPOL in (Fissel et al. 2016), and SCUPOL (Poidevin et al. 2013), our analysis of interferometric data requires us to characterize how the filtering alters the polarization quantities we use in our statistics. With this aim, we use a set of synthetically observed non-ideal MHD simulations computed with RAMSES (Teyssier, R. 2002;Fromang et al. 2006) that follow the gravitational collapse of cores whose range of initial mass and turbulence reproduce the main characteristics of the sources from our sample. The set consists of six simulations of collapsing cores (with total masses of 30, 60, and 100 M ). We perform radiative transfer on these models using the POLARIS code , which produces the Stokes I, Q, and U maps and assumes either that a constant fraction of the dust grains are perfectly aligned everywhere (perfect alignment, known as "PA" hereafter) or that paramagnetic grains are aligned via radiative torques, known as "RATs" hereafter (e.g., Lazarian 2007). Note that the hypothesis of perfect alignment is not physical, and we do not aim to reproduce or interpret the polarized dust emission from Class 0 envelopes as resulting from perfect alignment. However, while we recognize that an hypothesis of perfect dust alignment is not a physical , and with or without filtering). The histogram lines have been smoothed with a 1D-Gaussian kernel of a size of 0.2% in P frac and 2 • in S. In both panels, the shaded areas correspond to the mean of the uncertainty in P frac and S within each bin of the histogram . In the right panel, we do plot the errors in S, derived following Alina et al. (2016). We do not show uncertainties for the synthetically observed simulations, as they have not been filtered by the CASA simulator.
model but a phenomenological one, it has been suggested that the properties of dust polarization at the larger scales of the diffuse ISM (especially the results of S × P frac ) can be explained and reproduced with perfect alignment (Planck Collaboration et al. 2018;Seifried et al. 2020). In the first part of our discussion we aim to compare our results with those obtained at larger scales, and thus perfect alignment remains an interesting point of comparison with RATs, and is a useful benchmark to compare how different physical models of grain alignment affect the statistical properties of the polarized emission. In addition, a case where the grains are perfectly aligned is only taking into account the source-specific geometrical effects governing the resulting polarization maps, and thus is useful to understand where alignment drops or is suppressed. We present all the details of the simulations and the radiative transfer calculations in Appendix D. In order to produce realistic synthetic observations to compare with the ALMA datasets, we use the CASA simulator (with the typical ALMA uv-coverage of these observations) to implement the effects of interferometric filtering and atmospheric noise on the POLARIS synthetic emission maps. In Figure 7 we present the histograms of S and P frac (where all simulations have been merged) before and after filtering, assuming RATs or PA 3 . Note that if no spatial 3 Note that merging all the simulations does not change the result, as each simulation of the six we present in Appendix D sees their values of S and P frac increase. This increase is also seen in Figures D.3 and D.4. The effects of filtering in the case of the three massive simulations (Figure D.4) are marginal, most likely because these cores are very bright and exhibit magnetic fields that are on average more organised than the three less massive simulations.
filtering is applied to the synthetic maps, it is difficult to reproduce the rather high values of polarization fraction typically observed in ALMA observations of protostellar cores using models that only include grain alignment via RATs (despite the fact that we include relatively large grains in our calculations; see details in Appendix D); this was pointed out previously in Valdivia et al. (2019). We find that spatial filtering systematically causes the entire distributions of the dispersion of polarization angles S and the polarization fraction P frac to increase. These increases also translate into an increase in the mean values of S ×P frac (see Figure 8 for the evolution of S × P frac as a function of the column density, in different simulations and implementing grain alignment via both PA and RATs). In addition, one can see that the effect of filtering in Figure 8 seems to be stronger at low column densities. This makes sense because the low column density regions lie at large scales (within the envelope probed by our ALMA observations), where the power spectra in the Stokes I versus Stokes Q and U maps show large discrepancies (see Section 4.1.2). It is therefore possible that some of the high values of polarization fraction P frac found in the ALMA dust polarization observations of protostars may be related to the spatial filtering of the intrinsically different spatial distributions of the polarized versus total intensity emission.
While there is significant uncertainty in the reliability of the calculated values of P frac given our analysis of both filtering and power spectra, we note that the statistical behavior of P frac corresponds to what is predicted by theory and models. In addition, the distributions of P frac we present in Figure 1 peak at reasonable values (∼ 5%); the high values of P frac that we find in the ALMA observations are located in the tail of the P frac distribution. As a thought experiment, we apply the same kind of filtering that ALMA produces to the Planck maps by artificially placing them further away so each map has the angular size of the ALMA field of view at 870 µm. We then synthetically observe them with the CASA simulator using a combination of ALMA antenna configurations similar to those used to observe the sample of cores analyzed here. We find that this simple exercise indeed confirms the change of the power law index of the correlation between S and P frac : from an initial 0.94 before filtering to 0.48 after filtering (see Figure 9). This change is drastic, and the power law index we obtain is much smaller than the ones we obtain from the ALMA observations. This can be explained by the fact that the emission from the Planck observations corresponds to very diffuse regions, and thus the filtering is removes a significant part of the initial flux in the three Stokes maps, which yields a more dramatic effect of the spatial filtering on the recovered correlation between S and P frac .

On the dust grain alignment efficiency inside a Class 0 protostellar core
In this section we discuss how our statistical analysis of polarized dust emission properties has improved our understanding of the dust grain alignment process at work in Class 0 protostellar envelopes. We mainly focus on the interpretation of the evolution of S × P frac as a function of column density, both in ALMA observations (Sections 4. According to the analytical model developed for the ISM by the Planck collaboration (see Appendix A), the product S × P frac is a proxy for the maximum dust grain alignment efficiency P frac,max 4 , and is statistically independent of the magnetic field configuration. This value of P frac,max is influenced by a variety of parameters, such as the collisional de-alignment of grains by gas particles (which scales with density); the dust grain size, shape, and composition; and the local irradiation conditions. We stress that the absolute average values of S × P frac that we present here cannot be compared directly with the values derived from the Planck data because of possibly different physical origins of the turbulent component of the magnetic field at ISM versus core scales. However, if the turbulent component of the magnetic field is still on average isotropic at core scales, then, as S and P frac are inversely dependent on the disorganization level of the apparent magnetic field, it is reasonable to assume that S × P frac traces the intrinsic capability of dust grains to align themselves with the local magnetic field.
We show in Figure 6 that the product S × P frac obtained with all the ALMA dust polarization observations is remarkably constant as a function of column density in protostars, with an average value of 0.36 +0.10 −0.17 . Despite the increasing complexity of the magnetic field topology from core to disk scales, the drastically different local physical conditions (e.g., density, pressure, temperature, and irradiation conditions), the flat profile of the average S × P frac over two orders of magnitude in column density suggests that, within the statistical uncertainties reported in Figure 6, the grain alignment efficiency remains approximately constant throughout a protostellar envelope. This is reminiscent of the Planck results in star-forming clouds: it suggest that both in the ISM and in cores, the dust grain alignment mechanism(s) at work do not appear to be very sensitive to local physical conditions. We stress that in the range of column densities accessible to ALMA, the models implementing RATs (when all averaged together) show a decrease of a factor of two in S × P frac relative to what we see from the ALMA data (see Figure 8 left). In the sections that follow, we explore possible reasons (e.g., different local environmental conditions) for the discrepancy between our ALMA results and the models implementing grain alignment via RATs.

The effects of environmental conditions on dust grain alignment efficiency
Our findings that the average grain alignment efficiency does not strongly depend on the local column density in protostellar cores is, however, at odds with the expected behavior of grain alignment with respect to the quite inhomogeneous local conditions in cores. For example, the rise of gas pressure and density near the center of the protostellar core, which causes the gaseous damping timescale to decrease, is a crucial factor that theoretically leads to a loss of dust grain alignment efficiency (Reissl et al. 2020). Furthermore, observations have revealed that radiation, presumably caused by accretion processes near the central protostar, causes enhanced polarized emission along the cavity walls of bipolar outflows (Hull et al. 2017a;Maury et al. 2018;Le Gouellec et al. 2019;Hull et al. 2020). Finally, indications of larger dust grain size with respect to the ISM dust grain population have been found in embedded objects (Miotello et al. 2014;Valdivia et al. 2019;Galametz et al. 2019;Le Gouellec et al. 2019;Hull et al. 2020). This suggests that, in the context of the RAT alignment mechanism, these phenomena may counter-balance one another, thus precluding a significant variation of alignment efficiency as a function of column density. Then the constant trend of S × P frac could be due to averaging all the observations from our sample, whereas the individual protostars may have very different local conditions at a given normalized column densities because, e.g., their luminosity and absolute densities are different. In addition, note that the statistical weights (in terms of number of independent points, see Table 2) of each observation are very different. Therefore, averaging all the observations may cause observations with larger weights to mask the results of those with lower weights. Hence, to examine in more detail some of these physical processes that are thought to cause the efficiency of RATs to vary, we perform two separations in our datasets. First, we separate the "low" and "high" luminosity cores (i.e., high Observed cores High − luminosity core models implementing RAT High − luminosity core models implementing RAT and filtering Low − luminosity core models implementing RAT Low − luminosity core models implementing RAT and filtering Fig. 8: Observed distributions of the mean values of S × P frac as a function of the column density NH 2 (normalized by its maximum value, N H 2 ,peak ) of all the ALMA cores (triangles) and of all the models (crosses). Left: the four lines representing the simulations correspond to results from all the simulations merged together, using RATs or PA, both filtered and not filtered. Right: we focus on the simulations implementing RATs only, both filtered and not filtered, separating the three simulations with low protostellar accretion luminosity from the three with high accretion luminosity. The shaded areas represent ± the standard deviation of the Gaussian fit performed on each bin of points. The error bars correspond to these standard deviation values divided by the square root of the number of points in each bin. 8(N), B335, VLA 1623, L1448 IRS2) in order to investigate whether the strength of the central source of irradiation affects the dust grain alignment efficiency in the entire core. Note that the value of 10 L is arbitrary, and simply allow us to separate the sample into two bins in luminosity. Second, we separate the emission from the outflow cavities versus that from the envelope on each side of the outflow by taking the emission from inside and outside of the cone of the bipolar outflow (see Figure B.2 and Appendix B). While the magnetic field lines in outflow cavity walls are very organized (having small values of S), which contributes to the observed high values of polarization fraction, the enhancement of polarized emission in these regions seems to be linked with irradiation conditions favourable to dust grain alignment via RATs. Indeed, detections of CCH spectral line emission-a molecular tracer known to be sensitive to irradiation-toward outflow cavity walls in B335 (Imai et al. 2016) and in Serpens Emb 8(N) (Le Gouellec et al. 2019) support this hypothesis. On the contrary, the envelope emission not associated with the outflow is not expected to have such favorable irradiation conditions because of the larger amount of dense envelope material located at all depths along the photon propagation path through the envelope.
In Figure 10 we present, for each of these four cases, the evolution of S × P frac as a function of the normalized column density. The trends of S × P frac in outflow cavity versus envelope emission are very similar. Thus, despite very different irradiation conditions in the outflow cavities and the surrounding envelope, these differences seem insufficient to cause observable changes in the grain alignment efficiency between these two regions. Assuming grains are aligned via RATs, the grains embedded in the envelope thus must still receive amounts of anisotropic irradiation that are sufficient to align grains. However, we do see significant differences in the trends of S × P frac between the low-and high-luminosity cores. The distribution from the high-luminosity cores follows the constant trend from the two previous cases (i.e., outflow cavity and envelope emission), as well as the trend seen when all the datasets are merged together. On the other hand, the distribution from the low-luminosity cores shows a clear decrease in S × P frac as a function of column density, which indicates less efficient grain alignment in the highest density regions of these cores. As the product S × P frac is also dependent on the function f m (δ) (see Section 4.1.1), this may indicate that these low-luminosity cores are subject to different amounts of the gravo-turbulent motions responsible for the correlation between S and P frac . However, while this would modify the average values of S and P frac , it is unclear that this would result in such a strong decrease of S × P frac with the column density, and there are no obvious reasons why only the low-luminosity cores would be affected.
A possible explanation for this decrease of S × P frac toward low-luminosity cores is the amount of irradiation received by dust grains with respect to their position in the protostellar envelope. Indeed, the sub-sample of lowluminosity cores tend to have smaller envelope masses than the high-luminosity cores. Therefore, the optically thin regions of dust emission in the low-luminosity cores tend be closer to the protostellar embryo. However, as the irradiation emanating from the low-luminosity protostellar embryo at the center of the envelope is smaller relatively to the high-luminosity protostars, dust grains located close to the  Figure 4. Interferometric filtering degrades the quality of the S versus P frac correlation (smaller R 2 ) and affects its power law index α, which flattens from 0.94 before filtering to 0.48 after filtering.
peak in column density may be less efficiently aligned in the low-luminosity cores, which causes the drop we see of S × P frac toward these sources. On the contrary, the thermal dust emission emanating from the outer envelope of the lowluminosity cores, i.e., where 10 −2 < N H2 /N H2,peak < 10 −1 , may correspond to column densities that are low enough to still be permeated by the interstellar radiation field. This would increase dust grain alignment efficiency and could potentially justify the on-average higher S × P frac values of the low-versus the high-luminosity cores, in this range of normalized column density. In addition, the higher irradiation emanating from the central protostar of the high-luminosity cores may propagate far enough to maintain a relatively high grain alignment efficiency, even at larger envelope radii.
Finally, note that the overall larger column densities in more massive protostars (which is the case for the highluminosity cores of our sample) also leads to a larger amount of material being optically thick: it is possible that a decrease of S × P frac also happens in these cores, but within the inner envelope radii where dust emission becomes optically thick (and is thus hidden) at (sub)-mm wavelengths. . We examine the maps of integrated optical depth in the synthetic observations of our numerical protostellar models and find that, indeed, the central regions (i.e., the inner ∼ 100 au in the three lowluminosity simulations, and the inner ∼ 200 au of the highluminosity simulations) are optically thick. Nevertheless, the trends of S × P frac as a function of column density for the simulations implementing RATs show relatively flat profiles for both the low-and high-luminosity cases (Figure 8).  the column density NH 2 , which is normalized to its maximum value N H 2 ,peak , for all the cores, and for other cases we tried, including separating high-and low-luminosity cores, as well as outflow cavity walls versus envelope emission. The solid lines are the running means, and the shaded areas are ± the standard deviation of each bin of points. The error bars correspond to these standard deviation values divided by the square root of the number of points in each bin.

The role of protostellar luminosity in aligning grains via RATs in MHD models
Analyzing the statistics from the simulations using the mean of S × P frac estimator as a function of N H2 yields the result Article number, page 15 of 38 A&A proofs: manuscript no. STAPs_arxiv shown in Figure 8, where we compare in the left panel the observational distribution presented in Figure 6 with the distributions from all the simulations, implementing RATs or PA, both filtered and not filtered. Similarly, in the right panel of Figure 8 we compare the same observational distribution with the distributions from the simulations implementing RATs, separating the three simulations with more massive cores and hence higher luminosities from the three others with lower masses and lower luminosities.As we index the irradiation emanating from the central protostar directly to its mass, the three more massive simulations have a stronger radiation field in the core. We find that, in the case of RATs, S × P frac is on average higher in high-luminosity sources than in low-luminosity sources (Figure 8 right panel). This shows that the statistical estimator S × P frac is sensitive to dust grain alignment efficiency, as expected from RAT theory, which predicts that grain alignment efficiency scales with the strength of the local radiation field. Finally, note that in the case of perfect alignment, as all of the susceptible grains are aligned with the magnetic field, we see even higher average values of S × P frac (Figure 8 left panel). Note that the three massive simulations correspond to much higher-mass cores that those in our sample of Class 0 protostellar cores. Moreover, these simulations do not have initial turbulence, and are quite axisymmetric at the time steps we choose, which may be inconsistent with the requirement of the analytical model that the disorganized component of the magnetic field be isotropic in order to justify tracing dust grain alignment efficiency with S × P frac . We stress that, however, at all column densities probed in our sample of protostars, the observed S × P frac values are overall larger than those predicted by models with dust grains aligned via RATs, but are consistent with S × P frac values predicted by models where grains are perfectly aligned (see the left panel of Figure 8). The perfect alignment hypothesis allows us to estimate the typical S × P frac values produced by the combination of perfect local alignment and geometrical effects along the line of sight and in the plane of the sky, and in consequence suggests that the grain alignment efficiency in protostellar envelopes is higher than the efficiency produced by standard RATs alone. Only models implementing RATs in high-luminosity cores (see Figure 8, right panel) produce values that are marginally consistent with the observed values of S × P frac in our sample of protostars. Our results thus suggest that the efficiency of grain alignment via RATs does not match most observations, and highlight the importance of investigating the potentially key role of protostellar irradiation, in our future efforts to reproduce the observed S × P frac . We stress, however, that implementing different dust-grain properties (see Section 4.2.4) as well as different grain alignment mechanisms (see Section 4.2.5) could potentially allow models to approach the values of S × P frac seen in ALMA observations.

Different dust grain properties
Other potential origins for the differences in dust polarization properties in protostellar environments with respect to the ISM, are the different dust properties, such as, e.g., dust grain size, structure, or even composition. The plethora of dust polarization detections toward the densest regions of young protostellar cores indicate that dust grains are still aligned down to very small scales (∼ 100 au) close to the embedded protostar, where dust grain characteristics may not be well constrained.
Estimations of photon-penetration length scales in submillimeter wavelength ALMA observations of protostellar cores have revealed that, given the wavelength of the radiation impinging on the dust grains, the detected dust polarization should emanate from dust grains larger than the (sub-)micron-sized dust grain size expected in a typical ISM population (Le Gouellec et al. 2019;Hull et al. 2020). Radiative transfer studies assuming dust grains aligned by RATs of simulations of low-mass collapsing cores have shown that the typical amount of polarization detected in observations can be reproduced by simulations if the implemented maximum dust grain size exceeds 10 µm ). In addition, indications of such large grains have been found in multi-wavelength observations of protostellar envelopes in studies analyzing dust grain emissivities (Miotello et al. 2014;Galametz et al. 2019). Finally, while the typical elongation of dust grains in star-forming environments is unconstrained observationally, grain alignment may strongly depend on this parameter; further efforts that use dust models to produce predictive observational tests would help further constrain the effect of grain elongation.
Finally, grain alignment efficiency closer to the observed values traced by S × P frac in protostellar envelopes (see Section 4.2.3, where we find that the observed level of dust grain alignment efficiency is not reproduced by standard RATs in the radiative transfer calculations of our models) may be reached if we change the paramagnetic properties of the dust grains used in our radiative transfer calculations. Assuming RATs are the main mechanism aligning grains in protostellar environments, considering dust grains with super-paramagnetic inclusions allows RATs to align more grains (Hoang & Lazarian 2016). Note that this modified RAT theory was tested in models of the diffuse interstellar medium by Reissl et al. (2020), who found that RATs acting on super-paramagnetic grains produce values of grain alignment efficiency very similar to those obtained when grains are perfectly aligned.

Other grain alignment mechanisms
Despite a global agreement in the trends, the statistical properties of polarized dust emission seen in ALMA observations of protostellar cores cannot be fully reproduced by the synthetic observations of MHD simulations of young, collapsing protostellar cores. The distributions of S versus P frac from the simulations (see Figures D.3 and D.4) show clear trends, but they do not match those found from the ALMA observations ( Figure 4). As mentioned in Section 4.2.4, we lack detailed understanding of dust grain properties in the protostellar envelopes probed by ALMA observations. In addition, although we demonstrate the influence of irradiation on the efficiency of grain alignment via RATs in the S × P frac analysis of our models (see Section 4.2.3), we see in Section 4.2.2 that no major differences in dust grain alignment efficiency were found between two regions of the cores that should experience different irradiation conditions: the outflow cavity walls, and the regions in the envelope not associated with the outflow. Finally, the average values of S × P frac from ALMA observations do not seem to match the values obtained from models implementing RATs, but rather show a better match with perfect alignment.
One possibility to explain the aforementioned issues, assuming our protostellar MHD models accurately represent the environments of young protostars, is that we may not fully understand all of the mechanisms causing the linear polarization we detect, and that an additional mechanism(s) may be dominant over RATs when the latter are no longer effective. The dynamical context of some dust polarization observations, especially in the outflow cavity walls or accretion streamers, may favor other theories of grain alignment such as Mechanical Alignment Torques (MATs). Introduced in Hoang et al. (2018), the MAT theory describes the alignment of dust grains with respect to the magnetic field orientation via mechanical torques induced by supersonic gas-dust drift (in an outflow, for example). It is, however, not yet clear how one could identify the occurrence of this new dust grain alignment mechanism in the objects we study here. In addition, the dust grain size distribution may be affected by RAdiative Torque Disruption (RATD; Hoang et al. 2019;Hoang & Tram 2020), which predicts that the dust-grain size distribution will shift to smaller values due to the disruption of large aggregates that are spun-up to suprathermal rotation speeds and thus broken apart by radiative torques.
The statistical analysis of the dust polarization we perform in this article, thanks to the tools of S and P frac , is able to characterize the processes at work in the alignment of dust grains in the envelopes of Class 0 protostellar cores, such as the role played by the radiation field. However, additional methods must be developed in order to identify the potential occurrence of recently proposed grain alignment mechanisms.
1. We find a significant correlation between the dispersion of polarization angles S and polarization fraction P frac in the polarized dust emission from protostellar envelopes, with a resulting correlation of S ∝ P frac −0.79 . This correlation is sensitive to the morphology of the turbulent component of the magnetic field and to other intrinsic characteristics of the polarized emission. This correlation found in the ALMA cores has a smaller power law index than the correlation found at larger scales in the Planck observations of star-forming clouds, where they found S ∝ P frac −1 . This could be a consequence of the different nature of turbulence in Class 0 sources versus the ISM (i.e., due to gravitational infall, rotation, and outflowing motions); or due to interferometric filtering, which produces artificially high P frac in ALMA observations; as well as of the possibility that grain alignment varies with local conditions, which are significantly different between the star-forming molecular clouds and protostellar cores. Finally, our observations and their comparison to synthetic observations of protostellar models suggest that additional alignment mechanisms may be at work in protostars (see point 6). 2. We find that the flattening of the correlation between S and P frac in our ALMA results versus the larger-scale Planck results can be reconciled with the Planck analytical model if it is modified to include only one layer of randomly oriented magnetic field to represent the turbulent component of the field. This results in an overall less turbulent component of the apparent magnetic field compared with that produced by the multi-scale turbulence at work in the ISM. 3. The product S × P frac , which is sensitive to dust grain alignment efficiency, shows a constant profile as a function of column density in the sample of cores analyzed, with a constant value of 0.36 +0.10 −0.17 . This suggests that the grain alignment mechanism producing the polarisation observed at millimeter wavelengths, over 3 orders of magnitude in column density (from N H2 = 10 22 cm −2 to N H2 = 10 25 cm −2 ), may not depend strongly on the local conditions such as gas density and temperature. 4. We examine the statistical properties of polarized dust emission emanating from the outflow cavity walls versus the regions of the envelope not associated with the outflow. These regions are expected to experience drastically different irradiation conditions. We do not find any obvious difference in dust grain alignment efficiency between the two. 5. However, we find hints that, contrary to the highest luminosity cores in our sample, the lowest luminosity sources experience a decrease of their dust grain alignment efficiency at higher column densities. The environmental conditions in the central regions of the envelopes are indeed expected to disfavor the alignment of dust grains via the Radiative Alignment Torque (RAT) mechanism, as a result of the lower level of irradiation emanating from the central protostar of the low-luminosity cores relative to the high-luminosity cores. The density of the outer envelope of these low-luminosity cores may be tenuous enough to be permeated by the interstellar radiation field, thus increasing dust grain alignment efficiency with increasing radii. Finally, the higher irradiation emanating from the central protostar of the high-luminosity cores may propagate far enough to maintain a relatively high grain alignment efficiency, even at larger envelope radii. 6. We use synthetic observations of the polarized dust emission in a small sample of outputs from non-ideal MHD simulations of protostellar collapse. We apply the S × P frac analysis to these synthetic maps of polarized dust emission, assuming either grain alignment via Radiative Alignment Torques (RATs) or perfect alignment (PA; i.e., alignment of all susceptible grains), and we show that the statistical estimators used in our work seem to be sensitive to the overall efficiency of grain alignment. Furthermore, our S × P frac analysis of the simulations implementing RATs suggests that the average value of this estimator is sensitive to the radiation field strength in the core. Finally, the simulations with perfect alignment yield on average higher S ×P frac values than those implementing RATs. 7. When implementing RAT alignment in our radiative transfer calculations, we do not reproduce with our simulations the S versus P frac statistics obtained from the ALMA observations. This may suggest that the simulations are not fully adequate representations of the Class 0 protostellar envelopes in our observations, or that the S versus P frac correlation is not sensitive to the details of the physical mechanism(s) aligning the dust grains. The values of S × P frac obtained from the ALMA observations seems to lie among the values predicted by PA, and are significantly higher than those found in models including RATs alone, especially at high column density. This suggests that, to be able to reproduce the dust alignment efficiency found in cores, one needs either more efficient RATs than the classical RATs with paramagnetic grains, an extra alignment mechanism(s), or different irradiation conditions than those assumed in models. 8. Our results suggest that the continuum and polarized dust emission in the ALMA observations have different intrinsic spatial scales, which affects the statistics. We show that the differences in emitting power of the different Stokes parameters as a function of spatial scale can produce artificially high P frac , especially at large scales where Stokes I has on average less power with respect to Stokes Q and U . Finally, this work on synthetic observations suggests that interferometric filtering biases the values of S and P frac , causing artificially high values of both.
While the work we present here has shed light on the physics of dust grain alignment in Class 0 protostellar cores, many open questions remain about the details of the physical environment at envelope scales. Future investigations involving detailed comparisons of the observations of cores with those reproduced by simulations, alongside observations of chemical tracers associated with the polarized dust emission, will illuminate the role played by the local conditions in producing the polarization observed at small scales in Class 0 protostellar envelopes.  2018) developed an analytical model able to reproduce the phenomenological properties of polarized dust emission. They assumed the total emission arises from a small number N of independent layers, each of them emitting a fraction 1/N of the total intensity. The magnetic field was described as the sum of a uniform and an isotropic turbulent component. This model is based on a few essential parameters, including the maximum polarization fraction P frac,max (which will tell us about the intrinsic capability of the grains to align themselves with respect to the magnetic field), the ratio of the standard deviation of the turbulent magnetic field to the magnitude of the ordered field f m , and the spectral index α M of this turbulent component. At a given location, the analytical relationship between the dispersion of polarization angles S and the polarization fraction P frac was found to be the following: where φ is the polarization position angle, ∆φ i = φ − φ i , δ is the lag (as introduced in Section 3.1, the lag describes the surface over which the dispersion S is derived, and thus corresponds to the characteristic length scale at which we quantify the disorganization of polarization position angles), and Γ i is the inclination angle of the magnetic field − → B i in a given layer i with respect to the plane of the sky. The value of A is approximated as follows: such that we obtain: where factor f m (δ) represents the typical relative fluctuation of the magnetic field at the scales corresponding to the annulus between δ/2 and 3δ/2. This factor was defined as follows: where σ Bi (δ) is the fluctuation of the magnetic field − → B i . This function was modeled as follows: where ω is the full width half maximum of the spatial resolution of the observations. The parameter values used in Planck Collaboration et al. In the plots where we relate S and P frac , we plot this relation in red, using the analytical coefficient of 0.339. The Planck team found a coefficient of 0.31 from their observations. When considering the results of the Planck team, this analytical model yields a dispersion of polarization angles S that is proportional to P frac −1 . However, for our study, the specific case of N = 1 is relevant. This gives: In this specific case: P frac = P frac,max cos 2 Γ , (A.9) and thus we obtain: Consequently, in the specific case of N = 1, S is proportional to P frac −1/2 in the analytical model.

Appendix B: ALMA cores
In Table B.1 we list the details and the coordinates of each source in our sample. In Table B.2 we present the outflow properties of each source. In Figure B.2 we present the maps of Stokes I, polarized intensity P , and polarization angle dispersion S, for all of the ALMA observations. The white dotted lines show the separation between the outflow cavities and the envelope emission not associated with the outflow cavities. We characterize this separation in Figure  B.1, where e and c denote the thickness of the equatorial planes (separating the two lobes of the outflow) and the outflow cavity walls, respectively, as fit by eye using the polarized intensity maps. The outflow position angles and opening angles are taken from the literature. Figure B.3 shows the S vs. P frac correlations for all the datasets. (see Figure B.1).
Article number, page 21 of 38 A&A proofs: manuscript no. STAPs_arxiv   (Enoch et al. 2011;Kristensen et al. 2012). The values calculated for Serpens Emb 8 and Emb 8(N) include both sources together (Enoch et al. 2009(Enoch et al. , 2011    The parameters e and c have been determined by eye thanks to the previously published CO emission. Furthermore, IRAS 16293B lies within one of the outflow cones of IRAS 16293A; however, we consider it to be envelope emission. f VLA1623A/B do not have clearly identifiable outflow cavities in the polarized dust emission. Therefore we consider all the emission of this source to be coming from the envelope. However, we keep the parameters e and c (determined by eye thanks to the CO emission) in this Table and the white-dotted lines in Figure B.2 to simply indicate the location and shape of the outflow. g We report e and c for the two datasets separately, as they have very different angular resolutions. h We report values separately for the two protostars in this core, IRAS4A1 and IRAS4A2.

Appendix C: Power spectra as a function of spatial scale in the ALMA observations
The power spectra of the maps of the Class 0 sources observed with ALMA are shown in C.1. To produce these power spectra we perform 2D Fourier transforms of the Stokes I, Q, and U maps and take their normalized absolute magnitude, after which we calculate the azimuthal average to derive the power with respect to the spatial scale. We make an attempt to recover the missing flux at large spatial scales of Stokes I with respect to Stokes Q and U . To correct the Stokes I power spectrum across a given range of spatial scales, we scale up the flux in Stokes I using the differences between the Stokes I and the Stokes Q and U power spectra at those same spatial scales. The idea behind this correction is to solve the problem of high P frac values discussed in Sections 4.1.2 and 4.1.3. Unfortunately, this simplistic method is not robust; as the Stokes I signal at large scales is buried in noise, we do not properly recover the initial missing flux. Furthermore, this simple method creates artifacts in the Stokes I maps.
Note that we produced power spectra of the Stokes maps from the simulations in the same way that we do for the ALMA observations. We find the same discrepancies between the power in Stokes I and the power in Stokes Q and U at large scales. However, as our simulations do not reproduce rigorously the variety of morphologies we detect in the Class 0 ALMA observations, we do not show these additional plots.

Appendix D: MHD simulations and their synthetic observations
In order to characterize better the statistics we obtain from our analysis of ALMA dust polarization observations of Class 0 protostars, we perform synthetic observations of nonideal radiation-magneto-hydrodynamic (MHD) simulations of protostellar collapse, exploring the impact that a range of parameters-such as the dust grain alignment hypothesis, the initial mass and turbulence of the simulation, and the effect of interferometric filtering-have on the statistics from these simulations. We use six different setups for the simulations performed with the RAMSES code (Teyssier, R. 2002;Fromang et al. 2006;Commerçon et al. 2011;Masson et al. 2012), where sink particles are implemented (Krumholz et al. 2004;Bleuler & Teyssier 2014). Mignon-Risse et al. (2020) and Verliat et al., in preparation present in detail similar simulations; however, their simulations employ a novel radiative transfer method that the simulations we use here do not. Three of the simulations follow the collapse of magnetized, intermediatemass dense cores without initial turbulence, while the three others follow the collapse of weakly magnetized, low-mass cores with initial turbulence. The idea behind our analysis of these simulations is to choose different physical conditions to represent the variety of environments present in the observed ALMA cores. To this aim, we use as our models six simulation outputs with central stellar masses between 0.5 and 7 M that sample randomly a domain in initial mass, magnetic energy, turbulent energy, as is also the case with our observations. The details of models we use can be found in Table D.1. In this paper, we do not aim to reproduce or interpret the polarized dust emission from Class 0 envelopes as resulting from perfect alignment. With the goal to assess the statistical properties of polarization, one model already provides enough data points to make a statistical analysis: the inclusion of several models only allow to illustrate that there may be some local conditions (turbulence) affecting slightly the trends, and that irradiation due to the central protostar is key.
We perform radiative transfer calculations on these simulations using the POLARIS code , which calculate the local dust temperature and dust grain alignment efficiency of oblong dust grains with respect the magnetic field orientation following the RAT theory developed in Lazarian & Hoang (2007); Hoang & Lazarian (2014). In each run of POLARIS, we either choose to calculate the grain alignment of each dust grain via RATs, or we employ the Perfect alignment (PA) hypothesis, which assumes that all susceptible grains are aligned with their long axes perpendicular to the local magnetic field orientation. We derive the temperature of the central object from the luminosity of the blackbody, which is indexed to the mass of the sink following the empirical correlation from Weiss et al. (2004). We also include the interstellar radiation field using a value of G 0 = 1 (Mathis et al. 1983). Note that in order to compute the radiative transfer in a reasonable amount of time, we must delete the mass in the highest density cells surrounding the sink (in a ∼ 15 au diameter region). This may result in an overestimated radiation field in the core, as the photons will not be processed by the material we remove. We assume a gas-to-dust ratio of 100. The dust grain population is composed of 62.5% astronomical silicates and 37.5% graphite grains (Mathis et al. 1977); note that this composi-tion governs the ultimate number of aligned grains in the PA regime, as silicates can be aligned with the magnetic field much more easily than graphite/carbonaceous grains (Andersson et al. 2015, and references therein). The dust grains are oblate with an aspect ratio of 0.5 (Hildebrand & Dragovan 1995) and they follow a standard MRN-like distribution (Mathis et al. 1977) with cutoff sizes of a min = 2 nm and a max = 10 µm. We choose this latter value as the maximum grain size in POLARIS in light of recent work that has hinted at the presence of grains larger than the typical ∼ 0.1 µm ISM dust grains in Class 0 envelopes (e.g., Valdivia et al. 2019;Le Gouellec et al. 2019;Galametz et al. 2019;Hull et al. 2020). The radiation field resulting from the radiative transfer, impinging on the dust grains in the protostellar envelope, comprises low-energy submillimeter photons whose wavelength need to be comparable to the size of dust grains in order to efficiently align the grains via RATs.
We synthetically observe the MHD simulations with PO-LARIS at 870 µm, at a distance of 400 pc, in maps 8000 au in size with pixel sizes of 8 au. We observed each of the six simulations along two independent, orthogonal lines of sight. As a result, we analyse twelve different POLARIS synthetic observations, each of which was produced assuming grain alignment via either RATs or PA. Thus, this sample of twelve models and our fifteen ALMA datasets yield a similar numbers of cases to which we can apply our statistical analysis. From POLARIS we obtain the three Stokes parameters I, Q, and U , which we convolve with a 2D Gaussian kernel to smooth out the different resolutions of the cells based on local density, which is due to the use of Adaptive Mesh Refinement (AMR). We choose a pixel size of 8 au with PO-LARIS; however, within the 8000 au core, there are many different cell sizes, which degrades the spatial resolution in some regions of our radiative transfer maps. In order to have independent points while running our statistics, we smooth the resulting Stokes maps to 80 au resolution, which is the largest cell size in the central region of the synthetic observation. Beyond this central region, the AMR cell sizes are even larger than 80 au, but we compute the statistics within the central ∼ 1500 au zone, where the AMR cell size is smaller than 80 au. At this point, we obtain in this central zone a first set of "perfect" maps on which we calculate the same statistical estimators used to study the polarization properties of our ALMA observations; we denote these perfect synthetic observations, "without filtering." In addition, we use the CASA simulator to interferometrically filter the synthetically observed maps, mimicking ALMA observations. For each simulation, we combine synthetic observations from ALMA configurations C-3, C-5, and C-6, with an exposure time of 6000 s per antenna configuration. The resulting synthesized beams (resolution elements) of these filtered maps have an effective size of 80 au. After filtering the maps with the CASA simulator, we compute our statistics in the same way that we do with the ALMA data, using the threshold criterion of Stokes I explained in Section 3.1. We denote this latter set of results, "with filtering." Similar to Figure  B  All the simulations have a spatial resolution of a few au, and implement ambipolar diffusion. The initial density profile is ρ ∝ 1/(1+r 2 ) in runs I, II, III, and is uniform in runs IV, V, and VI. For each of them we select a time step at which the simulation exhibits compact features in density similar to those that we see in the ALMA observations, i.e., bright emission from the infalling envelope, disc-like structures, and bipolar outflow cavities. a The sink mass is that of the largest central sink. This core is fragmenting, and thus there are several other, smaller sinks. b The jet is implemented by hand, with a speed of 66% of the escape speed, and an opening angle of 30 • . The corresponding outflowing mass ejected by the sink is 1/3 of the mass accreted by the sink. We present in Figures D.3 and D.4 the distribution of the polarization fraction P frac as a function of the dispersion of polarization angles S in the synthetically observed maps, separating the three simulations that implement initial turbulence and have lower total mass from the three others that do not implement turbulence and have a much higher total mass. In each Figure, we plot P frac versus S using PA and RATs, before and after spatial filtering.
In Figure D.3, we see that both P frac and S are higher in the case of perfect alignment. In the case where we use the perfect alignment hypothesis, the detected polarized emission covers larger regions of the core than with RATs. Indeed, fewer grains are aligned if we assume RATs; this explains the distribution of polarization fraction, which is directly sensitive to the dust grain alignment efficiency, and is lower in the case of RATs. It also explains the lack of detection when assuming RATs and filtering the maps, as we add atmospheric noise to non-filtered maps that are only marginally polarized. The correlation between S and P frac seems to be the closest to the observational correlation presented in Figure 4 when we consider the lower-mass cases where the simulations have initial turbulence, have perfect grain alignment, and are spatially filtered.
The results from the statistics using the other set of three simulations (see Figure D.4), which have higher mass and no turbulence, behave differently. The observed S versus P frac correlations do not vary significantly whether we filter the maps or not, or whether we use RATs or perfect alignment. This can be explained by the fact that the central heating source is much hotter than in the lower-mass simulations: the simulations used in Figure D.4 have larger initial masses, and are synthetically observed at later times in terms of core evolution, which means that their sinks are more massive, being on the order of a few solar masses. In consequence, RATs appear to be so efficient that the statistics obtained from these simulations are very close to those obtained when we assume perfect alignment. The correlations fitted to the distributions do not vary significantly within these four sub-cases; however, we still notice that on average, the distributions from the perfect alignment cases tend to have larger values of S and P frac , which is the expected behavior.

Appendix E: Π Investigations
In the analytical model of Planck Collaboration et al. (2018), they demonstrate the S × P frac estimator can trace the P frac,max parameter, given some assumptions such as that the intensity maps should not vary strongly and there should be only small differences of polarization position angles ∆φ between adjacent cells, implying that tan ∆φ ≈ ∆φ and Q j Q − U j U P 2 . Finally, the assumption that S and S 2 behave the same is also made in the analytical model. These assumptions may not be valid in the emission maps of Class 0 protostellar cores, as we observe, for example, strong gradients in Stokes I maps. We present here a new estimator of P frac,max , called Π, the derivation of which does not require these assumptions. This new estimator will be investigated in detail by Guillet et al., in preparation. In order to derive the relation between Π, P frac,max , and f m (δ), we follow the same method presented in Appendix A and Appendix E of Planck Collaboration et al. (2018): Therefore, we have: In Figures E.1, E.2, and E.3 we show the comparisons between the results provided by Π and S × P frac . These plots show that the results are only marginally different, and thus we do not recompute our results using this new estimator.  (left) and Π (right) as a function of the column density NH 2 , normalized in each core by its maximum value N H 2 ,peak . The color scale represents number density of points in the plots. The solid black (red) lines and black (red) points represent the running mean of S × P frac (Π); the associated error bars are ± the standard deviation of each bin. To facilitate the visual comparison, the running means of both S × P frac and Π a re plotted in both panels.