Distortion of Magnetic Fields in the Dense Core CB81 (L1774, Pipe 42) in the Pipe Nebula

The detailed magnetic field structure of the starless dense core CB81 (L1774, Pipe 42) in the Pipe Nebula was determined based on near-infrared polarimetric observations of background stars to measure dichroically polarized light produced by magnetically aligned dust grains in the core. The magnetic fields pervading CB81 were mapped using 147 stars and axisymmetrically distorted hourglass-like fields were identified. On the basis of simple 2D and 3D magnetic field modeling, the magnetic inclination angles in the plane-of-sky and line-of-sight directions were determined to be 4° ± 8° and 20° ± 20°, respectively. The total magnetic field strength of CB81 was found to be 7.2 ± 2.3 μG. Taking into account the effects of thermal/turbulent pressure and magnetic fields, the critical mass of CB81 was calculated to be Mcr = 4.03 ± 0.40 M⊙, which is close to the observed core mass of Mcore = 3.37 ± 0.51 M⊙. We thus conclude that CB81 is in a condition close to the critical state. In addition, a spatial offset of 92″ was found between the center of the magnetic field geometry and the dust extinction distribution; this offset structure could not have been produced by self-gravity. The data also indicate a linear relationship between polarization and extinction up to AV ∼ 30 mag going toward the core center. This result confirms that near-infrared polarization can accurately trace the overall magnetic field structure of the core.


