An observational correlation between magnetic field, angular momentum and fragmentation in the envelopes of Class 0 protostars?

To assess the potential role of magnetic fields in regulating the envelope rotation and the fragmentation of Class 0 protostars, we carried out observations of the dust polarized emission at 0.87 mm with the SMA, in the envelopes of a large sample of 20 Class 0 protostars. We estimate the mean magnetic field orientation over the central 1000 au envelope scales and compared it to that of the protostellar outflow in order to study the relation between their misalignment and the kinematics of the circumstellar gas. We discover a strong relationship between the misalignment of the magnetic field orientation with the outflow and the amount of angular momentum observed at similar scales in the protostellar envelope, revealing a potential link between the kinetic and the magnetic energy at envelope scales. The relation could be driven by favored B misalignments in more dynamical envelopes or a dependence of the envelope dynamics with the large-scale B initial configuration. Comparing the trend with the presence of fragmentation, we observe that single sources are mostly associated with conditions of low angular momentum in the inner envelope and good alignment of the magnetic field with protostellar outflows, at intermediate scales. Our results suggest that the properties of the magnetic field in protostellar envelopes bear a tight relationship with the rotating-infalling gas directly involved in the star and disk formation: we find that it may not only influence the fragmentation of protostellar cores into multiple stellar systems, but also set the conditions establishing the pristine properties of planet-forming disks.


Introduction
In our Galaxy, a majority of stars are found in multiple stellar systems, and a significant fraction of solar-type stars will host planetary systems (Duchêne & Kraus 2013;Hsu et al. 2019). Most of the final stellar mass is collected during a short but vigorous accretion phase. During this so-called protostellar phase, the star forms at the center of an infalling-rotating core, concomitantly with a surrounding disk of gas in circular orbits around the star: while the star will inherit from the majority of the accreted mass, most of the angular momentum contained in the protostellar envelope is expected to be expelled or stored in the protostellar disk, which evolution will eventually lead to protoplanetary systems (Zhao et al. 2020). Class 0 objects are the youngest accreting protostars, surrounded by a dense envelope accreted onto the central protostellar embryo during a short (t < 5×10 4 yr) accretion phase (André et al. 2000;Evans et al. 2009).
Characterizing the dynamics of the gas and the physical processes of these youngest protostars is crucial to understand the efficiency of the star formation process, the global properties of stars in our Galaxy and what sets the conditions allowing disks and planets to form around them.
Magnetic fields (hereafter B) are ubiquitous in the Universe (Vallée 2004), and have been observed to permeate the interstellar material deep down into star-forming cores and protostellar environments (Girart et al. 2006;Hull & Zhang 2019). From a theoretical point of view, the presence of B in star-forming cores has been shown to significantly alter the dynamics of the gas participating in the building of stars during the accretion phase, influencing the resulting properties of these stars and associated circumstellar disks (Terebey et al. 1984;Wurster & Li 2018;Hennebelle & Inutsuka 2019). The mechanism for evacuating angular momentum from the infalling gas through magnetic torques applied by Alfvén waves is called 'magnetic brak-ing'. The initial global collapse, driven by gravity, drags the field lines, leading to an hourglass morphology of the field lines and amplifying the magnetic intensity. The magnetic braking seems to be particularly enhanced by the pinching of field lines that lengthens the magnetic lever arms and efficient at transporting the angular momentum from the inner envelope toward the outer one (Galli et al. 2006). As less angular momentum is transported towards the forming star, only small protoplanetary disks form while the star grows (Allen et al. 2003;Hennebelle & Fromang 2008;Masson et al. 2016;Hirano & Machida 2019).
The importance of the misalignment between magnetic field and rotation axis has been stressed by several authors (Ciardi & Hennebelle 2010;Joos et al. 2012;Gray et al. 2018). They found that in ideal magnetohydrodynamic (MHD) calculations, only small disks or even no disk form in the aligned configurations, when the field is strong enough; and it is comparatively much easier to form a disk in the misaligned case. Numerical simulations taking into account non-ideal MHD effects (such as ambipolar diffusion or Hall effect) were able to overcome the magnetic braking catastrophe, leading to the formation of disks similar to those observed Zhao et al. 2018;Wurster & Bate 2019). Studies by Hennebelle et al. (2020) or Wurster & Lewis (2020) seem to also predict that the misalignment of the magnetic field with the envelope rotation axis directly affects the protostellar disk formation, for instance leading to the formation of larger planet-forming disks in the misaligned cases investigated compared to the smaller disks observed in aligned cases. This is particularly clear if the field intensity is such that the mass-to-flux ratio is on the order of 10.
Another key prediction of magnetized models is that a strong, organized B partly alters the ability of the envelope to fragment, hence suggesting B is one of the regulating agents driving the birth of the multiple stellar systems commonly observed in the Galaxy (Hennebelle & Teyssier 2008).
The influence of the various characteristics of the B field (orientation with the collapse direction, strength) is, however, poorly quantified, observationally speaking. Few studies attempted to test the predicted relationship linking the B-field orientation in protostellar cores to the magnitude of the angular momentum of the gas responsible for disk properties and the formation of multiple stellar systems (see for instance the works from Chapman et al. 2013;Hull et al. 2013;Zhang et al. 2014). Because of the difficulty to trace B in these small embedded astrophysical structures, it has been difficult so far to reach the statistical significance that would allow to draw firm conclusions on the role of magnetic braking to form stars and disks (Yen et al. 2015a;Maury et al. 2018). In order to statistically investigate the B-field orientation, we carried a SubMillimeter Array (SMA) survey of 20 low-mass Class 0 protostars, using 0.87 mm polarized dust emission. Since asymmetric dust particles of the interstellar medium align themselves with their minor axis parallel to the B-field lines (Andersson et al. 2015), the polarized angle observed provides us with a robust proxy of the direction perpendicular to the magnetic field orientation. Class 0 objects were ideal to perform this analysis as most of the mass collapsing onto the central embryo still resides in the envelope, allowing us to trace the B orientation at envelope (1000-2000 au) scales. Galametz et al. (2018) was presenting results on a first subsample of 12 sources, focusing on the properties of polarization fractions and general alignment between B and the outflow at envelope scale. We reported the detection of linearly polarized dust emission in all the objects of the sample. By comparing the B orientation with that of the outflow axis, usually used as a proxy for the rotational axes of these systems, we noticed that at the scales traced in our analysis, the B-field lines were preferentially misaligned in sources where large equatorial velocity gradients were reported in the literature. A potential link between the envelope dynamics and the B orientation could be an additional signature that B has a strong impact of the collapse and fragmentation, as suggested as well by the analysis of massive cores by Zhang et al. (2014).
We complement the Galametz et al. (2018) observations with 8 additional Class 0 envelopes observed with the SMA at 0.87mm at comparable scales. The full SMA B measurements are combined with gas kinematics information obtained homogeneously from N 2 H + observations of velocity gradients in the envelopes, either from the Continuum And Lines in Young ProtoStellar Objects survey (CALYPSO; Maury et al. 2019;Gaudel et al. 2020, 7 sources) or observations published in the literature that we reprocessed if needed (12 sources, see Table  4). Our goal is to observationally test the theoretical predictions of the conditions needed for magnetic braking affecting the collapse and assess the potential role of B in regulating the matter infall and envelope rotation, the formation of disks and the fragmentation in multiple systems. Galametz et al. (2018) presented 12 low-mass Class 0 protostars observed in polarization at 345 GHz with the SMA. We complement this first subsample with 8 additional low-mass protostars. The 20 sources cover a wide range of protostellar properties, with isolated, binary, triple or quadruple systems forming in cores which masses range from 0.2 to 12M . Details of the full sample are provided in Table 1.

