Revealing the diverse magnetic field morphologies in Taurus dense cores with sensitive sub-millimeter polarimetry

We have obtained sensitive dust continuum polarization observations at 850 $\mu$m in the B213 region of Taurus using POL-2 on SCUBA-2 at the James Clerk Maxwell Telescope (JCMT), as part of the BISTRO (B-fields in STar-forming Region Observations) survey. These observations allow us to probe magnetic field (B-field) at high spatial resolution ($\sim$2000 au or $\sim$0.01 pc at 140 pc) in two protostellar cores (K04166 and K04169) and one prestellar core (Miz-8b) that lie within the B213 filament. Using the Davis-Chandrasekhar-Fermi method, we estimate the B-field strengths in K04166, K04169, and Miz-8b to be 38$\pm$14 $\mu$G, 44$\pm$16 $\mu$G, and 12$\pm$5 $\mu$G, respectively. These cores show distinct mean B-field orientations. B-field in K04166 is well ordered and aligned parallel to the orientations of the core minor axis, outflows, core rotation axis, and large-scale uniform B-field, in accordance with magnetically regulated star formation via ambipolar diffusion taking place in K04166. B-field in K04169 is found to be ordered but oriented nearly perpendicular to the core minor axis and large-scale B-field, and not well-correlated with other axes. In contrast, Miz-8b exhibits disordered B-field which show no preferred alignment with the core minor axis or large-scale field. We found that only one core, K04166, retains a memory of the large-scale uniform B-field. The other two cores, K04169 and Miz-8b, are decoupled from the large-scale field. Such a complex B-field configuration could be caused by gas inflow onto the filament, even in the presence of a substantial magnetic flux.


