© V. J. M. Le Gouellec et al. 2020

1 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 Int. XXXV. 2016; Soler 2019; Seifried et al. 2020), and it appears to 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 that is 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 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, that is 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 observationsof 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) have found a significant decrease in 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 alternatively as a possible lossof 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 as follows: 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 the 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 (Brauer et al. 2016 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 in the polarization fraction in observations of cores at a coarse angular resolution were partly explained by beam smearing, that is to say 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 Int. XIX 2015; Planck Collaboration Int. XX 2015). 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 XII (2020) analyzed these statistical estimators as a function of column density from the diffuse interstellar medium (ISM) to molecular clouds, probing column density up to 10²² cm⁻². They found no variation of dust grain alignment efficiency with varying conditions typical 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₂, 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, 2014), 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 and 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. 2017a,b, 2020; Maury et al. 2018; Sadavoy et al. 2018a,b; Kwon et al. 2019; Takahashi et al. 2019; Le Gouellec et al. 2019). 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 Sect. 2, we introduce the statistical tools we use to analyze 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 Sect. 3. Finally, we discuss the results we obtain regarding the dust grain alignment in young protostellar objects in Sect. 4, along with comparisons between the ALMA observations and synthetic observations of MHD simulations. We draw our conclusions in Sect. 5.

2 ALMA observations of Class 0 protostars

2.1 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 denote the polarized intensity P (defined as $P=\sqrt{Q^2+U^2}$, which we systematically debias, see Sect. 3.1), the polarization fraction ${\mathcal{P}}_{\text{frac}}$ (defined as ${\mathcal{P}}_{\text{frac}}=P/I$), and the polarization position angle ϕ (defined as $\varphi =0.5\arctan \frac{U}{Q}$).