The SMA dust polarization observations
The observations (taken in both compact and sub-compact configurations), data reduction and polarization maps of the first 12 low-mass Class 0 protostars observed are presented in Galametz et al. (2018). The observations of the 8 additional low-mass protostars were obtained with the SMA settled in compact configuration in the 345 GHz band (Project 2018B-S015, PI: A. Maury). The antennas used for each observation date are mentioned in Table 2. The polarimeter on the SMA makes use of quarter-wave plate (QWP) in order to convert the linear polarization into circular polarization. The antennas are switched between polarizations (QWP are rotated at various angles) in a coordinated temporal sequence to sample the various combinations of circular polarizations on each baseline. A variety of observational modes (single and dual receiver polarization modes) were used for the observations of Galametz et al. (2018). The dual-receiver full polarization mode, fully commissioned, was then the only mode used to observe the additional 8 targets. The new observations were also taken using the SMA Wideband Astronomical ROACH2 Machine (SWARM) rather than the Application-Specific Integrated Circuit (ASIC) correlator: the added bandwidth has helped increase the SMA sensitivity. The reader will find a detailed description of the SMA polarimeter system in Marrone (2006) and Marrone & Rao (2008). Frequent observations of various calibrators were taken between the target observations to ensure the future gain and polarization calibration. Flux calibrators (Callisto and Neptune) were also observed but will not be used when we performed the flux calibration of the observations (see § 2.3). a These are the sources whose SMA polarization observations are described in this paper. The polarization results for the rest of the sources are described in Galametz et al. (2018), Girart et al. (2014) and Rao et al. (2009). b The total globule mass is probably a factor of 3-5 larger (Stutz et al. 2008

Data reduction, self-calibration and flux calibration
We perform the data reduction on the raw visibilities using the IDL-based software MIR (Millimeter Interferometer Reduction). The calibration includes: an initial flagging of high system temperatures T sys and other wrong visibilities, a bandpass calibration, a correction of the cross-receiver delays, a gain and a flux calibration. The various calibrators observed for each of these steps and the list of antennae used for the observations are summarized in Table 2. Data are then exported to MIRIAD (Sault et al. 1995) 1 for additional processing (i.e. extra flagging) and in particular to perform the instrumental polarization calibration. Quasars were observed to calculate the leakage terms. The continuum data of the targets are used to perform an iterative self-calibration of the Stokes I visibilities. The process is repeated with deeper cleans and shorter intervals until its converges (no rms improvement). We finally use the Atacama Large Millimeter/submillimeter Array (ALMA) calibration source cat-  Galametz et al. (2018). b after self-calibration. c Mean polarization fraction defined as the unweighted ratio between the mean polarization over total flux. alogue 2 to gather the fluxes of quasars we also observed with the SMA (3C 84, 3C 454.3, 3C 273, 3C 279) at similar dates than those of our observations. These fluxes are compared to the SMA amplitudes in order to derive the multiplying factors to apply to the targets visibility amplitude in order to flux-calibrate the dataset. These scaling factors are reported for each date in Table 2. The visibilities covered by the observations range from 15 to 85 kλ.

Deriving the continuum and polarization maps
Stokes parameters are defined as with Q and U the linear polarization and V the circular polarization. We use a robust weighting of 0.5 to transform the visibility data into a dirty map (using the MIRIAD invert task). We produce cleaned images of the various Stokes parameters (using the MIRIAD clean task). Finally, maps of the polarized intensity (debiased) P i , polarization fraction p frac and polarization angle P.A. are produced (using the MIRIAD impol task) as follow: with σ Q,U the average rms of the Q and U maps. We apply a 5-σ cutoff on Stokes I and a 3-σ cutoff on Stokes Q and U in order to only discuss locations where polarized emission is robustly detected. The synthesized beams and rms of the various cleaned maps are provided in Table 3. Maps are produced with a pixel size of 0.6 . In appendix (see Table B.1), we test that this choice does not impact the mean B-field orientation we derive. The Stokes I dust continuum emission maps are shown in Fig. 1. A description of the morphology of this continuum emission as well as details on the source multiplicity are provided in appendix A. The polarized intensity and polarization fraction maps are shown in Fig. D.1. we present, to our knowledge, the 2 https://almascience.eso.org/sc/ first detections of polarized dust emission at envelope scales toward BHR7 and HH25.

Mean magnetic field orientation
The polarization angles are rotated by 90 • to obtain the magnetic field direction. The B vectors obtained for the 20 sources are overlaid on the Stokes I maps in Fig. 1. We note that the strength of SMA observations is two-fold. First, their interferometric nature allows to filter out the large-scale B field permeating the surrounding host cloud, and focus on the fields present in the inner protostellar envelopes.
Second, the modest spatial resolution of our observations allows to cancel out the more complex topology of the field occurring at small (< 500 au) scales, due to both the intense gravitational pull of the infalling material and the launching of protostellar outflows . They should also be less prone than high-angular resolution observations to selectively trace the B field in locations where dust grain alignment efficiency may be highly inhomogeneous, for instance along irradiated cavity walls located very close to the central protostellar objects (see for example LeGouellec et al. 2019). Note that since most Class 0 disks are only contributing at scales much smaller than the scales probed by our beam (less than 25% of Class 0 disks extend beyond 60 au; see Maury et al. 2019), dust polarization due to self-scattering is unlikely to contribute to the polarization observed at envelope scales with the SMA.
To trace the main direction of B at envelope scales, we extract the mean B-field orientation within the central 1000 au region of each source. To perform the calculation, we use the polarization angle and polarization angle error maps produced with our data reduction procedure within the idl/wmean function. The weighted mean is calculated as: with µ the mean position angle of the B field, x i the individual position angles detected within the central 1000 au region and σ i their associated errors. We report the magnetic field position angles in Table 4 and overlaid them (with red segments) on the Stokes I maps for the full sample in Fig. 2. The errors provided in Table 4 are the external uncertainties eu based on the spread 1. B-field vectors (derived from the polarization vectors assuming a 90 • rotation) overlaid as orange segments on the SMA 850 µm Stokes I continuum maps. Color scales are in Jy/beam. Contours at [-3,5,10,20,30,40,50,60,70,80,90,100] σ appear in blue. The filled ellipses on the lower left corner indicate the synthesized beam of the SMA maps.
of the values obtained multiplying the internal uncertainty iu by the square-root of the reduced chi-squared, with: and eu = iu This method is appropriate to calculate a mean B-field angle recovered at 1000-au scales in protostellar envelopes, from individual detections presenting a large dynamic range in Signal-to-Noise ratios (SNRs) -as are some of the polarization detections in the SMA map of each individual source -but also to propagate the individual errors and angle dispersions into an error on the mean value. While most protostars show small dispersions of their individual detections around the mean B field (18 sources out of 20 have dispersions <20 • ), we note that two protostars (IRAS16293A and Per-B1c) present large To quantify the impact of the pixel size used during the map making procedure on the mean B-field orientation, we re-derive the polarization maps for pixels equal to 0.6 (our fiducial pixel scale), 0.7 0.8 0.9 and 1.0 and re-estimate the mean magnetic field position angles (see Table B.1). We observe that the pixel size only has a minor influence of the position angles derived, with uncertainties on the orientation mostly dominated by the dispersion of the position angles within the 1000 au central region. We note that for HH211, the interval of values obtained by changing the pixel size ranges from 162 • to 175 • : the error in position angle may be closer to 10 • than the 5 • we report in Table 4 for this source.
Additionally, because other methods have been developed to calculate mean polarization angles, we propose in Appendix C a simple comparison between the mean angles obtained from our method and three other averaging methods used to analyze single-dish observations of dust polarized emission. Respectively, the three additional methods rely on (i) a simple averaging (no weighting), (ii) an averaging of the Stokes values before computing an angle (for instance used in the single-dish maps of clouds by Li et al. 2006), or (iii) on summing individual Stokes fluxes to remove the variations around a mean value (considering that Stokes Q and U are positive and negative, their sum cancels most variations and should converge towards the most widespread value in the map). This second test shows that the position angle of the mean B field computed with the different methods presents only small variations in the mean position angle, with a dispersion smaller than the error bars from our method reported in Table 4, which indicates that our measurement are robust envelope-scale values. One exception is HH212, for which the interval of values obtained using the various methods ranges from 51 to 59 • . As for HH211, the error on position angle may be closer to 10 • than the 4 • error reported in Table 4 for this source.
The magnetic field orientation of 12 sources of the sample is discussed in Galametz et al. (2018). We add here extra notes concerning the choices made for the current analysis. For IRAS4B, the B field is complex, with an average B-field direc-tion on the eastern part of 149±45 • and an average direction on the western part of 51±19 • . Both sides give a misalignment with the outflow of 30-50 • . In this analysis, we use a misalignment of 40 • for this source. The B-field orientation used for IRAS03292 is the average value detected in the off-centered region where B is detected by the SMA (see Galametz et al. 2018). For L1448-2A, the hourglass shape is resolved with the SMA: we thus used the central vector detected in the TADPOL survey (Hull et al. 2014). Finally, for IRAS16293-A, we estimated the mean magnetic field orientation in a smaller (3 in radius) aperture to avoid contamination by the companion source IRAS16293-B. The average value (173±40 • ) is consistent with the NS orientation also found in Rao et al. (2009), and the uncertainty is large.
The reader can find details on the magnetic field orientation of the extra 8 sources in appendix A. Most of the additional sources present misaligned configurations between the B and outflow position angles. For B1-bS, we note that we decide to use the average NE value for our analysis, as it seems to better trace the magnetic field at envelope scales connected with that traced on larger scales. For L483, the eastern line segment shows discrepancy with the global orientation of the western line segments : we have not included it in our analysis (see the discussion in Appendix A showing this choice does not affect however the correlation we find).
Finally, Table 4 also provides the outflow position angles retrieved from the literature (see references in the table). The outflow position angle uncertainties can vary from one source to another and depend on the outflow inclination and potential overlap of the red-shifted and blue-shifted components. We assume a conservative 10 • error on the outflow position angle for all sources.

Kinematic properties of the protostellar envelopes
The CALYPSO sources -For sources that are part of the CALYPSO sample (i.e. L1157, L1448C, L1448N-B, L1448-2A, IRAS4A, IRAS4B and SVS13-B), Gaudel et al. (2020) recently presented observations of the dense gas kinematics using C 18 O and N 2 H + measurements and derive specific angular momentum estimates throughout their collapsing protostellar envelopes from 50 au to 10000 au scales. In their analysis, velocity maps are produced from a combined Plateau de Bure (PdBI) + 30-m N 2 H + dataset. A hyperfine structure line profile is used to determine the velocity of the molecular line emission in order to produce a velocity map. Then velocity gradients are fitted over a 40 × 40 region around the sources and determined via the least-square minimization, with v grad = v 0 + a∆α + b∆β, with ∆α and ∆β the offsets with respect to the central source (Goodman et al. 1993). Serp SMM18 is also part of the CALYPSO sample but is not included in the sample studied in Gaudel et al. (2020). For this source, we determine the velocity gradient strength and position angle using the same 2D fitting technique. The velocity map is shown in appendix in Fig. E.1. CB230, HH211-mm, IRAS03282 and L483 -For these sources, 2D velocity gradients and position angles on similar scales have been estimated in Tobin et al. (2011) using the N 2 H + line and fitting a plane to the entire velocity field (with R > 10 for L483 and R of about 20-25 for the other sources). We use their velocity gradients and position angles in the following analysis (see their Table 10). L1157 is also part of their sample, with larger ve-locity gradient strength (from 3.5 to 10.5 km s −1 pc −1 ) than those derived in Gaudel et al. (2020) (0.8 km s −1 pc −1 ) and different position angles depending on the use of PdBI, Combined Array for Research in Millimeter-wave Astronomy (CARMA) or Very Large Array (VLA) data. We keep the values from Gaudel et al. (2020) in this paper, but remind the reader that the N 2 H + velocity field is extremely complex for this source, with a velocity gradient direction close to the outflow axis but whose symmetry is perturbed by a redshifted line emission southeast of the protostar, with potential contribution of the outflow cavity walls to the N 2 H + emission for this source (Tobin et al. 2011;Gaudel et al. 2020).
Per B1-bS -We calculate the 2D velocity gradients from the N 2 D + datacube obtained from N. Hirano and presented in Huang & Hirano (2013). Compared to its companion source B1-bN, the B1-bS line profiles are dominated by the V LS R = 6.3 km s −1 component. We fit this hyperfine structure with the hfs cube procedure (Estalella 2017), fitting the line when detected at a 5-σ level. The procedure returns a velocity map from which we determine the gradient on a 20 × 20 region around the source to avoid contamination by the companion. The map of the gas velocities in the envelope of B1-bS, that we use to extract the main velocity gradient in this source, is shown in Fig. E.1. Because of the smaller region fitted, the magnitude of the velocity gradient could potentially only provide an upper limit for this source, as velocity gradients tend to increase when probed at smaller physical radii, as shown for instance by Gaudel et al. (2020) (their Table 2).
Per B1-c -We calculate the 2D velocity gradients from the N 2 H + velocity map obtained from B. Matthews and presented in Matthews et al. (2006). As explained in their analysis, the spectral resolution was not sufficient to separate the hyperfine splitting components so the moment map was taken over an isolated line component of N 2 H + . The velocity map of this source is presented in Fig. E.1. We fit the velocity gradient in a 20 × 20 region around the source to encompass the full N 2 H + emission presented in Matthews et al. (2006) (their Fig. 8). As for Per B1-bS, the magnitude of the velocity gradient might provide an upper limit for this source because of the smaller region fitted.
BHR7-MMS -We calculate the 2D velocity gradients from both the H 2 CO and N 2 D + datacubes obtained from J. Tobin and presented in Tobin et al. (2018). We used Gaussian line profiles to model the H 2 CO emission and the hfs cube procedure (Estalella 2017), fitting the hyperfine structure of the N 2 D + line when detected at a 5-σ level. For H 2 CO, as there is no robust detection beyond a radius of 7 (equivalent to 3000 au for BHR7), we estimated the velocity gradient on a 14 × 14 region around the source and obtain a gradient of 40 km s −1 pc −1 . We note however that in the cold envelope, there might not be enough CO in the gas phase to form H 2 CO: a significant part of the H 2 CO emission could thus come from warmer gas belonging to the outflow, as suggested by the position angle (-19 • , thus aligned with the outflow) of the velocity gradient we derive. In that respect, the N 2 D + emission could better trace the kinematics of the gas in the envelope. The velocity map derived from the N 2 D + datacube is presented in Fig. E.1. Using this data over a 30 × 30 region (equivalent to the ±20 regions used for the Perseus sources), we obtain a weaker, yet still strong velocity gradient of ∼15 km s −1 pc −1 . We choose to use this measurement in the rest of the analysis.
B335 -B335 is a particular case, with very little rotation observed in the source. Menten et al. (1984) only observe, from  NH 3 observations, a velocity shift of a few 10 −2 km s −1 over half an arcmin scale, confirmed by the weak velocity gradient estimated by Caselli et al. (2002). From Saito et al. (1999) (their Fig  5), we derive a velocity gradient of 1 km s −1 pc −1 over the interval ±60 (equivalent to the ±20 regions used in Perseus) in the NS direction (i.e. perpendicular to the outflow). A recent study by Watson (2020), based on the reflection nebulosity of a nearby star, is suggesting that the distance to B335 could be 165 pc (we use 100 pc in this study). This further distance, though it does not impact much the mean B-field position angle derived for the source would lead to a velocity gradient about a factor of 2 lower than that used in this analysis. We note that there is still a lot of uncertainties in the position angle of the velocity gradient for this source, with a C 18 O velocity gradient tilted by 18 • with respect to the EW outflow direction (Yen et al. 2015a,b). We are not analyzing the position angle of the velocity gradient for this source. and SMA are both observing this large velocity gradient that might be partly contaminated by the various outflows emerging from this complex system. If CN observations are used (Antonio Hernández-Gómez's PhD thesis 3 , private communication), then the velocity gradient observed perpendicular to the main EW outflow decreases to 25 km s −1 pc −1 . We use this value in the paper.
HH25-MMS and HH797 -Finally, we could not find observations for HH25 that would allow us to estimate a velocity gradient. In spite of it being in the window of observations, the H 13 CO + line is unfortunately not detected in this object by SMA. For HH797, the complex velocity pattern derived from the C 18 O datacube from Palau et al. (2014) did not allow us to estimate a clean velocity gradient strength or direction. We decide to drop these two sources for the rest of the analysis.
The velocity gradient strengths (in km s −1 pc −1 ) and position angles used in this analysis are summarized in  In both cases, the nature of the extra source needs to be better investigated to be confirmed or infirmed as a companion. Finally, CB230-A's companion seems to host two NIR objects, thus could be part of a triple system (Massi et al. 2008). For this source, we also indicate both the velocity gradient derived by Tobin et al. (2011) from a 1D fit perpendicular to the outflow direction (empty symbol) and a 2D fit of the total velocity field (field symbol), since the two measurements lead to different values of the velocity gradient.
range. This range is consistent with that derived for a sample of 17 nearby protostellar systems by Tobin et al. (2011) (their Fig. 26 right). We note that velocity gradients aligned in the equatorial plane are commonly interpreted as envelope rotation. Recent analysis have revealed a more complex interpretation, with sources showing shifts or even reversal of the gas velocity gradients within envelopes. In some sources, the observed gradients might actually originate from on-going infall, or linked with turbulence (Gaudel et al. 2020). In all cases however, large velocity gradients are tracing more dynamical envelopes with larger kinetic energy.

Misaligned B fields associated to small angular momentum of the protostellar gas
In Galametz et al. (2018), we qualitatively noticed a larger occurrence of misalignment of B-field lines in sources where large velocity gradients were detected in the equatorial plane at thousands of au scales. The measurements of the envelope kinematics combined with B position angles measurements allow us to quantify the relationship. In Fig. 3, we show the projected angle between the B-field orientation and the outflow axis as a function of the velocity gradient (in km s −1 pc −1 ) that is used as a proxy to probe the gas dynamics in the surrounding envelope. Errors on the misalignment angles (x-axis) are the addition of the B-field orientation error quoted in Table 4 and that of the outflow position angle error. The colors used in the plot will be discussed in § 4.2.
We observe a reasonably good (R=0.68) positive correlation between the misalignment of B with respect to the outflow and the strength of the velocity gradient traced at envelope scales. This quantitatively demonstrates that a relation exists between the orientation of the magnetic field and the kinematic energy in envelopes. This is consistent with the result presented in Yen et al. (2015a) and based on a sample of 17 Class0/I protostars where no source with large specific angular momenta were found with a strongly aligned configuration. The correlation observed in Fig. 3 could have various interpretations. One interpretation could be that the misalignments of B at envelope scales are driven by the strong rotational or infall motions of the envelope while another interpretation could be that an aligned B field could have favored the smaller velocity gradients we observe. Protostars are color-coded as a function of their velocity gradient, with increasing gradients as the color darkens. Top right: Relation between the velocity gradient position angle and the outflow direction. We use the same convention for the stripes and lines. We note that the blue regions indicate sources whose velocity gradient might partly be tracing the outflow motion rather than the envelope motions. Bottom: Corresponding histograms of the misalignment of (left) the B field and (right) the outflow axis with respect to the velocity gradient position angle. Colors and lines delineate the same angle offsets as in the top panels.
If the misalignment of B with respect to the outflow at envelope scales would be driven by the envelope kinematics (i.e. if the initially-aligned B-field lines were to be twisted by the infalling or rotating matter), we should expect a relation between the position angle of B and that of the velocity gradient. We plot the velocity gradient position angle as a function of the mean magnetic field position angle in Fig. 4 (left). To ease the visualization, the light and darker stripes indicate where the projected angle between the two directions is smaller than 20 and 10 • respectively. The dashed line indicates a difference of 45 • between the two position angles. We color-code sources as a function of the velocity gradient strength. The bottom left panel shows the corresponding histogram of the angular difference between the position angle of B and the velocity gradient position angle.
We observe that however high (so potentially influential) the velocity gradient, we do not observe an alignment between the magnetic field direction and the velocity gradient position angle. This suggests that the B-field lines do not preferentially follow the direction of the matter infall and collapse dynamics. On the contrary, the sources seem to be scattered over the plot, as highlighted by the relatively flat corresponding histogram (Fig. 4  bottom left). This could indicate that the correlation observed in Fig. 3 has a more complex explanation than the matter infall or rotation misaligning the magnetic field lines. This could favor the second interpretation, namely that in sources showing an aligned configuration of B, magnetic braking could be more efficient at removing angular momentum, leading to the smaller velocity gradient we observe.
Reinforcing further this interpretation that the correlation may be due to more efficient magnetic braking at removing angular momentum in sources initially in an aligned configuration of B, we stress that most sources with small values of their envelope gas velocity gradient have also an envelope B field well aligned with the field observed in their surrounding environment. The large scales magnetic field lines probed around L1157 by Planck or probed at intermediate scales with SHARP (Stephens et al. 2013;Chapman et al. 2013) are consistent with the SMA B orientation. We note that results from optical polarimetry are also consistent, with a small angle offset <20 • (Sharma et al. 2020). For B335, the orientation of the near-IR polarization vectors seems to also fit with the orientation of the sub-mm polarization vectors (Bertrang et al. 2014) as well as with the east-west direction found with JCMT-POL by Yen et al. (2019). Finally, the B-field lines in NGC 1333 IRAS4A, HH211 or L1448-2A seem to be extremely well ordered at scales traced by the SHARP or SCUBA instruments (Hull et al. 2014, see) down to the SMA scales. B1-c does not follow this trend though, with a large scale B field oriented with a position angle of 35 • (Matthews & Wilson 2002), so nearly perpendicular to the outflow direction, compared to the SMA 100 • orientation.
The nature of the gas kinematics recovered at these envelope scales (i.e. whether the angular momentum originates from rotation, infall or even from turbulent motions inherited from the initial conditions and/or turbulence), is unknown. Yet, its amount seems intimately linked to the magnetic configuration in protostellar envelopes.

A hint of higher multiplicity in systems with misaligned magnetic fields ?
In order to investigate the impact of the magnetic field on the envelope fragmentation, we color-code the sources of Fig. 3 as a function of whether or not they are fragmented below 5000 au scales. We indicate in particular whether the source hosts a single, double, triple, and quadruple dust peak (using submm direct imaging). Details on the fragmentation for each individual source are provided in appendix A. We stress that very close multiplicity (namely below 100 au scales) is not considered in our analysis. We observe that sources standing as single objects (in red) seem to mostly reside in the bottom left corner of the plot, i.e. with relatively small velocity gradients of their surrounding envelopes and aligned magnetic field orientation with respect to the outflow axis. Uncertainties also remain on the nature of potential companions detected in HH211 and L483, hence these sources appear with two colors in Fig. 3. In order to assess whether the two populations (single versus multiple) belong to the same distribution, we apply a 2D Kolmogorov-Smirnov test. We use the two python scripts ndtest.py 4 and KS2D.py 5 to estimate the K-S statistics and p-values. These 2D testings are based on statistical methods developed by Peacock (1983) and Fasano & Franceschini (1987). Excluding the sources for which the exact multiplicity nature is unsure (i.e. HH211 and L483), both methods return the same low p-value of 0.13, indicating that the single and multiple source populations are likely not belonging to the same population. More sources would, however, be necessary to reinforce this statistical test. We also note that the two sources SVS13-B and IRAS4B have large scale companions (SVS13-A and IRAS4B2 located at 4200 au and 3200 au respectively from SVS13-B and IRAS4B1) but they do not seem to be fragmented below 3000 au scales, to our knowledge. These sources have moderate velocity gradients and B misalignment, so are both located in the bottom left quadrant of Fig. 3. If our separation 4 https://github.com/syrte/ndtest/blob/master/ndtest.py 5 https://github.com/Gabinou/2DKS/blob/master/KS2D.py criterion is whether or not fragmentation is observed above and below 3000 au scales (compared to 5000 au scales as before), the p-value of the previous K-S test drops to 0.04, so consistent as well with our conclusion that a dichotomy exists between the single versus multiple source population and that the magnetic field alignment might impact how the envelope fragments.
Observational studies already suggested that the magnetic field could influence the fragmentation rate at molecular clouds or filaments scales (for instance Teixeira et al. 2016;Koch et al. 2018). Our analysis at protostellar envelope scales seems to support the theoretical predictions that the magnetic field orientation in the envelope also plays a role in favoring or inhibiting the fragmentation processes of a dense protostellar core into multiple systems. We discuss further the predictions from MHD simulations and how they relate to our results in the following section.

Consistency with predictions from MHD simulations
In ideal MHD models, B fields are shown to strongly regulate the transport of angular momentum and hence modify the final properties of stars and disks (see Hennebelle & Inutsuka 2019;Wurster & Li 2018;Zhao et al. 2020, and references therein). The inclusion of more realistic non-ideal MHD effects changed this view, with finer effects which may play a crucial role regarding the ability of B to interplay with gas kinematics, such as ionization and dust properties (Zhao et al. 2018). Some studies have shown that B is, in particular, less efficient at transporting angular momentum when initially misaligned with the rotation of the system (Joos et al. 2012).
That a link might exist between the initial core-scale magnetic field orientation with respect to the rotation axis driving the small-scale outflow launching and the presence of fragmentation has been predicted by models of protostellar formation. Through the collapse of a dense rotating-infalling core and in the extreme case of an envelope rotation axis initially perpendicular to the magnetic field, Price & Bate (2007) have, for instance, showed that the magnetic tension seems to play an even more significant role in helping fragmentation. Inserting magnetic fields seems to also be a necessary condition to reproduce the fragmentation rates now observed in massive cores with ALMA, as suggested by Fontani et al. (2016). The observations provide information of the main field direction at envelope scales: in an aligned configuration, thus more organized magnetic field configuration, the magnetic pressure will be more efficient at stabilizing the envelope, reducing the rotationally-induced fragmentation at comparable scales (i.e below 5000 au scales). If fragmentation could be favored by less efficient magnetic braking processes, it could also be enhanced by 'turbulent' fragmentation, i.e. fragmentation linked with the increase of the kinematic energy of envelopes in misaligned configurations. Following our observational results, a study of non-ideal MHD models of protostellar collapse, to look for the roots of the correlation we report, was initiated and is currently carried out by the team.
How about observing a potential influence of the misalignment of the B-field orientation on the final protostellar disk sizes? If less angular momentum is transported from the envelope to the disk scales in the case of B-aligned configuration, as suggested by our study at envelope scales, one should expect to form smaller rotationally-supported disks. Observationally, the small size of disks in some Class 0 protostars has been interpreted as a consequence of an efficient magnetic braking potentially disrupting the disk formation (Maury et al. 2010;Yen et al. 2015a;Segura-Cox et al. 2018). MHD simulations from Hirano et al. (2020) recently confirm that in the later accretion phase, the smallest disk radius and mass are produced in alignment configuration cases. Using the CALYPSO sample, Maury et al. (2019) have showed that less than 25% of 26 Class 0 protostars may harbor large protostellar disks resolved at radii >60 au. Their results also favor a magnetized scenario for the disk formation. Unfortunately, most sources showing the presence of large disklike structures (L1527, SMM4, MM22, GF9-2) are not included in the current sample. The disk sizes are mostly unresolved at ∼50 au scales in the remaining sources of our sample overlapping with the CALYPSO sample. This prevents us from studying the impact of the B misalignment on the disk formation itself. Such studies are currently hampered by the need of very high spatial resolution to study those disks, with ALMA for instance. Recent ALMA results from Cox et al. (2018) have, for instance, analyzed the polarization angle dispersion and try to connect the signature/randomness of the B field with the disk/non-disk nature of their sources. The extremely complex morphology of magnetic fields in Class 0 protostars observed with ALMA at 50au typical disk scales may suggest a less dynamically relevant B field at these scales, as expected from non-ideal MHD models (for example with ambipolar diffusion leading to a weaker coupling of the magnetic field lines with the circumstellar gas; Mellon & Li 2009;Tsukamoto et al. 2015). We stress, however, that the characterization of B in disks remains problematic, as several mechanisms can contribute to the production of polarized dust emission at these scales (Kataoka et al. 2015;Cox et al. 2018).

Outflow contamination
Assessing the preferential direction and strength of the moving gas in protostellar envelopes is complicated by the presence of outflows (sometimes more than one in the case of IRAS16293-A for instance) that can contribute to the gradients observed, although our choice of mostly using a dense/cold gas tracer such as N 2 H + should largely limit the contamination. In Fig. 4 (right), we plot the velocity gradient position angle as a function of the outflow direction. The dashed line delineates a difference of 45 • between the two position angles. The light and darker stripes indicate when the projected angle between the two is smaller than 20 and 10 • respectively, thus regions with sources whose velocity gradient might be more related to the outflow dynamics than to the envelope kinematics we are analyzing in this study. The bottom right panel of Fig. 4 provides the histogram of the Outflow -Gradient misalignment. The dispersion in the Outflow -Gradient misalignment is consistent with the results from Tobin et al. (2011) (their Fig. 27). In the end, only 3 sources have a velocity gradient roughly aligned with the outflow direction; i.e. with ∆(P.A.) < 30 • (see Fig. 4 right).

Projection effects
The position and misalignment angles quoted in this analysis are projected in the plane of the sky and might differ from the real angles between the various directions analyzed (velocity gradients, B, outflows). In order to assess the uncertainties linked with projection effects, we perform a Monte-Carlo test. We develop a python script that randomly generates a 3D plane. This plane will be our plane of sky. We then generate 10000 random pairs of vectors in this 3D space. For each pair, vectors are projected onto the 3D plane (using the orthogonal basis of the plane) to estimate the two 'projected vectors' and thus 'observed angle' between the two. The 'real' and 'observed' angles are then compared. Figure 5 shows the density distribution of the observed versus real angles. The 2D density plot is generated with the hexbin python function (Seaborn package). For simplicity, we modify the color bar to indicate the probabilities of each observed-real angle pair. As one could expect, the projected angles are most of the time equal (i.e. on the diagonal) or lower (i.e. on the bottom right area) than the corresponding real angle. Qualitatively, the bright diagonal highlights that we are not making a fundamental mistake by taking the observed angle as a proxy for the real angle.
To quantify the effects, we separate the 'observed' angles in 18 bins of 5 • each and provide the mean 'real angle' in Table 5. The standard deviation for each observed angle bins are provided as uncertainties. We note that the 'true angle' distributions corresponding to each bin are not Gaussian. We observe that the largest discrepancies between the 'real' and 'observed' angles appear for 'observed' angles below 40 • , with large error bars. Consequently for Fig. 3, these projection effects could shift the sources with a B misalignment lower than 40 • to the right, realigning these sources along a more global linear correlation, as the projection effects do not strongly affect sources beyond a 40 • misalignment. Synthetic observations derived from MHD simulations are currently developed by the team, to complement these tests on projection effects (Valdivia et al. in prep.).

Dependency on tracers
One of the caveats of this analysis is also its dependence on the molecular line tracers chosen. Both C 18 O and N 2 H + datacubes are available for the CALYPSO sources. We decided to select N 2 H + as a tracer of the envelope kinematics. Indeed, several analysis showed it as a robust tracer of the outer envelopes, thus the scales we trace with the SMA, rather than of the central regions where this molecule is usually depleted and whose kinematics gas are then usually traced via C 18 O or H 2 D + (Bergin et al. 2002;Anderl et al. 2016;Gaudel et al. 2020;Maret et al. 2020). We used the NH 3 molecule for HH212: the joined analy-  Tobin et al. (2011) has confirmed that both lines are tracing similar physical conditions at the scales we are studying in this analysis. We also selected N 2 D + and CN as complementary tracers of cold dense gas. Thus, unless some strongly different chemical effects are at work in the different protostellar envelopes, all sources have been observed with homogeneous gas tracers. Variations in the opacity from one source to another could bias the physical scales at which the velocity gradients are probed, with the observed N 2 H + arising from projected shells or more external regions than what we assume. N 2 H + or N 2 D + could also be absent near the protostar, although this should not have a strong effect on the scales probed in this analysis. More analysis using other tracers of the gas dynamics would help test the robustness of our conclusions at various envelope radii and/or for gas subject to different local conditions.

Summary
We have carried out dust polarization observations at 0.87 mm for 20 Class 0 protostars with the SMA. We analyze the misalignment of the magnetic field orientation (derived over the central 1000 au scales) with the outflow orientation in order to compare this misalignment with the gas kinematics. The dynamics of the gas is traced via measurements of velocity gradients at envelope scales using molecular line emission. Our analysis provides for the first time a quantification of the striking correlation linking the misalignment of the magnetic field orientation at envelope scales and the angular momentum of the gas reservoir that is directly involved in the star formation process. Comparing the trend with the presence of multiple stellar systems, we also show that sources that tend to stand as single objects mostly reside in environments with weak velocity gradient and/or rather aligned B-field orientation compared to the outflow axis. Altogether, the observations yield towards a coherent picture for the role of magnetic field in forming stars and their protoplanetary disks: they suggest that strong B fields in an aligned configuration may be more efficient in regulating both the gas kinematics and the level of fragmentation during the early, embedded phases of star and disk formation. Our findings could be in line with theoretical expectation from state-of-the art magnetized models of star formation, which predict reduced angular momentum at smaller scales due to magnetic braking, although a thorough exploration of the physical causes behind the observed correlations should be explored in numerical models. Our results provide a strong observational confirmation of the cornerstone role of B to regulate the formation of stellar systems and settle the primordial conditions from which the future disk, star and planets will form. 2019), the B1-b system already shows strong variations in its polarization position angle, with a mean value of the magnetic field direction position angle toward B1-bS of about 30 • . This orientation is consistent with the average B-field orientation obtained with the SMA in the northeast region (26 • ). The western vectors are oriented at 111 • , so perpendicular to the eastern one and along the outflow direction. We decide to use the average NE value for this analysis, as it seems to better trace the magnetic field at envelope scales connected with that traced on larger scales.

B1-c
Continuum morphology at 850 µm -The large scale 2.7 and 3.3 mm continuum emission observed by Matthews et al. (2006) present extensions mostly along the outflow. On the contrary, the 850 µm SMA continuum is flattened in the direction perpendicular to the source outflow. This is consistent with the SMA 1.3 mm dust continuum image presented in Chen et al. (2013). The northern and southwestern plumes are also observed with ALMA at 870 µm (see Cox et al. 2018).
Fragmentation -The source stands as single source from the SMA scales (this work), to the VLA scales probed by the VLA/ALMA Nascent Disk and Multiplicity (VANDAM) survey (Tobin et al. 2016) down to the ALMA scales probed by Cox et al. (2018).
Magnetic field orientation -On the large scales observed by SCUBA (Matthews & Wilson 2002), B1-c has a polarization angle at 35 • , so a magnetic field nearly perpendicular to the outflow direction (estimated at -55 • ; Matthews et al. 2006). Matthews et al. (2008)

BHR7-MMS
Continuum morphology at 850 µm -The SMA 850 µm continuum emission is elongated in the north-south direction. By comparison, the SMA 1.3 mm continuum emission presented in Tobin et al. (2018) is much more compact and marginally extended in the east-west direction but was observed with the Very Extended Configuration, thus a resolution 3-4 times higher than that of the current analysis. Fragmentation -BHR7-MMS is an isolated dark cloud. The source does not seem to be fragmented at the intermediate scales probed with the SMA , and this work).
Recently taken but not yet published ALMA observations of the source should reveal the inner morphology of the source in the coming years.
Magnetic field orientation -To our knowledge, this is the first time that polarized dust emission is used to probe the magnetic field direction in this source. We observe that B is oriented east-west, perpendicular to the direction of the outflow. The northern and southern vectors are slighted tilted compared to the orientation of the central detection, which might be the signature of an hourglass morphology.

