Observations of Magnetic Fields Surrounding LkHα 101 Taken by the BISTRO Survey with JCMT-POL-2

JCMT is a 15 m diameter telescope situated near the summit of Maunakea. It is the largest single-dish telescope working at submillimeter wavelengths between 1.4 and 0.4 mm. The present work studies the polarized emission received from the LkHα 101 region at 850 μm. The beam size at this wavelength is 141. The JCMT detector, SCUBA-2, consists of four arrays of 40 × 32 bolometers each, covering a solid sky angle of ∼45' × 45'. The BISTRO program has been allocated two observation campaigns, called BISTRO-1 and BISTRO-2, each having 224 hours to map the regions towards Auriga, IC 5146, Ophiuchus, L1689B, Orion A&B, Perseus B1, NGC 1333, Serpens, Taurus B211/213, L1495, Serpens Aquila, M16, DR15, DR21, NGC 2264, NGC 6334, Mon R2, and Rosette. Recently, a third campaign, BISTRO-3, which focuses on mapping various massive clouds, some nearby prestellar cores, and the Galactic Center clouds, has been approved. General descriptions of the BISTRO survey and the measurement of magnetic fields have been reported for Orion A (Pattle et al. 2017; Ward-Thompson et al. 2017), Ophiuchus A/B/C (Kwon et al. 2018; Soam et al. 2018; Liu et al. 2019; Pattle et al. 2019), IC 5146 (Wang et al. 2019), Perseus B1 (Coudé et al. 2019), M16 (Pattle et al. 2018), and NGC 1333 (Doi et al. 2020).

LkHα 101 was the last region to be observed within the framework of the BISTRO-1 program. The data were taken over nine days between 2017 and 2019 with 21 visits and a total integration time of about 14 hours. The last two observations were made on 2019 January 8. Data were read and reduced using Starlink (Currie et al. 2014). In this analysis, we use gridded data with pixel size of 12'' × 12'', similar to the telescope beam. The Stokes parameter Q, U, and I time streams are reduced using pol2map, a POL-2-specific implementation of the iterative map-making procedure makemap (Chapin et al. 2013). The instrumental polarization (IP) is mainly caused by the wind blind; it is corrected for by using the IP model determined from POL-2 measurements with the wind blind in place from observations of Uranus (Friberg et al. 2018). In each pixel, the total intensity, I, the Stokes parameters Q, U, and their uncertainties, δI, δQ, and δU, are provided. The POL-2 data presented in this paper are available at https://doi.org/10.11570/20.0011.

Figure 1 (left) shows the location of the POL-2 observed field in the LkHα 101 area. The 850 μm intensity map, I map, is shown in the right panel of the figure. The center of the map is at R. A. = 4^h30^m15^s, $\mathrm{decl}.=35^\circ 17^{\prime} 05^{\prime\prime}$. Emission is observed to come from a region at the center of the map and a lane of dust in the southwestern direction. It matches the densest area of the region shown in the top-left panel of Figure 2 of Broekhoven-Fiene et al. (2018).

Zoom In Zoom Out Reset image size

Figure 2 shows the maps of the measured Stokes Q (upper left) and U (upper right) and the distributions of their uncertainties, δQ and δU (lower), which are used to calculate other polarization parameters for further analysis. The distributions of δQ and δU have the same mean and rms values of 1.1 mJy beam⁻¹ and 0.5 mJy beam⁻¹, respectively. The mean value of δQ and δU measured in the region is at the level expected for the BISTRO survey (∼3 mJy beam⁻¹ for the 4'' pixel map).

Zoom In Zoom Out Reset image size

The debiased polarized intensity, PI, is calculated using the following formula:

$\begin{eqnarray}&&{PI}=\sqrt{{Q}^{2}+{U}^{2}-\delta {{PI}}^{2}}\end{eqnarray} \tag{ 1 }$

(Montier et al. 2015a, 2015b), where the uncertainty on PI is

$\begin{eqnarray}&&\delta {PI}=\sqrt{\displaystyle \frac{{Q}^{2}\delta {Q}^{2}+{U}^{2}\delta {U}^{2}}{{Q}^{2}+{U}^{2}}}.\end{eqnarray} \tag{ 2 }$

The polarization angle is defined as

$\begin{eqnarray}&&\theta =0.5{\tan }^{-1}\left(\displaystyle \frac{U}{Q}\right)\end{eqnarray} \tag{ 3 }$

and its uncertainty

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

We note that the orientation of the magnetic field line is perpendicular to the polarization angle assuming that the polarization of the emission comes from elongated grains that interact with an underlying magnetic field. The angle of the magnetic field line is east of north, ranging from 0° to 180°. Finally, the polarization fraction, P, and its uncertainty, δP, are calculated as

