The JCMT BISTRO Survey: An 850/450$\mu$m Polarization Study of NGC 2071IR in OrionB

We present the results of simultaneous 450 $\mu$m and 850 $\mu$m polarization observations toward the massive star forming region NGC 2071IR, a target of the BISTRO (B-fields in Star-Forming Region Observations) Survey, using the POL-2 polarimeter and SCUBA-2 camera mounted on the James Clerk Maxwell Telescope. We find a pinched magnetic field morphology in the central dense core region, which could be due to a rotating toroidal disk-like structure and a bipolar outflow originating from the central young stellar object, IRS 3. Using the modified Davis-Chandrasekhar-Fermi method, we obtain a plane-of-sky magnetic field strength of 563$\pm$421 $\mu$G in the central $\sim$0.12 pc region from 850 $\mu$m polarization data. The corresponding magnetic energy density of 2.04$\times$10$^{-8}$ erg cm$^{-3}$ is comparable to the turbulent and gravitational energy densities in the region. We find that the magnetic field direction is very well aligned with the whole of the IRS 3 bipolar outflow structure. We find that the median value of polarization fractions, 3.0 \%, at 450 $\mu$m in the central 3 arcminute region, which is larger than the median value of 1.2 \% at 850 $\mu$m. The trend could be due to the better alignment of warmer dust in the strong radiation environment. We also find that polarization fractions decrease with intensity at both wavelengths, with slopes, determined by fitting a Rician noise model, of $0.59 \pm 0.03$ at 450 $\mu$m and $0.36 \pm 0.04$ at 850 $\mu$m, respectively. We think that the shallow slope at 850 $\mu$m is due to grain alignment at the center being assisted by strong radiation from the central young stellar objects.

the magnetic field direction is very well aligned with the whole of the IRS 3 bipolar outflow structure. We find that the median value of polarization fractions is 3.0% at 450 μm in the central 3′ region, which is larger than the median value of 1.2% at 850 μm. The trend could be due to the better alignment of warmer dust in the strong radiation environment. We also find that polarization fractions decrease with intensity at both wavelengths, with slopes, determined by fitting a Rician noise model of 0.59 ± 0.03 at 450 μm and 0.36 ± 0.04 at 850 μm, respectively. We think that the shallow slope at 850 μm is due to grain alignment at the center being assisted by strong radiation from the central young stellar objects.