Introduction
CB81 (cataloged by Clemens & Barvainis 1988, also referred to as Pipe 42 by , and as L1774 in the catalog of Lynds 1962) is a dense core in the Pipe Nebula, representing a nearby dark cloud complex having a filamentary shape located at a distance of = -+ d 130 58 24 pc (Lombardi et al. 2006, see also Alves & Franco 2007;Dzib et al. 2018). A search for young stellar objects (YSOs) was previously conducted using the Spitzer telescope at 24 and 70 μm, covering the entire surface of the Pipe Nebula, and no YSOs were reported toward CB81 (Forbrich et al. 2009). A search for (diskless) pre-main-sequence stars based on X-ray observations also resulted in no such objects being detected in the region of CB81 (Forbrich et al. 2010). Thus, CB81 can be considered a starless dense core, despite a report of the possible detection of a CO outflow (Yun & Clemens 1992 near an IRAS source toward CB81. NH 3 observations were conducted toward CB81 in radio wavelengths . The temperature of the core was reported to be T kin =12.4±0.4 K. The velocity dispersions, σ v , associated with the NH 3 (1,1) and NH 3 (2,2) lines, were determined to be 0.11 km s −1 and 0.09 km s −1 , respectively. Using a value of 0.1 km s −1 for the NH 3 velocity dispersion for the core, the calculated turbulent velocity dispersion, σ turb , is 0.0627 km s −1 . Franco et al. (2010) employed R-band linear polarimetry toward the CB81 region in the Pipe Nebula, and reported a polarization angle of θ R =158°.1±5°.35. The magnetic field direction in the CB81 region is roughly perpendicular to the global filament of the Pipe Nebula.
In addition to the studies described above, many surveys and statistical analyses have been performed toward the Pipe Nebula, including extinction map studies (Lombardi et al. 2006;Lada et al. 2008;Román-Zúñiga et al. 2010), dust emission studies using data acquired from the Herschel satellite (Peretto et al. 2012), radio wave molecular line studies (Onishi et al. 1999;Muench et al. 2007;Rathborne et al. 2009), polarimetric studies , and dust temperature studies (Forbrich et al. 2014;Hasenberger et al. 2018). However, there have only been a few research projects focusing specifically on CB81. CB81 is located in the "stem" region of the Pipe Nebula, and appears bright in far infrared maps generated by the Herschel satellite. This core is also evident in a previously reported dust extinction map (Román-Zúñiga et al. 2010), with a peak A V of greater than 20 mag. CB81 is the densest core in the "stem" region and is relatively isolated from neighboring cores (see Figure 4 of Román-Zúñiga et al. 2010). Moreover, the Pipe Nebula has a rich stellar backdrop consisting of galactic bulge stars. These characteristics make CB81 an ideal object for studying the magnetic fields of dense cores.
Since the search of the hourglass-shaped magnetic fields toward dense cores is rather limited (Kandori et al. 2017a; hereafter Paper I; Kandori et al. 2019), we chose CB81 for the target in our near-infrared (NIR) polarimetric surveys of dense cores. Increasing the number of cores with hourglass-shaped magnetic fields is important to take statistics of the total magnetic field strength, mass-to-flux ratio, and polarizationextinction relationship through three-dimensional (3D) modeling of the hourglass fields (Kandori et al. 2017b, hereafter Paper II;Kandori et al. 2018a, hereafter Paper III;Kandori et al. 2018c). These statistics are closely related to the study of the initial conditions of star formation, the core formation process, and the magnetically aligned dust in dense environments.
In the present study, wide-field background star polarimetry at NIR wavelengths was conducted for CB81 and the Stokes I data in the H and K s bands were used to map the dust extinction characteristics of CB81 as a means of determining the density structure of the core (Appendix). The plane-of-sky magnetic field structure was ascertained using several thousands of stars in and around the core radius, while removing contamination from the ambient field component. In addition, the total magnetic field strength of the core was estimated based on the Davis-Chandrasekhar-Fermi method (Davis 1951;Chandrasekhar & Fermi 1953) and the 3D magnetic field modeling of the core. Using the resulting magnetic field information, the kinematical stability of CB81 is discussed herein. The relationship between polarization and extinction (P-A) for CB81 was derived by subtracting ambient polarization and correcting for the depolarization effect caused by distortions of the magnetic field shape and the line-of-sight magnetic inclination angle. The resulting linear P-A relationship indicates that the observed polarization reflects the overall magnetic field structure in the core.

Observations and Data Reduction
We observed CB81 using the JHK s -simultaneous imaging camera SIRIUS (Nagayama et al. 2003) with the associated polarimetry mode SIRPOL (Kandori et al. 2006) on the IRSF 1.4 m telescope at the South African Astronomical Observatory (SAAO). IRSF/SIRPOL is a useful instrument for NIR polarization surveys because it simultaneously provides wide-field JHK s polarization images (7 7×7 7 with a scale of 0 45 pixel −1 ). The uncertainty in the polarization degree values resulting from sky variations during exposures is typically 0.3%. The uncertainty in the determination of the polarization angle origin (i.e., the correction angle) is less than 3° (Kandori et al. 2006;Kusune et al. 2015).
Observations of CB81 were conducted on the nights of 2017 June 15 and 24. A number of 15 s exposures were performed at four half-waveplate angles (in the sequence of 0°, 45°, 22°.5, and 67°.5) at 10 dithered positions (equivalent to one set). The total integration time was 3000 s (20 sets) per waveplate angle.
The data thus acquired were reduced using the Interactive Data Language (IDL) software. This reduction included flat-field correction with twilight flat frames, median sky subtraction, and frame combining after registration and was performed in the same manner described by Kandori et al. (2007). Note that the accuracy of the calibration based on twilight flat images is discussed in the Appendix. Point sources with peak intensities greater than 10σ above the local sky background could be detected on the Stokes I images, and aperture polarimetry was performed for a number of sources for each waveplate position angle image (I 0 , I 45 , I 22.5 , and I 67.5 ). In this manner, 5159, 9576, and 5895 sources were detected in the J, H, and K s bands, respectively. The aperture radius equaled the FWHM of the stars (2.8 pixels for the J, H, and K s bands), and the sky annulus was set to 10 pixels with a 5 pixel width. The limiting magnitudes were 18.5, 18.0, and 17.0 mag in the J, H, and K s bands, respectively. A relatively small aperture radius was used to suppress flux contamination from neighboring stars. We did not use psf-fitting photometry because the goodness of fit for each star on different waveplate angle images can cause systematic errors in polarization measurements. Sources with photometric errors greater than 0.1 mag were rejected.
The Stokes parameter for each star was derived using the relationships I=(I 0 +I 45 +I 22.5 +I 67.5 )/2, Q=I 0 −I 45 , and U=I 22.5 −I 67.5 . The polarization degree, P, and polarization angle, θ, were determined from the equation = + P Q U I 2 2 and q = U Q 0.5atan( ). Because P is a positive quantity, the resulting P values tend to be overestimated, especially in the case of low signal-to-noise ratio sources. We corrected for this bias using the equation & Kronberg 1974). In the present study, we discuss the results obtained working in the H band, for which dust extinction effects are less severe than those in the J band, and the polarization efficiency is greater than that in the K s band. Figure 1 presents the resulting polarization vector map for CB81 in the H band. Here, stars for which P H /δP H 5 are shown. CB81 is located slightly to the east from the center of the image, appearing as a region of dark obscuration. The white circle indicates the core radius (R=177″) determined from NIR extinction studies (see the Appendix).