$\begin{eqnarray}&&P( \% )=100\times \displaystyle \frac{{PI}}{I}\end{eqnarray} \tag{ 5 }$

and

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

In an attempt to select measurements to be retained for further analysis, we plot in Figure 3 correlations relating the uncertainties of the measured parameters (no cut has been applied). From the left panel of the figure, it is clear that δI and δPI are strongly correlated. When I is well measured, so is PI; the converse also holds. It is not the case for other parameters in which the correlations are weak (see Figure A1) except for the case of δP versus δθ (Figure 3, right), where we find that most of the good measurements are confined within a parabola shown as a red solid curve. The equation of the parabola is δP = −0.07δθ² + 1.4δθ. In fact, δP and δθ are related quantities; Serkowski (1962) gives δθ = 28.65°δP/P. In the following analyses, we use only the data obeying δP < −0.07δθ² + 1.4δθ (hereafter called the δP−δθ cut). Because the 850 μm polarized emission from the LkHα 101 region is weak, possibly the weakest in comparison with other regions surveyed by BISTRO, after several attempts, we found that the best compromise between keeping the highest quality data and having good statistics is to use the δP−δθ cut, which can still ensure the robustness of the results and conclusions of the current studies.

Zoom In Zoom Out Reset image size

Figure 4 shows the distributions of I/δI, PI/δPI, and P/δP for the observed region before and after the δP−δθ cut is applied. These signal-to-noise ratios (S/Ns) measure the quality of the measurements of the corresponding quantities I, PI, and P. Their dimensionless values (mean, rms) are (12.1, 19.8), (2.1, 1.9), and (1.7, 2.0) for I/δI, PI/δPI, and P/δP, respectively. When applying the δP−δθ cut, these numbers become (38.5, 22.7), (4.0, 2.4), and (3.9, 2.4): the δP−δθ cut eliminates the parts of the distributions having low S/N. Some more detailed information about the raw data set can be found in the Appendix.

Zoom In Zoom Out Reset image size

Figure 5 shows the map of the inferred B-fields of the LkHα 101 region keeping only pixels obeying the δP−δθ cut. The line segments are rotated 90° from the polarization vectors to follow the magnetic field in the plane of the sky. These are called half-vectors because there is an ambiguity in their directions. After application of the δP−δθ cut, the number of remaining line segments is 419. This map is a direct result of the observations and will be used for further analysis and discussions in the next sections.

Zoom In Zoom Out Reset image size

It was shown by Davis (1951) and Chandrasekhar & Fermi (1953) that turbulent motions generate irregular magnetic fields. Based on the analysis of the small-scale randomness of the magnetic field lines, assuming that the dispersion in the magnetic field angles is proportional to the Alfvén Mach number, the field strength can be estimated. This is called the Davis–Chandrasekhar–Fermi (DCF) method. A variant of the method has been proposed by Crutcher (2012); it gives an estimate of the magnitude of the magnetic field in the plane of the sky, B_POS, as

$\begin{eqnarray}&&{B}_{\mathrm{POS}}=Q\sqrt{4\pi \rho }\displaystyle \frac{{\sigma }_{V}}{{\sigma }_{\theta }}\approx 9.3\sqrt{n({{\rm{H}}}_{2})}\displaystyle \frac{{\rm{\Delta }}{\rm{V}}}{{\sigma }_{\theta }}\quad (\mu {\rm{G}}),\end{eqnarray} \tag{ 7 }$

where Q is the factor used to correct for the line-of-sight and beam-integration effects (Ostriker et al. 2001), ρ is the gas density, σ_V the one-dimensional nonthermal velocity dispersion in km s⁻¹, ΔV = 2.355σ_V, σ_θ is the dispersion of the polarization position angles about a mean B-field in degrees, and n(H₂) is the number density of molecular hydrogen in units of cm⁻³.

In the next sections, we evaluate the magnetic field angle dispersion, σ_θ, a measurable of BISTRO, using two different methods.

To estimate the dispersion of the polarization angle, we first use the "unsharp masking" method which was introduced by Pattle et al. (2017). The principle of the method is to look for the turbulent component of the magnetic field by removing a mean field using a boxcar filter. In practice, this is done by looping over all the pixels of the map, calculating the difference between the original polarization angle map and the intensity-weighted mean map consisting of a 3 × 3 pixel box centered on the considered pixel. The standard deviation of the distribution of deviation angles, Δθ = θ_meas − 〈θ〉, represents the turbulent component of the field and gives the angular dispersion of the region. This process is illustrated in Figure 6.

Zoom In Zoom Out Reset image size

