The JCMT BISTRO Survey: Revealing the Diverse Magnetic Field Morphologies in Taurus Dense Cores with Sensitive Submillimeter Polarimetry

According to the filamentary paradigm of star formation, low-mass stars predominantly form in dense cores that are distributed in a chain-like fashion along gravitationally unstable filamentary clouds (Hartmann 2002; André et al. 2014; Tafalla & Hacar 2015; Marsh et al. 2016). The 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 the B-field, gravity, and turbulence in the formation of cores and their collapse to form stars is still a subject of investigation.

Studies of the B-field on cloud scales with Planck 850 μm low-resolution (∼5' or ∼0.2 at 140 pc) polarization observations and optical and near-infrared (NIR) polarimetry of background stars have revealed that low-density gas striations are mostly aligned with the B-field, and high-density filamentary structures are oriented perpendicular to the B-field (Alves et al. 2008; Sugitani et al. 2010; Chapman et al. 2011; 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), it governs not only the contraction, stability, and collapse of the core (Mestel & Spitzer 1956; Mouschovias & Spitzer 1976) but also the properties of the circumstellar disk by helping to remove 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 a 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, the subcritical core initially 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 hourglass 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 the 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 bipolar outflow (Fiedler & Mouschovias 1992, 1993; Galli & Shu 1993; Mocz et al. 2017; Hull & Zhang 2019).

Evidence for magnetically regulated star formation through observations of a coherent B-field across orders of magnitude in size scale (e.g., Li et al. 2006, 2009; Hull et al. 2014) is not always the norm. A departure from coherency, especially at smaller scales, can occur in regions dominated by turbulence (e.g., Hull et al. 2017b), shocks from outflows (e.g., Hull et al. 2017a), 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 arising from gravitational collapse (Ching et al. 2017, 2018). These observations suggest that the very local environment can determine the morphology and role of the B-field.

We emphasize here that B-field observations of low-mass dense cores (i) formed out of a single natal filament, (ii) characterized by an ordered B-field at larger scales (subparsec to several parsecs; see Figure 1(a)), (iii) having signposts 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. Taurus B213 is one of these rare regions, making the B213 cores the ideal laboratories to understand the role of the B-field in the star formation process.

Zoom In Zoom Out Reset image size

We conduct sensitive dust polarization observations at 850 μm toward B213 as part of the B-fields In STar-forming Region Observations (BISTRO; Ward-Thompson et al. 2017) survey to resolve its B-field. BISTRO is a large program on the 15 m James Clerk Maxwell Telescope (JCMT), making use of its SCUBA-2 camera and POL-2 polarimeter. The B213 filament is nearby (distance ∼140 pc; Elias 1978), well studied, and part of the ∼10 pc filament LDN 1495, as shown in Figure 1(a). It is fragmented into a chain of cores that are in the early evolutionary stages of low-mass star formation (Figure 1(b)). These include three prestellar cores, namely, Miz-8b, Miz-2, and HGBS-1 (Mizuno et al. 1994; Marsh et al. 2016); two class 0/I protostellar cores, IRAS 04166+2706 and IRAS 04169+2702 (Ohashi et al. 1997; Tafalla et al. 2010; Takakuwa et al. 2018); and one evolved object, J04194148+2716070, classified as a class II T Tauri star (Davis et al. 2010). We hereafter refer to IRAS 04166+2706 and IRAS 04169+2702 as K01466 and K04169, respectively (see Kenyon et al. 1990, 1993), adopting the core nomenclature of Bracco et al. (2017).

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

The POL-2 observations of two fields in B213 were carried out as part of the JCMT BISTRO survey (JCMT project code M16AL004) between 2017 November 5 and 2019 January 8. 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 12' diameter, of which the central ∼7' represents usable coverage, 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 (τ₂₂₅) varying between 0.02 and 0.06. The total exposure time for the two fields is ∼28 hr (14 hr 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). ⁹¹ The final mosaicked maps, calibrated in millijanskys per beam, are produced from coadded Stokes I, Q, and U maps with a 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, U, and PI maps, binned to a pixel size of 12'', are ∼1.3, ∼0.9, ∼0.9, and ∼1.0 mJy beam⁻¹, respectively. Here PI represents the polarized intensity of the dust emission, debiased using the asymptotic estimator method; our PI map is shown in Figure A2. The instrumental polarization (IP) of POL-2 was corrected for using the "2019 August" 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. (2021).

3.1.  B-field on Small Scales

We present the data of 28 polarization measurements satisfying the following criteria: (i) the ratio of intensity to its uncertainty I/σ_I > 10 and (ii) the degree of polarization to its uncertainty P/σ_P > 3, where P = PI/I. These measurements are listed in Table 1. The resulting PI within the core boundaries (see Appendix A) ranges from ∼2 to ∼4 mJy beam⁻¹ with a median uncertainty in PI, σ_PI, of 0.64 mJy beam⁻¹. 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⁻¹), as well as weak polarized emission in comparison to the other regions studied by the BISTRO program (Ward-Thompson et al. 2017; Kwon et al. 2018; Soam et al. 2018; Coudé et al. 2019; Liu et al. 2019; Pattle et al. 2019; Wang et al. 2019; Doi et al. 2020).