Distortion of Magnetic Fields
Although the polarization distribution is relatively uniform outside the core radius, a kink or bent structure is evident in the polarization vectors inside the radius. The polarization vectors located outside the core are regarded as off-core polarizations generated by the ambient medium. As shown in Figures 2 and  3, the off-core vectors were found to have values of = P H,off  2.40% 0.60% and q =    152 . 9 8 . 9 H,off (representing the median value and one standard deviation). The off-core polarizations are relatively well ordered ( Figure 3) and have similar polarization degrees ( Figure 2). The polarization angle, q H,off , is consistent with the value obtained based on observations in the R band toward the CB81 region, q =    158 . 1 5 . 35 R (Franco et al. 2010). Following the same procedure described in our previous paper (Paper I), we fitted the off-core vectors on the sky plane in Q/I and U/I. The distributions of the Q/I and U/I values were subsequently modeled as , where x and y are the pixel coordinates, and A, B, and C are the fitting parameters. Note that we did not employ higher-order spatial fitting (e.g., second-order fitting). The F-test based on residuals of the firstand second-order fittings resulted in significance values of 58.6% for Q/I and 66.3% for U/I, both of which greatly exceed the 5% or 1% level, and these fittings cannot be distinguished. Comparing the plane fitting to the median value obtained from subtraction analysis (i.e., the zero order), the F-test gives significance values of 0.6% for Q/I and 9.3×10 −8 % for U/I. Therefore, it is evident that the plane-fitting process is suitable.
The off-core vectors generated in this work are provided in Figure 4. The regression vectors of the off-core polarization   components were subsequently subtracted from the original vectors in order to isolate the polarization vectors associated with CB81 ( Figure 5). Comparing Figures 1 and 5, the ordered fields going from the northwest to southeast disappear, and a distorted polarization distribution appears in and around the core. Following the subtraction, the polarization degree of the off-core vectors is successfully suppressed to give a value close to 0% (Figure 2), while the polarization angle becomes randomly distributed ( Figure 3). Therefore, the subtraction analysis was evidently successful.
The distance to the off-core stars was examined using the Gaia DR2 data set (Gaia Collaboration et al. 2016 and all the Gaia DR2 sources (with faint limiting magnitudes of G∼20 mag) distributed within 4′ from CB81 were found to be located in the background to the core. Thus, when determining the ambient off-core polarization vectors, contamination from foreground stars could be ignored.
Interestingly, in Figure 5, the distribution of the axisymmetrically distorted fields seems to have been offset based on the dust obscuration distribution. We thus treated the center of magnetic field geometry as a free parameter in the model fitting process presented in the subsequent parabolic model section.
It is notable that our analysis of the subtraction of the ambient off-core polarization component may look wrong, because the analysis may remove the very low polarization angle dispersion (i.e., a strong magnetic field component that laces both the surrounding diffuse interstellar medium and dense molecular cloud). However, this is not the case. The Stokes parameters observed toward the core can be written as = = + Q I q q q obs obs ISM core ( ) and = = + U I u u obs obs ISM ( ) u core , where q ISM and u ISM show the polarization component for the diffuse surrounding medium and q core and u core show the polarization component arisen in the core. The ambient subtraction analysis just isolates q core and u core , and this analysis does not affect the alignment of polarization vectors solely associated with the core.

Parabolic Model
The most probable configuration for the magnetic field lines, estimated using a parabolic function and its shift and rotation, is shown in Figure 6 (solid white lines). In this figure, 147 polarization vectors having P H 0.5% are included and the field of view is the same as the diameter of CB81 (354″) in the δ direction, with a value of 294″ in the α direction. The size of the field of view is therefore different in the α and δ directions. This is the case because the center of the distorted field is located slightly to the east relative to the center of obscuration, and the eastern end of the distorted field was not fully incorporated into the field of view during the observations. A comparison of the distorted field to the A V distribution is presented in Figure 7, in which the white and red crosses indicate the center of the distorted field and the centroid center of the A V distribution, respectively (see the Appendix). It is evident from this image that there is a displacement between the center of the magnetic field structure and the center of the mass distribution.
The best-fit parameters were determined to be θ mag =4°±8°a nd C=1.43(±0.68)×10 −5 pixel −2 (=7.06×10 −5 arcsec −2 = 6.4×10 −9 au −2 ) for the parabolic function = + y g gCx 2 , where g specifies magnetic field lines, θ mag is the position angle of the magnetic field direction (from north through east), and C determines the degree of curvature of the function. Note that we used the parabolic function in the 90°-rotated form so that, in the case that θ mag is 0°, the x-direction corresponds to the direction of decl.
A parabolic function was employed because this is the simplest means of approximating the analytically described hourglassshaped magnetic field model (Mestel 1966;Ewertowski & Basu 2013;Myers et al. 2018; see also Kandori et al. 2020, hereafter Paper VI, for a comparison of the parabolic function to the hourglass model). In a subsequent paper, we intend to use the Myers model to fit the same data and to discuss the core formation mechanism.
The observational error associated with each star was taken into account when calculating c dq i 2 , where n is the number of stars, x and y are the coordinates of these stars, θ obs and θ model denote the polarization angles from observations and from the model, and δθ i is the observational error) during the fitting procedure. The coordinate origin of the parabolic function is R.A.=17 h 22 m 47 86, decl.=−27°04′ 40 5 (J2000), as determined by searching for the minimum χ 2 for the variables x and y in and around the core. The center coordinate of the core (that is, the centroid center of the A V distribution; see the Appendix for the derivation of the A V map) is R.A.=17 h 22 m 41 42, decl.=−27°05′12 8 (J2000). The angular distance between these two centers is 92″, which corresponds to roughly half the core radius.
The uncertainties to determine the center coordinates were estimated to be 18 0 and 16 7 for the magnetic center and the extinction center, respectively. For the latter value, we used the difference of the centroid center between the A V distribution and the Herschel 500 μm emission map (André et al. 2010;Peretto et al. 2012;Roy et al. 2019) as the 1σ uncertainty. We discuss the origin of this displacement in Section 3.4.
The parabolic fitting was evidently satisfactory, since the standard deviation of the residual angles, θ res =θ obs −θ fit , was smaller when using the parabolic function (δθ res =20°.03± 0°. 78, see Figure 8) than in the case of a uniform field (δθ uni =30°.78±0°.69). The uncertainties associated with δθ res and δθ uni were estimated using the bootstrap method. In this process, a random number following a normal distribution with the same width as the observational error was added for each star, and fitting was performed using the parabolic function. This process was repeated 1000 times to ascertain the dispersion of the resulting δθ res distribution. Comparing δθ res and δθ uni , there is a statistically significant difference, suggesting that a distorted field is most likely present. This was also confirmed from an F-test with a significance of 3.3×10 −7 %, indicating that δθ res and δθ uni are statistically different. The intrinsic dispersion, dq dq dq = int res 2 err 2 1 2 ( ) , estimated using the parabolic fitting, was found to be 18°.65±0°.84 (0.3260 radian), where δθ err is the standard deviation of the observational error in the polarization measurements.
Assuming frozen-in magnetic fields, the intrinsic dispersion of the magnetic field direction, δθ int , can be attributed to perturbation of the Alfvén wave by turbulence. The strength of the plane-of-sky magnetic field (B pos ) can be estimated from the relationship pr s dq = B C 4 pos corr 1 2 turb int ( ) , where ρ and σ turb are the mean density of the core and the turbulent velocity dispersion (Davis 1951;Chandrasekhar & Fermi 1953), and C corr =0.5 is a correction factor suggested by theoretical studies (Ostriker et al. 2001;see also, Heitsch et al. 2001;Padoan et al. 2001;Heitsch 2005;Matsumoto et al. 2006). Using the mean density (ρ=3.93 (±0.79)×10 −20 g cm −3 ; see the Appendix) and turbulent velocity dispersion (σ turb = 0.0627 kms 1 ) from Rathborne et al. (2008) and the δθ int derived in this study, a weak magnetic field is obtained as the lower limit of the total field strength ( B | |): B pos =6.8 μG. For the purposes of error estimation, we assumed that the turbulent velocity dispersion value was accurate within 30%. Thus, the uncertainty in B pos was estimated to be 2.2 μG.