We apply this method to calculate the angle dispersion for the central region and the southwestern dust lane (hereafter called CR for the central region and DL for the dust lane). The CR is defined as having R < 120'' from the map center, and the DL is the region surrounded by an ellipse having a center at RA ∼ 4^h30^m47 and $\mathrm{DEC}\sim +35^\circ 14^{\prime} 47\buildrel{\prime\prime}\over{.} 8$, major × minor axes of 300'' × 108'', and a position angle of 135° (see Figure 1, right, and the white curves in Figure 9). The angle dispersion obtained by requiring ∣Δθ∣ < 90° (to avoid the effect of the ±180° ambiguity of the magnetic field lines) are σ_θ = 17.2° ± 0.4° and 171 ± 05 for the CR and the DL, respectively. To understand how the dispersion depends on the measurement uncertainties, we list in Table 1 the polarization angle dispersion as a function of cuts applied on δI, in addition to the δP−δθ cut. Exploring the δI dependence is sufficient because δI and δPI are strongly correlated (Figure 3 left). We note that δI ranges from ∼0.6 to ∼5.0 mJy beam⁻¹ (Figure A2).

We see from Table 1 that the angle dispersion is robust while applying an additional cut requiring higher S/N measurements of the total and polarized intensities. This is an indication that δP−δθ is a robust cut. Only in the case of the most rigorous cut requiring δI < 1 mJy beam⁻¹ for the DL does the number of remaining half-vectors drastically decrease from ∼190 to 12 with the angle dispersion decreasing from σ_θ ∼ 17° to 58.

The uncertainties of the polarization angle dispersion are also studied following the procedure of this same method (Pattle et al. 2017). Figure 7 displays the dependence of the angle dispersion on the maximum allowed uncertainties, $\delta {\theta }_{\max }$, for each boxcar filter. In practice, this means that the angle dispersion of the map is calculated by requiring that the maximum uncertainty of all the pixels in each 3 × 3 boxcar filter be smaller than $\delta {\theta }_{\max }$. The polarization angle dispersion is expected to increase with the maximum allowed uncertainty of the boxcar filters. Indeed, as can be seen from Figure 7 (left for the CR and right for the DL), the angle dispersion increases as $\delta {\theta }_{\max }$ increases. As was done by Soam et al. (2018) and Liu et al. (2019), from the angle dispersion for different $\delta {\theta }_{\max }$ (Figure 7), we calculate the mean angle dispersion ± uncertainties of σ_θ = 149 ± 27 and 147 ± 26 for the CR and the DL, respectively. These values are kept as the final polarization angle dispersion for further analysis. We note that the results are obtained excluding a first few bins of $\delta {\theta }_{\max }$, where the numbers of included pixels are smaller than 20 and with the application of the δP−δθ cut.

Zoom In Zoom Out Reset image size

Another method to estimate the dispersion of polarization position angles makes use of the so-called structure function (Hildebrand et al. 2009), which is defined as follows. We consider pairs of pixels, at locations x and x + l, each associated with a polarization position angle θ. For a given pixel separation l, the structure function is defined as the mean square angle between the polarization half-vectors of the pair:

$\begin{eqnarray}&&\langle {\rm{\Delta }}{\theta }^{2}(l)\rangle =\displaystyle \frac{1}{N(l)}\displaystyle \sum _{i=1}^{N(l)}{[\theta ({\boldsymbol{x}})-\theta ({\boldsymbol{x}}+{\boldsymbol{l}})]}^{2},\end{eqnarray} \tag{ 8 }$

where N(l) is the number of such pairs. We assume that the magnetic field B can be approximated by the sum of a large-scale structured field of mean amplitude B₀ and a turbulent component δB;  we also assume that the correlation length of the turbulent component is much smaller than the distance over which B varies significantly. For small separations l, we can write

$\begin{eqnarray}&&\langle {\rm{\Delta }}{\theta }^{2}(l)\rangle ={b}^{2}+{m}^{2}{l}^{2}+{\sigma }_{M}^{2}(l),\end{eqnarray} \tag{ 9 }$

where b is the rms contribution of the turbulent component and m measures the contribution of the gradient of B;  in addition, the contribution of the measurement uncertainties, σ_M(l), is included. All of these contributions add in quadrature. Ignoring the contribution of the large-scale structured field in Equation (9) and from the definition of the polarization angle dispersion, σ_θ, we have ${\sigma }_{\theta }^{2}={b}^{2}/2$.

Hildebrand et al. (2009) expressed the ratio of the turbulent field, δB, to the large-scale underlying field (i.e., mean field), B₀, as