Table 1. Polarization Data along with the Celestial Coordinates of the Pixels

R.A. (J2000)	Decl. (J2000)	I ± σI	Q ± σQ	U ± σU	PI ± σPI	P ± σP	${\theta }_{\mathrm{core},B}\pm {\sigma }_{{\theta }_{\mathrm{core},B}}$
(deg)	(deg)	(mJy beam−1)	(mJy beam−1)	(mJy beam−1)	(mJy beam−1)	(%)	(deg)
(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)
64.982471	27.157917	25.0 ± 0.6	−0.7 ± 0.6	−2.5 ± 0.7	2.5 ± 0.7	9.8 ± 2.8	37 ± 7
65.004946	27.161250	16.9 ± 0.7	1.1 ± 0.7	2.7 ± 0.7	2.8 ± 0.7	16.8 ± 4.4	124 ± 7
64.989963	27.161250	33.0 ± 0.7	0.3 ± 0.6	3.0 ± 0.6	3.0 ± 0.6	9.0 ± 1.9	132 ± 6
64.986217	27.161250	37.3 ± 0.7	−0.7 ± 0.6	3.1 ± 0.7	3.1 ± 0.7	8.3 ± 1.8	141 ± 5
65.004946	27.164583	35.7 ± 0.7	1.1 ± 0.7	2.7 ± 0.7	2.9 ± 0.7	8.1 ± 1.9	124 ± 7
65.001200	27.164583	78.6 ± 0.7	−0.2 ± 0.7	2.7 ± 0.7	2.6 ± 0.7	3.3 ± 0.9	137 ± 7
64.989963	27.164583	113.9 ± 0.9	2.3 ± 0.6	0.3 ± 0.6	2.2 ± 0.6	1.9 ± 0.5	94 ± 8
65.008696	27.167917	20.3 ± 0.6	−1.6 ± 0.7	1.6 ± 0.7	2.2 ± 0.7	10.8 ± 3.3	157 ± 8
64.993708	27.167917	318.4 ± 1.7	2.7 ± 0.7	0.2 ± 0.6	2.6 ± 0.7	0.8 ± 0.2	92 ± 7
64.989963	27.167917	140.9 ± 0.8	3.1 ± 0.6	1.6 ± 0.6	3.5 ± 0.6	2.5 ± 0.4	104 ± 5
64.986212	27.174583	45.4 ± 0.6	1.9 ± 0.6	0.9 ± 0.6	2.0 ± 0.6	4.4 ± 1.3	103 ± 8
64.967479	27.181247	12.2 ± 0.5	−1.5 ± 0.6	−1.6 ± 0.6	2.1 ± 0.6	17.6 ± 5.1	24 ± 8
64.971225	27.184581	15.8 ± 0.6	1.4 ± 0.6	−1.9 ± 0.7	2.3 ± 0.6	14.4 ± 4.0	63 ± 7
64.959979	27.191247	46.1 ± 0.6	0.2 ± 0.6	2.0 ± 0.6	1.9 ± 0.6	4.1 ± 1.3	132 ± 9
64.971221	27.194583	26.1 ± 0.6	0.2 ± 0.6	2.9 ± 0.6	2.9 ± 0.6	11.0 ± 2.5	133 ± 6
64.967475	27.194581	48.4 ± 0.6	−2.0 ± 0.6	0.4 ± 0.6	1.9 ± 0.6	4.0 ± 1.3	174 ± 9
64.963729	27.194581	58.6 ± 0.6	−1.5 ± 0.6	−1.8 ± 0.6	2.3 ± 0.6	3.8 ± 1.0	25 ± 7
64.933733	27.217906	54.7 ± 0.6	−0.8 ± 0.6	−2.4 ± 0.6	2.5 ± 0.6	4.6 ± 1.1	36 ± 6
64.926237	27.217903	25.0 ± 0.6	−0.6 ± 0.5	−2.7 ± 0.7	2.7 ± 0.7	10.8 ± 2.7	39 ± 6
64.926237	27.221236	72.0 ± 0.6	0.4 ± 0.6	−2.9 ± 0.6	2.8 ± 0.6	3.9 ± 0.8	49 ± 6
64.937479	27.227906	37.7 ± 0.6	1.8 ± 0.6	−1.9 ± 0.7	2.5 ± 0.6	6.7 ± 1.7	66 ± 7
64.933729	27.227906	113.4 ± 0.7	2.2 ± 0.6	−0.0 ± 0.7	2.1 ± 0.6	1.8 ± 0.6	90 ± 9
64.929979	27.227903	290.0 ± 1.4	1.2 ± 0.6	−2.3 ± 0.6	2.5 ± 0.6	0.9 ± 0.2	59 ± 7
64.926233	27.227903	253.9 ± 1.5	−1.8 ± 0.6	−2.2 ± 0.7	2.8 ± 0.6	1.1 ± 0.3	25 ± 6
64.933725	27.234572	14.6 ± 0.6	0.5 ± 0.7	−2.2 ± 0.7	2.2 ± 0.7	14.9 ± 4.5	51 ± 9
64.854983	27.247850	29.3 ± 1.1	3.6 ± 1.0	−2.0 ± 1.0	3.9 ± 1.0	13.4 ± 3.5	76 ± 7
64.907471	27.251225	22.1 ± 0.8	4.0 ± 0.8	−0.5 ± 0.8	3.9 ± 0.8	17.8 ± 3.8	86 ± 6
64.922458	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. R.A. and decl: celestial coordinates. I: total intensity. Q and U: Stokes parameters. PI: debiased polarized intensity. P: debiased degree of polarization. ${\theta }_{\mathrm{core},B}$: B-field orientation, determined by applying an offset of 90° to the polarization angle.

