JCMT BISTRO Survey: Magnetic Fields within the Hub-filament Structure in IC 5146

Our polarimetric continuum observations toward the IC 5146 dark cloud system were carried out between 2016 May and 2017 April. The observed field targeted the brightest HFS located at the eastern end of the IC 5146 main filament, as shown in Figure 1. We performed 20 sets of 40- minute observations toward the IC 5146 region with τ_{225 GHz} ranging from 0.04 to 0.07.

Zoom In Zoom Out Reset image size

The POL-2 observations were made using POL-2 DAISY scan mode (Friberg et al. 2016; P. Bastien et al. 2019, in preparation), producing a fully sampled circular region of 11 arcmin diameter. Within the DAISY map, the noise is uniform and lowest in the central 3 arcmin diameter region and increases to the edge of the map. The POL-2 data were simultaneously taken at 450 μm with a resolution of 9.6 arcsec and at 850 μm with a resolution of 14.1 arcsec. The 450 μm data are not reported in this paper since the 450 μm instrumental polarization (IP) model has been only recently commissioned.

The IC 5146 data were reduced in a three-stage process using pol2map,⁶⁶ a script recently added to the SCUBA-2 mapmaking routine smurf (Berry et al. 2005; Chapin et al. 2013).

In the first stage, the raw bolometer time streams for each observation are converted into separate Stokes Q, U, and I time streams using the process calcqu.

In the second stage, an initial Stokes I map is created from the I time stream from each observation using the iterative map-making routine makemap. For each reduction, areas of astrophysical emission are defined using a signal-to-noise ratio (S/N) based mask determined iteratively by makemap. Areas outside this masked region are set to zero after each iteration until the final iteration of makemap (see Mairs et al. 2015 for a detailed description of the role of masking in SCUBA-2 data reduction). Convergence is reached when successive iterations of the mapmaker produce pixel values that differ by <5% on average. Each map is compared to the first map in the sequence to determine a set of relative pointing corrections. The individual I maps are then coadded to produce an initial I map of the region.

In the third stage, the final Stokes I, Q, and U maps are created. The initial I map described above is used to generate a fixed S/N-based mask for all further iterations of makemap. The pointing corrections determined in the second stage are applied during the map-making process. In this stage, skyloop, a variant mode of makemap, is invoked. In this mode, rather than each observation being reduced consecutively as is the standard method, one iteration of the mapmaker is performed on each of the observations. At the end of each iteration, all the maps created are coadded. The coadded maps created after successive iterations are compared, and when these coadded maps on average vary by <5% between successive iterations, convergence is reached. Using skyloop typically improves the mapmaker's ability to recover faint extended structure, at the expense of additional memory usage and processing time. The mapmaker was run three times successively to produce the output I, Q, and U maps from their respective time streams. The Q and U data were corrected for IP using the final output I map and the latest IP model (2018 January) (Friberg et al. 2016, 2018).

In all pol2map instances of makemap, the polarized sky background is estimated by doing a principal component analysis of the I, Q, and U time streams to identify components that are common to multiple bolometers. In the first run of makemap (stage 2), the 50 most correlated components are removed at each iteration. In the second run (stage 3), 150 components are removed at each iteration, resulting in smaller changes in the map between iterations and lower noise in the final map.

The output I, Q, and U maps were calibrated in units of Jy beam⁻¹, using a flux conversion factor of 725 Jy pW⁻¹—the standard 850 μm SCUBA-2 flux conversion factor multiplied by 1.35 to account for additional losses due to POL-2 (Dempsey et al. 2013; Friberg et al. 2016).

Finally, a polarization vector catalog was created from the coadded Stokes I, Q, and U maps. To improve the sensitivity, we binned the coadded Stokes I, Q, and U maps into 12'' pixels, and the binned data reached rms noise levels of 1.1 mJy beam⁻¹ for Stokes Q and U.

We calculated the polarization fractions and orientations in the 12'' pixel map. We debiased the former with the asymptotic estimator (Wardle & Kronberg 1974) as

$$\begin{eqnarray}&&P=\displaystyle \frac{1}{I}\sqrt{\left({U}^{2}+{Q}^{2})-\displaystyle \frac{1}{2}(\delta {Q}^{2}+\delta {U}^{2}\right)},\end{eqnarray} \tag{ 1 }$$

where P is the debiased polarization percentage, and I, Q, U, δI, δQ, and δU are the Stokes I, Q, U, and their uncertainties. The uncertainty of polarization fraction was estimated using

$$\begin{eqnarray}&&\delta P=\sqrt{\displaystyle \frac{({Q}^{2}\delta {Q}^{2}+{U}^{2}\delta {U}^{2})}{{I}^{2}({Q}^{2}+{U}^{2})}+\displaystyle \frac{\delta {I}^{2}({Q}^{2}+{U}^{2})}{{I}^{4}}}.\end{eqnarray} \tag{ 2 }$$

The polarization position angle (PA) was calculated as

$$\begin{eqnarray}&&\mathrm{PA}=\displaystyle \frac{1}{2}{\tan }^{-1}\left(\displaystyle \frac{U}{Q}\right),\end{eqnarray} \tag{ 3 }$$

and its corresponding uncertainties were estimated using

$$\begin{eqnarray}&&\delta \mathrm{PA}=\displaystyle \frac{1}{2}\sqrt{\displaystyle \frac{({Q}^{2}\delta {U}^{2}+{U}^{2}\delta {Q}^{2})}{{\left({Q}^{2}+{U}^{2}\right)}^{2}}}.\end{eqnarray} \tag{ 4 }$$

The magnetic field orientations used in this paper were assumed to be PA + 90°.

The SCUBA-2 850 μm waveband covers the wavelength of the CO (J = 3−2) rotational line, and thus our measured continuum flux could be affected by CO line emission (e.g., Drabek et al. 2012; Coudé et al. 2016). Furthermore, the CO (J = 3−2) rotational line is known to be polarized via the Goldreich–Kylafis effect (Goldreich & Kylafis 1981, 1982). For example, the typical polarization fraction of CO (J = 3−2) could be ≲3% in dense clouds and outflows (Ching et al. 2016), calculated using the formulation in Deguchi & Watson (1984) and Cortes et al. (2005). The polarization angle of CO line is either parallel or perpendicular to the magnetic fields depending on optical depth and the relative angle between magnetic fields and gas velocity fields (Cortes et al. 2005). If a typical polarization fraction of 2% for CO (J = 3−2) is assumed, the polarized intensity from the line would be only 0.02%–0.14% of the total 850 μm flux, which is insignificant compared with the uncertainties of polarization, ≳0.2%–0.5%, in the central hub.