$\begin{eqnarray}&&\displaystyle \frac{\delta B}{{B}_{0}}=\displaystyle \frac{b}{\sqrt{2-{b}^{2}}},\end{eqnarray} \tag{ 10 }$

where $b=\sqrt{2}{\sigma }_{\theta }$. We note that Equations (9) and (10) hold only if the correlation length of the turbulent component satisfies δ < l (Hildebrand et al. 2009).

We calculate the structure functions for both the CR and the DL; the result is displayed in Figure 8. At large distances, the structure functions tend to the random field value of ∼52° (Serkowski 1962). A fit to Equation (9) for short distances, 12'' < l < 36'', gives b = 209 ± 68 and 239 ± 71, namely σ_θ = 148 ± 48 and 169 ± 50 for the CR and the DL, respectively. When applying additional cuts on δI, δI < 1 (δI < 1.5) mJy beam⁻¹, we obtain values of σ_θ = 143 ± 61 (145 ± 49) for the CR. For the DL, this cannot be done with the 1 mJy beam⁻¹ cut because too few half-vectors are retained. With the 1.5 mJy beam⁻¹ cut, we obtain σ_θ = 146 ± 31. All of the angle dispersions obtained when applying additional δI cuts are in agreement within one standard deviation with the values obtained when applying the δP−δθ cut only. This once again confirms that the choice of the δP−δθ cut is robust.

Zoom In Zoom Out Reset image size

The dust column density of the Auriga–California molecular cloud has been constructed by Harvey et al. (2013) using four Herschel wavebands at 160, 250, 350, and 500 μm. Using the KOSMA 3 m telescope to detect the CO(J = 2−1) and (J = 3−2) line emissions, Li et al. (2014) studied the morphology of the L1482 molecular filament of the Auriga–California molecular cloud, of which our studied region is part. The column and number densities of 23 identified clumps along L1482 were measured. Our CR, which fits in a circle of 120'' radius, encloses their Clump 10 (an ellipse of 228'' × 110'' major × minor axes); the DL, extending over 600'' × 216'', overlaps with their Clump 12 (an ellipse of 195'' × 68'' major × minor axes), but is significantly larger. We use the column density map (Figure 9) from Harvey et al. (2013) to calculate the average column densities and number density; we find N(H₂) = (0.96 ± 0.39) × 10²² cm⁻² and n(H₂) = 1.22 × 10⁴ cm⁻³ over the CR, and N(H₂) = (1.44 ± 0.53) × 10²² cm⁻² and n(H₂) = 1.25 × 10⁴ cm⁻³ over the DL. The number densities are calculated following the same strategy used by Li et al. (2014) in their Section 3.3.2. The masses of the CR and the DL are calculated using $M=\beta {m}_{{{\rm{H}}}_{2}}{N}_{\mathrm{total},{{\rm{H}}}_{2}}{(D{\rm{\Delta }})}^{2}$, where β = 1.39 is a factor that takes into account the contribution of He in addition to H₂ to the total mass, ${m}_{{{\rm{H}}}_{2}}$ is the mass of a hydrogen molecule, ${N}_{\mathrm{total},{{\rm{H}}}_{2}}$ is the total column density, and Δ = 14'' is the pixel size of Herschel data. For consistency, as was done by Harvey et al. (2013), we use the same distance to the LkHα 101 region of D = 450 pc to calculate the column densities. The radius of the CR is 120'', and the radius of the DL is assumed to be $R=\sqrt{({ab})/2},$ where a and b are the major and minor axes of the DL ellipse. Then, the number density of hydrogen molecules is obtained using ${n}_{H2}=3M/(4\pi {R}^{3}{m}_{{{\rm{H}}}_{2}})$. The relative uncertainties in the number densities are taken to be equal to those of the column densities (42% and 37% for the CR and the DL, respectively).

Zoom In Zoom Out Reset image size

The nonthermal velocity dispersion along the line of sight has been measured by the Purple Mountain Observatory (PMO) 13.7 m radio telescope using the ¹³CO(1−0) transition (Li et al. 2014). The results are σ_V = 0.52 km s⁻¹ (meaning ΔV = 2.355σ_V = 1.22 km s⁻¹) for Clump 10 and σ_V = 0.64 km s⁻¹ (meaning ΔV = 2.355σ_V = 1.52 km s⁻¹) for Clump 12.