3D Magnetic Field
The 3D magnetic field modeling was performed following the procedure described in a previous paper (Paper II, see also Paper VI). The 3D version of the simple parabolic function, j = + z r g g gCr , , ) in the cylindrical coordinate system (r, z, j) was used for the modeling of the core magnetic fields, where g specifies the magnetic field lines, C is the curvature of the lines, and j is the azimuth angle (measured in the plane 354 (0.19 0.23 pc) in the α and δ directions at a distance of 130 pc. The field size in the δ direction is equal to the diameter of the core. The background image for this figure is a magnified version of that in Figure 5. The white lines indicate the direction of the magnetic field inferred from parabolic fitting. The scale bar above the image indicates 5% polarization. perpendicular to r). Using this function, the magnetic field lines are axially symmetric around the r axis, and the value of z(r, j, g) is not affected by the parameter j.
The 3D model was virtually observed after rotating in the line-of-sight (γ mag ) and the plane-of-sky (θ mag ) directions. The analysis was performed in the same manner described in Section 3.1 of Paper VI (see also Sections 2 and 3.1 of Paper II). The resulting polarization vector maps for the 3D   . c 2 distribution of the polarization angles (c q 2 ). The best magnetic curvature parameter (C) was determined at each q inc . q =  0 inc and 90°c orrespond to the edge-on and pole-on geometries in the magnetic axis. parabolic model for various γ mag are shown in Figure 15 of Paper VI. The density distribution of the model core was calculated based on the Bonnor-Ebert sphere model (Ebert 1955;Bonnor 1956) with a solution parameter of ξ max =13.0 (see the Appendix). The density distribution of CB81 was shifted from the center of the magnetic field configuration, as described in Sections 3.1 and 3.2, so the Bonnor-Ebert density distribution in the 3D model was shifted by 92″ to be consistent with observations. The model map can be compared with observations for each parameter set (see Figure 15 of Paper VI). Since the polarization distributions in the model core differ from one another depending on the viewing angle toward the line of sight, χ 2 fitting of these distributions with the observational data can be used to restrict the line-of-sight magnetic inclination angle and 3D magnetic curvature. Figure 9 summarizes the distribution of c q = å - , where n is the number of stars, x and y show the coordinates of stars, θ obs and θ model denote the polarization angle from observations and the model, and δθ i is the observational error) calculated using the model and observed polarization angles. The optimal magnetic curvature parameter, C, was determined at each inclination angle γ mag to obtain c q 2 . From the polarization angle fitting, it is clear that the large inclination angle (g    30 40 mag -) is unlikely. Note that the large g mag model shows a radial polarization pattern that does not match observations. The distribution of c q 2 is relatively flat for the region g   20 mag . The minimization point for c q 2 is g =  20 mag with an uncertainty of 24°. Since c q 2 increases rapidly from g =  20 mag going toward larger inclination angles, the uncertainty in this direction should be less than 24°. We can therefore conclude that the g mag value is 20°, with an uncertainty of 20°. The magnetic curvature obtained at g =  20 mag is =´-C 2.5 10 4 arcsec 2 . Note that the c 2 values in Figure 9 are relatively large. This result can likely be attributed to the scatter of the polarization angles caused by Alfvén waves, which cannot be included in the observational error term in the calculation of c 2 . Figure 10 shows the best-fit 3D parabolic model along with the observed polarization vectors. The direction of the model vectors generally agrees with the observations. The standard deviation of the differences in the plane-of-sky polarization angles between the 3D model and the observations is 15°.24, which is less than the 2D fitting result, as well as that for a uniform field (30°.78). Figure 11 presents the same observational data as Figure 10 but with the background image processed using the line integral convolution (LIC) technique (Cabral & Leedom 1993). Here, we used the publicly available IDL code developed by Diego Falceta-Gonçalves. The direction of the LIC "texture" is parallel to the magnetic field direction, and the background image is based on the polarization degree of the model core.