The CO contamination in total intensity might also decrease the observed polarization fraction. Johnstone et al. (2017) calculated the fraction of CO (J = 3−2) line emission to the total flux in the JCMT 850 μm waveband toward several clumps in the IC 5146 system. The fraction of CO (J = 3−2) to total flux in our target region is mostly ∼1%–3%, but a higher fraction of ∼7% was found in the central hub. Hence, the CO contamination would reduce the measured polarization fraction by a factor of 1%–7%. Nevertheless, this effect is insignificant to our analysis because the S/Ns of our polarization detections are typically only ∼2–4.

We show the observed magnetic field orientations traced by POL-2 850 μm polarization, with pixel size of 12'', overlaid on the Stokes I map, with pixel size of 4'', in Figure 2. We selected the 139 vectors with I/δI > 10 to ensure that the selected data are associated with the target core-scale HFS. Montier et al. (2015a) suggest that the uncertainty in Stokes I may enhance the bias in polarization fraction for data with I/δI < 10, and thus our I/δI > 10 selection criterion could exclude these biased data. Among these I/δI > 10 data, 30 of them have 2 < P/δP < 3 and 42 of them have P/δP > 3. In order to better probe the magnetic fields, we further added the P/δP > 2 criterion to exclude the samples with higher uncertainties in PA, and the final selected samples have a maximum δPA of 127 and a mean δPA of 85. Figure 3 shows the zoom-in polarization map toward the HFS and our final selected samples. We note that the CO contamination in Stokes I only has insignificant effect on our sample selection. Assuming the CO contamination in total intensity is 7% everywhere, as the worst case, the number of P/δP > 2 vectors would only decrease to 68 from 72.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The Stokes I map shows a central massive clump in which three filaments intersect. The observed morphology is consistent with the typical HFS. The central massive clump hosts ∼80% of the total intensity within the system, and thus we recognize the central massive clump as the hub of the HFS. Three filaments are identified extending from the central hub to the north, east, and south. The magnetic field revealed by our polarization map seems to have small angular dispersion but also shows a change of orientations from the north to the south.