We also analyzed archival CO(3−2) JCMT/HARP⁸⁵ data (Buckle et al. 2009) in order to evaluate the velocity dispersion. Figure 10 (left) displays the velocity-integrated intensity map, and Figure 10 (right) displays the intensity-weighted mean Doppler velocity of the region observed by HARP. We note that when we superimpose the intensity map of the 850 μm emission collected by POL-2 (the blue contours in Figure 10 right) on the mean velocity map, we find that the location of the emission matches very well that of the redshifted arc of the cloud. It suggests a cloud emitting polarized light and moving away from us over a blueshifted background (Figure 10 right). A small clump in the northern part of the map is seen to account for most of the blueshifted emission beyond −2 km s⁻¹; its spectrum is shown in the right panel of Figure 11.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The left panels of Figure 11 show the integrated spectra for the CR (left) and the DL (center); they display a two-component structure. The fits of these spectra to two Gaussians give standard deviations of the blue- and redshifted components (σ_B, σ_R) = (0.90 ± 0.02, 0.56 ± 0.01) km s⁻¹ and (0.60 ± 0.03, 0.79 ± 0.02) km s⁻¹ with mean values of (−1.25, 0.89) km s⁻¹ and (−1.08, 0.57) km s⁻¹ for the CR and the DL, respectively. The average temperatures are 29.7 K for the CR and 20.8 K for the DL (Harvey et al. 2013). The corresponding velocity dispersion caused by thermal turbulence are only at per mil level and therefore negligible: the nonthermal FWHM line widths are equal to (ΔV_B, ΔV_R) = (2.12 ± 0.05, 1.32 ± 0.02) km s⁻¹ and (1.41 ± 0.07, 1.86 ± 0.05) km s⁻¹ for the CR and the DL, respectively. Because the redshifted parts of the cloud trace well the 850 μm polarized emission from the region, we use their velocity dispersion to calculate the B-field strength in the next section. However, we conservatively use uncertainties on velocity dispersion estimated from the combination of the three independent measurements added in quadrature: ¹³CO(1−0) emission (Li et al. 2014), HARP blueshifted, and HARP redshifted spectra instead of using only the uncertainties from the fits to the HARP redshifted spectra. The final velocity dispersions retained for further analysis are 1.32 ± 0.40 km s⁻¹ and 1.86 ± 0.19 km s⁻¹ for the CR and DL, respectively. Though the velocity dispersion obtained from CO(3−2) (redshifted part) and from ¹³CO(1−0) agree within ∼20%, which supports the use of the CO(3−2) line emission for calculating the gas velocity dispersion, we note that the CO(3−2) line available to us may not be optically thin. Therefore, the results regarding the magnetic field strength in the current paper are obtained under the assumption that the CO(3−2) line traces the observed dust volume.

The B-fields are measured over a region of ∼1.6 pc at a spatial resolution of ∼0.03 pc. With the irregular mass distribution and the presence of several protostar candidates and an early B star in the region, the magnetic field morphology of the observed region is expected to be complex. In fact, it is the case of the B-field patterns at the CR of the map, the highest emission region, which shows a drastic change of the plane-of-the-sky field directions (see Figures 5 and 13). The field lines are perpendicular to each other, running north–south and east–west. This is an indication of the existence of important field turbulence or of the magnetic field tangling in the dense region. The measured polarization angle dispersion of the order of σ_θ ∼ 15° also supports the existence of the field turbulence. However, it is identified from the measurement of the ratio of the turbulent field to the large-scale underlying field that the B-fields of the whole LkHα 101 region is still generally dominant over the turbulence (see Section 4.3). Going farther away from the center of the map where the density is lower, the B-fields tend to follow the periphery of the matter distribution. However, for the outermost parts, in particular, where the contour level (green curve in Figure 13) is highly curved, the B-fields are perpendicular to the curve.

In the DL, the main orientation of the B-fields has the tendency to follow the filamentary structure (northwest–southeast) of the DL.

A similar tendency of the field running along the filamentary structure is also found in the low-density clumps enclosed by the green contours in Figure 13 scattered around the CR and the DL. This trend is better seen with the elongated clumps.