Magnetic Properties of the Core
Using the determined magnetic inclination angle, g mag , of 20°, the total magnetic field strength of CB81 was found to be  Although CB81 was found to be magnetically supercritical, this does not necessarily indicate that the core was in an unstable state. The critical mass of CB81, taking into account both magnetic and thermal/turbulent support effects equal to ) K cm −3 (see the Appendix). Therefore, we conclude that CB81 is in a condition not far from the critical state. This critical (or marginally stable) state of CB81 is consistent with the starless characteristics of the core.
The relative importance of magnetic fields with regard to the support of the core was investigated using the ratios of thermal and turbulent energy to the magnetic energy, b º C V 3 s 2 A 2 and b s s where C s , s turb , and V A denote the isothermal sound speed at 12.4 K, the turbulent velocity dispersion, and the Alfvén velocity. These ratios were found to be b~12.4 and b~1.12 turb . Thermal support is therefore dominant in CB81, with a small contribution from static magnetic fields and turbulence. The magnetic field is weak and the turbulence appears to have largely dissipated.
The A V distribution (Figure 15; see the Appendix) of CB81 is almost circular from the center to the boundary, although there is a filamentary, tail-like structure extending toward the northwest. The plane-of-sky magnetic axis (PA=4°) of CB81 is roughly perpendicular to this filamentary structure, and the circular shape of CB81 is consistent with its weak magnetic field strength. While a magnetohydrostatic configuration can produce a flattened inner-density distribution whose major axis is perpendicular to the direction of magnetic fields (Tomisaka et al. 1988), the resulting density distribution can be circular if M M mag BE is small. The magnetic field strength of CB81 is weak, but the polarization vectors shown in Figures 6 and 7 are relatively well-aligned along the distorted field, and the magnetic fields appear to be strong enough to ensure alignment. Li et al. (2015) conducted magnetohydrodynamic (MHD) simulations and concluded that, in the slightly magnetically supercritical case (l = 1.62), magnetic fields are well-aligned. In contrast, in the very magnetically supercritical case (l = 16.2), magnetic field lines are highly tangled by the turbulent motions, and this does not match observations. Comparing these two scenarios, CB81 is closer to the former case involving moderately magnetically supercritical conditions (l = 4.04), and the relatively good alignment of magnetic fields in CB81 is consistent with the theoretical simulation. Figure 11. Same data as presented in Figure 10, but with the background image generated using the line integral convolution (LIC) technique (Cabral & Leedom 1993). The direction of the LIC "texture" is parallel to the direction of the magnetic fields and the background image is based on the polarization degree of the model core.
The origin of the spatial offset of  92 between the center of the magnetic field configuration and the A V distribution is presently uncertain. This offset may result from the initial core formation conditions. Because the magnetic field around the core was found to be relatively uniform, it is reasonable to assume that that the initial magnetic field distribution during core formation was also uniform. The density of the medium surrounding CB81 is likely~10 3 cm −3 (P T ext kin for the core; see the Appendix), and the magnetic fields should be frozen in the medium at this density. In this case, the initial mass distribution would not be expected to be uniform, to account for the observed offset hourglass-shaped magnetic fields. In the absence of external forces such as turbulence or shocks, the center of mass for the medium and the center of the distorted magnetic fields should be in the same position, because the mass will move toward the common center of gravity and the frozen-in magnetic fields will be dragged by the medium. Thus, turbulence and/or shocks are required to migrate the medium toward a position different from the center of gravity. In this simple scenario, an external force such as turbulence or shocks plays a role in core formation, and a core formation scenario driven solely by self-gravity can be rejected. Therefore, the existence of offset hourglass-shaped fields may restrict the formation scenario of dense cores.