To compare the magnetic fields in the observed HFS with the large-scale magnetic fields shown in Figure 1, we plot a histogram of the PAs of the local magnetic fields from POL-2 and WLE17 data within our field of view (diameter of 11') in Figure 4 with a bin size of 10° that is close to our mean δPA of 85.

Zoom In Zoom Out Reset image size

The PA histogram of our data shows two major components separated by a dip at 10°. The PA > 10° component has a mean PA of 37° and a PA dispersion of 15°, which is similar to the large-scale magnetic fields (28° ± 21°). In contrast, the PA < 10° component has a mean PA of −27° and a PA dispersion of 27°. The PA difference of 64° between the two components is much greater than the PA dispersion for large-scale magnetic fields as well as for our mean observational uncertainty (85), suggesting that the observed magnetic field morphology is significantly different from the large-scale magnetic fields.

Figure 5 shows the locations of these two components. To represent the major axis of the main filament, we plotted the yellow dashed line that shows the direction across the intensity peaks of the two clumps along the parsec-scale filament. This major axis has an orientation of −73° and roughly separates the spatial distribution of the two magnetic field orientation components. Within the HFS, the red and green components tend to be distributed in the northern and southern half of the HFS; the tendency, however, is reversed in the western clump. In addition, the orientation perpendicular to the main filament (17°) is also close to the dip between the two components. These features favor the possibility that the magnetic fields in the HFS are curved along the main filament. In contrast, the WLE17 data only show a major peak (28°) in the PA histogram, which is roughly perpendicular to the large-scale filament but with a ≈10° offset in PA. A minor component peaked at ≈−75° is also shown by our data; however, this component is diffusely distributed over the area, and thus more vectors are needed to reveal these structures.

Zoom In Zoom Out Reset image size

To investigate whether our polarization data trace the dust grains in high-extinction regions, we plot the 850 μm emission polarization fraction P_emit versus A_V in Figure 6. To reveal the complete P_emit–A_V distribution, this figure includes all the data with I/δI > 10, and the data points are color coded based on their S/N of P_emit.

Zoom In Zoom Out Reset image size

To estimate the A_V, we calculated the τ_{850 μm} from the observed 850 μm intensity using I_{850 μm} = τ_{850 μm}B(T_dust), assuming that the dust emission is optically thin at 850 μm. We used the dust temperature B(T_dust) derived in Arzoumanian et al. (2011) via fitting the Herschel data at five wavelengths with a modified blackbody function assuming a dust emissivity index of 2. The dust temperature map and the Stokes I map were both resampled on a 12'' grid to match our polarization catalog. The τ_{850 μm} was converted to A_V using the R_V = 3.1 extinction curve in Weingartner & Draine (2001). We note that the extinction curve may vary at dense regions due to grain growth. If R_V changes from 3.1 to 5.5 within the observed regions, we would underestimate the P_emit versus A_V slope by 10%.

The P_emit are equivalent to the extinction polarization percentages divided by the optical depth (P_ext/τ_λ) in the optically thin case (Andersson et al. 2015) and so are proportional to the PE (defined as P_ext/A_V). Thus, the observed P_emit versus A_V slope is equivalent to the PE versus A_V slope.

We further plotted in Figure 6 the PE versus A_V revealed by WLE17 optical and infrared polarization data to show how PE varies in low A_V regions. The PE at 850 μm are in the form of P_ext,850/τ₈₅₀, and the PE obtained in the H-band are represented by P_ext,H/A_V. Thus, a scaling factor $\tfrac{{P}_{\mathrm{ext},850}}{{P}_{\mathrm{ext},H}}\cdot \tfrac{{A}_{V}}{{\tau }_{850}}$ is required to convert the PE at the two wavelengths, which is determined by the unknown dust properties (Andersson et al. 2015). Via matching the fitting results of PE versus A_V relation in the H-band (WLE17) and in the 850 μm bands (described in Section 3.3) at A_V = 20 mag, we found a scaling factor of 48.3, which we used to match the two data sets. This scaling factor is not a universal constant, as discussed by Jones et al. (2015), and varies with physical conditions in different clouds.

To determine the P_emit versus A_V slope, the conventional approach is to apply a least-squares power law fit to data selected by an S/N cut in P_emit. Following this approach, we fitted the P/δP > 2 and >3 data with a power-law function. The best-fit functions are shown in Figure 6 by blue and green dashed lines, and the best-fit power-law indices are −1.02 ± 0.02 and −1.08 ± 0.02, respectively. Nevertheless, because Figure 6 shows that the P_emit–A_V distribution is significantly truncated by the S/N cut, and also because the best-fit trends are very similar to the truncated boundary, it raises doubts about whether or not the fitting is biased by the sample selection.

We investigated how the sample selection biased on P/δP affects the fitting to the P_emit versus A_V distribution by performing Monte Carlo simulations of data sets with underlying P_emit versus A_V function and randomly generated measurement errors in the Appendix. We found that the fitted power index would be dominated by the S/N cut and approaches −1 rather than the true underlying value, if the P_emit versus A_V distribution is significantly truncated by the applied P/δP selection criteria. Hence, to obtain an unbiased power index, it is recommended to include the low P/δP data, so that a complete probability density function (PDF) of P_emit can be recovered. Nevertheless, the use of low P/δP data would break the Gaussian PDF assumption (Wardle & Kronberg 1974; Vaillancourt 2006) required for least-squares fit, and therefore we turn to using a Bayesian approach to fit the observed P_emit versus A_V distribution.

We used a Bayesian approach to apply a non-Gaussian PDF and fit the observed P_emit versus A_V trend with a power-law model. The Bayesian statistical framework provides a model fitting tool based on probability theory (see detailed introduction in Jaynes & Bretthorst 2003). The general form of the Bayesian inference is

$$\begin{eqnarray}&&P(\theta | D)=\displaystyle \frac{P(\theta )P(D| \theta )}{P(D)}\end{eqnarray} \tag{ 5 }$$

or

$$\begin{eqnarray}&&\mathrm{Posterior}=\displaystyle \frac{\mathrm{Prior}\times \mathrm{Likelihood}}{\mathrm{Evidence}},\end{eqnarray} \tag{ 6 }$$

where D is the observed data and θ represents the model parameters. The posterior $P(\theta | D)$ describes the probability of the model parameters matching the given data, which is what we are interested in. The evidence P(D) is the probability of obtaining the data, which mainly serves as a normalization factor for the posterior. The prior P(θ) serves as the initial guessed probability of the model parameters based on our prior knowledge. The likelihood $P(D| \theta )$ describes how likely it is for a given model parameter set to match the observed data.

Assuming the measurements in Stokes Q and U have similar and Gaussian distributed noise, the PDF of the observed polarization fraction has been known to follow the Rice distribution (Rice 1945; Wardle & Kronberg 1974; Simmons & Stewart 1985; Quinn 2012)

$$\begin{eqnarray}&&F(P| {P}_{0})=\displaystyle \frac{P}{{\sigma }_{P}^{2}}\exp \left[-\displaystyle \frac{{P}^{2}+{P}_{0}^{2}}{2{{\sigma }_{P}}^{2}}\right]{I}_{0}\left(\displaystyle \frac{{{PP}}_{0}}{{\sigma }_{P}^{2}}\right),\end{eqnarray} \tag{ 7 }$$

where P is the observed polarization fraction, P₀ is the real polarization fraction, σ_P is the uncertainty in polarization fraction, and I₀ is the zeroth-order modified Bessel function. The likelihood function of polarization measurements is defined as

$$\begin{eqnarray}&&L({P}_{0})=\displaystyle \prod _{i=1}^{n}F({P}_{n}| {P}_{0}),\end{eqnarray} \tag{ 8 }$$

where P_n represents the nth measurement.

To perform the fit to the P_emit versus A_V trend using a Bayesian approach, we assumed the following power-law model such that

$$\begin{eqnarray}&&{P}_{0}=\beta {A}_{V}^{-\alpha },\end{eqnarray} \tag{ 9 }$$

where α and β are the free model parameters and A_V is the observed visual extinction. The uncertainty in the polarization fraction is

$$\begin{eqnarray}&&{\sigma }_{P}={\sigma }_{{QU}}/I.\end{eqnarray} \tag{ 10 }$$

Here, the I is the observed total intensity, and the σ_QU is a free model parameter describing the dispersion in Stokes Q and U, which has contributions from both the instrumental uncertainty and the intrinsic dispersion due to source properties. We then simply used uniform priors within reasonable limits as

$$\begin{eqnarray}\begin{array}{rcl}P(\alpha ) & = & \left\{\begin{array}{cc}\mathrm{Uniform} & 0\lt \alpha \lt 2\\ 0 & \mathrm{otherwise}\end{array}\right.\\ P(\beta ) & = & \left\{\begin{array}{cc}\mathrm{Uniform} & 0\lt \beta \lt 400\\ 0 & \mathrm{otherwise}\end{array}\right.\\ P({\sigma }_{{QU}}) & = & \left\{\begin{array}{cc}\mathrm{Uniform} & 0\lt {\sigma }_{{QU}}\lt 10\\ 0 & \mathrm{otherwise}.\end{array}\right.\end{array}\end{eqnarray} \tag{ 11 }$$

The Bayesian model fitting was performed with the Python Package PyMC3 (Salvatier et al. 2016) via a Markov Chain Monte Carlo method using the Metropolis–Hastings sampling algorithm. The 12 arcsec pixel data were used for the fitting to ensure that each measurement is independent.

The derived posterior of each model parameter is shown in Figure 7, and the 95% highest posterior density (HPD) interval of each parameter is plotted to represent the uncertainty. The 95%, 68%, and 50% confidence regions (CRs) predicted by the posterior distribution in P_emit versus A_V space are shown in Figure 8, assuming a dust temperature of 13 K. Since the error distribution of P is asymmetric and also varies with P/σ_P and A_V, the predicted P versus A_V is not simply a straight line on a logarithmic scale, even though the input model is a power law. Almost all the data points are well within the 95% CRs predicted by the posterior.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The derived α has a mean value of 0.56 with a 95% confidence interval from 0.22 to 0.83. The α derived by the Bayesian method is less steep than the α ≈ 1.0 derived from the conventional approach, confirming that the conventional method is biased (see Figure 6). The α range of 0.22–0.83 includes the index of 0.25 ± 0.06 obtained from near-infrared polarization data in A_V of 3–20 mag regions (see Wang et al. 2017), and thus no significant difference in PE was found between the A_V < 20 mag and A_V = 20–300 mag regions (Figure 6). In addition, the fitted σ_QU of 1.78 mJy beam⁻¹ is greater than our estimated instrumental noise of 1.1 mJy beam⁻¹, which indicates a significant intrinsic dispersion in PE.

The value of α smaller than unity indicates that the extinction polarization fraction (P_ext = τP_emit) increases with column density. Since the extinction polarization fraction is defined as tanh(Δτ) (Jones 1989), where Δτ is the differential optical depth between two polarization directions, the increase of extinction polarization fraction indicates a higher amount of aligned dust grains. Hence, our results suggest that the dust grains in the IC 5146 dense regions can still be aligned with magnetic fields. The α of ∼0.5 is also predicted by the simulations based on RATs theory (e.g., Whittet et al. 2008) in low- density regions, where the radiation field is sufficiently strong to align dust grains.

Three possibilities could explain why the dust grains within these dense regions can still be aligned. First, because our target is an active star-forming region, the embedded young stars could be the sources of radiation needed to align the dust grains in dense regions. Second, WLE17 found that the dust grains in IC 5146 have significantly grown from the diffuse interstellar medium. These large dust grains could be aligned more efficiently by the radiation with longer wavelengths, which can penetrate the dense regions (Lazarian & Hoang 2007a; Hoang & Lazarian 2009). Third, the mechanical torques due to infalling gas and outflows in the star-forming regions could possibly align the dust grains in the absence of a radiation field (Lazarian & Hoang 2007b; Hoang et al. 2018). These possibilities will be further investigated in upcoming BISTRO papers probing the PE in different environments.

To investigate whether magnetic fields influence the clump fragmentation within the IC 5146 HFS, we examined the correlation between the magnetic field orientations and the clump morphologies. Based on the JCMT 850 μm Gould Belt Survey data, Johnstone et al. (2017) identified eight submillimeter clumps within the regions where we had polarization detections (Figure 9). To represent the orientation of each clump, we used our total intensity map and performed a 2D Gaussian fit to each clump to find the PA of its FWHM major axis. The 2D Gaussian fit had typical orientation uncertainties of ∼10°. The obtained clump orientations are listed in Table 1 and plotted over the polarization map in Figure 9.

Zoom In Zoom Out Reset image size

Table 1.  Geometric and Polarization Properties of the Clumps

IDa	Total Massa	Major Axis	Minor Axis	Reff	Clump Orientation	B-field PApeakb	σNT	${\alpha }_{\mathrm{vir}}=\tfrac{{{\rm{M}}}_{\mathrm{vir}}}{{M}_{\mathrm{clump}}}$ c	αvir,B
	(M⊙)	(arcsec)	(arcsec)	(arcsec)	(deg)	(deg)	(km s−1)
42	11 ± 3	51 ± 2	33 ± 2	43 ± 2	14 ± 10	134.8 ± 28.5	0.25 ± 0.01	0.9 ± 0.2	⋯
43	2.0 ± 0.5	29 ± 1	16 ± 1	22 ± 1	4 ± 10	100.1 ± 6.6	0.12 ± 0.01	1.4 ± 0.3	⋯
45	0.77 ± 0.18	37 ± 3	16 ± 1	24 ± 2	174 ± 10	0.7 ± 9.5	0.18 ± 0.02	5.1 ± 1.0	⋯
46	6.4 ± 1.6	26 ± 2	23 ± 1	24 ± 2	25 ± 10	51.5 ± 19.7	0.21 ± 0.03	0.7 ± 0.1	⋯
47	85 ± 20	40 ± 2	36 ± 2	38 ± 2	86 ± 10	21.0 ± 2.7	0.36 ± 0.12	0.2 ± 0.1	0.2–1.0d
48	6.0 ± 1.5	35 ± 3	25 ± 1	30 ± 3	164 ± 10	−4.1 ± 1.6	0.14 ± 0.01	0.7 ± 0.1	⋯
50	0.97 ± 0.24	29 ± 1	18 ± 1	23 ± 1	9 ± 10	13.6 ± 11.4	⋯	⋯	⋯
52	7.6 ± 1.9	31 ± 2	17 ± 1	23 ± 2	135 ± 10	64.7 ± 40.7	0.29 ± 0.01	0.9 ± 0.2	⋯
53	1.7 ± 0.4	34 ± 3	18 ± 1	24 ± 3	140 ± 10	−19.7 ± 10.1	⋯	⋯	⋯

Notes.

^aThe clumps ID and mass were obtained from Johnstone et al. (2017), but the masses were scaled to a distance of 813 ± 106 pc. ^bMean magnetic field orientation averaged using the 3 × 3 pixels at the intensity peaks. ^cThe virial masses of clumps were calculated considering the support from thermal pressure and turbulence. ^dIf the inclination angle of the magnetic field with respect to the line of sight is >15°.

Download table as:  ASCIITypeset image

To estimate the local magnetic field orientation within each clump, we calculated the mean magnetic field orientations by averaging the data within 3 × 3 pixels at the clump intensity peaks. The size of 3 × 3 pixels is comparable to the typical radius of these clumps (≈20''–40'', see Table 1), so the average represents the mean magnetic field orientation over the dense center of these clumps. In addition, The 3 × 3 pixels also provide an estimation of orientation dispersions, which were used as the uncertainties of the averaged orientations; if only one vector was obtained for a given clump, the instrumental uncertainty would be used. To handle the ±180° ambiguity, the mean and dispersion of the magnetic field PAs were calculated in a new coordinate system, where the PA dispersion was minimized, and the calculated results were converted back to the standard coordinate system.

The derived local magnetic field orientations versus the clump orientations are plotted in Figure 10. The comparison between the clump axis and the magnetic field orientations in the clump is limited by the small number of statistics. Nevertheless, there appears to be a tendency that the observed clumps are likely either parallel or perpendicular to the mean magnetic field orientation within ∼20°, as shown in Figure 10. The upcoming BISTRO data would provide a much bigger sample set from various systems to statistically confirm this tendency. If the tendency is real, it would suggest that the magnetic fields are a key factor in the fragmentation of these clumps. We note that clumps 43 and 52 only contain two polarization vectors, which are almost perpendicular to each other, and thus the mean magnetic field orientations are not meaningful for these two cases. Because the orientation of these two clumps is still parallel to one of the vectors and perpendicular to the other, these two clumps are still consistent with the tendency.

Zoom In Zoom Out Reset image size

We further plot in Figure 10 the mean large-scale magnetic field orientation, 28°, obtained from WLE17. Only clump 47 has a magnetic field orientation similar to the large-scale magnetic fields within 20°. All other clumps are aligned either parallel or perpendicular with the local magnetic field and show no significant correlation with the large-scale magnetic fields. Hence, these clumps are more likely formed after the local magnetic fields were distorted by the process of the clump formation.

The Davis–Chandrasekhar–Fermi (DCF) method (Davis 1951; Chandrasekhar & Fermi 1953) is commonly used to evaluate the magnetic field strength from dust polarizations. The DCF method assumes that turbulent kinematic energy and the turbulent magnetic energy are in equipartition, and hence the magnetic field strength can be estimated using

$$\begin{eqnarray}&&{B}_{\mathrm{pos}}=Q\,\sqrt{4\pi \rho }\displaystyle \frac{{\sigma }_{v}}{\delta {\phi }_{\mathrm{intrinsic}}},\end{eqnarray} \tag{ 12 }$$

where δϕ is the intrinsic angular dispersion of the magnetic fields, σ_v is the velocity dispersion, ρ is the gas density, and Q is a factor accounting for the complicated magnetic field and inhomogeneous density structure. Ostriker et al. (2001) suggested that Q = 0.5 yields a good estimation of magnetic field strength on the plane of sky if magnetic field angular dispersion is <25°.

3.5.1. Magnetic Field Angular Dispersion

We used the 12'' pixel polarization data to calculate the magnetic field angular dispersion to ensure that all vectors we used are independent measurements. To avoid small number statistics (fewer than 10 vectors), we only perform the angular dispersion estimation using the polarization vectors (45 vectors) within the central hub (clump 47) (Figure 9).

The DCF method requires an estimation of magnetic field distortion caused by turbulence, and the underlying magnetic field geometry might bias the estimation. Thus, we calculated the magnetic field angular dispersion in a local area to avoid the angular dispersion contribution from the large-scale nonuniform magnetic field geometry. Specifically, we selected each 24'' × 24'' box (i.e., the width of two independent beams) and calculated the mean and the corresponding sum of squared differences (SSD = ${\sum }_{i=1}^{n}{(\bar{\mathrm{PA}}-{\mathrm{PA}}_{i})}^{2}$) of the PA using the, at most, nine vectors within each box. The SSD from all boxes were averaged with inverse-variance weighting, and the square root of the mean SSD was taken as the observed angular dispersion. Next, the mean instrumental uncertainties (δϕ_ins) of 80 were removed from the observed angular dispersion (δϕ_obs) to obtain the intrinsic dispersion (δϕ_intrinsic) via

$$\begin{eqnarray}&&\delta {\phi }_{\mathrm{intrinsic}}^{2}=\delta {\phi }_{\mathrm{obs}}^{2}-\delta {\phi }_{\mathrm{ins}}^{2}.\end{eqnarray} \tag{ 13 }$$

With these corrections, the calculated δϕ_intrinsic for clump 47 is 174 ± 06.

3.5.2. Velocity Dispersion

To estimate the velocity dispersion, we used the C¹⁸O (J = 3−2) spectral data taken by Graham (2008) with the JCMT HARP receiver (Buckle et al. 2009). CO and its isotopomers are well mixed with H₂ and are commonly used to trace the gas kinematics. The C¹⁸O (J = 3−2) line, in particular, is expected to trace gas with volume density up to ∼10⁵ cm⁻³ (e.g., Di Francesco et al. 2007), which is comparable to the densities of our target field. In addition, the C¹⁸O (J = 3−2) line is likely optically thin in this field (Graham 2008), so it traces the kinematics of all the gas in the clump. Therefore, we assumed that the C¹⁸O (J = 3−2) line width can well represent the gas velocity dispersion in our observing regions.

The C¹⁸O data reveal, at least, three velocity components within the central hub, which peaked at 3.7, 4.1, and 4.5 km s⁻¹. Because the three components have very similar velocities, and also multiple YSOs have been identified in the central hub by Harvey et al. (2008), the multiple components are possibly structures within the hub, instead of foreground or background components. We performed a multicomponent Gaussian fit to estimate the C¹⁸O (J = 3−2) line width, using the python package PySpecKit (Ginsburg & Mirocha 2011). We only accepted the fitted Gaussian components with amplitudes >5σ. To estimate the thermal velocity dispersion, we adopted a gas kinematic temperature (T_kin) of 10 ± 1 K, which is the same as the excitation temperature estimated from the ¹³CO (J = 3−2) line in Graham (2008), leading to $\sqrt{\tfrac{{k}_{B}{T}_{\mathrm{kin}}}{{m}_{{{\rm{C}}}^{18}{\rm{O}}}}}=0.05\pm 0.01$ km s⁻¹. The thermal velocity dispersions were then removed from the fitted line widths to obtain the nonthermal velocity dispersions via

$$\begin{eqnarray}&&{\sigma }_{\mathrm{NT}}^{2}={\sigma }_{\mathrm{obs}}^{2}-\displaystyle \frac{{k}_{B}{T}_{\mathrm{kin}}}{{m}_{{{\rm{C}}}^{18}{\rm{O}}}},\end{eqnarray} \tag{ 14 }$$

where σ_NT is the nonthermal velocity dispersion, σ_obs is the observed C¹⁸O Gaussian line width, and ${m}_{{{\rm{C}}}^{18}{\rm{O}}}$ is the molecular weight. The inverse-variance weighted mean of the σ_obs and σ_NT of all velocity components over the central hub was 0.37 and 0.36 km s⁻¹, respectively, and the dispersion of σ_obs of 0.12 km s⁻¹ among the central hub was used as its uncertainty.

3.5.3. Volume Density

3.5.4. Magnetic Field Strength and Mass-to-flux Ratio

Using Equation (12) and the quantities estimated above (Table 2), the plane-of-sky magnetic field strength (B_pos) is estimated to be 0.5 ± 0.2 mG. To evaluate the relative importance of magnetic fields and gravity in the central hub, we calculate the mass-to-flux critical ratio via

$$\begin{eqnarray}&&{\lambda }_{\mathrm{obs}}=\displaystyle \frac{{\left(M/{\rm{\Phi }}\right)}_{\mathrm{obs}}}{{\left(M/{\rm{\Phi }}\right)}_{\mathrm{cri}}},\end{eqnarray} \tag{ 15 }$$

where the observed mass-to-flux ratio is

$$\begin{eqnarray}&&{\left(M/{\rm{\Phi }}\right)}_{\mathrm{obs}}=\displaystyle \frac{\mu {m}_{{\rm{H}}}{N}_{{{\rm{H}}}_{2}}}{{B}_{\mathrm{pos}}},\end{eqnarray} \tag{ 16 }$$

where μ = 2.8 is the mean molecular weight per H₂ molecule and the (M/Φ)_cri is the critical mass-to-flux ratio defined as

$$\begin{eqnarray}&&{\left(M/{\rm{\Phi }}\right)}_{\mathrm{cri}}=\displaystyle \frac{1}{2\pi \sqrt{G}}\end{eqnarray} \tag{ 17 }$$

(Nakano & Nakamura 1978). Due to the unknown inclination of the clumps, the observed mass-to-flux ratio λ_obs is also only an upper limit. Crutcher (2004) suggested that a statistically average factor of one-third could be used to estimate the real mass-to-flux ratio accounting for the random inclinations for an oblate spheroid core, flattened perpendicular to the orientation of the magnetic field. Since we have shown that the central clump is elongated with its major axis perpendicular to the local magnetic field, we adopt a factor of one-third to estimate the mass-to-flux ratio via

$$\begin{eqnarray}&&\lambda =\displaystyle \frac{1}{3}{\lambda }_{\mathrm{obs}}.\end{eqnarray} \tag{ 18 }$$

Table 2.  Derived Magnetic Field Strength of Clump 47

ID	σNT	δϕ	${n}_{{{\rm{H}}}_{2}}$	Bpos	λ
	(km s−1)	(deg)	(cm−3)	(mG)
47	0.36 ± 0.12	17.4 ± 0.6	(9.8 ± 2.4) × 105	0.5 ± 0.2	1.3 ± 0.4

Note. The magnetic field strength estimated using the DCF method. The variables σ_NT, δϕ, ${n}_{{{\rm{H}}}_{2}}$, B_pos, and λ represent the velocity dispersion, magnetic field angular dispersion, H₂ volume density, plane of sky magnetic field strength, and mass-to-flux ratio, respectively.

Download table as:  ASCIITypeset image

The estimated mass-to-flux ratio for the central clump is 1.3 ± 0.4. The DCF method often tends to overestimate the magnetic field strength, because the effect of integration over the telescope beam and along the line of sight might smooth out part of the magnetic field structure, resulting in an underestimated angular dispersion (Heitsch et al. 2001; Ostriker et al. 2001; Crutcher 2012). In addition, our target region has a complicated velocity structure, and therefore the measured velocity dispersion might have contributions from gas accretion or contraction motions instead of only isotropic turbulence, also leading to an overestimated magnetic strength. Hence, our estimation of the mass-to-flux ratio only represents a lower limit. A mass-to-flux ratio ≳1.0 suggests that the central hub is supercritical and that magnetic fields and gravity are comparably important at subparcsec scales.

3.5.5. Angular Dispersion Function

Hildebrand et al. (2009) developed an alternative method to improve the DCF method using a polarization angular dispersion function (ADF) to accurately extract the turbulent component from the polarization data. Houde et al. (2009) further generalized the ADF method by including the effect of signal integration along the thickness of the clouds and over the telescope beam. In this section, we test whether the magnetic field strength estimated using the ADF method is significantly different from our estimation in Section 3.5.4.

The ADF method assumes that the magnetic fields in clouds are combinations of ordered large-scale component B₀ and turbulent component B_t, and the ratio of these two components determines the intrinsic polarization angular dispersion that

$$\begin{eqnarray}&&\delta {\phi }_{\mathrm{intrinsic}}={\left[\displaystyle \frac{\langle {B}_{t}^{2}\rangle }{\langle {B}_{0}^{2}\rangle }\right]}^{\tfrac{1}{2}},\end{eqnarray} \tag{ 19 }$$

where $\langle ...\rangle$ denotes an average. Hence, the DCF equation (Equation (12)) can be written as

$$\begin{eqnarray}&&{B}_{\mathrm{pos}}=\sqrt{4\pi \rho }\,{\sigma }_{v}{\left[\displaystyle \frac{\langle {B}_{t}^{2}\rangle }{\langle {B}_{0}^{2}\rangle }\right]}^{-\tfrac{1}{2}}.\end{eqnarray} \tag{ 20 }$$

The detailed derivation given by Hildebrand et al. (2009) and Houde et al. (2009) shows that the ratio of turbulent to magnetic energy can be estimated from the ADF using the following equation:

$$\begin{eqnarray}&&1-\langle \cos [{\rm{\Delta }}{\rm{\Phi }}({\ell })]\rangle \simeq \displaystyle \frac{1}{N}\displaystyle \frac{\langle {B}_{t}^{2}\rangle }{\langle {B}_{0}^{2}\rangle }(1-{e}^{-{{\ell }}^{2}/2({\delta }^{2}+2{W}^{2})})+a{{\ell }}^{2},\end{eqnarray} \tag{ 21 }$$

where ΔΦ(ℓ) is the difference in the polarization angle measured at two positions separated by a distance ℓ. The quantities δ and a are unknown parameters, representing the turbulent correlation length and the first-order Taylor expansion of the large-scale magnetic field structure. The quantity W is the telescope beam radius, which is 6.2 arcsec at 850 μm. The quantity N is the number of turbulent cells observed along the line of sight and within the telescope beam and can be estimated from

$$\begin{eqnarray}&&N={\rm{\Delta }}^{\prime} \displaystyle \frac{{\delta }^{2}+2{W}^{2}}{\sqrt{2\pi }{\delta }^{3}},\end{eqnarray} \tag{ 22 }$$

where Δ' is the effective cloud thickness, which is assumed to be the clump effective radius. Via fitting the above equation to the observed $1-\langle \cos [{\rm{\Delta }}{\rm{\Phi }}({\ell })]\rangle$ versus ℓ distribution, the three unknown parameters δ, $\tfrac{\langle {B}_{t}^{2}\rangle }{\langle {B}_{0}^{2}\rangle }$, and a can be derived.

We applied the ADF method to estimate the magnetic field strength in clump 47 using the same selected polarization vectors as in Section 3.5.1. We calculated the $\cos [{\rm{\Delta }}{\rm{\Phi }}({\ell })]$ and ℓ from each pair of the polarization vectors within clump 47, and the results are averaged in bins of width ℓ = 12 arcsec to estimate the ADF $1-\langle \cos [{\rm{\Delta }}{\rm{\Phi }}({\ell })]\rangle$ versus ℓ. The calculated ADF is plotted in Figure 11. We fitted the observed ADF using Equations (21) and (22), and the best-fitting parameters are shown in Table 3. The obtained turbulent to magnetic energy ratio $\tfrac{\langle {B}_{t}^{2}\rangle }{\langle {B}_{0}^{2}\rangle }$ is 0.27 ± 0.03, suggesting that the turbulent magnetic field component is weaker than the ordered large-scale magnetic field. With the previously derived gas velocity and the volume density (Table 2), the magnetic field strength is estimated to be 0.4 ± 0.2 mG, which is consistent with our estimation using the DCF method (0.5 ± 0.2 mG) within the uncertainties.

Zoom In Zoom Out Reset image size

Table 3.  Magnetic Field Strength of Clump 47 Using the ADF Method

Fit Result			Derived Quantities
δ	$\tfrac{\langle {B}_{t}^{2}\rangle }{\langle {B}_{0}^{2}\rangle }$	a	N	Bpos
(arcsec)		(arcsec−2)		(mG)
15.1 ± 3.9	0.27 ± 0.03	(1.5 ± 0.3) × 10−5	1.73 ± 0.4	0.4 ± 0.2

Download table as:  ASCIITypeset image

3.5.6. Alfvénic Mach Number

The turbulent Alfvénic Mach number (M_A) describes the relative importance of magnetic fields and turbulence, and hence it is a key factor in most cloud evolution models (e.g., Padoan et al. 2001; Nakamura & Li 2008). In the sub-Alfvénic case (M_A ≤ 1), magnetic fields are strong enough to regulate turbulence and cause an organized magnetic field and cloud structure. In the super-Alfvénic case (M_A > 1), the turbulence can efficiently change the magnetic field morphology, and the magnetic field morphology is expected to be random.

The Alfvénic Mach number can be estimated from the angular dispersion of the magnetic field if the same assumptions as for the DCF method are made. In doing so, the definition of the Alfvénic Mach number

$$\begin{eqnarray}&&{M}_{{\rm{A}}}=\displaystyle \frac{{\sigma }_{\mathrm{NT}}}{{V}_{{\rm{A}}}}=\displaystyle \frac{{\sigma }_{\mathrm{NT}}\sqrt{4\pi \rho }}{B}\end{eqnarray} \tag{ 23 }$$

can be combined with the equation of the DCF method (Equation (12)), yielding

$$\begin{eqnarray}&&{M}_{{\rm{A}}}=\displaystyle \frac{\delta \phi \cdot \sin \theta }{Q},\end{eqnarray} \tag{ 24 }$$

where θ is the inclination of the magnetic fields, with respect to the line of sight, so that B_pos = $B\sin \theta$. For Q = 0.5, the obtained magnetic field angular dispersion 17° corresponds to M_A of $0.6\sin \theta$, and hence the central hub is likely sub-Alfvénic.

In this section, we use the virial analysis to investigate whether thermal pressure, turbulence, and magnetic fields are sufficient to support clumps against gravity. If a clump with uniform density is supported by only thermal pressure and turbulence, the virial mass (M_vir) is

$$\begin{eqnarray}&&{M}_{\mathrm{vir}}=\displaystyle \frac{5{R}_{\mathrm{eff}}}{G}({\sigma }_{\mathrm{NT}}^{2}+{c}_{s}^{2})\end{eqnarray} \tag{ 25 }$$

(Bertoldi & McKee 1992; Pillai et al. 2011; Liu et al. 2018), where R_eff is the geometric mean of major and minor radius and c_s = 0.19 km s⁻¹ is the sound speed for a kinematic temperature of 10 K and mean molecular weight. Virial mass is the maximum mass of a stable clump with the support from kinematic and thermal energy. Hence, a clump mass greater than the virial mass, or a virial parameter (α_vir = M_vir/M_clump) less than unity, indicates that the clump is gravitationally unstable. We calculate the virial parameter of each clump and list the results in Table 1. Except for clumps 43 and 45, most of the clumps have α_vir less than unity, suggesting that thermal pressure and turbulence are insufficient to support the clump against gravity. Hence, these clumps require additional support from magnetic fields to stop gravitational collapse.

If support from magnetic fields is taken into account, the virial mass of a clump becomes

$$\begin{eqnarray}&&{M}_{\mathrm{vir}}^{B}=\displaystyle \frac{5{R}_{\mathrm{eff}}}{G}\left({\sigma }_{\mathrm{NT}}^{2}+{c}_{s}^{2}+\displaystyle \frac{{V}_{{\rm{A}}}^{2}}{6}\right)\end{eqnarray} \tag{ 26 }$$

(Bertoldi & McKee 1992; Pillai et al. 2011; Liu et al. 2018), where the additional term $\tfrac{{V}_{{\rm{A}}}^{2}}{6}$ stands for the support from magnetic field pressure. We have estimated an Alfvénic Mach number of $0.6\sin \theta$ for clump 47 in Section 3.5, which corresponds to an Alfvénic velocity of $0.64/\sin \theta$. With support from magnetic fields, the α_vir of clump 47 becomes 0.2–1.0 for a θ of 15°–90° and greater than unity if θ < 15°. Hence, if the direction of magnetic field is not very close to the line of sight, clump 47 is likely gravitationally unstable, which is consistent with the existence of YSOs in the central hub (Harvey et al. 2008). In addition, the presence of YSOs in clump 47 indicates a density structure more complicated than our simple assumption, which could further reduce the virial mass (Sanhueza et al. 2017), and thus this clump could be even more unstable than our estimation.

In Section 3, we show that the magnetic field orientation around the HFS has two major components, tending to be distributed in the northern and southern part of the system. The two components can be explained by either a curved magnetic field or a foreground/background component overlaid on a uniform magnetic field. Nevertheless, because the C¹⁸O (J = 3−2) spectral data taken by Graham (2008) show that all components in the HFS are within a narrow velocities range (∼3–5 km s⁻¹), the first possibility is favored, unless the foreground/background component coincidentally has a velocity very similar to the HFS.

The curved magnetic field could be originated by an uniform large-scale magnetic field dragged by the contraction of the large-scale main filament. The dragging along the major axis of the large-scale filament would cause the single peak in the large-scale magnetic field to broaden and split into two peaks, and thus the center of the splitting shown in the histogram (∼15°) is consistent with the orientation perpendicular to the large-scale filament. In addition, the spatial distribution tendency of the two components can also be explained, because the contraction along the major axis would lead to an axisymmetric pattern. The feature of parsec-scale magnetic fields being perpendicular to filaments but modified by core collapsing at a smaller scale has been found in other filamentary clouds, such as Serpens South (Sugitani et al. 2011), Orion A (Pattle et al. 2017), and W51 (Koch et al. 2018).

Supercritical main filaments are expected to fragment along their major axis and trigger star formation (e.g., André et al. 2010, 2014; Pon et al. 2011; Miettinen 2012; Clarke et al. 2016), which could be a possible origin of the observed HFS. The main filament connected to our observed HFS has a supercritical mass per unit length (152 M_⊙ pc⁻¹, Arzoumanian et al. 2011), and the submillimeter clumps identified in Johnstone et al. (2017) also indicate that some filament fragmentation has already taken place.

Some theoretical work suggests that fragmentation of filaments with aspect ratios greater than 5 tends to first begin at their ends, where the edge-driven collapsing mode is more efficient than the homologous collapse mode over the whole filament (Pon et al. 2011, 2012). In contrast, the centralized collapse mode is more important in shorter filaments with high initial density perturbations or no magnetic support (Seifried & Walch 2015). The edge-driven collapsing mode is consistent with the HFSs found at the end of the main filament in the IC 5146 cloud. In addition, Graham (2008) found the gas velocity within the main filament increases from the center to both ends, based on ¹³CO and C¹⁸O line observations. The velocity gradient suggests that the gas within the main filament is likely fragmented toward the two massive HFSs at the ends.

The center-to-ends filament fragmentation picture might seem inconsistent with the observed magnetic field morphology, which shows a pattern of end-to-center contraction. The curved magnetic field morphology, however, might be shaped at an early evolutionary stage, when the filament was still contracting and accumulating mass until its density was sufficiently high to trigger fragmentation. In addition, the global end-to-center contraction and the local center-to-end fragmentation could be occurring simultaneously but at different scales, as suggested by hierarchical gravitational fragmentation models (Gómez & Vázquez-Semadeni 2014; Gomez et al. 2018).

To explore the magnetic field morphology within collapsing clouds, Gomez et al. (2018) simulated molecular clouds undergoing global, multiscale gravitational collapse. In this simulation, the magnetic fields would be dragged by the gravitational contraction but eventually reach a stationary state in which the ram pressure of the flow balances the magnetic tension. Hence, the model predicts a random magnetic field morphology on parsec scales and a U-shaped magnetic field within the filaments following the equation

$$\begin{eqnarray}&&{\left(\displaystyle \frac{{v}_{l}}{{v}_{{\rm{A}}}}\right)}^{2}=2\sin (2\alpha ),\end{eqnarray} \tag{ 27 }$$

where v_l is the gas velocity along the filaments, v_A is the Alfvénic velocity and the α is the angle between the magnetic field line and the direction perpendicular to the filament, illustrated in Figure 12. Although the predicted large-scale random magnetic field morphology is inconsistent with the uniform magnetic fields shown by the WLE17 data, a U-shaped magnetic field within the filaments has been observed in our POL-2 data, suggesting that this model might be important when the filaments become dense enough.

Zoom In Zoom Out Reset image size

The observed α is ∼30°, estimated by the two components shown in the PA histogram (Figure 4), so v_l/v_A is expected to be 1.3 by Equation (27). We assume that the v_l is approximately the velocity difference along the filament around the central hub. The line-of-sight C¹⁸O centroid velocity of the central hub (clump 47) is ∼4.1 km s⁻¹, and the western clump 42 has a centroid velocity of ∼3.8 km s⁻¹. Hence, the velocity difference along the filament between clumps 42 and 47 is $0.3/\cos \phi$ km s⁻¹, where ϕ is the inclination angle of the filament with respect to the line of sight. With the Alfvénic velocity of $0.64/\sin \theta$ estimated in Section 3.5, the observed v_l/v_A is $\sim 0.5\sin \theta /\cos \phi$. Due to the unknown inclination angle, we can only speculate that the model expectation would be correct if the filament is nearly perpendicular to the line of sight (ϕ > 67°).

Based on our observed magnetic field morphologies, we propose a three-stage scenario to explain the origin of the observed HFS, illustrated in Figure 13. In the first stage, the large-scale magnetically subcritical filaments are first formed with dynamically important magnetic fields as described in strong magnetic field filament formation models (e.g., Nakamura & Li 2008; Van Loo et al. 2014), and these filaments appear perpendicular to a uniform large-scale magnetic field, as revealed by the WLE17 data. In the second stage, the large-scale filaments accumulate mass via accretion along magnetic field lines or filament mergers (e.g., Li et al. 2010; André et al. 2014) and eventually become magnetically and thermally supercritical. The contraction of supercritical filaments would bend the uniform primordial magnetic fields, similar to the case in Orion A (Pattle et al. 2017). In the third stage, the dense clumps within filaments, often at the ends of filaments, would tend to fragment along magnetic fields and form second-generation filaments with hub-filament morphologies, because density perturbations parallel to the magnetic fields grow faster than those perpendicular to the fields (e.g., Nagai et al. 1998; Van Loo et al. 2014). The collapse of the cores within the second-generation filaments is also regulated by the bent magnetic fields, and so the cores are oriented either parallel or perpendicular to local magnetic fields, as shown in Figure 10, and are less correlated with the primordial magnetic field.

Zoom In Zoom Out Reset image size

Stars form predominantly from high column density filaments (André et al. 2010). Although most filaments are either oriented parallel or perpendicular to the large-scale magnetic fields (Li et al. 2014; Planck Collaboration et al. 2016), only a few young stars have been observed with hourglass magnetic field morphologies, which favor a star formation scenario where the core collapse is regulated by strong magnetic fields (e.g., Girart et al. 2006; Rao et al. 2009; Tang et al. 2009). As a counterexample, ALMA polarization observations toward the embedded source Ser-emb 8 show chaotic magnetic fields (Hull et al. 2017), indicating that this star was formed under weak magnetic field conditions. This difference poses the question of how physical scales and environments generally determine the role of magnetic fields in star formation.

To address the role of magnetic fields in star formation, the SMA polarization survey toward massive cores (Zhang et al. 2014) revealed that magnetic fields on the core scale (0.1–0.01 pc) are mostly either parallel or perpendicular to the magnetic fields on the parsec scales. Li et al. (2015) further analyzed the magnetic field morphologies in NGC 6334 on the 100–0.01 pc scales and found that local magnetic fields on all these scales are either parallel or perpendicular to the local cloud elongation. Both these results suggest that magnetic fields are dynamically important during the collapse and fragmentation of clouds, possibly guiding the contraction of filaments and cores. Koch et al. (2014) further used a large sample set (50 sources) to examine the bimodal distribution of the relative orientation between the magnetic fields and the density structures and found that the distribution is more scattered than those in previous surveys, although a bimodal distribution cannot be ruled out.

In Section 3.4, we find a tendency toward clumps within the observed HFS having orientations parallel or perpendicular to the local magnetic fields (at the 0.1–0.01 pc scale). The local magnetic fields in many of these clumps, however, have orientations 30°–60° different from the parsec-scale magnetic field, which is inconsistent with the findings of Zhang et al. (2014) and Li et al. (2015). The inconsistent cases are mainly clumps within the extending filaments, which follow the orientation of the curved magnetic fields (see Section 4.1). These clumps are much fainter than those in the central hub, which possibly explains why they were missed in previous surveys. Nevertheless, because we still found the orientations of these clumps to be well coupled with the host filaments and the local magnetic fields, our results support the idea that magnetic fields are important in regulating the core and filament collapse on spatial scales of 0.01–0.1 pc.