Download table as:  ASCIITypeset image

Assuming a distance to Taurus of 140 pc, our observations allow us to delineate the 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 (see 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; see 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 (PA) of the core B-field, ${\bar{\theta }}_{\mathrm{core},B}$, using uncertainties in polarization angle as weights. These values are given in Table 2. The ${\bar{\theta }}_{\mathrm{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 ${\bar{\theta }}_{\mathrm{core},B}$ and the large-scale mean B-field orientation (${\theta }_{B}^{\mathrm{largescale}};$ see Appendix B) based on multiwavelength polarimetry. Also listed are the offset between ${\bar{\theta }}_{\mathrm{core},B}$ and the PA of each core's major axis (${\theta }_{\mathrm{core}};$ see Appendix C).

Zoom In Zoom Out Reset image size

Table 2. Various Parameters for K04166, K04169, and Miz-8b

No.	Parameter	K04166	K04169	Miz-8b
	Weighted Mean B-field Orientation along with Various Offset PAs
1	No. of B-field segments	8	10	4
2	Weighted mean B-field orientation (${\bar{\theta }}_{\mathrm{core},B}$ $\pm {\epsilon }_{{\bar{\theta }}_{{\rm{c}}{\rm{o}}{\rm{r}}{\rm{e}},B}}$; deg) a , b	48 ± 2	121 ± 2	158 ± 4
3	Angular dispersion (δθ ; deg)	18 ± 4	20 ± 3	35 ± 7
4	Standard deviation in ${\theta }_{\mathrm{core},B}$(deg) b	20	33	35
5	PA of core major axis (${\theta }_{\mathrm{core}};$ deg) c	127 ± 20	126 ± 16	119 ± 30
6	$| {\bar{\theta }}_{\mathrm{core},B}-{\theta }_{B}^{\mathrm{largescale}}|$ (deg) d	19	88	51
7	$| {\bar{\theta }}_{\mathrm{core},B}-{\theta }_{\mathrm{core}}|$ (deg)	79	5	39
	Core Dimensions, Mass, Column and Number Densities, etc.
1	Semimajor axis (a) (pc)	0.032 ± 0.005	0.036 ± 0.006	0.029 ± 0.010
2	Semiminor axis (b) (pc)	0.024 ± 0.003	0.025 ± 0.005	0.023 ± 0.005
3	Effective radius (Reff) (pc) e	0.028 ± 0.003	0.030 ± 0.004	0.026 ± 0.005
4	Median dust temperature (Td) (K) f	12.2 ± 0.4	12.2 ± 0.8	10.9 ± 0.2
5	Integrated flux (Fν ) (mJy)	1564	1969	516
6	Mass (M) (M☉)	0.55 ± 0.28	0.69 ± 0.36	0.22 ± 0.11
7	Column density (N(H2)) × 1021 (cm−2)	10 ± 6	11 ± 6	5 ± 3
8	Number density (n(H2)) × 104 (cm−3)	9 ± 5	9 ± 6	5 ± 4
	B-field Strength and Magnetic and Turbulent Pressures, etc.
1	B-field strength (using DCF method) (μG)	38 ± 14	44 ± 16	12 ± 5
2	B-field pressure (PB ; ×10−10 dyn cm−2)	0.6 ± 0.4	0.8 ± 0.6	0.06 ± 0.05
3	Turbulent pressure (Pturb; ×10−10 dyn cm−2)	0.4 ± 0.3	0.7 ± 0.5	0.2 ± 0.1
4	PB /Pturb	1.3 ± 1.3	1.0 ± 1.0	0.3 ± 0.4
5	Mass-to-flux ratio criticality (λ = (M/ϕ)/(M/ϕ)cri)	0.7 ± 0.4	1.4 ± 1.0	3 ± 2
6	Alfvén velocity (VA; km s−1)	0.17 ± 0.09	0.21 ± 0.12	0.07 ± 0.04
7	Alfvénic Mach number (MA)	1.1 ± 0.6	1.2 ± 0.7	2 ± 1
	Energy Parameters
1	p	0.27	0.27	0.27
2	Rotational energy Erot (×1041 erg)	0.02	0.1	0.01
3	Magnetic energy Emag (×1041 erg)	2	3	0.1
4	Erot/Emag	0.01	0.04	0.05
	Various PAs
1	Core minor axis PA (deg)	∼37	∼36	∼29
2	Outflow PA (deg) g	∼33	∼58	⋯
3	Core θG (deg; from W to N) h	−169	−144	93
4	Core θG (deg; from N to E) i	101	126	3
5	PA of the core rotation axis (deg; from N to E) i	11	36	93

Notes.

^a While estimating ${\bar{\theta }}_{\mathrm{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 30°, parallel to the large-scale B-field (with a PA of 29°). ^b We also cross-check the ${\bar{\theta }}_{\mathrm{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. ^c Here ${\theta }_{\mathrm{core}}$ is the PA of the major axis of the core obtained by fitting an ellipse to the 13 mJy beam⁻¹ POL-2 Stokes I contours of the cores. ^d Here ${\theta }_{B}^{\mathrm{largescale}}$ is the mean large-scale mean B-field orientation (29°; Table B1). ^e R_eff= $\sqrt{{ab}}$. ^f Based on the HGBS column density map. ^g The PA of the outflows is determined based on the midline that passes through the center of the bipolar cones. The outflow data are from Tokuda et al. (2020; see also Ohashi et al. 1997). ^h The core θ_G is based on the total velocity gradient derived from N₂H⁺ data (Punanova et al. 2018). The negative sign corresponds to the angle in the clockwise direction from W to S (see Table B.1 of Punanova et al. 2018). ⁱ The core θ_G and PA of the core rotation are perpendicular to each other.

Download table as:  ASCIITypeset image

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 halfway 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, appears to be rather complex.

Furthermore, we observe a good alignment between the core B-field (∼48°) and outflows (∼33°) in K04166 (Figures 2 and 3(d)), consistent with studies by Davidson et al. (2011) and Chapman et al. (2011). In contrast, we see a misalignment between the core mean B-field (∼121°) and the outflows (∼58°) in K04169 (Figures 2 and 3(e)), in accordance with studies by Hull et al. (2013, 2014), Hull & Zhang (2019), and Yen et al. (2021). Despite the fact that these two cores lie within ∼0.25 pc of each other, they exhibit different B-field/outflow orientations.

Zoom In Zoom Out Reset image size

3.2.  B-field Strength

Based on the assumption that turbulence-induced Alfvén waves can distort B-field orientations, the plane-of-the-sky component of the B-field strength (B_pos) can be estimated using the Davis–Chandrasekhar–Fermi relation (DCF method; Davis 1951; Chandrasekhar & Fermi 1953),

$$\begin{eqnarray}&&{B}_{\mathrm{pos}}\simeq Q\sqrt{4\pi \,\mu {m}_{{\rm{H}}}\,{n}_{{{\rm{H}}}_{2}}}\left(\displaystyle \frac{{\delta }_{{V}_{\mathrm{NT}}}}{{\delta }_{\theta }}\right),\end{eqnarray} \tag{ 1 }$$

where ${n}_{{{\rm{H}}}_{2}}$ is the gas number density, ${\delta }_{{V}_{\mathrm{NT}}}$ is the nonthermal gas velocity dispersion, and δ_θ is the dispersion in polarization angles about the mean B-field orientation. Here 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). This suggests that without this correction factor, the 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 (${\delta }_{{V}_{\mathrm{LSR}}}$) from archival N₂H⁺ (1–0) data (Punanova et al. 2018) ⁹² obtained using the IRAM 30 m telescope. The spatial and velocity resolutions of the N₂H⁺ data are 265 and 0.063 km s⁻¹, respectively. The thermal contributions to the observed velocity dispersions (${\delta }_{{V}_{T}}$) 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 (${\delta }_{{V}_{\mathrm{LSR}}}$) to obtain nonthermal velocity dispersions (${\delta }_{{V}_{\mathrm{NT}}}$). The angular dispersion in the B-field is calculated using the relation for the inverse variance–weighted standard deviation of the 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 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 of ∼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 low-mass star-forming regions (Crutcher et al. 2004; Chapman et al. 2011; Crutcher 2012) and 2 orders of magnitude less than the ∼1 mG values estimated in massive star-forming regions (e.g., 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²/8π and P_turb= $\rho {\delta }_{\mathrm{NT}}^{2}$, respectively. Second, we estimate the Alfvénic Mach number using the relation M_A = $\sqrt{3}\left(\tfrac{{\delta }_{\mathrm{NT}}}{{V}_{{\rm{A}}}}\right)$, where Alfvén velocity ${V}_{{\rm{A}}}=\tfrac{{B}_{\mathrm{los}}}{\sqrt{4\pi \,\rho }}$ (where $\rho ={n}_{{{\rm{H}}}_{2}}\,\mu \,{m}_{{\rm{H}}}$). Third, we use the mass-to-magnetic flux ratio to infer how important the B-field is 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.

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.

4.1. K04166 and K04169

4.2. Miz-8b

We have performed deep dust polarization observations toward 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-fields on large scales and (ii) quiescent physical conditions and (iii) being 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 DCF method, we estimate the B-field strengths in K04166, K04169, and Miz-8b to be 38 ± 14, 44 ± 16, 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 of Miz-8b.
• 4.  
  Based on the correlation between the PAs of the core-scale B-field, 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 the B-field, consistent with the paradigm of low-mass star formation 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 the B-field in the star formation process and the subsequent properties of the forming systems.

We thank the referee for constructive suggestions that have improved the content and flow of this paper. This work is supported by National Natural Science Foundation of China (NSFC) grant Nos. 11988101, 11725313, and U1931117 and the International Partnership Program of Chinese Academy of Sciences grant No. 114A11KYSB20160008. This work is also supported by Special Funding for Advanced Users, budgeted and administrated by the Center for Astronomical Mega-Science, Chinese Academy of Sciences (CAMS). C.E. thanks Tokuda Kazuki for kindly providing us with the ALMA data on the protostellar outflows. C.E. thanks Sihan Jiao and Yuehui Ma for help with the analyses. K.Q. acknowledges support from the National Natural Science Foundation of China (NSFC) through grant U1731237. C.W.L. is supported by the Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Education, Science and Technology (NRF-2019R1A2C1010851). N.O. is supported by Ministry of Science and Technology (MOST) grant MOST 109-2112-M-001-051 in Taiwan. C.L.H.H. acknowledges the support of the NAOJ Fellowship and JSPS KAKENHI grants 18K13586 and 20K14527. W.K. was supported by the New Faculty Startup Fund from Seoul National University. D.J. is supported by the National Research Council of Canada and a Natural Sciences and Engineering Research Council of Canada (NSERC) Discovery Grant. A.S. acknowledges support from the NSF through grant AST-1715876. K.H.K is supported by the center for Women In Science, Engineering and Technology (WISET) grant funded by the Ministry of Science and ICT (MSIT) under the program for returners into R&D (WISET-2019-288; WISET-2020-247). P.N.D and N.B.N are funded by Vietnam National Foundation for Science and Technology Development (NAFOSTED) under grant number 103.99-2019.368. Y.W.T acknowledges the grant MOST 108-2112-M-001-004-MY2.

The James Clerk Maxwell Telescope is operated by the East Asian Observatory on behalf of the National Astronomical Observatory of Japan, the Academia Sinica Institute of Astronomy and Astrophysics, the Korea Astronomy and Space Science Institute, and the Center for Astronomical Mega-Science. Additional funding support is provided by the Science and Technology Facilities Council of the United Kingdom and participating universities in the United Kingdom, Canada, and Ireland.

Facility: JCMT. -

Figure A1 plots PI 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 PI can be seen up to I ∼ 100 mJy beam⁻¹, beyond which PI remains approximately constant, although there exist fewer data points.

Zoom In Zoom Out Reset image size

To extract the reliable data from our POL-2 measurements of the B213 cores, we adopt the selection criteria I/σ_I > 10 and P/σ_P > 3, which yield 28 polarization measurements (black filled circles in Figure A1). The resulting median σ_PI is 0.64 mJy beam⁻¹. The PI values lie in the range 1.88–3.94 mJy beam⁻¹, with a median of 2.60 mJy beam⁻¹, whereas I ranges from 13 to 318 mJy beam⁻¹, as shown in Figure A1. The P values lie in the range 0.82%–17.8% with a minimum of ∼1%, median of ∼7%, and standard deviation of ∼5%. Above the 3σ level in PI, a clear detection of PI (yellow contours) within the core boundaries determined from the Stokes I map (red contour at 13 mJy beam⁻¹) can be seen in the PI map, as shown in Figure A2.

Zoom In Zoom Out Reset image size

In order to compare the core-scale B-field (see 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, ⁹³ 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 (${\theta }_{B}^{\mathrm{largescale}}\pm \sigma$) 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 an 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 (see 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 (${\theta }_{B}^{\mathrm{largescale}}$) with a mean orientation of 29° spans spatial scales from ∼0.2 pc (∼5' resolution of Planck) to ∼2.4 pc (1° area around B213).

Table B1. Mean B-field Orientations, Determined from Optical, NIR, and Low-resolution Submillimeter (Planck/850 μm) Polarization Observations

Wavelength	Diameter (arcmin)	No. of Stars/Segments	${\theta }_{B}^{\mathrm{largescale}}\pm \sigma$ (deg)	Offset PA (deg)
(1)	(2)	(3)	(4)	(5)
Optical a	60	15	29 ± 14	104
NIR	60	42	37 ± 17	96
Submillimeter (Planck b )	60	445	29 ± 17	104

Notes. Here ${\theta }_{B}^{\mathrm{largescale}}\pm \sigma$ are the mean and standard deviation values resulting from a Gaussian distribution fitted to the data. Offset PA is the difference in angle between the PA of the B213 filament (∼133°) and the large-scale mean B-field (${\theta }_{B}^{\mathrm{largescale}}$; column (5)).

^a Two measurements with significant deviation in either P or θ are excluded from the optical data. ^b Pixels with values <0.008 K_CMB have been excluded from the Planck data in order to prevent randomization of our inferred B-field direction by measurements dominated by noise.

Download table as:  ASCIITypeset image

To estimate various energy and pressure terms for the cores, we extract their masses and 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⁻¹) of each core. The resulting core dimensions (a = semimajor and b = semiminor), effective radius (${R}_{\mathrm{eff}}=\sqrt{{ab}}$), and PA in degrees east of north are given in Table 2.

The integrated fluxes (F_ν ) and median dust temperatures (T_d; from the Herschel Gould Belt Survey (HGBS) temperature map) over the core are used to estimate core masses using the relation (Hildebrand 1983)

$$\begin{eqnarray}&&M=\displaystyle \frac{{F}_{\nu }\,{D}^{2}}{{B}_{\nu }({T}_{{\rm{d}}}){\kappa }_{\nu }},\end{eqnarray} \tag{ C1 }$$

where D = 140 pc is the distance of B213, κ_ν = 0.0125 cm² g⁻¹ (e.g., Johnstone et al. 2017) is the dust mass opacity, and B_ν (T_d) is the 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 SCUBA-2 (Dempsey et al. 2013), and a 50% uncertainty in dust mass opacity (e.g., Roy et al. 2014).

The column and number densities of the cores are estimated using the following relations:

$$\begin{eqnarray}&&N({{\rm{H}}}_{2})=\displaystyle \frac{M}{\mu \,{m}_{{\rm{H}}}\,\pi \,{R}_{\mathrm{eff}}^{2}}\end{eqnarray} \tag{ C2 }$$

and

$$\begin{eqnarray}&&n({{\rm{H}}}_{2})=\displaystyle \frac{3\,M}{4\,\mu \,{m}_{{\rm{H}}}\,\pi \,{R}_{\mathrm{eff}}^{3}}.\end{eqnarray} \tag{ C3 }$$

Estimated masses and column and number densities and their corresponding uncertainties are given along with T_d and F_ν values in Table 2.

To infer the importance of the 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),

$$\begin{eqnarray}&&\lambda =\frac{(M/\phi )}{(M/\phi {)}_{\mathrm{crit}}}=7.6\,{N}_{\parallel }({{\rm{H}}}_{2})/{B}_{\mathrm{tot}},\end{eqnarray} \tag{ D1 }$$

where N_∥(H₂) is the mean column density (N(H₂)_POL2) in units of 10²¹ cm⁻² 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/$\sqrt{(4{\pi }^{2}G)}$ (Nakano & Nakamura 1978), corresponds to the stability criterion for an isothermal gaseous layer threaded by a perpendicular B-field. A cloud region with (M/ϕ) > (M/ϕ_crit), i.e., λ > 1, will collapse under its own gravity, 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 the B-field. Taking the mean N(H₂)= N_∥(H₂) as (10 ± 6) × 10²¹, (11 ± 6) × 10²¹, and (5 ± 3) × 10²¹ cm⁻² and B = B_tot as 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 (i) the projection effects between N_∥(H₂)/B_tot and the actual measured N(H₂)/B_∥ (where B_∥ is the plane-of-the-sky B-field strength), (ii) the B-field being perpendicular to the core elongation in the case of an oblate spheroid or parallel to the core elongation in the case of a prolate spheroid, and (iii) the assumption that the 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 the mean B-field is parallel to the major axis (Planck Collaboration et al. 2016a; see their Appendix D.4 ⁹⁴ ). No correction was applied on the λ value of Miz-8b because of the misalignment between the 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.

By assuming that the cores are uniform density spheres, we measure the rotational and magnetic energies using the following relations (see Wurster & Lewis 2020):

$$\begin{eqnarray}&&{E}_{\mathrm{rot}}=\displaystyle \frac{{{pMR}}^{2}{{\rm{\Omega }}}^{2}}{5}\end{eqnarray} \tag{ E1 }$$

and

$$\begin{eqnarray}&&{E}_{\mathrm{mag}}=\displaystyle \frac{{B}^{2}V}{8\pi }.\end{eqnarray} \tag{ E2 }$$

In Equation (E1), the correction factor, p = $\tfrac{2(3-{ \mathcal A })}{3(5-{ \mathcal A })}$= 0.27 (where ${ \mathcal A }$ is the power index in the density distribution of the form $\rho \propto {r}^{-{ \mathcal A }}$ and we consider ${ \mathcal 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 = $\sqrt{{ab}}$ (where a = semimajor axis and b = semiminor axis), and the volume of the core V = (4/3)$\pi {R}_{\mathrm{eff}}^{3}$. Here Ω is the angular velocity or magnitude of the velocity gradient of the core measured from the N₂H⁺ data and is found to be 2.05 ± 0.02 km s⁻¹ pc⁻¹ for K04166, 3.86 ± 0.04 km s⁻¹ pc⁻¹ for K04169, and 1.88 ± 0.02 km s⁻¹ pc⁻¹ for Miz-8b (Punanova et al. 2018; see their Table B.2). Here M is the mass of the cores (see Appendix C). The derived energy values and their ratios are given in Table 2.

We model the observed line-of-sight centroid velocities of N₂H⁺ (V_LSR, km s⁻¹) and the corresponding offset length scales in sky coordinates (R.A (Δ_α in pc) and decl. (Δ_δ in pc)) around each pixel in terms of velocity gradients in R.A (${{\rm{\nabla }}}_{{V}_{\alpha }}$, km s⁻¹ pc⁻¹) and decl. (${{\rm{\nabla }}}_{{V}_{\delta }}$, km s⁻¹ pc⁻¹) and the constant systematic velocity of that reference pixel (V₀, km s⁻¹) using a first-degree bivariate polynomial of the form (Goodman et al. 1993; Henshaw et al. 2016; Sokolov et al. 2018)

$$\begin{eqnarray}&&{V}_{\mathrm{LSR}}={V}_{0}+{{\rm{\nabla }}}_{{V}_{\alpha }}{{\rm{\Delta }}}_{\alpha }+{{\rm{\nabla }}}_{{V}_{\delta }}{{\rm{\Delta }}}_{\delta }.\end{eqnarray} \tag{ F1 }$$

We have used the IDL algorithm mpfit to perform weighted, nonlinear, minimum χ² fitting to constrain the velocity gradients ${{\rm{\nabla }}}_{{V}_{\alpha }}$ and ${{\rm{\nabla }}}_{{V}_{\delta }}$ and their corresponding uncertainties. These are further used to derive the magnitude (${ \mathcal G }$) and direction (${{\rm{\Theta }}}_{{ \mathcal G }}$) of the velocity gradients using the following relations:

$$\begin{eqnarray}&&{ \mathcal G }\equiv | {{\rm{\nabla }}}_{V}| =\sqrt{{{\rm{\nabla }}}_{{V}_{\alpha }}^{2}+{{\rm{\nabla }}}_{{V}_{\delta }}^{2}}\end{eqnarray} \tag{ F2 }$$

and

$$\begin{eqnarray}&&{{\rm{\Theta }}}_{{ \mathcal G }}=\arctan \left(\displaystyle \frac{{{\rm{\nabla }}}_{{V}_{\alpha }}}{{{\rm{\nabla }}}_{{V}_{\delta }}}\right).\end{eqnarray} \tag{ F3 }$$

We considered at least six adjacent pixels lying within the beam size of the IRAM 30 m telescope for N₂H⁺ (1–0), 265, around each pixel when the fitting was performed. In addition, we estimated the uncertainties in the velocity gradients 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 velocity gradients superimposed on the POL-2 Stokes I maps of K04166, K04169, and Miz-8b.