Prestellar cores, protostars, and then low-mass stars are believed to form in filaments (André et al. 2014). This paradigm is supported by simulations (e.g., Inoue & Inutsuka 2012; Soler et al. 2013). B-fields help funnel matter onto the filaments. Because the DL is a subcritical filament (see Section 4.3), the overall magnetic field direction parallel to its structure is in agreement with the popular picture of the B-field evolution in star-forming regions. Indeed, at large scale, B-fields are typically perpendicular to the main structure of filaments (e.g., Matthews et al. 2014; Planck Collaboration et al. 2015). At the scale of the core and filament size, it is found roughly that magnetic field runs perpendicular to a filament when the filament is gravitationally supercritical, but parallel when it is subcritical (Palmeirim et al. 2013; Ward-Thompson et al. 2017). However, we note that criticality is not the only parameter to decide the configurations of field versus structure. Planck and BLASTPol⁸⁶ data show a parallel-to-perpendicular transition at visual extinction A_V ∼ 3 mag (Planck Collaboration et al. 2016; Soler et al. 2017). In addition, using SOFIA data in Serpens South, Pillai et al. (2020) found another transition from perpendicular back to parallel at A_V ≳ 21 mag. Our DL has an average column density of 1.44 × 10²² cm⁻² (Table 2), which corresponds to A_V ∼ 15.3 mag. The visual extinction is obtained using the standard conversion factor between column densities and visual extinction N(H₂) = 9.4 × 10²⁰ cm⁻² A_V mag (Bohlin et al. 1978). Comparing with what is found in Serpens South, the DL lies in the region where the field is perpendicular to the filament (see Figure 3 of Pillai et al. 2020) and close to the region where the median relative orientation of the field and filament crosses 45°. This suggests that the perpendicular-to-parallel transition may vary depending not only on visual extinction but also on other parameters of a cloud. The CR is an example of this. It is also subcritical but the fields are complex with the presence of protostar candidates and an early B star. More statistics will help to settle this issue.

For the DCF method to be applicable, Ostriker et al. (2001) suggested that the polarization angle dispersion should be smaller than ∼25°, which is the case of the present data set. Using the H₂ number densities and the nonthermal velocity dispersion along the line of sight given in Sections 3.4 and 3.5, Equation (7) gives the magnitude of the magnetic field in the plane of the sky, B_POS. The uncertainties of B_POS are propagated using the following relation:

$\begin{eqnarray}&&\displaystyle \frac{\delta {B}_{\mathrm{POS}}}{{B}_{\mathrm{POS}}}=\sqrt{{\left(\displaystyle \frac{1}{2}\displaystyle \frac{\delta n({{\rm{H}}}_{2})}{n({{\rm{H}}}_{2})}\right)}^{2}+{\left(\displaystyle \frac{\delta {\rm{\Delta }}V}{{\rm{\Delta }}V}\right)}^{2}+{\left(\displaystyle \frac{\delta {\sigma }_{\theta }}{{\sigma }_{\theta }}\right)}^{2}},\end{eqnarray} \tag{ 11 }$

where δn(H₂), δΔV, and δσ_θ are the uncertainties of n(H₂), ΔV, and σ_θ, respectively.

The magnetic field strengths obtained from two methods are listed in Table 2 for the CR and the DL. As shown, the two methods yield a similar magnetic field strength of B_POS ∼ 91 μG for the CR and B_POS ∼ 138 μG for the DL.

It is very interesting to note that the mean measured field strength of ∼115 μG is very close to the value of 100 μG adopted by Zhang et al. (2020) to explain the observed fragmentation length scale (core spacing) of star-forming filaments in the X-shape nebula of the California molecular cloud.

The relative importance of gravity to magnetic fields is usually described by the mass-to-flux ratio, M/Φ. In the units of the critical value, the mass-to-flux ratio is given by the formula from Crutcher (2004),

$\begin{eqnarray}&&\lambda =\displaystyle \frac{{(M/{\rm{\Phi }})}_{\mathrm{observed}}}{{(M/{\rm{\Phi }})}_{\mathrm{critical}}}=7.6\times {10}^{-21}\displaystyle \frac{N({{\rm{H}}}_{2})}{{B}_{\mathrm{POS}}},\end{eqnarray} \tag{ 12 }$

where ${(M/{\rm{\Phi }})}_{\mathrm{critical}}=1/(2\pi \sqrt{G})$, G is the gravitational constant, N(H₂) is the gas column density measured in cm⁻², and B_POS is the strength in μG. As B_POS is the magnetic field component in the plane of the sky, a factor of 3 is introduced to correct for geometrical biases (Crutcher 2004).

Plugging B_POS and N(H₂) obtained for the CR and the DL into Equation (12), we obtain λ for these regions, which are listed in Table 2.

The left panel of Figure 12 shows the measured values of B_POS as a function of the column density N(H₂) for the different regions using the results from Table 2 together with those of previous studies (Pattle et al. 2017; Soam et al. 2018; Coudé et al. 2019; Soam et al. 2019; Wang et al. 2019) surveyed by POL-2. The dashed line corresponds to the separation between subcritical and supercritical regions. The line is obtained from Equation (12), setting λ = 1, where we have an equal contribution of mass and magnetic flux. Compared to other regions, LkHα 101 has particularly low values of N(H₂). The Auriga–California region lies well above the dashed line, which is subcritical.

Zoom In Zoom Out Reset image size

The right panel of Figure 12 shows the variation of B_POS with the mass-to-flux ratio, λ, for the different regions in the left panel. Compared to the other regions, the Auriga–California has a rather low value of λ and lies in the subcritical regime.