Polarization-Extinction Relationship
As described in a previous paper (Paper III), the relationship between the dichroic polarization of dust (P) and extinction (A) in dark clouds is important. This relationship affects the observational interpretation of the interstellar polarization angle, which in turn is closely related to the plane-of-sky magnetic field direction, and also the investigation and restriction of dust grain alignment models. The former point is essential to guarantee that our polarization observations can trace the dust alignment and magnetic fields deep inside the CB81 core, while the latter point is important for the comparison of observations with dust alignment theories.
The observed polarizations toward CB81 represent superpositions of the polarization from the core and from the ambient medium surrounding the core. The polarization arising from the ambient off-core medium should therefore be subtracted from the observed polarization in order to isolate the polarization associated with the core. As described in Section 3.3, distorted magnetic fields surrounding the core can also produce a depolarization effect. Furthermore, the line-ofsight inclination angle in a magnetic axis can affect observed plane-of-sky polarizations. These effects should be corrected in order to accurately ascertain the relationship between polarization and extinction in dense cores. Figure 12(a) shows the observed P H versus -H K s relationship with no correction (original data). The polarization is seen to generally increase with increases in the extinction value, up to -H K 0.6 s mag, but exhibits a highly scattered distribution above -H K 0.6 s mag. The correlation coefficient of the relationship is only 0.49. The observed polarization is the superposition of the polarization from the core itself and ambient polarization. The off-core stars located outside the core radius, >  R 177 , can be used to estimate off-core polarization vectors to be subtracted from the observed polarizations toward the core (see Section 3.1), and a relatively linear P H versus -H K s relationship was obtained after subtracting ambient polarization components as shown in Figure 12(b). The slope of the relationship is 1.90±0.07% mag 1 . CB81 is associated with inclined distorted magnetic fields, as discussed in Sections 3.1-3.3, which provide depolarization effects inside the core. Based on the known 3D magnetic field structure, a depolarization correction factor was estimated, as shown in Figure 13. This figure was created simply by dividing the polarization degree map of the distorted field model with edge-on geometry (g =  0 mag ) by the inclined distorted field model (g =  20 mag ). The correction factor map can thus simultaneously correct for both the depolarization and lineof-sight inclination angles. In Figure 13, the factors distributed around the equatorial plane are less than unity (indicating depolarization). This is due to the crossing of the polarization vectors at the front and back sides of the core along the line of sight (see the explanatory illustration in Figure 7 of Kataoka et al. 2012). Because the center of mass and the center of the distorted magnetic fields are different in CB81, this slightly breaks the symmetry of the distribution of depolarization correction factors. Figure 12(c) shows the depolarization and inclination-corrected P-A relationship obtained by dividing the Figure 12(b) relationship by the correction factor map ( Figure 13). In Figure 12(c), the P-A relationship is steeper than that in Figure 12(b), reflecting the correction. The slope of the relationship is 3.16±0.09% mag 1 . The correlation coefficients obtained in these analyses were 0.82 (ambient subtraction, Figure 12(b)) and 0.81 (depolarization and inclination correction, Figure 12(c)). Both are higher than the value of 0.49 determined for the original relationship.
Comparing Figures 12(a) and (c), the changes in both scatter and slope are dramatic. A relatively linear relationship between polarization and dust extinction in the range up toÃ 30 V mag was obtained. This result confirms that our NIR polarimetric observations accurately trace the overall polarization (magnetic field) structure of CB81. Figure 14 presents the P A H V versus A V diagram. The corrected P H data as shown in Figure 12   fit to the data points with The linear relationship in Figure 12(c) is reflected in the shallow slope of the linear fitting result. The dotted line shows the fitting of the entire data set using the power law µ a P A A H V V , giving an α index of−0.47±0.10. The power-law fitting appears to be identical to the linear fitting for > A 5 V mag. The dotted-dashed line presents the observational upper limit reported by Jones (1989). The relationship was calculated based on the equation ) and the parameter η is set to 0.875 (Jones 1989). t K denotes the optical depth in the K band and » P A 0.62 The polarization efficiency of CB81 is about half that of FeSt 1-457 (Paper VI). CB81 and FeSt 1-457 are associated with the same dark cloud complex (Pipe Nebula), and this difference in efficiency implies that the dust properties and/or the radiation environment (based on radiative torque grain alignment theory: Dolginov & Mitrofanov 1976;Draine & Weingartner 1996, 1997Lazarian & Hoang 2007) can vary even in the same dark cloud complex.