INTRODUCTION
According to the filamentary paradigm of starformation, low-mass stars predominantly form in dense cores which are distributed in a chain-like fashion along gravitationally unstable filamentary clouds (Hartmann 2002;Tafalla & Hacar 2015;Marsh et al. 2016;André et al. 2014). Magnetic field (B-field) is important at all scales during this process (Shu et al. 1987;McKee & Ostriker 2007;Crutcher 2012;Ward-Thompson et al. 2020). Nevertheless, the interplay between B-field, gravity, and turbulence in the formation of cores and their collapse to form stars is still a subject of investigation.
Studies of B-field on cloud scales with Planck 850 µm low-resolution (∼5 arcmin or ∼0.2 at 140 pc) polarization observations and optical and near-infrared (NIR) polarimetry of background stars have revealed that lowdensity gas striations are mostly aligned to B-field and high-density filamentary structures are oriented perpendicular to B-field (Chapman et al. 2011;Alves et al. 2008;Sugitani et al. 2010;Planck Collaboration et al. 2016a;Wang et al. 2020). These observations imply that material can accumulate along field lines and aid in the assembly of dense structures perpendicular to the B-field as a result of gravitational collapse and/or converging flows (see Ballesteros-Paredes et al. 1999a;Hartmann et al. 2001;Soler & Hennebelle 2017).
If the large-scale, uniform B-field is inherited down to core scale (< 0.1 pc), they govern not only the contraction, stability, and collapse of the core (Mestel & Spitzer 1956;Mouschovias & Spitzer 1976), but also the properties of circumstellar disk by help removing angular momentum via magnetic braking (Mouschovias 1991;Allen et al. 2003;Li et al. 2014). According to the theory of isolated, low-mass star-formation via ambipolar diffusion (Mouschovias 1991;Mouschovias et al. 2006), the gravitational collapse of a dense core is regulated by strong, ordered B-field such that the core preferentially contracts along field lines. As a result, the core acquires an oblate-like structure over 10 000 au scales. After gaining sufficient mass via B-field-mediated contraction, initially the subcritical core becomes supercritical and eventually collapses under its own gravity. At this stage, the flux-freezing condition will no longer be valid due to efficient neutral-ion decoupling. As a result of this ambipolar diffusion, the B-field will acquire an hour-glass morphology on protostellar envelope scales, < 1000 au (e.g., Galli & Shu 1993;Girart et al. 2006;Stephens et al. 2013). This model predicts a positive correlation between angle of the mean B-field and that of the minor axes of the filament and core, and the axes of both pseudodisk symmetry and of the bipolar outflow (Fiedler & Mouschovias 1992Galli & Shu 1993;Mocz et al. 2017;Hull & Zhang 2019).
Evidence for magnetically regulated star-formation through observations of coherent B-field across orders of magnitude in size scale (e.g., Li et al. 2006, Li et al. 2009, Hull et al. 2014 is not always the norm. Departure from coherency, especially at smaller scales, can occur in regions dominated by turbulence (e.g., Hull et al. 2017a), shocks from outflows (e.g., Hull et al. 2017b), gravity-driven gas flows (e.g., Pillai et al. 2020), stellar feedback driven by expanding ionization fronts from H ii regions (Arthur et al. 2011;Pattle et al. 2018;Eswaraiah et al. 2020), or gas dynamics arise from gravitational collapse (Ching et al. 2017(Ching et al. , 2018. These observations suggest that the very local environment can determine the morphology and role of B-field. We emphasize here that B-field observations of lowmass dense cores (i) formed out of a single natal filament, (ii) characterized with ordered B-field at larger scales (sub-pc to several pc; see Figure 1(a)), (iii) having signpost of accretion flows (Palmeirim et al. 2013;Shimajiri et al. 2019), and (iv) hosting pristine physical conditions, unaffected by any disruption by strong stellar feedback are sparse. The Taurus B213 is one of these rare regions, making the B213 cores the ideal laboratories to understand the role of B-field in the star formation process.
Here, for the first time, we resolve the B-field in the three cores of B213 on 0.01 -0.1 pc spatial scales. In this letter our key aims are to examine whether (i) Bfield at scales < 0.1 pc are coherent with, or decoupled from, the uniform large-scale B-field, and (ii) the paradigm of magnetically regulated, isolated low-mass star-formation holds in these cores. This paper is organized as follows: Section 2 describes the observations and data reduction. Sections 3 and 4 present the results and discussion, respectively; and Section 5 summarizes our main findings.  Figure 1. (a) The overall structure of the B213 filament is shown using a high-resolution (18. 2) column density map taken from the Herschel Gould Belt Survey (HGBS) archive (Palmeirim et al. 2013). B213 is a part of the 10-pc-long large-scale filament LDN 1495 in the Taurus molecular cloud. Measurements of the large-scale B-field orientation are overlaid: optical polarimetry measurements are shown as magenta segments (Heyer et al. 1987;Heiles 2000;Goodman et al. 1990); NIR polarimetry measurements, as cyan segments (Goodman et al. 1992;Chapman et al. 2011); and Planck 850 µm dust emission polarimetry measurements, as yellow segments. (b) Our inferred core-scale B-field geometry (red segments) superimposed on the area of our total intensity (Stokes I) map which contains the fragmented cores of B213. The extent of the map is shown as a white box in panel (a). Plotted segments correspond to a detection at a minimum of 3σ in polarization fraction and 10σ in total-intensity. Note that all the segments are shown with equal lengths to better display the B-field morphology. The white contour marks a column density of 1×10 22 cm −2 (Palmeirim et al. 2013), as measured in the Herschel data, which outlines the structure of the parental B213 filament, which has fragmented into the cores shown. Red contours drawn in both panels mark 10σ total intensity (Stokes I) values, where 1σ = 1.3 mJy beam −1 is the rms noise. Note that the apparent >10σ intensity values seen in low-column-density regions in panel (a) result from low exposure times at the edges of the POL-2 map, and so delineate the extent of our mapped area. Reference scale and beam size (∼14. 1) are shown.
and 2019 January 08. The two fields, shown in Figure  A2, have a center-to-center angular separation of ∼5 . The fields were each observed 20 times using the POL-2 DAISY mapping mode (Holland et al. 2013;Friberg et al. 2016). This mode results in maps with a 12arcminute diameter, of which the central ∼7 arcmin represents usable coverage, and so these two pointings represent a tightly-spaced mosaic. The observations were made in JCMT weather bands 1 and 2, with 225 GHz atmospheric opacity (τ 225 ) varying between 0.02 and 0.06. The total exposure time for the two fields is ∼28 hrs (14 hrs in each of the two overlapping fields), resulting in one of the deepest observations yet made by the BISTRO Survey. The 850 µm POL-2 data were reduced using the pol2map routine recently added to SMURF (Berry et al. 2005;Chapin et al. 2013) 1 . The final mosaiced maps, 1 http://starlink.eao.hawaii.edu/docs/sun258.htx/sun258ss73.html calibrated in mJy beam −1 , are produced from co-added Stokes I, Q, and U maps with pixel size of 4 , while the final debiased polarization vector catalog is binned to 12 to achieve better sensitivity. The rms noise values in our Stokes I, Q, and U , and P I maps, binned to a pixel size of 12 , are ∼1.3, ∼0.9, ∼0.9, and ∼1.0 mJy beam −1 , respectively. P I represents polarized intensity of the dust emission, debiased using the asymptotic estimator method; our P I map is shown in Figure A2. The instrumental polarization (IP) of POL-2 was corrected for using the 'August 2019' IP model (Friberg et al. 2018). The POL-2 data reduction process is described in detail by Doi et al. (2020) and Pattle et al. (2020).

B-field on small scales
We present the data of 28 polarization measurements satisfying the criteria: (i) the ratio of intensity to its uncertainty I/σ I > 10 and (ii) degree of polarization to its uncertainty P/σ P > 3, where P = P I/I. These measurements are listed in Table 1. The resulting P I within the core boundaries (see Appendix A) ranges from ∼2 to ∼4 mJy beam −1 with a median uncertainty in polarized intensity, σ P I , of 0.64 mJy beam −1 . The polarization fraction ranges from ∼0.8 to ∼18% with a median value of ∼7%. The B213 cores are characterized by weak dust emission (∼12 -318 mJy beam −1 ) as well as weak polarized emission in comparison to the other regions studied regions by the BISTRO program Kwon et al. 2018;Soam et al. 2018;Wang et al. 2019;Coudé et al. 2019;Liu et al. 2019;Pattle et al. 2019;Doi et al. 2020).
Assuming a distance to Taurus of 140 pc, our observations allow us to delineate B-field in B213 on scales ranging from ∼2000 au (∼0.01 pc) up to ∼0.25 pc, the length over which the cores K04166, Miz-8b and K04169 are distributed. The resulting B-field geometry, based on the 28 polarization measurements (cf., Table 1), is shown in Figure 1(b). Since the three cores, T-Tauri, Miz-2, and HGBS-1, have only a single measurement each (and also because the background noise dominates at the locations of these cores -cf., Appendix A), we exclude them from further analysis and discussion. The overall B-field morphology appears to be uniform within K04166 and K04169, but the mean field directions are offset by ∼90 • from one another. In contrast, the B-field morphology in Miz-8b is complex.
We compute the weighted mean position angle of the core B-field,θ core,B , using uncertainties in polarization angle as weights. These values are given in Table 2. Thē θ core,B along with the low-resolution B-field morphology based on Planck 850 µm polarization data, is shown in Figure 2. Table 2 lists the offset betweenθ core,B and the large scale mean B-field orientation (θ largescale B ; see Appendix B) based on multi-wavelength polarimetry. Also listed are the the offset betweenθ core,B and the position angle of each core's major axis (θ core ; see Appendix C).
Interestingly, we see completely different B-field geometry in each of the three cores. The B-field in K04166 lies roughly parallel to the large-scale field (or perpendicular to the filament), while that in K04169 lies roughly perpendicular to the large-scale field (or roughly parallel to the filament). The field direction in Miz-8b lies roughly half-way between the other two, albeit with a larger standard deviation in B-field orientations (35 • ; see Table 2). Hence, we see that the core-scale B-field, in a set of cores spanning ∼ 6 , or ∼0.25 pc, appear to be rather complex.

B-field strength
Based on the assumption that turbulence-induced Alfvén waves can distort B-field orientations, the planeof-sky component of the B-field strength (B pos ) can be estimated using the Davis-Chandrasekhar-Fermi relation (Davis 1951;Chandrasekhar & Fermi 1953, DCF method): where n H2 is the gas number density, δ VNT is the nonthermal gas velocity dispersion, and δ θ is the dispersion in polarization angles about the mean B-field orientation. Q is a factor accounting for line-of-sight and beam dilution effects, which we take as 0.5 based on studies using synthetic polarization maps generated from numerically simulated clouds (Ostriker et al. 2001), which suggest that without this correction factor, DCF-measured B-field strength is overestimated by a factor of 2 when the angular dispersion in the B-field is 25 • .
As illustrated in Appendix C, we have used the 850µm Stokes I map to extract core dimensions, column and number densities, and masses. To quantify the nonthermal velocity dispersion induced by the turbulence, we estimated the average velocity dispersion (δ VLSR ) from archival N 2 H + (1 -0) data (Punanova et al. 2018) 2 , obtained using the IRAM 30-m telescope. The spatial and velocity resolutions of the N 2 H + data are 26. 5 and 0.063 km s −1 , respectively. The thermal contributions to the observed velocity dispersions (δ VT ) are estimated (based on the mean dust temperatures of the cores given in Table 2). These components are quadratically subtracted from the observed velocity dispersions (δ VLSR ) to obtain non-thermal velocity dispersions (δ VNT ). The angular dispersion in the B-field is calculated using the relation for the inverse-variance-weighted standard deviation of B-field (e.g, Wang et al. 2020). These estimated parameters are listed in Table 2.
Using equation 1 and the parameters listed above, the B-field strength is estimated to be 38±14 µG for  Figure 1(b) but now we only show the mean B-field orientation in the three cores K04166, Miz-8b, and K04169, based on the weighted mean of the position angles which we measure. The blue and red dashed arrows denote the protostellar outflows (lengths are not to scale) emanated from K04166 and K04169. The red contour around each core is drawn at I = 13 mJy beam −1 , corresponding to 10σ in total intensity. The large-scale B-field morphology, as determined from the over-sampled Planck 850 µm polarization data (pixel size 1 ), is shown as yellow segments. The white contour is as described in Figure 1 K04166, 44±16 µG for K04169, and 12±5 µG for Miz-8b. Since the majority of the B-field segments in K04166 and K04169 are confined to the core radii ∼20 -50 , the B-field strengths in these cores are mainly valid to the core-envelopes. Further, we caution here that the B-field strength of Miz-8b could be highly uncertain because of the limited number of B-field segments, and hence the larger angular dispersion, used in the DCF method. The current estimations are similar to the B-field strengths of ∼10 -100 µG estimated in relatively unperturbed lowmass star-forming regions (Crutcher et al. 2004;Chapman et al. 2011;Crutcher 2012), and two orders of magnitude less than the ∼1 mG values estimated in massive star forming regions (eg., Curran & Chrysostomou 2007;Hildebrand et al. 2009;Pattle et al. 2017;Liu et al. 2020).
We can use our estimated B-field strength to infer the dynamic state and physical properties of the cores (see Table 2). First, we estimate the magnetic and turbulent pressures using the relations P B = B 2 /8π and P turb = ρδ NT 2 , respectively. Second, we estimate the Alfvénic Mach number using the relation . Third, we use the mass-tomagnetic flux ratio to infer how important is the B-field in comparison to gravity. We measure the mass-to-flux ratio in units of the critical value, as described in Appendix D. Finally, we estimate the rotational energy of each core to determine how rotation may influence the B-field in Appendix E. The derived energy values, along with all other parameters, are listed in Table 2.

DISCUSSION
Since the two protostellar cores, K04166 and K04169, are at a similar evolutionary stage and share similar characteristics (see Table 2), we discuss their energy parameters and gas kinematics with reference to the differences in B-field morphology in Section 4.1. These aspects for the prestellar core Miz-8b are addressed in section 4.2.

K04166 and K04169
Magnetic-to-turbulent pressure ratio is seen to be ∼1 in both cores (see Table 2). This suggests that B-field and turbulence are near equilibrium with each other. Equivalently, the Alfvénic Mach number (∼1) suggests that turbulent motions are trans-Alfvénic. Therefore turbulent motions are not dominant over, and so do not shape the morphology of, the B-field in these cores. The mean mass-to-flux ratio criticality of the cores, λ, is found to be ∼1, suggesting that the core-envelopes may be magnetically critical and marginally supported by B-field. The ratio of rotational to magnetic energy (cf., Appendix E), E rot /E mag 1, which infers that the core rotational energy is too weak to alter the B-field orientation.
Our analysis indicates that there is an equipartition among magnetic, turbulent, and gravitational energies in the core-envelopes of K04166 and K04169. Then the question arises as to why the mean B-field orientation in the two cores are different from each other. We use the morphological correspondence between N 2 H + velocity gradients and B-field, as shown in Figures 3(a) & (b), to shed light on this.
The velocity field in K04166 as inferred from velocity gradient map is well defined, fairly uniform, and is almost perpendicular to the B-field segments. This could be interpreted as bulk core rotation, with the angular momentum (or core rotation-axis with PA ∼11 • ) being parallel to the B-field direction. In addition, the out-flow is well-collimated and exhibits extremely high velocity (EHV) components (Wang et al. 2014), suggesting a possible role of B-field in channeling the outflow and transporting energy and angular momentum away from the rotating circumstellar disk. The position angles of the core (and filament) minor axis (∼37 • ), core rotation-axis (∼11 • ), and bipolar outflows (∼33 • ) are all roughly aligned with bothθ core,B (48 • ) and θ largescale B (29 • ) to within ∼30 • , as shown in Figure 3(d). This strong geometrical correspondence suggests that B-field, which is inherited from the large-scale uniform B-field, have played a significant role in core evolution by allowing gas contraction along field lines to form the core, subsequently governing its collapse via ambipolar diffusion, and finally collimating the outflows. These signatures are in accordance with paradigm of low-mass star-formation process driven by ambipolar diffusion in K04166. However, we could not trace an hour-glass morphology in the inner core (radii < 20 or < 2800 au) due to limited resolution (14. 1 ∼ 2000 au).
On the other hand, the velocity gradient map in K04169 appears to be rather complex, which displays two converging flows, from the northeast and the southwest (Figure 3(b)). Counter-rotation between the disk and the envelope in K04169 is also reported (Takakuwa et al. 2018). We see that core mean B-field (θ core,B ∼121 • ) is nearly aligned parallel to the mean orientation of velocity gradient (θ G ∼126 • ; cf., Table 2). We suggest that this complex gas flows might have altered the B-field from being parallel to the core minor axis in the earlier stage to current perpendicular configuration in K04169. This might have also caused the misalignment of outflows (PA ∼58 • ), core rotation-axis (PA ∼36 • ), and core minor axis (PA ∼36 • ) withθ core,B as shown in Figure 3(e) (see Table 2 for more details).

Miz-8b
Unlike those in K04166 and K04169, Miz-8b has a disordered B-field. It has a magnetic-to-turbulent pressure ratio of ∼0.3 and a Alfvénic Mach number of ∼2 (cf., Table 2), which suggests that turbulence is super-Alfvénic and dominates over B-field. Cores formed in a weakly magnetized, turbulent cloud would have a chaotic Bfield configuration, because of the dominance of turbulent eddies over structural dynamics and field lines (Stone et al. 1998;Ballesteros-Paredes et al. 1999b;Mac Low & Klessen 2004;Li et al. 2014). We suggest that B-field in Miz-8b are complex because of the dominance of turbulent flows (Figure 3(c)). As a result B-field is decoupled from large-scale ordered field (Figure 3(f)). Since the inferred B-field strength in Miz-8b is weaker (12±5 µG), it will not support the core against grav- The integrated velocity ranges are −10 to 0 and 10 to 20 km s −1 for the blue-and red-shifted outflows of K04166, and −5 to 4 and 11 to 17 km s −1 for the blue-and red-shifted outflows of K04169. The lowest and subsequent contour steps are 0.4 and 1.2 K km s −1 , respectively. The angular resolution of the ALMA data is 6.8 × 6.5 .
ity as the mass-to-flux ratio is found to be supercritical (λ ∼3).

SUMMARY
We have performed deep dust polarization observations towards the Taurus B213 filament at 850 µm using SCUBA-2 and POL-2 on the JCMT as part of the BISTRO survey. We successfully detected polarized signal in, and studied in detail the B-field of, two protostellar cores (K04166 and K04169) and one prestellar core (Miz-8b) on scales from 2000 au to 0.25 pc.The main findings of this work are as follows: 1. Despite having (i) ordered B-field on large scales, (ii) quiescent physical conditions, and (iii) formed out of the same natal filament, the three B213 cores exhibit diverse magnetic field properties.
2. Among the three cores, only one, K04166, retains a memory of the large-scale B-field, with a field orientation consistent with those seen on larger scales. The other two cores appear to have decoupled from the large-scale field.
3. Using the Davis-Chandrasekhar-Fermi method, we estimate the B-field strengths in K04166, K04169, and Miz-8b to be 38±14 µG, 44±16 µG, and 12±5 µG, respectively. The associated magnetic energies are in equipartition with both turbulent and gravitational energies in the core envelopes of K04166 and K04169, while being much smaller than the turbulent energy in the core Miz-8b.
4. Based on the correlation between the position angle of the core-scale B-field and the large-scale field, the minor axis of the core, outflows, and the core rotation axis, we suggest that the formation and evolution of K04166 are regulated by B-field, consistent with the paradigm of low-mass starformation via ambipolar diffusion. However, as revealed by their complex velocity fields, the evolution of the other two cores, K04169 and Miz-8b, could be regulated by converging accretion flows and turbulent motions respectively.
We suggest that cores formed in a magnetically regulated molecular cloud may not necessarily retain a memory of the large-scale B-field of the cloud in which they form. Instead, localized differences in gas kinematics, which probably arise due to gas inflows onto the filament, can affect the role of B-field in the star formation process, and the subsequent properties of the forming systems.  27.267900 91.9 ± 1.3 1.3 ± 0.9 3.7 ± 0.9 3.8 ± 0.9 4.2 ± 1.0 125 ± 7 Note-RA and Dec: Celestial coordinates.

Note-I: Total intensity.
Note-Q and U : Stokes parameters.

Note-P : Debiased degree of polarization
Note-θcore,B: B-field orientation, determined by applying an offset of 90 • to the polarization angle.

APPENDIX
A. POLARIZATION PROPERTIES: DETECTION OF WEAKLY POLARIZED DUST EMISSION Figure A1 plots P I versus I for each core, using the selection criterion I/σ I > 10 (gray filled circles). In at least three cores, K04166, Miz-8b, and K04169, and also in the plot showing all of the cores, a slowly increasing trend in Column density (N (H2)) ×10 21 (cm −2 ) 10 ± 6 11 ± 6 5±3 8 Number density (n(H2)) ×10 4 (cm −3 ) 9 ± 5 9 ± 6 5±4 B-field strength, and magnetic and turbulent pressures, etc. Note-a While estimatingθcore,B, one B-field segment associated with K04169, with PA ∼ 37 • , was ignored as it belongs to an another condensation to the south of K04169. Similarly, two segments associated with Miz-8b, with PAs of 24 • and 63 • , were ignored as they fall in the core boundary and do not represent the B-field orientation in Miz-8b. These ignored segments are, within the 30 • , parallel to the large-scale B-field (with PA of 29 • ) .

Note-b:
We also cross-check theθcore,B and corresponding standard deviation values against the circular mean and circular standard deviation values (Doi et al. 2020, see their Appendix C), which are estimated to be 51±19 • , 121±21 • , and 158±35 • for K04166, K04169, and Miz-8b, respectively. These values are quite consistent with the quoted weighted mean and standard deviations.
Note-c: θcore is the position angle of the major axis of the core obtained by fitting an ellipse to the 13 mJy beam −1 POL-2 Stokes I contours of the cores.

B
is the mean large scale mean B-field orientation (29 • ; Table B1).
Note-f : Based on the Herschel Gould Belt Survey (HGBS) column density map.
Note-g: The PA of outflows is determined based on the midline that passes through the center of the bipolar cones. The outflow data is from (Tokuda et al. 2020, see also Ohashi et al. 1997).
Note-h: Core θG is based on the total velocity gradient derived from N2H + data (Punanova et al. 2018). Negative sign corresponds to the angle in clock-wise direction from W to S (see Table B.1 of Punanova et al. 2018).
Note-i, j: Core θG and PA of the core rotation are perpendicular to each other. Note-a: Two measurements with significant deviation in either P or θ are excluded from optical data.
Note-b: Pixels with values < 0.008 KCMB have been excluded from the Planck data, in order to prevent randomisation of our inferred B-field direction by measurements dominated by noise.
Note-θ largescale B ± σ are the mean and standard deviation values resulting from a Gaussian distribution fitted to the data.
Note-Offset PA is difference in angle between the PA of B213 filament (∼133 • ) and the large scale mean B-field (θ largescale B ) (column 5). P I can be seen up to I ∼100 mJy beam −1 , beyond which P I remains approximately constant, although there exist fewer data points.
To extract the reliable data from our POL-2 measurements of the B213 cores, we adopt the selection criterion I/σ I > 10 and P/σ P > 3, which yield 28 polarization measurements (black circles in Figure A1). The resulting median σ P I is 0.64 mJy beam −1 . The P I values lie in the range 1.88 and 3.94 mJy beam −1 , with a median of 2.60 mJy beam −1 , whereas I ranges from 13 to 318 mJy beam −1 as shown in Figure A1. The P values lie in the range 0.82 -17.8% with a minimum ∼1%, median ∼7% and standard deviation ∼5%. Above the 3σ level in P I, a clear detection of polarized intensity (yellow contours) within the core boundaries determined from Stokes I map (red contour at 13 mJy beam −1 ) can be seen in the PI map, as shown in Figure A2. In order to compare core-scale B-field (c.f. Section 3.1) with that in the large-scale, low-density surrounding region, we make use of archival optical, NIR, and Planck/850 µm low-resolution dust polarization data 3 ; and the B-field morphologies inferred from these data sets are shown in Figure 1. We select the data within 1 • of B213; the resulting values of the Gaussian mean and standard deviation in B-field orientation (θ largescale B ± σ) are given in Table B1. There exist a significant number of optical/NIR B-field segments around B213; however, NIR polarization measurements are confined to an area to the west of B213 and therefore may not reveal the local B-field of B213. Visual inspection suggests that optical-and Planck-inferred B-field is ordered. This is confirmed by their mean B-field orientations which are respectively found to be 29±14 • and 29±17 • . These values are nearly identical, with slightly different standard deviations (cf., Table B1), whereas NIR polarization data show a curved morphology, which follows the compressed and curved shell of LDN 1495, with a slightly different mean B-field orientation of 37±17 • . Therefore, to delineate the mean B-field in and around B213, we select the optical and Planck polarization data within 1 • of B213, which yielded a mean orientation of 29 • . This large-scale, coherent B-field (θ largescale B ) with a mean orientation of 29 • , span spatial scales from ∼0.2 pc (∼5 resolution of Planck) to ∼2.4 pc (1 • area around B213).

C. GEOMETRIES, EFFECTIVE RADII, MASSES, AND COLUMN AND NUMBER DENSITIES OF B213 CORES
To estimate various energy and pressure terms for the cores, we extract their masses, column, and number densities from the POL-2 Stokes I map. For this, core dimensions are obtained by fitting the ellipse function mpfitellipse.pro from the Marquardt library to the 10σ Stokes I contours (13 mJy beam −1 ) of each core. The resulting core dimensions Figure A1. Polarized intensity (P I) versus intensity (I) plots for each core, and all of the cores combined. The name of the core is stated in each panel. Gray filled circles represent the data satisfying the criterion I/σI > 10, whereas the black filled circles denote those satisfying both the criteria I/σI > 10 and P/σP ≥ 3. The dashed line represents the median σP I = 0.64 mJy/beam determined from the filled black points. Figure A2. Debiased polarized intensity (PI) map produced using our POL-2 Stokes Q and U maps of the B213 region. Non-smoothed PI contours are drawn at [2, 3, 4]×σP I , where σP I is the rms noise, ∼1 mJy beam −1 (estimated using the pixels in a signal-free region of PI map). Red contour corresponds to a POL-2 total intensity of 13 mJy beam −1 . The P I is nearly zero in the area surrounding the cores but within the cores themselves a clear detection can be seen. Cyan dashed circles mark the areas with diameters of 3 and 7 around the central positions of the two observed fields. Polarization measurements within the smaller circles as well as the common area covered by both larger circles should be useful. Therefore, the measurements of the three cores T-Tauri, Miz-2, and HGBS-1 may not be reliable due to the dominance of background noise at their locations. Each of the six cores are labelled.
(a = semi-major and b = semi-minor), effective radius (R eff = √ ab), and position angle in degrees east of north are given in Table 2.
The integrated fluxes (F ν ) and median dust temperatures (T d ; from the HGBS temperature map) over the core are used to estimate core masses using the relation (Hildebrand 1983): where D = 140 pc, is the distance of B213, κ ν = 0.0125 cm 2 g −1 (e.g., Johnstone et al. 2017) is the dust mass opacity, and B ν (T d ) is Planck function for a blackbody at temperature T d . The uncertainty in mass is estimated by propagating the standard deviation in T d , 10% of the value of F ν as the flux calibration uncertainty of SCUBA2 (Dempsey et al. 2013), and a 50% uncertainty in dust mass opacity (eg., Roy et al. 2014). The column and number densities of the cores are estimated using the following relations: and n(H 2 ) = 3 M 4 µ m H π R 3 eff . (C3) Estimated masses and column and number densities and their corresponding uncertainties are given along with T d and F ν values in Table 2.

D. MASS-TO-FLUX RATIO CRITICALITY
To infer the importance of B-field with respect to the gravity, we estimate the mass-to-magnetic flux ratio in units of the critical value (hereafter mass-to-flux ratio criticality) using the following relation (Crutcher et al. 2004;Chapman et al. 2011), where N (H 2 ) is the mean column density (N(H 2 ) POL2 ), in units of 10 21 cm −2 , along the magnetic flux tube and B tot is the total B-field strength in µG. The critical mass-to-flux ratio, (M/φ) crit = 1/ (4π 2 G) (Nakano & Nakamura 1978), corresponds to the stability criterion for an isothermal gaseous layer threaded by perpendicular B-field. A cloud region with (M/φ) > (M/φ crit ), i.e. λ > 1, will collapse under its own gravity, and so such a cloud is considered to be supercritical. A cloud with µ < 1 will be in a subcritical state because of the significant support rendered by B-field. Taking the mean N (H 2 ) = N (H 2 ) as (10±6)×10 21 , (11±6)×10 21 cm −2 , and (5±3)×10 21 cm −2 and B = B tot a s 38±17, 48±22, and 12±5 µG, we estimate µ values of 2±1, 2±1, and 3±2 for K04166, K04169, and Miz-8b, respectively. However, considering (a) the projection effects between N (H 2 )/B tot and the actual measured N(H 2 )/B (where B is the plane-of-the-sky B-field strength), (ii) B-field being perpendicular to the core elongation in the case of oblate spheroid, or parallel to the core elongation in the case of a prolate spheroid, and and (iii) the assumption that B-field is randomly oriented with respect to the line of sight, the actual value of µ becomes (1/3)λ obs for K04166 as the mean B-field is perpendicular to the core major axis; and (3/4)λ obs for K04169 as mean B-field is parallel to the major axis (Planck Collaboration et al. 2016a, see their Appendix D.4 4 ). No correction was applied on the λ value of Miz-8b because of the misalignment between mean B-field and core major axis. Therefore, the resulting λ values are 0.7±0.5, 1.4±1.0, and 3±1, which are given in Table 2.

E. RATIO OF MAGNETIC-TO-ROTATIONAL ENERGY
By assuming that the cores are uniform density spheres and we measure the rotational and magnetic energies using the following relations (see Wurster & Lewis 2020): and In Equation E5, the correction factor, p = 2(3−A) 3(5−A) = 0.27 (where A is the power-index in the density distribution of the form ρ ∝ r −A , here we we consider A = 1.6), accounts for the density distribution in the sphere (see Xu et al. 2020 for more details). We use the effective radii R = R eff = √ ab (where a = semi-major axis and b = semi-minor axis), and the volume of the core V = (4/3)πR 3 eff . Ω is the angular velocity or magnitude of velocity gradient of the core, measured from the N 2 H + data, and is found to be 2.05±0.02 km s −1 pc −1 for K04166, 3.86±0.04 km s −1 pc −1 for K04169, and 1.88±0.02 km s −1 pc −1 for Miz-8b (Punanova et al. 2018, see their Table B.2). M is the mass of the cores (cf., Appendix C).The derived energy values and their ratio are given in Table 2.

F. MORPHOLOGICAL CORRELATION BETWEEN B-FIELD AND THE GRADIENTS OF VELOCITY
We model the observed line-of-sight centroid velocities of N 2 H + (V LSR , km s −1 ) and the corresponding offset length scales in sky coordinates (R.A (∆ α in pc) and Dec (∆ δ in pc)) around each pixel in terms of velocity gradients in R.A (∇ Vα , km s −1 pc −1 ) and Dec (∇ V δ , km s −1 pc −1 ) and constant systematic velocity of that reference pixel (V 0 , km s −1 ) using the a first-degree bivariate polynomial of the form (Goodman et al. 1993;Henshaw et al. 2016;Sokolov et al. 2018) We have used the IDL algorithm mpfit to perform weighted, non-linear, minimum chi-square fitting to constrain the velocity gradients ∇ Vα and ∇ V δ , and their corresponding uncertainties. These are further used to derive the magnitude (G) and direction (Θ G ) of the velocity gradients using the following relations and We considered at least six adjacent pixels, lying within the beam size of IRAM 30-m telescope for N 2 H + (1 -0), 26. 5, around each pixel when fitting was performed. In addition, we estimate the uncertainties in the VGs using equation 2 of Punanova et al. (2018). These were used as weights while performing the weighted fits. The top panels of Figure 3 show VGs superimposed on the POL-2 Stokes I maps of K04166, K04169, and Miz-8b.