Being subcritical at the same time as being the densest region of Auriga–California may help to explain the very low star formation efficiency in comparison with that of the OMC as discussed in Section 1.

However, the star formation efficiency of a cloud depends on several parameters other than λ, such as matter distribution, evolutionary stage, and turbulence. More detailed studies are required to understand why the star formation efficiency in Auriga–California is lower than that of the OMC. Moreover, the measured ratio of the turbulent component to the large-scale component of the magnetic field is δB/B₀ ∼ 0.3, suggesting that the effect of B-fields is dominant over the turbulence in the region.

We now analyze the spatial variation of the polarization fraction within Auriga–California and explore grain alignment physics.

Figure 13 displays the map of the inferred magnetic line segments similar to Figure 5 but the length of the line segments is now proportional to the polarization fraction. It is clearly seen from the figure that the degree of polarization is higher in the more diffuse regions, but it drops significantly in the dense CR with maximum emission intensity (yellow contour). To see explicitly how the polarization fraction, P, changes with the total intensity I, in Figure 14 we show the variation of P(%) with I (the δP–δθ cut has been applied to the data). The observed polarization fraction tends to decrease with increasing total intensity, which is usually referred to as a polarization hole. A power-law fit of the form P ∝ I^−α gives the power index α = 0.82 ± 0.03. The uncertainty in α is obtained from the fit and does not include systematic errors. The value of α found for the LkHα 101 region is in the expected range for molecular clouds, between 0.5 and 1. In other regions surveyed by BISTRO, the estimated values are α = 0.8 for ρ Ophiuchus A (Kwon et al. 2018), 0.9 for ρ Ophiuchus B (Soam et al. 2018), 1.0 for ρ Ophiuchus C (Liu et al. 2019), and 0.9 for Perseus B1 (Coudé et al. 2019). Using a different approach working with non-debiased data, Pattle et al. (2019) found α = 0.34 for Oph A and α = 0.6 − 0.7 for Oph B and C, which are significantly smaller than those obtained by the previous authors.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We note that α is widely used as an indicator of dust grain alignment efficiency. One expects to have α = 0 for constant efficiency of grain alignment and α = 1 for grain alignment that only occurs in the outer layer of the cloud, with complete loss of grain alignment inside the cloud (Whittet et al. 2008). Various observations report the loss of grain alignment (i.e., α ∼ 1) at a large visual extinction of A_V ∼ 20 toward starless cores (Alves et al. 2014; Jones et al. 2014; Santos et al. 2019), which is consistent with the prediction by RAT theory (Hoang et al. 2020). Therefore, the best-fit value α = 0.82 here reveals that grain alignment still occurs inside the cloud, but with decreasing efficiency. This is consistent with the prediction from Radiative Torque (RAT) alignment theory (Lazarian & Hoang 2007; Hoang & Lazarian 2016) that grains inside the cloud can still be aligned due to stellar radiation from the LkHα 101 star.

Because LkHα 101 is the only early B star in the region, we display in Figure 15 the radial dependence of the polarization fraction averaged over 12'' wide rings centered on the star. There is evidence for a rapid decrease of the polarization fraction at distances from the star smaller than ∼150''. LkHα 101 is 15 times more massive than the Sun with a luminosity of 8 × 10³ L_⊙ (Herbig et al. 2004). The angular distance from LkHα 101 at which the mean energy density in empty space is equal to that of the interstellar radiation field, 2.19 × 10⁻¹² erg cm⁻³ (Draine 2010), is ∼3300''. This distance is significantly larger than the distance where we find P to decrease.

Zoom In Zoom Out Reset image size

To better understand the decreasing feature of P, in particular, the drastic decrease at r < 150'', we study the behavior of the polarization fraction P, column density N(H₂), and dust temperature T_dust as a function of distance r from LkHα 101. The two latter quantities are taken from Herschel (Harvey et al. 2013). The upper-left panel of Figure 16 shows the dust temperature map overlaid by 850 μm emission. Figure 16 (upper right) shows the observational data (symbols) and our power-law fit, P = ar^b, with two slopes of b = 0.86 ± 0.04 for r < 130'' and b = 0.35 ± 0.05 for r > 130''. Here, the radius r ∼ 130'' is chosen to be similar to the location of the separation between the CR and the DL (the valley in N(H₂) versus r at r ∼ 130'' in Figure 16, lower right). For the outer region (r > 130''), the polarization degree decreases slowly with decreasing r, which implies a slow decline in the grain alignment efficiency. However, the polarization degree decreases rapidly when approaching the location of LkHα 101 for r < 130'', whereas T_dust increases as expected from stronger heating by the star (Figure 16, lower left), except for only two data points close to the star at r < 25'' with decreased T_dust. Note that LkHα 101 is located at the location of the highest column density (Figure 16, lower right; for the column density map, see Figure 9).