Summary and Conclusion
The present study ascertained the detailed magnetic field structure of the starless dense core CB81 (L1774, Pipe 42) based on NIR polarimetric observations of background stars to measure dichroically polarized light produced by aligned dust grains in the core. After subtracting ambient polarization components, the magnetic fields pervading CB81 were mapped using 147 stars, and axisymmetrically distorted hourglass-like magnetic fields were identified. On the basis of simple 2D and 3D magnetic field modeling, magnetic inclination angles relative to the plane-of-sky and line-of-sight directions were determined to be    4 8 and    20 20 , respectively. Using these angles and the Davis-Chandrasekhar-Fermi method, the total magnetic field strength of CB81 was found to be Although the resulting M cr value is slightly greater than the core mass, M core (marginally stable), we conclude that CB81 is close to the critical state, wherẽ M M cr core . A spatial offset of 92″ was determined between the center of the magnetic field configuration and the A V distribution, possibly due to the initial conditions related to the formation of the core. Because the magnetic field around the core is fairly uniform, a simple interpretation is that the initial mass distribution during core formation was not uniform but rather biased, and this configuration was subsequently compressed by turbulence or shocks to create the observed offset hourglass-shaped magnetic fields. Self-gravity cannot drive such a process, so the existence of these offset hourglassshaped magnetic fields suggests a structure that restricts the core formation process. We obtained a linear relationship between the polarization and extinction values, up toÃ 30 V mag toward the stars with the deepest obscuration. The slope value obtained from this relationship, 3.16±0.09% mag −1 , is relatively small. This result can possibly be ascribed to a lack of a large grain population or the absence of a strong radiation field, or a combination of both effects. The linear relationship indicates that the observed polarizations reflect the overall magnetic field structure of the core. Further theoretical and observational studies are desirable with regard to explain the dust alignment in the dense core environment.
We are grateful to the staff of SAAO for their kind help during the observations. We wish to thank Tetsuo Nishino, Chie Nagashima, and Noboru Ebizuka for their support in the development of SIRPOL, its calibration, and its stable operation with the IRSF telescope. The IRSF/SIRPOL project was initiated and supported by Nagoya University, National Astronomical Observatory of Japan, and the University of Tokyo in collaboration with the South African Astronomical Observatory under the financial support of Grants-in-Aid for Scientific Research on Priority Area (A) No. 10147207 and No. 10147214, and Grants-in-Aid No. 13573001 and No. 16340061 of the Ministry of Education, Culture, Sports, Science, and Technology of Japan. R.K., M.T., N.K., K.T. (Kohji Tomisaka), and M.S. also acknowledge support by additional Grants-in-Aid No. 16077101, 16077204, 16340061, 21740147, 26800111, 16K13791, 15K05032, 16K05303, 19K03922.

A.1. Dust Extinction Map and the Bonnor-Ebert Model Fitting
Dust extinction (i.e., A V ) measurement, particularly at NIR wavelengths, is one of the most straightforward methods for revealing the density structure within dense dark clouds and cores. There are two approaches to these measurements: the star count method (e.g., Wolf 1923;Dickman 1978;Cambrésy 1999;Dobashi et al. 2005) based on determining the stellar density distribution and the NIR color excess (NICE) method (Lada et al. 1994;NICER by Lombardi & Alves 2001;NICEST by Lombardi 2009;XNICER by Lombardi 2018) based on the measurement of the stellar color excess. Considering the stellar densities at the H and K s bands in the densest region of CB81, we elected to use the latter approach, because the dust extinction toward the core center was not severe and we could detect many reddened background stars through the core.
The extinction map was produced following the procedure described in Kandori et al. (2005). Assuming that the population of stars across the CB81 field is invariable, we were able to use the mean stellar color in the reference field, á -ñ H K s ref , as the zero-point for color excess in the CB81 core. This reference field is situated at the northwest and southeast edges of the observation field. The color excess distribution of the core was obtained by subtracting á -ñ H K s ref from the observed H−K s color of stars. Since the extent of the analyzed field was relatively small (∼8′×8′), our assumption of an invariant stellar population is considered to be reasonable. The color excess distribution was determined by arranging circular cells (30″ in diameter) spaced 9″ apart on the observation field. The mean stellar color in each cell was calculated and the color excess was derived according to the ) is the color index of the ith star in a cell, N is the number of stars in a cell, and á -ñ H K s ref is the mean color of stars in the reference field. The H−K s color excess was converted to A V using =´-A E 21.7 V HK s (Nishiyama et al. 2008), following which the map was smoothed with an FWHM=15″ Gaussian filter. The resulting A V map is shown in Figure 15. The typical uncertainty associated with A V in the reference field (Ã 0 V mag) was estimated to be ∼0.5 mag. The position of the core center based on the centroid of the A V distribution is R.A.=17 h 22 m 41 42, decl.=−27°05′12 8 (J2000).
Following the procedures detailed in Sections 4.1 and 4.2 of Kandori et al. (2005), we produced a radial column density profile for CB81 and conducted the fitting using the Bonnor-Ebert model (Ebert 1955;Bonnor 1956). Note that we manually masked extinction features likely to be unrelated to CB81 (i.e., a diffuse streaming feature to the northwest in Figure 15) and masked regions were ignored in the circular averaging procedure. Figure 16 presents the results of the Bonnor-Ebert fitting of CB81. Here, the dots and error bars represent the average N(H 2 ) values at each annulus at  4. 5 intervals and the root mean square (rms) deviation of data points in each annulus, respectively. The solid line denotes the radial column density profile of the best-fit Bonnor-Ebert sphere convolved with the beam used in the A V measurements.
The dotted-dashed line indicates the Bonnor-Ebert model profile before the convolution.