Introduction
Magnetic fields are expected to influence the process of star formation, but their exact role at each evolutionary stage is not yet understood. The well-ordered magnetic field structures, which are observed from molecular cloud scales to protostellar scales, suggest that magnetic fields are an important part of the star formation process (e.g., Orion molecular cloud, Li et al. 2009; Taurus molecular cloud, Chapman et al. 2011;Orion A region, Pattle et al. 2018;NGC 1333IRAS4A, Girart et al. 2006. Theoretical studies predict that magnetic fields affect core collapse, disk, and outflow formation, multiplicity of protostars, and the star formation rate (Shu et al. 1987(Shu et al. , 2004Allen et al. 2003;Nakamura & Li 2005;Price & Bate 2007;Kudoh & Basu 2008;Machida et al. 2011Machida et al. , 2014. Dust polarization observations at submillimeter wavelengths are a good way to study dust properties as well as magnetic field properties in star-forming regions, most of which are heavily embedded in dense molecular clouds (Davis 1951;Goodman et al. 1995;Lazarian 2003). The observed polarization signature is attributed to aspherical dust grains that have their long axes aligned perpendicular to the local magnetic field direction (see reviews, e.g., Andersson et al. 2015;Lazarian et al. 2015;Pattle & Fissel 2019). Aligned dust grains result in the polarization of background starlight at optical/NIR wavelengths (Hall 1949;Hiltner 1949) and polarized thermal dust emission at far-IR/submillimeter wavelengths (Hildebrand 1989). The BISTRO (B-fields In STar-forming Region Observations) Survey project ), using the POL-2 polarimeter (Bastien et al. 2005a(Bastien et al. , 2005bFriberg et al. 2016) on the Submillimetre Common-User Bolometer Array 2 (SCUBA-2) camera (Holland et al. 2013) mounted on the James Clerk Maxwell Telescope (JCMT), provides dust polarization maps toward nearby star-forming regions at 450 μm and 850 μm with angular resolutions of 9 6 and 14 1, respectively (Dempsey et al. 2013). This project provides opportunities for the study of dust properties and magnetic field morphology around star-forming regions at linear scales of 2000 ∼ 5500 au. These observations therefore well fill the gap between the few-pc-scale magnetic field structures revealed by optical/infrared observations (e.g., Tamura et al. 2007;Li et al. 2009;Chapman et al. 2011) and the few-au-scale magnetic field structures seen in submillimeter/millimeter interferometry observations (e.g., Girart et al. 2006;Kwon et al. 2019).
Comparison of the magnetic field strength and structure with outflow properties is important in order to understand the star formation process. On the relatively large dense core scale of ∼0.1 pc, the rotational axis of a core (the outflow direction) is theoretically expected to be aligned with the magnetic field direction in strongly magnetized molecular cloud environments (Fiedler & Mouschovias 1993;Galli & Shu 1993; Allen et al. 2003). On the smaller star formation scale of ∼1000 au, simulations show that both a magnetic field and rotation are essential for launching, collimating, and stabilizing the jet/outflow (Gardiner et al. 2000;O'Sullivan & Ray 2000;Pudritz & Banerjee 2005). Without the inclusion of a magnetic field in their simulations, no stabilized collimated jet/outflow can be made.
Statistical studies of the alignment between outflow and magnetic field directions on both large and small scales have been undertaken, but have not yet produced a conclusive result. For example, Davidson et al. (2011) find good alignment between the magnetic field and outflow directions in the young stellar objects of L1527 and IC348-SMM2. Hull et al. (2013) find no tendency for magnetic field and outflow directions to align, instead finding perpendicular or random orientation of magnetic fields with respect to outflow directions based on observations of 16 low-mass protostars . Another recent study toward 62 low-mass young stellar objects has shown that magnetic field directions are misaligned with outflow directions by 50°± 15°in three-dimensional space, rather than being randomly oriented with respect to one other (Yen et al. 2021).
Magnetic fields have also been studied toward several outflows using molecular line polarization observations. Their results show that magnetic field directions are either parallel or perpendicular to outflow directions. For example, Glenn et al. (1997) found that the position angles of HCO + (1−0) line polarization vectors are 47°± 5°, which is about 20°offset with respect to the DR21 outflow direction. Cortes et al. (2006) found that 12 CO (2−1) line polarization vectors show position angles of 31°.5 ± 8°.7 and 48°.5 ± 7°. 9 in the redshifted velocity components and −50°.2 ± 9°. 2 and 53°.6 ± 8°. 9 in the blueshifted velocity components of the NGC 2071IR outflows. The derived polarization vectors are either parallel or perpendicular to the outflow direction of ∼40°-50°. Lee et al. (2018a) found that the position angles of SiO (8−7) line polarization vectors are mainly parallel to the HH 211 jet axis. Ching et al. (2016) found that 12 CO (3−2) line polarization vectors have a 20°o ffset relative to the NGC1333 IRAS 4A outflow direction. However, there is an ambiguity in the determination of a magnetic field direction due to the Goldreich-Kylafis Effect (Goldreich & Kylafis 1981). Molecular line polarization can be either parallel or perpendicular to the magnetic field direction, depending on the angles between the magnetic field and the line of sight. Meanwhile, it is quite straightforward to study the magnetic field morphology using dust polarized emission at submillimeter wavelengths, since the magnetic field is perpendicular to the polarization vector.
One of the BISTRO targets, NGC 2071IR, has a massive core with several near-infrared objects and nebulae in its central region (Walther et al. 1993;Tamura et al. 2007; Walther & Gaballe 2019).
The distance to our target is 417 ± 5 pc, adopted from the measured distance to NGC 2068 using Gaia data (Kounkel et al. 2018). Within the region, there are many outflows including a powerful northeast-southwest bipolar outflow. IRS 3 (∼1 M e ; Trinidad et al. 2009) has been identified as the launching source of this large-scale outflow, which extends over ∼0.6 pc (Bally 1982;Snell et al. 1984;Eislöffel 2000).
A polarization study toward NGC 2071IR has previously been performed by  using the SCUPOL polarimeter for SCUBA on the JCMT. They found a pinched magnetic field morphology in the central region, which has also been seen in several other collapsing star-forming regions (e.g., Girart et al. 2006;Attard et al. 2009;Rao et al. 2009;Hull et al. 2014;Cox et al. 2018;Lee et al. 2018b;Maury et al. 2018;Sadavoy et al. 2018;Kwon et al. 2019). However,  note that the pinched magnetic field structure of the NGC 2071IR core might be different from other hourglass shapes which are produced via a dragged magnetic field by the infalling material, since the core does not show any internal flattened shape in their continuum map. They derive a magnetic field strength of 56 μG in the dense core region using the Davis-Chandrasekhar-Fermi method (Davis 1951;Chandrasekhar & Fermi 1953). Their study also shows that the polarization fraction decreases with intensity with a power-law slope of −0.79. A polarization study toward the IRS 3 bipolar outflow has been also undertaken by Cortes et al. (2006) using the Berkeley-Illinois-Maryland Association (BIMA) array with a spatial resolution of about 4″. They found that magnetic field vectors inferred from 12 CO (2−1) line polarization observations are parallel to outflow directions.
Thanks to the high sensitivity of POL-2, we can investigate the magnetic field structure in detail not only within the central region but also across the whole area of NGC 2071IR over which the bipolar outflow spreads. Simultaneously obtained 450 μm and 850 μm continuum polarization data also provide an opportunity to study dust grain properties.
The paper is organized as follows. In Section 2, we describe our observations and the reduction of the JCMT POL-2 polarization data. In Section 3, we present the results of our 450 μm and 850 μm continuum polarization observations. We derive the magnetic field strength from polarization angles measured at 450 μm and 850 μm using the modified Davis-Chandrasekhar-Fermi method, and study grain properties using 450 μm and 850 μm polarization fractions. We also investigate the pinched magnetic field morphology and the alignment of the magnetic field direction with respect to the outflow. We finally compare the magnetic energy with the turbulent, gravitational, and outflow energies. In Section 4, we summarize our results. In the Appendix, we describe how we derive IRS 3 bipolar outflow kinetic energy density using HARP 12 CO (3−2) and 13 CO (3−2) molecular line data.

Observations and Data Reduction
Simultaneous 450 μm and 850 μm dust polarization observations were performed toward NGC 2071IR, with a reference position of R. A. = 05 h 47 m 05 040, decl. = 00°21′51 7 (J2000) using POL-2/SCUBA-2 on the JCMT. The NGC 2071IR region was observed twenty times between 2016 September 8 and 2017 November 11 in Band 2 weather (0.05 < τ 225 GHz < 0.08) for a total on-source time of 13.3 hr, using the POL-2-DAISY observing mode which produces maps with a diameter of 12′ (Friberg et al. 2016). The spatial resolutions are 14 1 (corresponding to ∼5800 au at a distance of 417 pc) at 850 μm and 9 6 (∼4000 au) at 450 μm. The achieved rms noise level in the 850 μm Stokes I map is ∼2.0 mJy beam −1 on 12″ pixels within the central 3′ diameter area of constant exposure time. The rms noise level in the 850 μm Stokes Q and U maps is ∼0.8 mJy beam −1 . The achieved rms noise level in 450 μm Stokes I is ∼6.5 mJy beam −1 on 12″ pixels, again within the central 3′ diameter area. The rms noise levels in 450 μm Stokes Q and U are ∼2.9 mJy beam −1 . We note that the derived noises are associated with the instrument, observing technique (Bastien et al. 2005a(Bastien et al. , 2005bHolland et al. 2013;Friberg et al. 2016) and data reduction process (Berry et al. 2005;Chapin et al. 2013).
The POL-2 data reduction process uses the pol2map script recently added to SMURF (Berry et al. 2005;Chapin et al. 2013). The raw timestream data are first converted into separate Stokes I, Q, and U timestreams using the process calcqu. We then produce an initial coadded I map from a set of Stokes I maps created from the Stokes I timestreams of each observation using the iterative mapmaker makemap.
Final I, Q, and U maps are produced using the skyloop task, which is a script that runs makemap simultaneously on the full set of 20 observations in order to find a solution that minimizes residuals across the full set of maps. The initial Stokes I map is used to define a "mask" identifying areas containing astrophysical signals; see Mairs et al. (2015) for a detailed description of the role of masking in SCUBA-2 data reduction. The final coadded set of maps is corrected to account for any small differences in pointing between the input observations. The final I map is used to estimate the Q and U signal caused by instrumental polarization (IP) at each point on the sky, using the "2019 August" IP model provided by the observatory. 86 The IP is approximately 1.5% and 2.3% at 450 μm and 850 μm, respectively, with a small dependence on elevation. The IP has been subtracted from all observations. We also remove 12 CO (3−2) emission contamination from the final 850 μm I map following the process described by Parsons et al. (2018), since NGC 2071IR has a strong CO molecular bipolar outflow covering the whole region (Drabek et al. 2012). For subtracting the CO emission, we use archival JCMT HARP (Heterodyne Array Receiver Program) 12 CO (3−2) data (project code MJLSG11) observed on 2007 November 25. We also use 13 CO (3−2) and C 18 O (3−2) molecular emission lines from these archival data to analyze the IRS 3 bipolar outflow.
The final I, Q, and U maps, with a pixel size of 4″, are combined using the task pol2stack to produce an output polarization vector catalog. All maps are calibrated in units of Jy beam −1 , using flux conversion factors (FCFs) of 725 mJy pW −1 at 850 μm and 962 mJy pW −1 at 450 μm (Friberg et al. 2016). Both FCF values take into account the additional losses of a factor of 1.35 at 850 μm and 1.96 at 450 μm due to the insertion of the POL-2 polarimeter into the SCUBA-2 light path. For 450 μm data, we convolve the 450 μm maps at each 4″ pixel with a 14 1 Gaussian beam using the SMOOTH450 parameter in the pol2map task. This convolution procedure converts the beam resolution of the original 450 μm map, 9 6, to a resolution of 14 1, which allows us to make a direct comparison between the 450 μm and 850 μm data. We bin every 3 × 3 pixels with the 4″ size into a pixel with a 12″ size to improve the signal-to-noise ratio. All further analyses in this paper have been undertaken using these final maps with a pixel size of 12″.
The polarization angle θ and polarization fraction P values used in this paper follow the conventional definitions of q = U Q 0.5 arctan( ) and = + P Q U I 2 2 0.5 ( ) , respectively. Polarization angle is measured from the north, increasing eastward. The debiased polarization fraction is where δI, δQ, and δU are measurement errors in those maps. We use the debiased data for all magnetic field studies, while we use the nondebiased data when using a Rician noise model to study the relationship between polarization fraction and total intensity. More detail on the data reduction process is given by Pattle et al. (2021). Figure 1 shows the final POL-2 850 μm polarization vector map, with vector selection criteria of (I/ δI) 10 (white vectors) and (I/δI) 20 (blue vectors). Figure 2 shows the magnetic field vector map, with each 850 μm polarization vector rotated by 90°, with vector selection criteria of (I/δI) 10 and (P/δP) 3. Figure 3 shows the 450 μm polarization vector map (left panel) and magnetic field vector map (right panel) with vector selection criteria of (I/δI) 50 and (P/δP) 3. Figure 4 shows, for the central region, a comparison of polarization angles and fractions at 450 μm and 850 μm, with a pixel size of 12″. We note that the polarization angles at the two wavelengths show a good agreement in the central dense region.
We take 220 GHz continuum and C 18 O (2−1) molecular line emission data from the Submillimeter Array (SMA) archive (Project code 2009B-S004) to study the kinematics in the central region. The SMA observation was carried out using the compact configuration, with projected baselines ranging from 8 to 76 m. The imaging process is undertaken using the Miriad software. The synthesized map using the natural weighting has a beam size of 3 4 × 2 8, with a position angle of −5°. Three continuum emission peaks are detected, with peak intensities of ∼0.22 Jy for IRS 1, ∼0.09 Jy for IRS 2, and ∼0.37 Jy for IRS 3.

Magnetic Field Strength Using the DCF Method
Knowledge of magnetic field strength is crucial in order to determine whether magnetic fields are important in star formation processes. We derive the plane-of-sky magnetic field strength (B pos ) from the 450 μm and 850 μm polarization data using the modified Davis-Chandrasekhar-Fermi method (DCF) (Davis 1951;Chandrasekhar & Fermi 1953) provided by Crutcher et al. (2004). Their equation is 86 is the mean molecular weight per hydrogen molecule (Kirk et al. 2013); σ ν is the velocity dispersion in km s −1 ; n s D = n 8 ln 2 is the FWHM velocity dispersion in km s −1 ; σ θ is the polarization angle dispersion in degrees; and Q is a factor to account for the unresolved complex magnetic field and density structure within the beam size, which we take to be 0.5 as suggested by Ostriker et al. (2001) (see also Pattle et al. 2017).
To determine the magnetic field strength, we first derive a polarization angle dispersion of 10°.7 ± 4°.9 within the central 2 5 × 2 5 region (marked as a blue box in Figure 2) using the unsharp masking method developed by Pattle et al. (2017). Following their procedure, we calculate a mean polarization angle in a 3 × 3 pixel box (corresponding to 36″ × 36″) using the relation q = U Q 0.5 arctan (¯¯), where the barred quantities, U Q and¯, represent the average values of these quantities within the box. We then calculate an angle dispersion of q q i i |¯| within the central 2 5 × 2 5 region using the equation where θ i is the polarization angle in the ith 12″ × 12″ pixel, and q ī is the mean polarization angle in the 3 × 3 pixel box centered on the ith pixel. The angle dispersion uncertainty of 4°.9 is from the median value of angle uncertainties of the selected polarization vectors. The nonthermal velocity dispersion, σ ν ∼ 0.45 km s −1 , is derived from the C 18 O (3−2) molecular line velocity dispersion obtained with HARP in the central region, σ ν , C 18 O ∼ 0.46 km s −1 , using the equation where k is the Boltzmann constant. The adopted temperature is 20.1 ± 1.8 K, estimated using Herschel data (Könyves et al. 2020). Therefore, the FWHM corresponding to the nonthermal velocity dispersion is n s D =ñ 8 ln 2 1.07 km s −1 . This value is almost identical to the FWHM of NH 3 (1,1) molecular line measured by Takano et al. (1986).
We next derive magnetic field strengths toward the same central point using two differently sized regions. One is for the central 60″ × 60″ region, representing a linear size of 0.12 pc, and the other is for the central 100″ × 100″ region, with a linear size of 0.2 pc. In the first compact core region, we obtain a column density of N(H 2 ) = 1.08(±0.62) × 10 23 cm −2 based on the measured total flux of 31.05 ± 3.42 Jy at 850 μm and adopting a dust opacity range of κ 850 μm = 0.018 ± 0.009 cm 2 g −1 , assuming a gas-to-dust ratio of 100 (Ossenkopf & Henning 1994). The derived average volume density and total mass in this region are n(H 2 ) = 3.67(±2.12) × 10 5 cm −3 and M = 36.3 ± 21.0 M e , respectively, assuming a uniform-density cylindrical shape with a length of L = 0.12 pc. We then derive a plane-of-sky magnetic field strength of 563 ± 421 μG using Equation (1).
In the central 100″ × 100″ area, representing a volume 4.5 larger than that of the 60″ × 60″ area, we obtain a column density of N(H 2 ) = 4.61(±2.67) × 10 22 cm −2 based on the measured total flux of 36.9 ± 4.1 Jy at 850 μm. The derived average volume density and total mass are N(H 2 ) = 9.41(±5.34) × 10 4 cm −3 and M = 43.1 ± 24.9 M e , respectively, assuming a uniform-density cylindrical shape with a length of L = 0.2 pc. The estimated average plane-of-sky magnetic field strength is 285 ± 212 μG.
We adopt B pos = 563 ± 421 μG as a representative magnetic field strength in NGC 2071IR, since the central 60″ × 60″ region holds most of the mass of NGC 2071IR (see Table 1). The 100″ × 100″ region is only 20% more massive than the central 60″ × 60″ region, despite being 4.5 times larger in volume. We additionally derive an angular dispersion of σ θ = 14°.0 ± 5°. 0 in the 450 μm polarization data, again using the unsharp masking method, within the central 2 5 × 2 5 region. The plane-of-sky magnetic field strength is therefore calculated to be about 429 ± 278 μG using the 450 μm  polarization data. Since the rms noises of our Q and U measurements at 850 μm are about three times smaller than those at 450 μm, we take B pos = 563 ± 421 μG measured at 850 μm as a representative magnetic field strength in NGC 2071IR.
Our derived magnetic field strength, 56 μG, is about ten times stronger than the strength derived in . Their smaller magnetic field strength is because, compared with our values, they used a larger derived polarization angle dispersion of 33°and a smaller volume density of N(H 2 ) ∼ 10 4 cm −3 , based on the critical density of the NH 3 (1,1) molecular line. Our adopted unsharp masking method for the polarization angle dispersion measurement ) provides a smaller value than the one measured with the method in . This is because the unsharp method enables us to measure the angle dispersion along a curved mean field direction, while the method used in  assumes a uniform mean field direction. Our derived volume density of 3.67(±2.12) × 10 5 cm −3 is also more likely to be accurate than the density value assumed from ammonia detection.
In order to compare our polarization data obtained using POL-2 with those of  using SCUPOL, we measure polarization angle dispersion and mean polarization angle from POL-2 data using the same method used by them. We derive a polarization angle dispersion of σ θ = 38°a nd a mean polarization angle of q =  33 , respectively. These obtained values are almost the same as theirs, which are σ θ = 33°and q =  34 , respectively. This result confirms that both data sets are consistent with each other at the central strong emission region.

Polarization Properties
We investigate the debiased and nondebiased polarization fractions as functions of total intensities at both 850 μm and 450 μm using two different fitting methods, the single power law, and the Rician-mean model. We used the following power-law model are the observed polarization fraction and its uncertainty, respectively, P i m is the polarization fraction calculated from the power-law model or the Rician-mean model, and N is the number of data. By minimizing χ 2 , we find the best-fitting parameters a s P and QU for the power-law and Rician-mean models. Table 2 lists the fitting results. Figure 5 shows the debiased and nondebiased polarization fractions as a function of total intensity at 850 μm within the central 3′ area. In the left panel of Figure 5, we fit the debiased polarization fraction as a function of the normalized total intensity with the power law. We obtain best-fit parameters of α = 0.35 ± 0.03 and P σQU = 0.09 ± 0.02 using 158 polarization vectors selected using the criterion of (I/δI) 10, where σ QU = 0.80 mJy beam −1 . For 96 polarization vectors selected using the criterion of (P/δP) 3, the best-fit parameters for the power-law relation are a =  =  s P 0.43 0.04 and 0.17 QU 0.05. The right-hand panel of Figure 5 shows the Rician-mean model fitting to the nondebiased polarization fraction as a function of total intensity for 166 data points without any criterion condition. We obtain the best Rician-mean model fitting parameters a =  =  s P 0.36 0.04 and 0.10 0.03 QU . The power-law index, 0.35 ± 0.03, fitted to debiased data points with the selection criterion based on the total intensity ((I/δI) 10), is almost the same as the slope value, 0.36 ± 0.04, obtained from the Rician-mean model fitting to the whole nondebiased data set (see Table 2). However, the slope of the power law, 0.43 ± 0.04, fitted to the data points (filled circles in the left panel of Figure 5) selected with the (P/ δP) 3 criterion, is steeper than the slope from the Ricianmean model. This is because the selection criterion based on the polarization fraction introduces a systematic bias by excluding data points with lower polarization fractions as shown in the left-hand panel (see also Figure 2 in . Figure 6 shows the debiased and nondebiased polarization fractions as a function of total intensity at 450 μm within the central 3′ area. The left-hand panel of Figure 6 shows a single power-law fit to the debiased polarization fraction as a function of the normalized total intensity with σ QU = 2.93 mJy beam −1 . We obtain best-fit parameters of a =  = s P 0.57 0.03 and QU  0.49 0.11 for the power-law relation using 150 polarization vectors selected using the criterion of (I/δI) 10. For 105 polarization vectors selected using the criterion of (P/δP) 3, the best-fit parameters for the power-law relation are a =  0.63 =  s P 0.03 and 0.82 0.17 QU . The right-hand panel of Figure 6 shows the Rician-mean model fitting to the full set of 164 data points without any selection criterion. The obtained best Ricianmean model fitting parameters are a =  = s P 0.59 0.03 and QU  0.55 0.12 (see Table 2). The power-law index, 0.57 ± 0.03, fitted to debiased data points with the selection criterion based on the total intensity, is almost the same as the slope value, 0.59 ± 0.03, obtained from the Rician-mean model fitting to the the Rician-mean model gives a best estimation of the α parameter from data without any selection criteria on nondebiased polarization fractions. We demonstrate that, from the comparisons of the slopes from the power-law and Rician-mean models, the slopes of the power-law model estimated from selected data, based on the polarization fraction criterion, have steeper slopes than those of the Rician-mean model.
We find that the polarization fraction decreases with intensity at both 450 μm and 850 μm. This trend of decreasing polarization fraction has often been shown in the dense regions of low-mass/massive star-forming molecular cores, such as Bok globules CB 54 and DC 253-1.6 (Henning et al. 2001), dense cores within the dark cloud Barnard 1 in Perseus molecular cloud (Matthews & Wilson 2002), massive molecular cores of W51 e1 and e2 (Lai et al. 2001), and giant . We note that negative error bars for some data points are not shown due to the logarithmic scale of the vertical axis. Some error bars are too small to be seen compared to the symbol size. , respectively. We note that negative error bars for some data points are not shown due to the logarithmic scale of the vertical axis. Some error bars are too small to be seen compared to the symbol size. molecular cloud complex Vela C (Fissel et al. 2016). The suggested possible explanations for this trend are due to: (1) the loss of the grain alignment because of the grain growth (Bethell et al. 2007;Brauer et al. 2016) or scattering and absorption of photons at high density (King et al. 2018(King et al. , 2019; (2) the optical depth effect (Dowell 1997;Hildebrand et al. 1999); or (3) the unresolved complex B-field structures along the line of sight. It is, however, beyond the scope of this paper to constrain them using only these data.
We note that the Rician-mean model fitting slope of 0.36 ± 0.04 at 850 μm is much shallower than the slope of 0.59 ± 0.03 at 450 μm. We test whether it is due to the optical depth effect at 450 μm and 850 μm. The optical depth at 850 μm is smaller than that at 450 μm by a factor of The Rician-mean model fitting slope of 0.36 ± 0.04 obtained from the 850 μm polarization data is almost the same as the slope of 0.33 toward Oph A derived by . They suggest that the shallower slope of 0.33 in Oph A, compared to a slope of 0.76 toward both Oph B and C, could be due to strong radiation from the B spectral type star S1, which may help maintain grain alignment in the higher optical depth region of Oph A. As is the case in Oph A, NGC 2071IR shows better grain alignment than do Oph B and C. The central three infrared young stellar objects, IRS 1, 2, and 3 in NGC 2071IR could be radiation sources, which are assisting grain alignment.
We find that the median value of polarization fractions is 3.0% at 450 μm with the condition of (P/δP) 3 in the central 3′ region, which is larger than the median value of 1.2% at 850 μm. The ratio of the median polarization fractions at 850 and 450 μm, P 850 /P 450 , is 0.37. The trend of a smaller polarization fraction at a larger wavelength has been also observed in M17, shown by Zeng et al. (2013). They interpreted that it is due to the better alignment of warmer dust in the strong radiation environment. In Figure 7, we visualize the polarization fraction distributions at 450 μm and 850 μm and the polarization ratio of P 850 /P 450 . It is clearly shown that the polarization fractions at both wavelengths decrease as the intensity increases. Combining the result of a shallow Rician-mean model fitting slope of 0.36 at 850 μm, and the higher polarization fraction at 450 μm than 850 μm (P 850 < P 450 ), we suggest that the central strong radiation from YSOs are indeed assisting the grain alignment suggested by the radiative alignment torques (RAT) theory (Lazarian & Hoang 2007).
Figure 7(c) shows the distribution of P 850 /P 450 values, which are <1, slightly increasing toward the center. It is because the dust alignment persists more efficiently at 850 μm compared to 450 μm, as presented by their Rician-mean model fitting slopes of 0.36 at 850 μm and 0.59 at 450 μm. We note that there might be some changes in dust grain properties (e.g., grain growth) in the center where the P 850 /P 450 ratio shows a small increase. Morphology   Figures 1 and 2 show the 850 μm polarization vector map and magnetic field vector map (rotated by 90°from the polarization vector map) of NGC 2071IR, respectively. We confirm the pinched magnetic field structure along a northeast and southwest direction around the center, as shown in the blue box of Figure 2. This pinched structure was also seen by . Pinched magnetic field structures have typically been interpreted to have been formed by the ambipolar diffusion process in dense star-forming cores (e.g., Girart et al. 2006;Attard et al. 2009;Rao et al. 2009;Hull et al. 2014;Cox et al. 2018;Lee et al. 2018b;Maury et al. 2018;Sadavoy et al. 2018;Kwon et al. 2019). Therefore, one interpretation of this pinched magnetic field morphology could be magnetically regulated core collapse. However, observational evidence of infalling material in this system is required to confirm this hypothesis.

Pinched Magnetic Field
Another possible interpretation is that the magnetic fields are shaped by a rotating disklike structure and bipolar outflows from the IRS 3 young stellar object, which could create an apparent hourglass morphology (see, e.g., the schematic diagram of L1448 IRS 2 shown in Figure 3 of Kwon et al. (2019)). In Figure 8(a), we show the JCMT HARP 12 CO (3−2) and C 18 O (3−2) molecular line emission maps. We find that C 18 O (3−2) emission is distributed perpendicular to the direction of the IRS 3 bipolar outflow traced by the 12 CO (3−2) line emission (black contours). In addition, we find a slight velocity gradient from northwest to southeast.
Careful examination of the central part of Figure 8(b) shows that the peak emission of the redshifted C 18 O (3−2) velocity component (red contours) is located northwest of the peak emission position of the blueshifted velocity component (blue contours). In Figure 8(c), the SMA C 18 O (2−1) molecular line emission over the blueshifted velocity range of 7.5 ∼ 9.4 km s −1 and the redshifted velocity range of 9.4 ∼ 11.0 km s −1 , respectively. The blue contour levels are 2.5, 2.7, 2.9, 3.1, 3.3, 3.5, 3.7, and 3.9 K km s −1 , and the red contour levels are 2.5, 2.9, 3.3, 3.7, 4.1, and 4.5 K km s −1 . Black contours represent the integrated 12 CO (3−2) emission, which traces a bipolar outflow from IRS 3. The northeastern part is the integrated emission over the blueshifted velocity range of −45.3 ∼ −19.9 km s −1 , and the southwestern part is the integrated emission over the redshifted velocity range of 37.6 ∼ 63.0 km s −1 . Black contour levels are 3σ, 10σ, 30σ, and 50σ, where σ = 0.7 K km s −1 . The blue contour levels are 3.6σ, 4.4σ, 5.2σ, 6.0σ, and 6.8σ, and red contour levels are 3.6σ, 7.2σ, 10.8σ, and 14.4σ, where σ = 0.25 Jy beam −1 km s −1 . The three 220 GHz SMA continuum emission peaks, IRS 1, 2, and 3, are marked as small circles in all three figures. emission map of the very central region shows a clear velocity gradient from the northwest redshifted velocity component to the southeast blueshifted velocity component. This gives us a hint that there is a rotating disklike structure perpendicular to the bipolar outflow direction. In Figure 9, we present the magnetic field vector map, overlaying a cartoon showing our suggestion of the rotating disklike structure and bipolar outflow. The shaded blue and red colors schematically represent outflow cavity walls and the shaded gray color represents a disklike structure. We suggest that the pinched magnetic field morphology is formed by the outflow cavity walls together with the rotating disklike structure at the center.

Alignment Between Outflow and Magnetic Field Directions
As discussed in Section 2, there is no conclusive consensus regarding the alignment between outflow and magnetic field directions. We consider NGC 2071IR to be an excellent environment in which to study the alignment issue. Thanks to the high sensitivity of our POL-2 map compared to the previous SCUPOL map presented by , we can investigate the magnetic field morphology in the outflow of IRS 3 in detail.
In Figure 10, we show the magnetic field vectors along with HARP 12 CO (3−2) and 13 CO (3−2) molecular line emission, which trace the bipolar outflow of IRS 3 well. In both panels, all of the magnetic field vectors shown are selected using the criteria of (I/δI) 10 and (P/δP) 5, derived from the 850 μm dust polarization data. In the left-hand panel, we present the high-velocity components of the 12 CO (3−2) molecular line emission, showing the blueshifted velocity range of −45.3 ∼ −19.9 km s −1 (blue contours) and the redshifted velocity range of 37.6 ∼ 63.0 km s −1 (red contours). We do not use the low-velocity component of 12 CO (3−2) to present the outflow, since the low-velocity 12 CO (3−2) emission is contaminated by the ambient molecular cloud. In the right-hand panel, we present the low-velocity components of the 13 CO (3−2) molecular line emission, showing the blueshifted velocity range of −1.1 ∼ 5.1 km s −1 (blue contours) and the redshifted velocity range of 13.4 ∼ 18.2 km s −1 (red contours).
In Figure 10, we show that the overall magnetic field direction is well aligned with the direction of the bipolar outflow. We also see the small-scale directional variation of magnetic field vectors over the region, like the curving shape in the blueshifted outflow component seen in both the 12 CO (3−2) and the 13 CO (3−2) emission maps. The directions of the magnetic field vectors are well aligned with the cavity walls shown in the low-velocity 13 CO (3−2) emission map in the right-hand panel of Figure 10. Due to this, we suggest that the magnetic field morphology of NGC 2071IR is formed by the strong influence of the outflow.

Comparison of Energy Densities
We compare the magnetic energy density with the turbulent, gravitational, and outflow energy densities of the NGC 2071IR star-forming region. First, we calculate the magnetic energy density using the equation where the three-dimensional magnetic field strength, B, is estimated from the plane-of-sky magnetic field strength as » p B B 4 pos assuming that, statistically, the three-dimensional magnetic field direction is isotropically distributed (Crutcher et al. 2004). The derived magnetic energy density is u B = 2.04( ± 3.05)×10 −8 erg cm −3 using B pos = 563 ± 421 μG. Second, we find the turbulent kinetic energy density, = u turb 3 2 ρ = v 3.01 turb 2 (±1.73) × 10 −8 erg cm −3 , adopting v turb = 1.07 km s −1 and ρ = μm H N(H 2 ) = 1.75(±1.01) × 10 −18 g cm −3 . We adopt a volume density of N(H 2 ) = 3.67(±2.12) × 10 5 cm −3 , which is a value derived in the central 60″ × 60″ area. Third, we derive the gravitational energy density assuming a uniform-density sphere with a radius of 30″. The estimated gravitational energy density is u g = 4πGρ 2 R 2 /5 = 1.81(±1.51) × 10 −8 erg cm −3 , where G is the gravitational constant. Finally, we obtain the kinetic energy density of the outflow, u outflow = 2.33(±0.32) × 10 −8 erg cm −3 , using HARP 12 CO (3−2) and 13 CO (3−2) emission lines. We describe how we obtain the energy of the bipolar outflow using CO emission lines in the Appendix.
We find that the magnetic, turbulent kinetic, gravitational, and outflow kinetic energy densities are all comparable to each other in the NGC 2071IR star-forming region (see Table 3). We note, however, that the kinetic energy density of the outflow is a lower limit since we count only the radial velocity of the outflow. We do not have any information about the inclination of the outflow. Therefore, NGC 2071IR is in a state of energy equipartition now, even though the gravitational energy might have been the largest term at an earlier epoch in order to have formed stars in the central region.

Summary
We summarize the results of our 450 μm and 850 μm polarization observations toward the massive star-forming region NGC 2071IR as follows: Figure 9. Magnetic field vectors derived from 850 μm polarization data using selection criteria of (I/δI) 10 and (P/δP) 3. The blue and red contours, which are the same as those in Figure 8, represent the HARP C 18 O (3−2) emission. We overlay a cartoon showing our suggestion of a rotating disklike structure and outflow. The shaded blue and red colors schematically represent outflow cavity walls, and the shaded gray color represents a disklike structure.
1. We derive a plane-of-sky magnetic field strength of B pos = 563 ± 421 μG in the central 60″ × 60″ region from polarization angle dispersion of 10°.7 ± 4°. 9 measured at 850 μm polarization data using the modified Davis-Chandrasekhar-Fermi method (Davis 1951;Chandrasekhar & Fermi 1953) provided by Crutcher et al. (2004). We also derive B pos = 429 ± 278 μG with the measurement of polarization angle dispersion, 14°.0 ± 5°.0, at 450 μm polarization data. 2. We find that the median value of polarization fractions is 3.0% at 450 μm in the central 3′ region, which is larger than the median value of 1.2% at 850 μm. The trend could be due to the better alignment of warmer dust in the strong radiation environment (Zeng et al. 2013). We also find that polarization fractions decrease with intensity at both wavelengths, with slopes, determined by fitting a Rician noise model of 0.59 ± 0.03 at 450 μm and 0.36 ± 0.04 at 850 μm, respectively. The obtained slope at 850 μm is similar to the result in Oph A, 0.33, derived by . We suggest that the dust grains of NGC 2071IR are well aligned, as in Oph A, with the support of radiation from the central IRS 1, 2, and 3 young stellar objects. 3. We confirm the pinched magnetic field morphology in the central area. This can be explained if the magnetic field is shaped by the rotating disklike structure and the bipolar outflow of IRS 3.
The James Clerk Maxwell Telescope is operated by the East Asian Observatory on behalf of The National Astronomical Observatory of Japan, Academia Sinica Institute of Astronomy and Astrophysics, the Korea Astronomy and Space Science Institute, and the Center for Astronomical Mega-Science (as well as the National Key R&D Program of China with No. 2017YFA0402700). 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. Additional funds for the construction of SCUBA-2 were provided by the Canada Foundation for Innovation. D.J. is supported by the National Research Council of Canada and by a Natural Sciences and Engineering Research Council of Canada (NSERC) Discovery Grant. A.S. acknowledges the financial support from NSF through grant AST-1715876. W.K. is supported by the New Faculty Startup Fund from Seoul National University. F.K. and L.F. acknowledge

Appendix
We describe how we derive the IRS 3 bipolar outflow kinetic energy density using HARP 12 CO (3−2) and 13 CO (3−2) molecular line data.
First, we estimate the approximate optical depth of the 12 CO (3−2) and 13 CO (3−2) molecular lines from the assumed isotopic line ratio of X = where T A * is the corrected antenna temperature for atmospheric attenuation, scattering, and spillover (see Figure A1). This gives t~4.8 CO 12 and t~0.08 CO

13
, with a ratio of R ∼ 13, in the blueshifted outflow, in the velocity range −0.5 ∼ 8.5 km s −1 . In the redshifted outflow, the optical depths are t~3.1 CO 12 and t~0.05 CO 13 , with a ratio of R ∼ 19, in the velocity range 10.5 ∼ 18.5 km s −1 . We do not include the central velocity range of 8.5 ∼ 10.5 km s −1 in the optical depth estimation, since 12 CO (3−2) is contaminated by ambient molecular emission.