HH25-MMS
Continuum morphology at 850 µm -HH25-MMS is part of a string of embedded young stellar objects (the HH24-26 complex; Bontemps et al. 1995;Gibb & Davis 1998). SMA 870 µm observations by Chen et al. (2013) has shown that three distinct sources composed the system (SMM1, SMM2 and SMM3) aligned in the north-south direction. If SMM1 could be the driving source of the HH25 outflow, the SMA polarization data focuses on the SMM2 source. Its SMA continuum seems to be extended in the east-west direction.
Fragmentation -As mentioned previously, SMM2 is part of a 3-source system, with separations of 13 and 11 with SMM1 and SMM3 respectively. On larger scales, HH25-MMS also has a Class I companion, HH26IR, located 1.5 in the southwest.
Magnetic field orientation -Elongated in the north-south direction, the HH25-MMS system is separated in three distinct sources (SMM1, SMM2 and SMM3) resolved by Chen et al. (2013). The SMA observation focuses on SMM2. This is the first known detection of polarized dust emission toward this source. We, however, do not detect polarization toward the very central region but in its western extension. The magnetic field at this position is tilted in the WSW-ENE direction (74 • ), i.e. mainly perpendicular to the projected line connecting the three sources. In this analysis, we do not compare the B-field orientation with the known outflow of the region (revealed by a VLA 3.6 cm survey) since there is strong evidence that the outflow originates from SMM1, not SMM2 (Bontemps et al. 1995;Chen et al. 2013).