In the diffuse ISM and molecular clouds, Planck Collaboration Int. XIX (2015), Planck Collaboration Int. XX (2015), and Planck Collaboration XII (2020) found a correlation between the local dispersion of the polarization position angle and the polarization fraction. 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 $\mathcal{S}$, which quantifies the local (non)-uniformity of the polarization angle, is defined as follows: $$S(r,\delta )=\sqrt{\frac{1}{N}{\sum }_{i=1}^N{[\varphi (r+{\delta }_i)-\varphi (r)]}^2},$$(1)

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 XII 2020). Planck Collaboration XII (2020) developed an analytical model (briefly described in Appendix A) that relates the two quantities $\mathcal{S}$ and ${\mathcal{P}}_{\text{frac}}$. They found, among other results, that $S\propto {\mathcal{P}}_{\text{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 $\mathcal{S}\times {\mathcal{P}}_{\text{frac}}$ – which is a proxy for grain alignment efficiency – as a function of the column density and the dust temperature, Planck Collaboration XII (2020) 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 $\mathcal{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 $\mathcal{S}$, the more disorganized the apparent magnetic field. $\mathcal{S}$ will saturate at the value of $\pi /\sqrt{N}$, as (ϕ = 0°) ≡ (ϕ = 180°) in a polarization map. In this study, 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 $\mathcal{S}$ will be larger. In contrast, a uniform magnetic field will have dispersion $\mathcal{S}$ close to zero, regardless of the line of sight. These facts limit the capability of $\mathcal{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 ${\mathcal{P}}_{\text{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 $\mathcal{S}$ reaches high values, P may drop below the detection limit.

The polarization fraction ${\mathcal{P}}_{\text{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 ${\mathcal{P}}_{\text{frac}}$ as seen by the observer. Consequently, assuming an isotropic turbulent component of the magnetic field, and given the caveats listed above, $\mathcal{S}$ and ${\mathcal{P}}_{\text{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 $\mathcal{S}$ and ${\mathcal{P}}_{\text{frac}}$ in Class 0 protostars. In this objective, we gathered ALMA observations of polarized dust emission taken toward these type of objects in the recent years.

2.2 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 observationstoward nearby, low- and 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 NGC 1333 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 Sect. 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 values of polarization fraction ${\mathcal{P}}_{\text{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 Fig. 1 the distribution of ${\mathcal{P}}_{\text{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 ${\mathcal{P}}_{\text{frac}}$, and ϕ (in radians) are showed as shaded areas, and are calculated as follows: $${\sigma }_{{\mathcal{P}}_{\text{frac}}}={\mathcal{P}}_{\text{frac}}\sqrt{{\left(\frac{{\sigma }_P}{P}\right)}^2+{\left(\frac{{\sigma }_I}{I}\right)}^2},$$(2) $${\sigma }_{\varphi }=\frac{1}{2}\frac{{\sigma }_P}{P}.$$(3)

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: $$N_{{\text{H}}_2}=\frac{S_{\nu }{d}^2}{A{\mu }_{{\text{H}}_2}{m}_{\text{H}}{\kappa }_{\nu }{B}_{\nu }(T_{\text{d}}(r))},$$(4)

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 observationsfor 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, ${\mu }_{{\text{H}}_2}$ is the mean molecular weight per hydrogen molecule (${\mu }_{{\text{H}}_2}=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): $$T_{\text{d}}(r)=38\text{K}{\left(\frac{r}{100\text{au}}\right)}^{-0.4}{\left(\frac{L_{\text{int}}}{1L_{\odot }}\right)}^{0.2},$$(5)

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 internal luminosities are known from the Herschel Gould Belt survey (André et al. 2010) and were used in Maury et al. (2019), Maret et al. (2020) and Belloche et al. (2020), for others we derived the internal luminosities using archival fluxes at 70 μm from Herschel PACS using the relation from Dunham et al. (2008) (see Table B.1). We find that the ALMA observations are sensitive to material in the inner envelope with typical column densities ~ 10²²−10²⁵ cm⁻². 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 self-scattering 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,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 inboth our observations at 1.3 and 0.87 mm. We also exclude these pixels from our analysis. We note that although Kwon et al. (2019) and Takahashi et al. (2019) disfavored 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.

Fig. 1 Histograms of polarization fraction ${\mathcal{P}}_{\text{frac}}$ (left) and polarization position angle ϕ (right) for all the sources of our sample. The histogram lines have been smoothed with a 1D-Gaussian kernel of a size of 0.2% in ${\mathcal{P}}_{\text{frac}}$ and 2° in ϕ. The shaded areas correspond to the mean of the uncertainty in the values of ${\mathcal{P}}_{\text{frac}}$ and ϕ within eachbin of the histograms. Among the two datasets we have toward OMC3 MMS6, “OMC3 MMS6 H” and “OMC3 MMS6 L” denote the high and low resolution observations, respectively.

3 Analysis of polarized dust emission in Class 0 protostellar cores

3.1 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 $\mathcal{S}$ in each pixel i of the Stokes maps, with respect to each of its n = 8 nearest neighboring pixels j as follows: $$S{(\delta )}_i=\sqrt{\frac{1}{n}{\sum }_{j=1}^n{\left[\frac{1}{2}\text{arctan}\frac{Q(j)U(i)-U(j)Q(i)}{Q(i)Q(j)+U(i)U(j)}\right]}^2}.$$(6)

Considering the sampling described above, the equivalent δ parameter (see also Eq. (1)) is approximately 1∕2 of a beam width (comparable to the value chosen in Planck Collaboration XII 2020). We note that the measured value of $\mathcal{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 $\mathcal{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 $\mathcal{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 Sect. 2.1, the way the dispersion $\mathcal{S}$ is derived causes the distribution to saturate, as a completely random distribution of polarization angles will produce a maximum value of $\mathcal{S}$ of $\pi /\sqrt{n}~{63}^{°}$ (as we chose n = 8).

While Planck Collaboration XII (2020) 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). We note that 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 XII (2020) to compute a pixel-selection criterion in order to test appropriately the correlation between $\mathcal{S}$ and ${\mathcal{P}}_{\text{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 Fig. 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 ${\langle P/{\sigma }_P\rangle }_I\ge 5$. We then take all points lying above this cutoff in Stokes I and form the sample to which we apply our statistical method. However, 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 Fig. 2, the values of $\mathcal{S}\times {\mathcal{P}}_{\text{frac}}$ diverge at Stokes I values lower than the aforementioned threshold, as a consequence of the noise-bias of $\mathcal{S}$ and ${\mathcal{P}}_{\text{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 Fig. 3 we show the maps of Stokes I, polarized intensity P, polarization fraction ${\mathcal{P}}_{\text{frac}}$, and dispersion of polarization angles $\mathcal{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 Fig. B.2.

3.2 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 Fig. B.3 we show he distribution of the dispersion of the polarization angles $\mathcal{S}$ as a function of the polarization fraction ${\mathcal{P}}_{\text{frac}}$ in the 15 maps probing the dust emission in protostellar cores. In these plots, the running mean of ${\mathcal{P}}_{\text{frac}}$ (shown as black points and line) shows the average trend and evolution of $\mathcal{S}$ with ${\mathcal{P}}_{\text{frac}}$. In particular, one can clearly see the area of the distribution affected by the saturation effect of $\mathcal{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 colored based on their Stokes I value. The relationship between $\mathcal{S}$ and ${\mathcal{P}}_{\text{frac}}$ observed byPlanck 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 ${\mathcal{P}}_{\text{frac}}$ associated with low dispersion in polarization angles $\mathcal{S}$ in regions of faint Stokes I values. Conversely, we see high $\mathcal{S}$ and low ${\mathcal{P}}_{\text{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 $\mathcal{S}\propto {\mathcal{P}}_{\text{frac}}{}^{-\alpha }$ 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 $\mathcal{S}$ as a function of ${\mathcal{P}}_{\text{frac}}$, along with the linear fit previously described, which is defined by the two parameters α and f such that $S=f/{\mathcal{P}}_{\text{frac}}{}^{\alpha }$. We note that at low values of ${\mathcal{P}}_{\text{frac}}$, the distribution of $\mathcal{S}$ flattens because of the saturation of $\mathcal{S}$ for high dispersion values. This is an artifact arising from the definition of $\mathcal{S}$, and thus these points should be excluded from the linear fit. To do so, we establish a threshold in ${\mathcal{P}}_{\text{frac}}$ of 1.3%, indicated by the vertical dot-dashed gray line in Fig. 4, which denotes the ${\mathcal{P}}_{\text{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\propto {\mathcal{P}}_{\text{frac}}{}^{-1}$.

Merging all the ALMA data does not enable us to investigate the relation between $\mathcal{S}$ and ${\mathcal{P}}_{\text{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 in Sect. 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 $\mathcal{S}$, ${\mathcal{P}}_{\text{frac}}$, and $\mathcal{S}\times {\mathcal{P}}_{\text{frac}}$ as a function of the local column density $N_{{\text{H}}_2}$ 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)¹. We notice a global trend in the merged data shown in Fig. 6 of increasing $\mathcal{S}$ and decreasing ${\mathcal{P}}_{\text{frac}}$ with increasing column density. Figure B.4 presents the results of linear regressions on the trends of $\mathcal{S}$ and ${\mathcal{P}}_{\text{frac}}$ as a function of $N_{{\text{H}}_2}$ in individual cores: the resulting polar-law indices and R-squared values suggest a significant decrease of ${\mathcal{P}}_{\text{frac}}$ with increasing $N_{{\text{H}}_2}$ across the sample, and hints of increasing $\mathcal{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).

Fig. 2 Illustration of how we compute the pixel-selection cutoff in Stokes I in the 1.3 mm observation of the B335 core. Top: evolution of the S/N of the polarized intensity, i.e., P∕σP, as a functionof 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 Sect. 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: $\mathcal{S}\times {\mathcal{P}}_{\text{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. We note that, as expected, one see that decreasing the resolution, and thus increasing the spatial length of the lag, causes on the dispersion $\mathcal{S}$ to increase as well, on average (see Sect. 2.1).

4 Discussion

4.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 $\mathcal{S}$ observed in the plane of the sky and the polarization fraction ${\mathcal{P}}_{\text{frac}}$ measured in each line of sight in these 11 protostellar envelopes. However, the $\mathcal{S}\propto {\mathcal{P}}_{\text{frac}}{}^{-0.79}$ relationship we find at core scales is shallower than the $\mathcal{S}\propto {\mathcal{P}}_{\text{frac}}{}^{-1}$ relationship found at larger scales in the Planck observations of star-forming clouds (Planck Collaboration XII 2020). Moreover, we obtain on average higher values of $\mathcal{S}$ and ${\mathcal{P}}_{\text{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 (Sect. 4.1.1); we then address the different intrinsic spatial scales of total intensity versus polarized emission (Sect. 4.1.2), and how interferometric filtering may affect observed polarization properties (Sect. 4.1.3).

Fig. 3 Maps from the 1.3 mm observations of the B335 core. Top left: total intensity (Stokes I) thermal dustemission in color scale. The emission is shown from 3σI where σI is the rms noise in the Stokes I map. Top right: polarized intensity P in colorscale.Bottom left: polarization fraction ${\mathcal{P}}_{\text{frac}}$ in colorscale,shown where I > 3σI. Bottom right: dispersion $\mathcal{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 Sect. 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°.

4.1.1 The physics governing the disorganized component of the magnetic field

The correlation between $\mathcal{S}$ and ${\mathcal{P}}_{\text{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 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 Sect. 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 $\mathcal{S}$ and ${\mathcal{P}}_{\text{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 Fig. 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 XII (2020), 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 Fig. 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². 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 (Fig. 4). For typical low-mass 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 star-forming cores scales (Friesen et al. 2017; Keown et al. 2017). In addition, within these cores, it is expected that the disorganized 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, we note that an adaption of the Planck Collaboration XII (2020) 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 $\mathcal{S}\propto {\mathcal{P}}_{\text{frac}}{}^{-0.5}$. It is therefore possible that the flatter correlation between $\mathcal{S}$ and ${\mathcal{P}}_{\text{frac}}$ observed in cores is due to a smaller number of contributing layers along the line of sight, resulting in an overall less disorganized component of the apparent magnetic field compared with that produced by the multiscale turbulence at work in the ISM.

The two left panels of Fig. 6 show the evolution of $\mathcal{S}$ and ${\mathcal{P}}_{\text{frac}}$ as a function of the normalized envelope column density $N_{{\text{H}}_2}$ in the envelopes. It seems that $\mathcal{S}$ and ${\mathcal{P}}_{\text{frac}}$ show opposite trends, which are the result of an increase in the fluctuations in the apparent magnetic field with increasing column density. We note that a similar trend in ${\mathcal{P}}_{\text{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_{{\text{H}}_2}>{10}^{24}$ cm⁻²). 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, 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 theangular 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 $\mathcal{S}$. Nevertheless, if the magnetic field is highly complex at smaller scales than the ones we probe, then indeed ${\mathcal{P}}_{\text{frac}}$ drops and conversely $\mathcal{S}$ rises toward its highest values.

Fig. 4 Distributions of the dispersion of polarization position angles $\mathcal{S}$ as a functionof the polarization fraction ${\mathcal{P}}_{\text{frac}}$ from all of the datasets merged together. The points were selected according the method developed in Sect. 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 ${\mathcal{P}}_{\text{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$\mathcal{S}$ in the derivation of the linear fit by applying a threshold in polarization fraction, indicated by the vertical dot-dashed gray line. This threshold denotes the ${\mathcal{P}}_{\text{frac}}$ level beyond which the distribution is linear. The solid red line corresponds to the Planck correlation from Planck Collaboration XII (2020), 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 ourALMA 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_{\text{frac}}^{\alpha }$. 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 subsamples of points. We calculate the uncertainties in f and α as standard deviations of Gaussian fits to those histograms.

Fig. 5 Distributions of the dispersion of polarization position angles $\mathcal{S}$ (top left), the polarization fraction ${\mathcal{P}}_{\text{frac}}$ (bottom left), and $\mathcal{S}\times {\mathcal{P}}_{\text{frac}}$ (right), as a function of the column density $N_{{\text{H}}_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 $\mathcal{S}$, ${\mathcal{P}}_{\text{frac}}$, and $\mathcal{S}\times {\mathcal{P}}_{\text{frac}}$; the associated black error bars are ± the standard deviation of each bin.

Fig. 6 Same as Fig. 5, but where the column density in each core has been normalized to the column density peak $N_{{\text{H}}_2,\text{peak}}$. The horizontal dot-dashed gray line indicates the mean of all the $\mathcal{S}\times {\mathcal{P}}_{\text{frac}}$ values, ${0.36}_{-0.17}^{+0.10}$. These uncertainty values are represented by the shaded area, which spans the range between the first and third quartiles of the $\mathcal{S}\times {\mathcal{P}}_{\text{frac}}$ distribution.

4.1.2 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 toward 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 Fig. 6 and the core-by-core analysis of Fig. 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 Sect. 4.1.3 and the simulation and radiative transfer presented in Appendix D), this may contribute to the shallower correlation we find in the $\mathcal{S}$ versus ${\mathcal{P}}_{\text{frac}}$ relation from the ALMA observations with respect to the trend obtained with large-scale cloud observations: indeed, a decrease in the highest values of ${\mathcal{P}}_{\text{frac}}$ would result in a relation closer to $\mathcal{S}\propto {\mathcal{P}}_{\text{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.

4.1.3 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 uv -coverage of the dataset. In contrast to the statistical analysis of dust polarization performed with single dish instruments such as Planck (Planck Collaboration Int. XIX 2015; Planck Collaboration Int. XX 2015; Planck Collaboration XII 2020), 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 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 (Reissl et al. 2016), which produces the Stokes I, Q, and U maps and assumeseither 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). We 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 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 $\mathcal{S}\times {\mathcal{P}}_{\text{frac}}$) can be explained and reproduced with perfect alignment (Planck Collaboration XII 2020; 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 Fig. 7 we present the histograms of $\mathcal{S}$ and ${\mathcal{P}}_{\text{frac}}$ (where all simulations have been merged) before and after filtering, assuming RATs or PA³. We note that if no spatial 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 $\mathcal{S}$ and the polarization fraction ${\mathcal{P}}_{\text{frac}}$ to increase. These increases also translate into an increase in the mean values of $\mathcal{S}\times {\mathcal{P}}_{\text{frac}}$ (see Fig. 8 for the evolution of $\mathcal{S}\times {\mathcal{P}}_{\text{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 Fig. 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 Sect. 4.1.2). It is therefore possible that some of the high values of polarization fraction ${\mathcal{P}}_{\text{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 ${\mathcal{P}}_{\text{frac}}$ given our analysis of both filtering and power spectra, we note that the statistical behavior of ${\mathcal{P}}_{\text{frac}}$ corresponds to what is predicted by theory and models. In addition, the distributions of ${\mathcal{P}}_{\text{frac}}$ we present in Fig. 1 peak at reasonable values (~5%); the high values of ${\mathcal{P}}_{\text{frac}}$ that we find in the ALMA observations are located in the tail of the ${\mathcal{P}}_{\text{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 $\mathcal{S}$ and ${\mathcal{P}}_{\text{frac}}$: from an initial 0.94 before filtering to 0.48 after filtering (see Fig. 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 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 $\mathcal{S}$ and ${\mathcal{P}}_{\text{frac}}$.

Fig. 7 Histograms of polarization fraction ${\mathcal{P}}_{\text{frac}}$ (left) and dispersion of polarization angles $\mathcal{S}$ (right) for all ALMA cores (dashed lines), and for all simulations (solid lines; with grain alignment via perfect alignment [PA] or radiative torques [RATs], and with or without filtering). The histogram lines have been smoothed with a 1D-Gaussian kernel of a size of 0.2% in ${\mathcal{P}}_{\text{frac}}$ and 2° in $\mathcal{S}$. In both panels, the shaded areas correspond to the mean of the uncertainty in ${\mathcal{P}}_{\text{frac}}$ and $\mathcal{S}$ within each bin of the histogram. Right panel: we do plot the errors in $\mathcal{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.

4.2 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 $\mathcal{S}\times {\mathcal{P}}_{\text{frac}}$ as a function of column density, both in ALMA observations (Sects. 4.2.1 and 4.2.2) and synthetic observationsof MHD models of protostellar evolution (Sect. 4.2.3). Finally, we investigate how the statistical properties of the polarized emission from Class 0 protostellar cores may be explained by different dust grain characteristics (Sect. 4.2.4) or additional grain alignment mechanisms (Sect. 4.2.5).

4.2.1 $\mathcal{S}\times {\mathcal{P}}_{\text{frac}}$ suggests no strong radial dependence of the average dust grain alignment efficiency in protostellar envelopes

According to the analytical model developed for the ISM by the Planck collaboration (see Appendix A), the product $\mathcal{S}\times {\mathcal{P}}_{\text{frac}}$ is a proxy for the maximum dust grain alignment efficiency ${\mathcal{P}}_{\text{frac}}{}_{,\text{max}}$⁴, and is statistically independent of the magnetic field configuration. This value of ${\mathcal{P}}_{\text{frac}}{}_{,\text{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 $\mathcal{S}\times {\mathcal{P}}_{\text{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 $\mathcal{S}$ and ${\mathcal{P}}_{\text{frac}}$ are inversely dependent on the disorganization level of the apparent magnetic field, it is reasonable to assume that $\mathcal{S}\times {\mathcal{P}}_{\text{frac}}$ traces the intrinsic capability of dust grains to align themselves with the local magnetic field.

We show in Fig. 6 that the product $\mathcal{S}\times {\mathcal{P}}_{\text{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.17}^{+0.10}$. 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 $\mathcal{S}\times {\mathcal{P}}_{\text{frac}}$ over two orders of magnitude in column density suggests that, within the statistical uncertainties reported in Fig. 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 $\mathcal{S}\times {\mathcal{P}}_{\text{frac}}$ relative to what we see from the ALMA data (see Fig. 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.

Fig. 8 Observed distributions of the mean values of $\mathcal{S}\times {\mathcal{P}}_{\text{frac}}$ as a functionof the column density $N_{{\text{H}}_2}$ (normalized by its maximum value, $N_{{\text{H}}_2,\text{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 onthe 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.

Fig. 9 Distributions of the dispersion of polarization position angles $\mathcal{S}$ as a functionof the polarization fraction ${\mathcal{P}}_{\text{frac}}$ from star formation regions observed by Planck without filtering (left), and with filtering (right). Same as Fig. 4. Interferometric filtering degrades the quality of the $\mathcal{S}$ versus ${\mathcal{P}}_{\text{frac}}$ correlation (smaller R2) and affects its power law index α, which flattens from 0.94 before filtering to 0.48 after filtering.

4.2.2 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. 2017b, 2020; Maury et al. 2018; Le Gouellec et al. 2019). 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; Hull et al. 2020; Le Gouellec et al. 2019). 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 $\mathcal{S}\times {\mathcal{P}}_{\text{frac}}$ could be due to averagingall the observations from our sample, whereas the individual protostars may have very different local conditions at a given normalized column densities because, for example, their luminosity and absolute densities are different. In addition, we note that the stat- istical 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 twoseparations in our datasets. First, we separate the “low” and “high” luminosity cores (i.e., high [bolometric luminosity >10 L_⊙]: NGC1333 IRAS4A, OMC3 MMS6, Serpens Emb 6, IRAS 16293, BHR71; low [bolometric luminosity <10 L_⊙]: Serpens Emb 8, Serpens Emb 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. We 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 Figs. B.1 and B.2 and Appendix B). While the magnetic field lines in outflow cavity walls are very organized (having small values of $\mathcal{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 favorable 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 Fig. 10 we present, for each of these four cases, the evolution of $\mathcal{S}\times {\mathcal{P}}_{\text{frac}}$ as a function of the normalized column density. The trends of $\mathcal{S}\times {\mathcal{P}}_{\text{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 $\mathcal{S}\times {\mathcal{P}}_{\text{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 $\mathcal{S}\times {\mathcal{P}}_{\text{frac}}$ as a function of column density, which indicates less efficient grain alignment in the highest density regions of these cores. As the product $\mathcal{S}\times {\mathcal{P}}_{\text{frac}}$ is also dependent on the function f_m(δ) (see Sect. 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 $\mathcal{S}$ and ${\mathcal{P}}_{\text{frac}}$. However, while this would modify the average values of $\mathcal{S}$ and ${\mathcal{P}}_{\text{frac}}$, it is unclear that this would result in such a strong decrease of $\mathcal{S}\times {\mathcal{P}}_{\text{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 $\mathcal{S}\times {\mathcal{P}}_{\text{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 subsample of low-luminosity 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 peak in column density may be less efficiently aligned in the low-luminosity cores, which causes the drop we see of $\mathcal{S}\times {\mathcal{P}}_{\text{frac}}$ toward these sources. On the contrary, the thermal dust emission emanating from the outer envelope of the low-luminosity cores (i.e., where 10⁻² < $N_{{\text{H}}_2}/N_{{\text{H}}_2,\text{peak}}$ < 10⁻¹) 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 $\mathcal{S}\times {\mathcal{P}}_{\text{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, we note that the overall larger column densities in more massive protostars (which is the case for the high-luminosity cores of our sample) also leads to a larger amount of material being optically thick: it is possible that a decrease of $\mathcal{S}\times {\mathcal{P}}_{\text{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. For example, the central region ~ 200 au of Serpens Emb 6 (Le Gouellec et al. 2019) was found to be likely optically thick at 870 μm; Takahashi et al. (2019) also reported an optically thick ~200 au central region in the high-luminosity source OMC3 MMS6. However, Le Gouellec et al. (2019) reported optically thin central regions in the low-luminosity cores Serpens Emb 8 and Serpens Emb 8(N). 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 low-luminosity simulations, and the inner ~200 au of the high-luminosity simulations) are optically thick. Nevertheless, the trends of $\mathcal{S}\times {\mathcal{P}}_{\text{frac}}$ as a function of column density for the simulations implementing RATs show relatively flat profiles for both the low- and high-luminosity cases (Fig. 8).

Fig. 10 Observed distributions of $\mathcal{S}\times {\mathcal{P}}_{\text{frac}}$ as a functionof the column density $N_{{\text{H}}_2}$, which is normalized to its maximum value $N_{{\text{H}}_2,\text{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.

4.2.3 The role of protostellar luminosity in aligning grains via RATs in MHD models

Analyzing the statistics from the simulations using the mean of $\mathcal{S}\times {\mathcal{P}}_{\text{frac}}$ estimator as a function of $N_{{\text{H}}_2}$ yields the result shown in Fig. 8, where we compare in the left panel the observational distribution presented inFig. 6 with the distributions from all the simulations, implementing RATs or PA, both filtered and not filtered. Similarly, in the right panel of Fig. 8 we compare the same observational distribution with the distributions from the simulations implementing RATs, separating the three simulationswith 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, $\mathcal{S}\times {\mathcal{P}}_{\text{frac}}$ is on average higher in high-luminosity sources than in low-luminosity sources (Fig. 8 right panel). This shows that the statistical estimator $\mathcal{S}\times {\mathcal{P}}_{\text{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, we 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 $\mathcal{S}\times {\mathcal{P}}_{\text{frac}}$ (Fig. 8 left panel).

We 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 disorganizedcomponent of the magnetic field be isotropic in order to justify tracing dust grain alignment efficiency with $\mathcal{S}\times {\mathcal{P}}_{\text{frac}}$. We stress that, however, at all column densities probed in our sample of protostars, the observed $\mathcal{S}\times {\mathcal{P}}_{\text{frac}}$ values are overall larger than those predicted by models with dust grains aligned via RATs, but are consistent with $\mathcal{S}\times {\mathcal{P}}_{\text{frac}}$ values predicted by models where grains are perfectly aligned (see the left panel of Fig. 8). The perfect alignment hypothesis allows us to estimate the typical $\mathcal{S}\times {\mathcal{P}}_{\text{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 (see also Kuffmeier et al. 2020 for a similar result). Only models implementing RATs in high-luminosity cores (see Fig. 8, right panel) produce values that are marginally consistent with the observed values of $\mathcal{S}\times {\mathcal{P}}_{\text{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 $\mathcal{S}\times {\mathcal{P}}_{\text{frac}}$. We stress, however, that implementing different dust-grain properties (see Sect. 4.2.4) as well as different grain alignment mechanisms (see Sect. 4.2.5) could potentially allow models to approach the values of $\mathcal{S}\times {\mathcal{P}}_{\text{frac}}$ seen in ALMA observations.

4.2.4 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 (Hull et al. 2020; Le Gouellec et al. 2019). 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 (Valdivia et al. 2019). In addition, indications of such large grains have been found in multiwavelength 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 $\mathcal{S}\times {\mathcal{P}}_{\text{frac}}$ in protostellar envelopes (see Sect. 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). 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.

4.2.5 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 cannotbe fully reproduced by the synthetic observations of MHD simulations of young, collapsing protostellar cores. The distributions of $\mathcal{S}$ versus ${\mathcal{P}}_{\text{frac}}$ from the simulations (see Figs. D.3 and D.4) show clear trends, but they do not match those found from the ALMA observations (Fig. 4). As mentioned in Sect. 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 $\mathcal{S}\times {\mathcal{P}}_{\text{frac}}$ analysis of our models (see Sect. 4.2.3), we see in Sect. 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 $\mathcal{S}\times {\mathcal{P}}_{\text{frac}}$ from ALMAobservations 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 $\mathcal{S}$ and ${\mathcal{P}}_{\text{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.

5 Conclusions

We perform a statistical analysis of the polarized dust emission emanating from a sample of fifteen ALMA observations toward eleven Class 0 protostellar cores (namely: Serpens Emb 6, Serpens Emb 8, Serpens Emb 8(N), BHR71 IRS1, BHR71 IRS2, B335, IRAS 16293, VLA 1623, L1448 IRS2, OMC3 MMS6, and NGC 1333 IRAS4A), at wavelengths ranging from 870 μm to 3 mm. The conclusions we draw are as follows.

• 1.

  We find a significant correlation between the dispersion of polarization angles $\mathcal{S}$ and polarization fraction ${\mathcal{P}}_{\text{frac}}$ in the polarized dust emission from protostellar envelopes, with a resulting correlation of $\mathcal{S}\propto {\mathcal{P}}_{\text{frac}}{}^{-0.79}$. This correlation is sensitive to the morphology of the disorganized 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\propto {\mathcal{P}}_{\text{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 ${\mathcal{P}}_{\text{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 $\mathcal{S}$ and ${\mathcal{P}}_{\text{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 disorganized component of the field. This results in an overall less disorganized component of the apparent magnetic field compared with that produced by the multiscale turbulence at work in the ISM.

• 3.

  The product $\mathcal{S}\times {\mathcal{P}}_{\text{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.17}^{+0.10}$. This suggests that the grain alignment mechanism producing the polarization observed at millimeter wavelengths, over 3 orders of magnitude in column density (from $N_{{\text{H}}_2}={10}^{22}$ cm⁻² to $N_{{\text{H}}_2}={10}^{25}$ cm⁻²), 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 proto- star of the low-luminosity cores relative to the high-lumino- sity 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 $\mathcal{S}\times {\mathcal{P}}_{\text{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 $\mathcal{S}\times {\mathcal{P}}_{\text{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 $\mathcal{S}\times {\mathcal{P}}_{\text{frac}}$ values than those implementing RATs.

• 7.

  When implementing RAT alignment in our radiative transfer calculations, we do not reproduce with our simulations the $\mathcal{S}$ versus ${\mathcal{P}}_{\text{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 $\mathcal{S}$ versus ${\mathcal{P}}_{\text{frac}}$ correlation is not sensitive to the details of the physical mechanism(s) aligning the dust grains. The values of $\mathcal{S}\times {\mathcal{P}}_{\text{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 ${\mathcal{P}}_{\text{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 $\mathcal{S}$ and ${\mathcal{P}}_{\text{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 comparisonsof 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.

Acknowledgments

V.J.M.L.G. acknowledges the support of the ESO Studentship Program. A.J.M. acknowledges support from the Joint ALMA Observatory Visitor Program, and ESO Visitor program. J.M.G. is supported by the Spanish grant AYA2017-84390-C2-R (AEI/FEDER, UE). C.L.H.H. acknowledges the support of both the NAOJ Fellowship as well as JSPS KAKENHI grants 18K13586 and 20K14527. This work has benefited from the support of the European Research Council under the Horizon 2020 Framework Programme (Starting Grant MagneticYSOs with grant agreement no. 679937). V.J.M.L.G. thanks Maud Galametz and Javier Abril for the precious help they offered for the statistical analysis. This paper makes use of the following ALMA data: ADS/JAO.ALMA# 2013.1.00726.S, ADS/JAO.ALMA#2017.1.00655.S, ADS/JAO.ALMA#2013.1. 01380.S, ADS/JAO.ALMA#2017.1.01392.S, ADS/JAO.ALMA#2018.1.01873.S, ADS/JAO.ALMA#2015.1.01112.S, ADS/JAO.ALMA#2016.1.00604.S, ADS/ JAO.ALMA#2015.1.00341.S, ADS/JAO.ALMA#2016.1.01089.S, ADS/JAO. ALMA#2015.1.00546.S. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), MOST and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. Facilities: ALMA. Software: APLpy, an open-source plotting package for Python hosted at http://aplpy.github.com (Robitaille & Bressert 2012), CASA (McMullin et al. 2007), Astropy (Astropy Collaboration 2018).

Appendix A Planck analytical model

Planck Collaboration XII (2020) 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 ${\mathcal{P}}_{\text{frac},\text{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 $\mathcal{S}$ and the polarization fraction ${\mathcal{P}}_{\text{frac}}$ was found to be the following: $${\mathcal{S}}^2(\delta )=\frac{f_{\text{m}}^2(\delta )}{3N}\frac{{\mathcal{P}}_{\text{frac},\text{max}}^2}{{\mathcal{P}}_{\text{frac}}^2}\mathcal{A},$$(A.1)

with $$\mathcal{A}={\sum }_{i=1}^N\left({\sin }^{2}2\Delta {\varphi }_i{\sin }^2{\Gamma }_i+{\cos }^{2}2\Delta {\varphi }_i{\cos }^2{\Gamma }_i\right),$$(A.2)

where ϕ is the polarization position angle, Δϕ_i = ϕ − ϕ_i, δ is the lag (as introduced in Sect. 3.1, the lag describes the surface over which the dispersion $\mathcal{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 $\overrightarrow{B_i}$ in a given layer i with respect to the plane of the sky. The value of $\mathcal{A}$ is approximated as follows: $${\langle \mathcal{A}\rangle }_{{\mathcal{P}}_{\text{frac}}}~1/\sqrt{2},$$(A.3)

such that we obtain: $${\langle \mathcal{S}\left(\delta \right)\rangle }_{{\mathcal{P}}_{\text{frac}}}\approx \frac{f_{\text{m}}(\delta )}{\sqrt{6N}}\frac{{\mathcal{P}}_{\text{frac},\text{max}}}{{\mathcal{P}}_{\text{frac}}},$$(A.4)

where the 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: $$f_{\text{m}}(\delta )=\frac{{\sigma }_{B_i}(\delta )}{B_i},$$(A.5)

where ${\sigma }_{B_i}(\delta )$ is the fluctuation of the magnetic field $\overrightarrow{B_i}$. This function was modeled as follows: $$f_{\text{m}}(\delta )=0.164f_{\text{M}}{(\frac{\omega }{{160}^{\prime }})}^{-1-{\alpha }_{\text{M}}/2},$$(A.6)

where ω is the full width half maximum of the spatial resolution of the observations. The parameter values used in Planck Collaboration XII (2020) are: α_M = −2.36, f_M = 0.9, N = 7, and P_frac,max = 0.26. Using these values yielded the following analytical relation: $${\langle \mathcal{S}\left(\delta \right)\rangle }_{{\mathcal{P}}_{\text{frac}}}=\frac{0.339}{{\mathcal{P}}_{\text{frac}}}{(\frac{\omega }{{160}^{\prime }})}^{\mathrm{0.18.}}$$(A.7)

In the plots where we relate $\mathcal{S}$ and ${\mathcal{P}}_{\text{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 $\mathcal{S}$ that is proportional to ${\mathcal{P}}_{\text{frac}}{}^{-1}$. However, for our study, the specific case of N = 1 is relevant. This gives: $${\mathcal{S}}^2(\delta )=\frac{f_{\text{m}}^2(\delta )}{3}\frac{{\mathcal{P}}_{\text{frac},\text{max}}^2}{{\mathcal{P}}_{\text{frac}}^2}{\cos }^2\Gamma .$$(A.8)

In this specific case: $${\mathcal{P}}_{\text{frac}}={\mathcal{P}}_{\text{frac},\text{max}}{\cos }^2\Gamma ,$$(A.9)

and thus we obtain: $$\mathcal{S}(\delta )=f_{\text{m}}(\delta )\sqrt{\frac{{\mathcal{P}}_{\text{frac},\text{max}}}{3{\mathcal{P}}_{\text{frac}}}}.$$(A.10)

Consequently, in the specific case of N = 1, $\mathcal{S}$ is proportional to ${\mathcal{P}}_{\text{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 Fig. B.2 we present the maps of Stokes I, polarized intensity P, and polarization angle dispersion $\mathcal{S}$, for all of theALMA 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 Fig. 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.

Fig. B.1 Scheme of the separation we have performed in our polarized intensity maps. The red and yellow areas are defined as the outflow cavity and the envelope emission not associated with the outflow cavities, respectively.

Fig. B.2 Left: total intensity (Stokes I); center: polarized intensity P; right: dispersion of polarization angles $\mathcal{S}$. The solid white contour represents the threshold in Stokes I used to select the data, as described in Sect. 3.1.The white dotted lines denote the cone of the bipolar outflow that we use to separate the statistics in the outflow cavities versus those in the regions of the envelope not associated with the outflow cavities. (see Fig. B.1).

Fig. B.3 Distributions of the dispersion of polarization position angles $\mathcal{S}$ with respect to the polarization fraction ${\mathcal{P}}_{\text{frac}}$ for each of the sources of our sample. Same as Fig. 4, except that here the color scale represents the value of Stokes I of the corresponding point, in Jy str−1.

Fig. B.4 Distributions of the dispersion of polarization position angles $\mathcal{S}$ (blue) and the polarization fraction ${\mathcal{P}}_{\text{frac}}$ (red) as a function of the normalized column density $N_{{\text{H}}_2}/N_{{\text{H}}_2,\text{peak}}$, for all the cores. The solid lines and points represent the running mean of $\mathcal{S}$ and ${\mathcal{P}}_{\text{frac}}$; the associated shaded areas are ± the standard deviation of each bin. The dashed lines are the linear fits, which are linear regressions done in the logarithmic space, whose slopes and uncertainties are calculated the same way as in Fig. 4.

Figure B.3 shows the $\mathcal{S}$ vs. ${\mathcal{P}}_{\text{frac}}$ correlations for all the datasets.

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 Fig. 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 ${\mathcal{P}}_{\text{frac}}$ values discussed in Sects. 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.

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.

Fig. C.1 Normalized power spectra of the three Stokes maps I, Q, and U, toward all the cores of the sample as a function of spatial scale, in au, ranging from the beam size up the maximum recoverable scale. Solid Lines: Stokes I. Dashed and dot-dashed lines: Stokes Q and U, respectively.Each color corresponds to an angular resolution. Red is the original angular resolution. Blue is four times lower resolution (in terms of beam area).

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 non-ideal 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 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 prep.) 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, intermediate-mass 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 (Reissl et al. 2016), which calculate the local dust temperature and dust grain alignment efficiency of oblong dust grains with respect to the magnetic field orientation following the RAT theory developed in Lazarian & Hoang (2007) and 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₀ = 1 (Mathis et al. 1983). 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); we note that this composition 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 or 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; Galametz et al. 2019; Hull et al. 2020; Le Gouellec et al. 2019). 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 POLARIS 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 analyze 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 POLARIS; 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 inSect. 3.1. We denote this latter set of results, “with filtering.” Similar to Fig. B.2, we present in Figs. D.1 and D.2 the maps of Stokes I, polarized intensity P, and dispersion of polarization angles $\mathcal{S}$ for all the projections of the MHD simulations we used, before and after filtering. In these Figures, we only show the perfect alignment case, as the emission in these maps is clearer than in the case of alignment via RATs.

We present in Figs. D.3 and D.4 the distribution of the polarization fraction ${\mathcal{P}}_{\text{frac}}$ as a function of the dispersion of polarization angles $\mathcal{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 ${\mathcal{P}}_{\text{frac}}$ versus $\mathcal{S}$ using PA and RATs, before and after spatial filtering.

Fig. D.1 Same as Fig. B.2, but for the simulations. On the left the cores are filtered, on the right they are not.

Fig. D.2 Same as Fig. B.2, but for the simulations. On the left the cores are filtered, on the right they are not.

In Fig. D.3, we see that both ${\mathcal{P}}_{\text{frac}}$ and $\mathcal{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 $\mathcal{S}$ and ${\mathcal{P}}_{\text{frac}}$ seems to be the closest to the observational correlation presented in Fig. 4 when we consider the lower-mass cases where the simulations have initial turbulence, have perfect grain alignment, and are spatially filtered.

Fig. D.3 Statistics toward all the low-luminosity simulations. Same as Fig. 4. Top-left: perfect alignment before filtering. Top-right: perfect alignment after filtering. Bottom-left: RATs before filtering. Bottom-right: RATs after filtering.

The results from the statistics using the other set of three simulations (see Fig. D.4), which have higher mass and no turbulence, behave differently. The observed $\mathcal{S}$ versus ${\mathcal{P}}_{\text{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 Fig. 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 subcases; however, we still notice that on average, the distributions from the perfect alignment cases tend to have larger values of $\mathcal{S}$ and ${\mathcal{P}}_{\text{frac}}$, which is the expected behavior.

Fig. D.4 Same as Fig. D.3, for all of the high-luminosity simulations.

Appendix E Π Investigations

In the analytical model of Planck Collaboration XII (2020), they demonstrate the $\mathcal{S}\times {\mathcal{P}}_{\text{frac}}$ estimator can trace the ${\mathcal{P}}_{\text{frac},\text{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_jQ − U_jU ≃ P². Finally, the assumption that ⟨S⟩ and $\sqrt{\langle S^2\rangle }$ 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 ${\mathcal{P}}_{\text{frac},\text{max}}$, called Π, the derivation of which does not require these assumptions. This new estimator will be investigated in detail by Guillet et al. (in prep.). In order to derive the relation between Π, ${\mathcal{P}}_{\text{frac},\text{max}}$, and f_m (δ), we follow the same method presented in Appendix A and Appendix E of Planck Collaboration XII (2020): $$\Pi =\frac{1}{2P}\sqrt{\frac{1}{N}{\sum }_j{\left(\frac{Q_{j}U-Q{U}_j}{I_j}\right)}^2}.$$(E.1)

Therefore, we have: $$\Pi ={\mathcal{P}}_{\text{frac},\text{max}}\frac{f_{\text{m}}(\delta )}{\sqrt{6N}}.$$(E.2)

In Figs. E.1–E.3 we show the comparisons between the results provided by Π and $\mathcal{S}\times {\mathcal{P}}_{\text{frac}}$. These plots show that the results are only marginally different, and thus we do not recompute our results using this new estimator.

Fig. E.1 Distributions of $\mathcal{S}\times {\mathcal{P}}_{\text{frac}}$ (left) and Π (right) as a function of the column density $N_{{\text{H}}_2}$, normalized in each core by its maximum value $N_{{\text{H}}_2,\text{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 $\mathcal{S}\times {\mathcal{P}}_{\text{frac}}$ (Π); the associated error bars are ± the standard deviation of each bin. To facilitate the visual comparison, the running means of both $\mathcal{S}\times {\mathcal{P}}_{\text{frac}}$ and Π a re plotted in both panels.

Fig. E.2 Same as Fig. E.1, but comparing the mean (shown in both panels of Fig. E.1) with the median.

Fig. E.3 Left: Stokes I; center: Π; right: $\mathcal{S}\times {\mathcal{P}}_{\text{frac}}$. The white contour represents the threshold in Stokes I used to select the data.

References

All Tables

All Figures

Fig. 1 Histograms of polarization fraction ${\mathcal{P}}_{\text{frac}}$ (left) and polarization position angle ϕ (right) for all the sources of our sample. The histogram lines have been smoothed with a 1D-Gaussian kernel of a size of 0.2% in ${\mathcal{P}}_{\text{frac}}$ and 2° in ϕ. The shaded areas correspond to the mean of the uncertainty in the values of ${\mathcal{P}}_{\text{frac}}$ and ϕ within eachbin of the histograms. Among the two datasets we have toward OMC3 MMS6, “OMC3 MMS6 H” and “OMC3 MMS6 L” denote the high and low resolution observations, respectively.

In the text

Fig. 2 Illustration of how we compute the pixel-selection cutoff in Stokes I in the 1.3 mm observation of the B335 core. Top: evolution of the S/N of the polarized intensity, i.e., P∕σP, as a functionof 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 Sect. 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: $\mathcal{S}\times {\mathcal{P}}_{\text{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. We note that, as expected, one see that decreasing the resolution, and thus increasing the spatial length of the lag, causes on the dispersion $\mathcal{S}$ to increase as well, on average (see Sect. 2.1).

In the text

Fig. 3 Maps from the 1.3 mm observations of the B335 core. Top left: total intensity (Stokes I) thermal dustemission in color scale. The emission is shown from 3σI where σI is the rms noise in the Stokes I map. Top right: polarized intensity P in colorscale.Bottom left: polarization fraction ${\mathcal{P}}_{\text{frac}}$ in colorscale,shown where I > 3σI. Bottom right: dispersion $\mathcal{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 Sect. 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°.

In the text

Fig. 4 Distributions of the dispersion of polarization position angles $\mathcal{S}$ as a functionof the polarization fraction ${\mathcal{P}}_{\text{frac}}$ from all of the datasets merged together. The points were selected according the method developed in Sect. 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 ${\mathcal{P}}_{\text{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$\mathcal{S}$ in the derivation of the linear fit by applying a threshold in polarization fraction, indicated by the vertical dot-dashed gray line. This threshold denotes the ${\mathcal{P}}_{\text{frac}}$ level beyond which the distribution is linear. The solid red line corresponds to the Planck correlation from Planck Collaboration XII (2020), 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 ourALMA 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_{\text{frac}}^{\alpha }$. 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 subsamples of points. We calculate the uncertainties in f and α as standard deviations of Gaussian fits to those histograms.

In the text

Fig. 5 Distributions of the dispersion of polarization position angles $\mathcal{S}$ (top left), the polarization fraction ${\mathcal{P}}_{\text{frac}}$ (bottom left), and $\mathcal{S}\times {\mathcal{P}}_{\text{frac}}$ (right), as a function of the column density $N_{{\text{H}}_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 $\mathcal{S}$, ${\mathcal{P}}_{\text{frac}}$, and $\mathcal{S}\times {\mathcal{P}}_{\text{frac}}$; the associated black error bars are ± the standard deviation of each bin.

In the text

Fig. 6 Same as Fig. 5, but where the column density in each core has been normalized to the column density peak $N_{{\text{H}}_2,\text{peak}}$. The horizontal dot-dashed gray line indicates the mean of all the $\mathcal{S}\times {\mathcal{P}}_{\text{frac}}$ values, ${0.36}_{-0.17}^{+0.10}$. These uncertainty values are represented by the shaded area, which spans the range between the first and third quartiles of the $\mathcal{S}\times {\mathcal{P}}_{\text{frac}}$ distribution.

In the text

Fig. 7 Histograms of polarization fraction ${\mathcal{P}}_{\text{frac}}$ (left) and dispersion of polarization angles $\mathcal{S}$ (right) for all ALMA cores (dashed lines), and for all simulations (solid lines; with grain alignment via perfect alignment [PA] or radiative torques [RATs], and with or without filtering). The histogram lines have been smoothed with a 1D-Gaussian kernel of a size of 0.2% in ${\mathcal{P}}_{\text{frac}}$ and 2° in $\mathcal{S}$. In both panels, the shaded areas correspond to the mean of the uncertainty in ${\mathcal{P}}_{\text{frac}}$ and $\mathcal{S}$ within each bin of the histogram. Right panel: we do plot the errors in $\mathcal{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.

In the text

Fig. 8 Observed distributions of the mean values of $\mathcal{S}\times {\mathcal{P}}_{\text{frac}}$ as a functionof the column density $N_{{\text{H}}_2}$ (normalized by its maximum value, $N_{{\text{H}}_2,\text{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 onthe 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.

In the text

Fig. 9 Distributions of the dispersion of polarization position angles $\mathcal{S}$ as a functionof the polarization fraction ${\mathcal{P}}_{\text{frac}}$ from star formation regions observed by Planck without filtering (left), and with filtering (right). Same as Fig. 4. Interferometric filtering degrades the quality of the $\mathcal{S}$ versus ${\mathcal{P}}_{\text{frac}}$ correlation (smaller R2) and affects its power law index α, which flattens from 0.94 before filtering to 0.48 after filtering.

In the text

Fig. 10 Observed distributions of $\mathcal{S}\times {\mathcal{P}}_{\text{frac}}$ as a functionof the column density $N_{{\text{H}}_2}$, which is normalized to its maximum value $N_{{\text{H}}_2,\text{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.

In the text

	Fig. B.1 Scheme of the separation we have performed in our polarized intensity maps. The red and yellow areas are defined as the outflow cavity and the envelope emission not associated with the outflow cavities, respectively.
In the text

Fig. B.2 Left: total intensity (Stokes I); center: polarized intensity P; right: dispersion of polarization angles $\mathcal{S}$. The solid white contour represents the threshold in Stokes I used to select the data, as described in Sect. 3.1.The white dotted lines denote the cone of the bipolar outflow that we use to separate the statistics in the outflow cavities versus those in the regions of the envelope not associated with the outflow cavities. (see Fig. B.1).

In the text

	Fig. B.3 Distributions of the dispersion of polarization position angles $\mathcal{S}$ with respect to the polarization fraction ${\mathcal{P}}_{\text{frac}}$ for each of the sources of our sample. Same as Fig. 4, except that here the color scale represents the value of Stokes I of the corresponding point, in Jy str−1.
In the text

Fig. B.4 Distributions of the dispersion of polarization position angles $\mathcal{S}$ (blue) and the polarization fraction ${\mathcal{P}}_{\text{frac}}$ (red) as a function of the normalized column density $N_{{\text{H}}_2}/N_{{\text{H}}_2,\text{peak}}$, for all the cores. The solid lines and points represent the running mean of $\mathcal{S}$ and ${\mathcal{P}}_{\text{frac}}$; the associated shaded areas are ± the standard deviation of each bin. The dashed lines are the linear fits, which are linear regressions done in the logarithmic space, whose slopes and uncertainties are calculated the same way as in Fig. 4.

In the text

Fig. C.1 Normalized power spectra of the three Stokes maps I, Q, and U, toward all the cores of the sample as a function of spatial scale, in au, ranging from the beam size up the maximum recoverable scale. Solid Lines: Stokes I. Dashed and dot-dashed lines: Stokes Q and U, respectively.Each color corresponds to an angular resolution. Red is the original angular resolution. Blue is four times lower resolution (in terms of beam area).

In the text

	Fig. D.1 Same as Fig. B.2, but for the simulations. On the left the cores are filtered, on the right they are not.
In the text

	Fig. D.2 Same as Fig. B.2, but for the simulations. On the left the cores are filtered, on the right they are not.
In the text

	Fig. D.3 Statistics toward all the low-luminosity simulations. Same as Fig. 4. Top-left: perfect alignment before filtering. Top-right: perfect alignment after filtering. Bottom-left: RATs before filtering. Bottom-right: RATs after filtering.
In the text

	Fig. D.4 Same as Fig. D.3, for all of the high-luminosity simulations.
In the text

Fig. E.1 Distributions of $\mathcal{S}\times {\mathcal{P}}_{\text{frac}}$ (left) and Π (right) as a function of the column density $N_{{\text{H}}_2}$, normalized in each core by its maximum value $N_{{\text{H}}_2,\text{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 $\mathcal{S}\times {\mathcal{P}}_{\text{frac}}$ (Π); the associated error bars are ± the standard deviation of each bin. To facilitate the visual comparison, the running means of both $\mathcal{S}\times {\mathcal{P}}_{\text{frac}}$ and Π a re plotted in both panels.

In the text

	Fig. E.2 Same as Fig. E.1, but comparing the mean (shown in both panels of Fig. E.1) with the median.
In the text

	Fig. E.3 Left: Stokes I; center: Π; right: $\mathcal{S}\times {\mathcal{P}}_{\text{frac}}$. The white contour represents the threshold in Stokes I used to select the data.
In the text