A.2. Physical Properties of CB81
The physical properties we obtained for the CB81 core include a radius, R, of 23,000±1000 au (    177 7 ),  K, where T kin is the kinematic temperature measured using NH 3 lines and T turb is the temperature calculated using a turbulent velocity dispersion of 0.0627 km s −1   (Lai et al. 2003), the consistency of T values obtained from the Bonnor-Ebert fitting and from the independent radio measurements indicates that the distance measured by Lombardi et al. (2006) is reasonably accurate. The starless dense core CB81 is now known to be a Bonnor-Ebert core with x~13.0 max . Because its critical value exceeds x = 6.5 max , CB81 is considered to be unstable with respect to Bonnor-Ebert equilibrium. As shown by Kandori et al. (2005), a series of solutions indicating an unstable Bonnor-Ebert equilibrium is indistinguishable from the density structure evolution of a collapsing sphere. Therefore, CB81 could be undergoing gravitational collapse if the core does not contain additional supporting force acting against gravity. A survey looking for gas infalling motion using the HCN ( = -J 1 0) molecular line failed to detect gas inward motion toward CB81 because the emission from the core is too faint (Afonso et al. 1998). To confirm the line-of-sight gas motion toward CB81, more sensitive radio observations are needed. It is obvious that the CB81 core has additional support from magnetic fields. With that support, the core can be in kinematically nearly critical condition. Kandori et al. (2005) suggested that the majority of starless cores show a nearly critical Bonnor-Ebert solution. Because CB81 is starless and has a steep density structure, as indicated by x~13.0 max , it can serve as an important object for the study of star formation by illustrating conditions just before or after the onset of gravitational collapse.

A.3. Accuracy of Flat-fielding
Here, we discuss the accuracy of flat-fielding associated with calibration of the SIRPOL data. At present, SIRPOL data are calibrated using a flat image based on data acquired toward the twilight sky (i.e., twilight flat, see Figure 17). Twilight sky images with different sky count levels can be used to generate pairs of images, following which differential images generated from each pair are combined and normalized to obtain the final twilight flat image. During the exposures to acquire these twilight flat images, the angle of the half-waveplate is fixed at 67°.5.
There are two issues that may affect the accuracy of the SIRPOL flat image: (1) the sky is generally polarized, and this may affect the accuracy of SIRPOL flat, and (2) it is also uncertain whether the twilight flat image created based on the fixed waveplate angle (67°.5) is applicable to the calibration of other waveplate angle images.
To ascertain these issues, we employed sky images obtained during the observations of CB81. A total of 10 sets of sky images (that is, 10 exposures repeated 10 times) were used for each waveplate angle image, and Figure 18 shows the Stokes I image of the median sky in the H band. Here, all the sky images generated at the four waveplate angles are combined. The vertical gap in the sky count at x=512 pixels is termed the "reset anomaly pattern." The median sky count of the image is 1408 ADU. Since the sky count in the H band is due to the OH airglow lines, we believe that the sky was not polarized. The effect of the moon (illumination=66%, distance=74°o n 2017 June 15 and dark on 2017 June 24) can be used to set the upper limit for the analysis, as described below.
We evaluated the fraction of polarized light using the median sky images of I 00 , I 22.5 , I 45 , and I 67.5 , where I angle is the image acquired with the waveplate angles of 0°, 22°.5, 45°, and 67°.5. If the SIRPOL flat image can generate an artificial pattern on calibrated images, the effect obtained from the flat image (i.e., an artificial polarization pattern) will be observed in the H-band median sky image. To remove small-scale variations, the median sky images were smoothed using a median filter with a width of 9 pixels.
The white vectors in Figure 18 show the polarization vector distribution of the median sky. During calculations of these polarization vectors, a uniform component was subtracted. Specifically, = -Q Q Q median corr ( ) and = U corr -U U median( ) were used because the subtraction analysis of uniform vector components is employed in our basic procedure for detecting hourglass-shaped magnetic fields. The values of Q median( ) and U median( ) are 2.07 ADU and 2.26 ADU, respectively. This subtraction analysis suppresses the resulting polarization degree by∼50%. Figure 19 presents a histogram of the polarization degrees, for which the median and standard deviation values are 0.13%±0.08%. This value is less than the polarization degree threshold for the analysis of an hourglass-shaped field (  P 0.5% H ). Moreover, the polarization angle distribution does not resemble an hourglass. Therefore, we conclude that the use of the present SIRPOL flat image does not affect the validity of our conclusion.