HH211-mm
Continuum morphology at 850 µm -As observed at 230 GHz by Gueth & Guilloteau (1999), the HH211-mm submm continuum emission is compact. Also resolved with the SMA, the continuum is slightly elongated in the south-west direction, so perpendicular to the outflow axis. Fragmentation -Higher-resolution SMA observations also performed at 850 µm have revealed that HH211-mm hosts two sources separated by 0.3 (Lee et al. 2009). The first, SMM1 is the protostar from which the collimated outflow originates while SMM2, in its southwest, is responsible for the southwestern extension we observe in our analysis. Modelling the jet wiggle, Lee et al. (2010) also suggested that the HH211-SMM1 source itself could be a proto-binary source with a separation of ∼5 au. Using ALMA data, Lee et al. (2019) have however recently questioned the nature of SMM2 as a secondary source. The companion source is also not detected as part of the VANDAM Perseus survey (Tobin et al. 2016).
Magnetic field orientation -The central field lines traced by the SMA observations have a north-south direction, i.e roughly perpendicular to the outflow direction. This average 175 • orientation is consistent with the results from the TADPOL survey from Hull et al. (2014) and the SCUBA-POL results from Matthews et al. (2009). This is also consistent with the central, northeast and northwest field lines traced at a 0.6 resolution with the SMA by Lee et al. (2014) down to the ALMA scales presented in Lee et al. (2018). Further from the center, the B-field lines seem to re-align in the direction of the outflow axis (see the eastern and western vectors).