Zoom In Zoom Out Reset image size

In theory, the rapid decrease of P for r < 130'' can arise from (1) significant loss of grain alignment, and (2) strong variation of the magnetic field. According to the popular RAT alignment theory (Lazarian & Hoang 2007; Hoang & Lazarian 2016), the loss of grain alignment is induced by the decrease of the incident radiation field that can align grains and/or the increase of the gas density that enhances grain randomization. In our situation, grains are subject to increasing radiation flux from the LkHα 101 star when r decreases, and the gas density of n(H₂) ∼ 10⁴ cm⁻³ in the CR is not very high. As a result, the degree of grain alignment is expected to increase with increasing local radiation energy density described by T_dust (Hoang et al. 2020), which would result in the increase of dust polarization (Lee et al. 2020). Therefore, the rapid decline of P for r < 130'' (Figure 16, upper right), or T_dust > 25 K (Figure 16, lower left), may challenge the popular theory of grain alignment based on RATs.

Recently, Hoang et al. (2019) suggested that RATs from an intense radiation field can spin grains up to extremely fast rotation. As a result, centrifugal stress can exceed the maximum tensile strength of grain material, resulting in the disruption of large grains into smaller fragments. A detailed modeling of grain disruption toward a dense cloud with an embedded source is presented in Hoang et al. (2020). Because such large grains dominate dust polarization at far-IR/submillimeter, the degree of dust polarization is found to decrease with increasing local radiation energy density (Lee et al. 2020). For the hydrogen density of n(H₂) ∼ 10⁴ cm⁻³ listed in Table 2, the numerical modeling in Lee et al. (2020) implies that the polarization degree first increases with grain temperature due to the increase of the radiation flux, then it decreases when T_dust exceeds ∼25 K. This can explain the decrease of P at small r (or high T_dust) observed in Figure 16. We note that previous observations by Planck (Guillet et al. 2018) and SOFIA/HAWC+ (Tram et al. 2021) also report the decrease of P when T_dust exceeds some value, which were explained by means of rotational disruption by RATs. A detailed modeling to understand the dependence of P versus r is beyond the scope of this paper.

The tangling of magnetic fields is usually invoked to explain the decrease of the polarization fraction, P, with the emission intensity, I, or column density (usually referred to as a polarization hole; see Pattle & Fissel 2019 for a review). However, there is no quantitative study that addresses the role of field tangling in causing the polarization holes at the scales of JCMT observations (or at smaller scales, i.e., CARMA⁸⁷ /SMA⁸⁸ and ALMA⁸⁹ ). The fact that we observe relatively ordered B-field in the region supports the RATD effect but cannot rule out the role of the field tangling. Single-dish versus interferometric observations of class 0 protostellar cores frequently show ordered fields on multiple spatial scales, even if the fields have very different morphologies on large (JCMT ∼ 10,000 au resolution) and small (ALMA ∼ 100 au resolution) scales. JCMT, CARMA, and ALMA maps of Ser-emb 8 (Hull et al. 2017a) and Serpens SMM1 (Hull et al. 2017b) show ordered fields on all scales. The case of SMM1 is particularly striking, where a well-ordered east–west field at JCMT scales (at least on the outskirts of the location of SMM1) is seen, and then a very well-ordered field north–south at CARMA scales, and then a highly complex field with multiple plane-of-sky components at ALMA scales. There are several other examples (see Sadavoy et al. 2018a, 2018b on IRAS 16293 and VLA 1623).

Analyses of the Planck polarization data (in particular, Planck Collaboration et al. 2015 and Planck Collaboration 2020) found that the polarization hole effect can be attributed entirely to turbulent tangling of the magnetic field along the line of sight and that the dust grain alignment efficiency is constant across a wide range of column densities. However, those spatial scales tend to be significantly larger than what we are dealing with in BISTRO observations. Similar conclusions are also reached in work performing similar statistical analyses of ALMA data, at spatial scales closer to the JCMT scales (Le Gouellec et al. 2020).

In summary, polarization holes have been observed in many star-forming regions (see Pattle & Fissel 2019). In the absence of field tangling, polarization holes observed toward protostars are inconsistent with the RAT alignment theory but can be explained by the joint effect of grain alignment and rotational disruption by RATs. There are several possible solutions that are being explored to explain the polarization holes, including RATD, magnetic field tangling. Detailed modeling of dust polarization taking into account grain alignment, disruption, and realistic magnetic fields is required to understand the origins of the polarization hole.