HH212
Continuum morphology at 850 µm -The SMA observation reveals a disk-like shape in the direction perpendicular to the outflow, consistent with the continuum emission detected in Wiseman et al. (2001) and observed at 0.9 mm with ALMA (Codella et al. 2014).
Fragmentation -SMA observations at higher resolution that those presented in this paper had suggested the presence of one, if not several, faint sources located around the main body of the HH212-mms source (Lee et al. 2008;Chen et al. 2013). However, more recent ALMA observations by Codella et al. (2014) taken at a similar resolution (0. 5 resolution) do not resolve, nor detect these sources. These fainter sources detected with the SMA are probably detections within HH212 flattened envelope. The VANDAM team, using VLA observations, derived a Toomre's Q parameter for this source consistent with a marginally unstable disk but do not detect multiplicity in the source (Tobin et al. 2020). We thus consider the source as a single source until further analysis confirm or infirm the presence of additional dust peaks.
Magnetic field orientation -The field lines in the central plane are oriented in a NE-SW direction, close to the outflow axis of the source. The southern vectors seem to be following the outflow cavity walls. We note that B traced with ALMA are, on the contrary, perpendicular to the outflow direction (Lee et al. 2018). The authors suggest that this orientation is possibly due to dust polarization arising from self-scattering at these small scales.

L483
Continuum morphology at 850 µm -Our SMA observations reveal a dust continuum elongated in the direction perpendicular to the outflow direction. The southwestern extension is consistent with that already observed by Park et al. (2000) at 3.4 mm.
Fragmentation -L483 seems to be an isolated dense core down to the 2 × 1. 5 scales traced by the SMA observations, as already suggested by Jørgensen (2004) and Chen et al. (2013). Recent 1.2mm ALMA observation resolves the southwestern extension, with a separate continuum source detected at a 3-σ level (Oya et al. 2017). Its nature, however, is not discussed.
Magnetic field orientation -Polarized emission is not detected toward the center of the source but in the southwestern extension (2 consistent detections). Its direction is perpendicular to the SE-NW outflow direction of the source and also perpendicular to the mean B-field orientation at core scales (93 • ) traced at 350 µm using SHARP by Chapman et al. (2013). We also have a detection in the eastern part of the source, although associated with weaker continuum emission (hence its high associated polarization fraction), with a B position angle of 90 • . We do not take this vector into account in our calculation of the magnetic field orientation because of its discrepancy with the other two detections. If we do take the 3 vectors into account, the weighted mean B position angle is 50 • , leading to a misalignment of 55 • between the magnetic field and the outflow orientation. This would put L483 in the center of the correlation observed in Fig. 3 and reinforce the general correlation we observe.

Serpens South MM18
Continuum morphology at 850 µm -The SMA continuum emission is extended in the western and southern part of the source, an extension consistent with that found from the PdBI by Maury et al. (2019).
Fragmentation - Maury et al. (2011) found that SerpS-MM18 is separated in two sources, MM18a being the primary protostar followed by a weaker secondary source MM18b 10 in its southwest, a result confirmed by the observations of the Serpens South complex by Plunkett et al. (2018). Magnetic field orientation -H and Ks-band polarization measurements have shown that the large-scale magnetic field is globally well ordered perpendicular to the main Serpens South filament, thus in an East-West direction (Sugitani et al. 2011). Zooming on SSMM18, we see that the magnetic field lines are also oriented in this EW direction, with lines perpendicular to the outflow axis. The divergence of the eastern vectors could be a signature of an hourglass morphology. associated errors, and emanating from rather homogeneous local conditions and large variations of their associated errors. As a reminder, the reference values used in this analysis are reported in Table 4. The values of the B position angles estimated with the different averaging methods are extremely consistent with each other (with a standard deviation of 5 • between the test and the nominal values), which is not surprising given the small area over which the calculation is performed. These tests show that the main component of B in the protostellar envelopes discussed in this study are, considering the intrinsic limitation of the data in hands, robust envelope-scale values and our handling of the polarization data is meaningful.  [-3,5,10,20,30,40,50,60,70,80,90,100] σ are overlaid in blue. The filled ellipses on the lower left corner indicate the synthesized beam of the SMA maps. Their sizes are reported in Table 3. Second row: Polarization fraction map.