Distortion of Magnetic Fields in BHR 71

The magnetic field structure of a star-forming Bok globule BHR 71 was determined based on near-infrared polarimetric observations of background stars. The magnetic field in BHR 71 was mapped from 25 stars. By using a simple 2D parabolic function, the plane-of-sky magnetic axis of the core was found to be θmag = 125° ± 11°. The plane-of-sky mean magnetic field strength of BHR 71 was found to be Bpos = 8.8–15.0 μG, indicating that the BHR 71 core is magnetically supercritical with λ = 1.44–2.43. Taking into account the effect of thermal/turbulent pressure and the plane-of-sky magnetic field component, the critical mass of BHR 71 was Mcr = 14.5–18.7 M⊙, which is consistent with the observed core mass of Mcore ≈ 14.7 M⊙. We conclude that BHR 71 is in a condition close to a kinematically critical state, and the magnetic field direction lies close to the plane of sky. Since BHR 71 is a star-forming core, a significantly subcritical condition (i.e., the magnetic field direction deviating from the plane of sky) is unlikely, and collapsed from a condition close to a kinematically critical state. There are two possible scenarios to explain the curved magnetic fields of BHR 71, one is an hourglass-like field structure due to mass accumulation and the other is the Inoue & Fukui mechanism, which proposes the interaction of the core with a shock wave to create curved magnetic fields wrapping around the core.


Introduction
Determining the magnetic fields in and around dense molecular cloud cores is important for many fields. For example, knowing the details of magnetic support against self-gravity can allow an assessment of the kinematical stability of dense cores, which is closely related to the initial physical conditions of star formation (e.g., Crutcher 2012; Kandori et al. 2017a, hereafter Paper I). If the line-of-sight magnetic inclination angle (γ mag ) is known, the total magnetic field strength can be obtained, allowing the magnetic criticality and kinematical stability of each core to be evaluated precisely. If a dense core is associated with hourglassshaped magnetic fields, specific polarization patterns can appear because of the depolarization effect of the inclined distorted magnetic field structure, and γ mag can be estimated through a simple model fitting (Kandori et al. 2017b, hereafter Paper II; see also Kataoka et al. 2012;Kandori et al. 2020d;Paper VI). To date, four dense cores with hourglass-shaped magnetic field structures have been identified (FeSt 1-457: Paper I, Barnard 68: Kandori et al. 2020a, Barnard 335: Kandori et al. 2020b, CB81 in Pipe Nebula: Kandori et al. 2020c). One core has been explained using the Inoue-Fukui mechanism (Inoue & Fukui 2013), showing a bending magnetic field structure (Corona Australis SL42 core: Kandori et al. 2020e).
Determining the depolarization pattern as well as γ mag is also important to correct the correlation between polarization and extinction , hereafter Paper III, see also Paper VI) and to compare the polarization data taken in different wavelengths (Kandori et al. 2018a, hereafter Paper V). Using the Davis-Chandrasekhar-Fermi technique in a modified form can help determine the radial distribution of the mass-tomagnetic flux ratio in a core (Kandori et al. 2018b, hereafter Paper IV).
Comparing observed hourglass-shaped field data and a theoretical hourglass field model assuming flux freezing (Paper VI; see also Mestel 1966;Ewertowski & Basu 2013;Myers et al. 2018) can help determine the core formation parameters, namely the initial density ρ 0 , initial radius R 0 (i.e., core formation radius), and initial magnetic field strength B 0 . These quantities are essential for discussing the formation and evolution of molecular clouds and cores.
As part of our magnetic field survey of dense cores, we investigated the BHR 71 core (Bourke et al. 1995a) located near the Coalsack region in the southern sky using near-infrared (NIR) polarimetry. BHR 71 is a Bok globule (Bok & Reilly 1947) and has several designations, such as Sa 136 (Sandqvist 1977), which is an extension of the catalog by Sandqvist & Lindroos (1976), and DC 297.7-2.8 (Hartley et al. 1986). In this paper, we use the name BHR 71 following Bourke et al. (1995a).
The distance to the BHR 71 core is often assumed to be the same as the distance to the Coalsack cloud, which is located in the vicinity of BHR 71. The distance to the Coalsack is generally taken to be 200 pc (Seidensticker & Schmidt-Kaler 1989;Straizys et al. 1994), but it may be as close as 150±30 pc (Corradi et al. 1997). Voirin et al. (2018) estimated a distance of 176±7 pc based on the Gaia data, and this is close to the values measured by Rodgers (1960) and Franco (1989) of 174±18 pc and 180±26 pc, respectively. We use Voirin's value of 176 pc in this paper.
NH 3 observations were made by Bourke et al. (1995b), and the kinematic temperature of BHR 71 was estimated to be 13 K. NH 3 (1,1) mapping observations with the Parkes 64 m radio telescope gave a full width at half maximum (FWHM) line width of ∼0.5 km s −1 over the core . The turbulent velocity dispersion σ turb of the core is thus 0.20 km s −1 .
The radius and mass of BHR 71 at a distance of 176 pc based on the Herschel spectrophotometry data were 0.28 pc and 14.7 M e , respectively (Yang et al. 2017). However, these quantities depend on the (molecular) probe used in the observations. For example, at the same distance, the radius and mass of the core measured in NH 3 (1,1) were 0.07 pc and 2.3 M e , respectively (Bourke et al. 1997). In this study, we use the value provided by Yang et al. (2017). Targon et al. (2011) conducted optical (R C band) polarimetry around the BHR 71 core. The obtained polarization angle (measured from north to east) was 102°, indicating that the direction of the global magnetic fields is roughly parallel to the elongation axis of BHR 71 (∼130°, see Figure 1, showing the elongation of the opaque region of BHR 71).
BHR 71 is a binary protostellar core, consisting of IRS1 and IRS2 separated by ∼15″ or ∼2600 au at a distance of 176 pc (Bourke et al. 1997;Bourke 2001). IRS1 is a Class 0 protostar (Bourke 2001;Green et al. 2013), but IRS2 may be more evolved (Bourke 2001). These protostars have misaligned outflows (Bourke 2001;Parise et al. 2006). Furthermore, Tobin et al. (2019) reported that the rotation of these protostars measured in C 18 O ( = -J 1 0) on scales of <1000 au around each source is in the opposite direction. These results may be evidence these protostars formed by turbulent fragmentation , and is consistent with a theoretical simulation of the formation of binary/multiple sources in turbulent environments (Offner et al. 2016).
Recently, Hull et al. (2019) reported the high-resolution dust emission polarimetry with ALMA toward IRS1 and IRS2, in which IRS1 appears to be an hourglass-shaped magnetic fields, and IRS2 exhibits a magnetic field that has been affected by its bipolar outflow.
In the present study, wide-field background star polarimetry at NIR wavelengths was conducted for BHR 71. The planeof-sky magnetic field structure was revealed using stars in and around the core. The plane-of-sky magnetic field strength was estimated based on the Davis-Chandrasekhar-Fermi method (Davis 1951;Chandrasekhar & Fermi 1953). Using the resulting magnetic field information, the kinematical stability and the origin of the distorted magnetic field structure in BHR 71 are discussed.

Observations
The observations of BHR 71 were conducted using the JHK s -simultaneous near-infrared polarimeter SIRPOL (Kandori et al. 2006; see also Nagayama et al. 2003 for camera) on the IRSF 1.4 m telescope in the South African Astronomical Observatory (SAAO). IRSF/SIRPOL provides a large field of view ( ¢´¢ 7.7 7.7 with a scale of 0 45 pixel −1 ), so that nearby dense cloud cores can be covered with a single pointing. Since SIRPOL is a single beam polarimeter, sky changes during the exposures at different half-waveplate angles. Typical uncertainty of polarization degree due to sky variation in a photometric night is about 0.3%. The instrumental polarization over the field of view is less than 0.3% (Kandori et al. 2020c). The accuracy of the zeropoint angle of the polarimeter is less than 3° (Kandori et al. 2006;Kusune et al. 2015). A polarimetric standard star, RCrA#88 Whittet et al. 1992), was observed on 2017 July 13, and we obtained the results (P H =2.82%±0.07% and θ H =91°.9±0°.9) which are consistent with the values in the literature.
We observed BHR 71 on the nights of 2017 June 13-16. Exposures of 10 s were performed at four half-waveplate angles (in the sequence of 0°, 45°, 22°.5, and 67°.5) at 10 dithered positions (one set). The total integration time was 4500 s (30 sets) per waveplate angle. Sky frames were observed before or after the observations of object frames.
Data were reduced using the Interactive Data Language software. The reduction procedures include dark subtraction, flat-field correction, median sky subtraction, and frame combining after registration (see, e.g., Kandori et al. 2007).
Point sources having peak intensities greater than 10σ above local sky background were cataloged on the Stokes I images. Aperture polarimetry was performed for these sources on the images take at four waveplate angles (I 0°, I 45°, I 22°. 5 , and I 67°. 5 ). Sources with photometric error greater than 0.1 mag were removed from the list, and 3442 sources were detected in the H band. The aperture radius was the same as the FWHM of stars (3.5 pixel). The sky radius and the width of sky annulus were set to 10 and 5 pixels, respectively. The limiting magnitudes were 18.0 mag in the H band. The relatively small aperture size was used to suppress the flux contamination from neighboring bright stars. A technique suitable for the photometry in a crowded field is the point-spread function fitting photometry. We did not employ this technique, because the different goodness of fit for each star on different waveplate angle images can cause systematic errors in polarimetric measurements. The Stokes parameters for each star were derived with the relationship 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 equations = P + Q U I 2 2 and θ=0.5 atan (U/Q), respectively. Because P is a positive quantity, the resulting P values tend to be overestimated especially for low S/N sources. We corrected for the bias in P measurements using the equation & Kronberg 1974). In this paper, we focus on the results in the H band, in which dust extinction is less severe in the J band and polarization efficiency is greater than in the K s band. Figure 1 shows the observed polarization vectors (yellow vectors) on the Stokes I image in the H band. To draw the vectors, we used 140 stars with P H /δP H 4. The BHR 71 core appears as a dark obscuration around the center of the image, and the coneshaped illumination of the outflow cavity wall by the protostar IRS1 can be seen. The shape of the dust obscuration of BHR 71 is elongated with the long axis toward the direction of ∼130°. There are many relatively small polarization vectors around the southeast corner and the southwest corner of the image, flowing from northwest to southeast. Reflecting this, the histogram of θ H for all sources has a peak fitted with a Gaussian function with a peak angle of 111°±4° (Figure 2), which is consistent with the optical polarimetry results (Targon et al. 2011).

Distortion of Magnetic Fields
The most striking feature in Figure 1 is the steeply curved polarization vector pattern mostly located in the northern side of the BHR 71 core. Strong polarization was associated with the curved field structure. If we isolate the strong polarization vectors in the field of view, only the observation data associated with the curved field component remain. Furthermore, they are all inside the radius of BHR 71. These vectors appear to be different from the θ H =111°-component, and are likely associated with the BHR 71 core.
To determine the basic physical properties of the BHR 71 core, we used the FIR data taken by the Herschel satellite. Sadavoy et al. (2018) provided an intensity-corrected Herschel map of BHR 71, and a map at an optical depth of 353 GHz (t 353GHz ) was included in their products. The map has 28×28 pixels with a pixel scale of 14″ pixel −1 and resolution of 36 3. The molecular hydrogen column density can be obtained from is the molecular weight per hydrogen molecule (assuming cosmic abundance ratios), and m H is the hydrogen atom mass. Using the above equation we converted the t 353GHz map to a molecular hydrogen column density map (see background image in Figures 9 and 12). The centroid center of the core was determined to be (R.A., decl.)=(12 01 36. 13, h m s - ¢  65 08 46. 7, J2000). The elongation of the core was determined based on an ellipse fit to the region ( ) > N H 10 2 22 cm −2 , giving ≈133°as the direction of its major axis.
To analyze the NIR polarimetry data further, we need to separate the global θ H =111°-component (ambient component) and the polarization component arisen from the BHR 71 core. While one simple way is to just ignore the weak polarization vectors and to use strong polarization vectors, we decided to use the method of ambient vector field subtraction following the methods described in Paper I. Note that the former method can provide results that are consistent with the latter method.
As shown in Figure 2, the ambient component of the polarization vectors is relatively well aligned. Spatial linear plane fitting of the component was thus conducted. To exclude the positions passing through the dense core region, the stars located within 300 pixels ( ¢¢ » 135 24000 au) of the core center were not used in the fitting. The distributions of the Q/I and U/I values were independently modeled as ( , where x and y are the pixel coordinates and A, B, and C are the fitting parameters. Figures 3 and 4 show the histograms of P H or θ H of  stars before (solid line) and after (dotted-dashed line) the subtraction. The results did not change the polarization vectors dramatically. However, after subtraction, the number of low P H vectors increased in Figure 3, and θ H became roughly randomly distributed in Figure 4. Therefore, the subtraction analysis works somewhat satisfactorily. Figure 5 shows the distribution of estimated ambient polarization vectors. Most of the vectors are as strong as P H ∼1%, distributed mainly in the southern half of the image. Figure 6 shows the polarization vectors after the subtraction of ambient vectors. The small polarization vectors distributed around the southeast corner and the southwest corner of the image disappeared, and the curved polarization vectors associated with the dense part of the core were reasonably well extracted. Figure 7 shows the relationship between H and K s color (i.e., A V ) and polarization degree toward the background stars with P H taken after subtraction of the ambient component. Stars with P H /δP H 4 were plotted. Though the scatter in the diagram is not small, the distribution of data points is positive, showing that the observations trace the dust polarization (i.e., magnetic field direction) inside a cold and dense environment. The slope of the relationship P H /(H−K s ) was found to be 2.52±0.04% mag −1 .

Parabolic Model
Since the global direction of the magnetic fields ( Figure 2) and the curved field associated with BHR 71 (Figure 6) are clear, we tried to fit the magnetic field structure of the core using a simple parabolic function, = + y g gCx 2 , where g specifies each magnetic field line and C determines the curvature of the parabola. We used the function in a 90°rotated form, so that the x-axis is toward 0°in position angle. We calculated the parabolic function for various C, various spatial positions, and various rotation angles to find the parameter that minimizes where n is the number of stars, x and y are the coordinates of the i-th star, q obs and θ model are the polarization angles from the observations and the model, and δθ i is the observational error. We used 25 stars with P H /δP H 4 and P H 1.5%. The latter threshold was set so as to exclude the stars affected by the ambient subtraction analysis. In Figure 3, most of the ambient sources show P H  1.5%. Figure 8 shows the result of the parabolic fit (white lines). Figure 9 is the same as Figure 8 but the column density map based on the Herschel data was used as a background image. The prominent curved magnetic fields located in the northeast part of the image were relatively well fitted with the parabolic function, whereas there is a small number of vectors in the southwest part of the image. In Figures 8 and 9, the red plus   . Figure 10 shows a histogram of θ res for all the stars. We subtracted the effect from the observational error using   H are plotted. The white lines indicate the direction of the magnetic field inferred from parabolic fitting. The red plus sign shows the core center determined on the column density map based on Herschel. The blue plus sign shows the center of the curved magnetic fields (i.e., parabolic function). The scale bar above the image indicates 5% polarization. Figure 9. Same as Figure 8, but with the column density map used as the background image. Figure 10. Histogram of residuals for the observed polarization angles after subtraction of the angles obtained by parabolic fitting (θ res ). All stars were used to make the histogram. (Davis 1951;Chandrasekhar & Fermi 1953)  where ρ is the mean density (1.2×10 −20 g cm −3 from Yang et al. 2017), σ turb is the turbulent velocity dispersion (0.20 km s −1 from Tobin et al. 2019), δθ int is 14°.6-24°.6 (0.255-0.429 radian, this study), and C corr =0.5 is the correction factor from theory (Ostriker et al. 2001, see also Heitsch et al. 2001Heitsch et al. , 2005Padoan et al. 2001;Matsumoto et al. 2006). From the above equation, we obtained a relatively weak magnetic field strength of B pos =8.8-15.0 μG. Note that the estimated field strength is the averaged value for the whole core.

Magnetic Properties of the Core and the Origin of the Curved Magnetic Fields
In this section, we first discuss the plane-of-sky magnetic properties of BHR 71, and then discuss the origin of the curved magnetic fields. Note that we did not conduct a threedimensional (3D) analysis of the curved field, because the number of stars is not sufficient.
The magnetic support of BHR 71 against gravity can be investigated using the parameter obs critical , which represents the ratio of the observed mass-to-magnetic flux to a critical value, ( ) p -G 2 1 2 1 , suggested by theory (Mestel & Spitzer 1956;Nakano & Nakamura 1978). We determined a value of -l = 1.44 2.43, indicating a magnetically supercritical state. The magnetic critical mass of the core, M mag =6.0-10.2 M e , is lower than the observed core mass of M core ≈14.7 M e (Yang et al. 2017). Note that this does not necessarily imply a dynamical collapse of BHR 71, because there are additional thermal and turbulent supports.
The critical mass of BHR 71, taking into account both magnetic and thermal/turbulent support effects, can be written as  M cr + M M mag BE (Tomisaka et al. 1988;McKee 1989), where M BE is the Bonnor-Ebert mass (Ebert 1955;Bonnor 1956). We obtained -= M 14.5 18.7 cr M e with M BE 8.5 M e . The Bonnor-Ebert mass was estimated using a kinematic temperature T kin of 13 K (Bourke et al. 1995b), a turbulent velocity dispersion of 0.2 km s −1 , and an assumed external pressure of 6×10 4 K cm −3 . The assumed external pressure is equivalent to T eff ×ρ, where ρ is the mean density of the core and T eff is the sum of the kinematic temperature T kin and the turbulence equivalent temperature = T turb m s = k 10.9 p turb 2 K. The mean molecular weight per free particle μ p is set to 2.33 and k is the Boltzmann constant. The obtained critical mass of -= M 14.5 18.7 cr M e is not far from the core mass M core ≈14.7 M e , suggesting a nearly critical state.
Since BHR 71 is a star-forming core, a stable subcritical state (i.e., large M cr ) for this core is unlikely. We found that M cr ∼M core , based on a plane-of-sky magnetic field strength B pos . This means that the line-of-sight magnetic field inclination angle cannot deviate from the plane of sky. If it did deviate, we would observe a large total field strength, which would lead to a large M cr . Thus, we conclude that the magnetic fields pervading BHR 71 lie near the sky plane, and the core started its collapse from a state near the kinematically critical condition.
Finally we propose a possible scenario for the origin of the curved magnetic fields in BHR 71. As shown in Figures 11 and  12, the magnetic fields pervading BHR 71 can be explained by a single curved field structure. The field lines (white lines) in Figures 11 and 12 are the same as those in Figures 8 and 9, but the southern field components are removed. The best model to explain this is the Inoue & Fukui (2013) mechanism, describing a shock wave propagating from the southwest in the direction of the position angle of 35°(perpendicular to θ mag ), sweeping the magnetic fields, such that the magnetic fields wrap around the core to create the curved magnetic field structure. Although we cannot specify the origin of the shock, BHR 71 is located in the Lower Centaurus-Crux association as a subgroup of the Scorpius-Centaurus association. The curved magnetic fields may be a remnant of the interaction of the core with the past shock wave. Note that the hourglass shape in Figures 8 and 9 is still possible, and thus we have two scenarios to explain the origin of the curved magnetic fields in BHR 71. To determine which of these the scenarios applies to, large-scale radio   Figure 11, but with the column density map used as the background image.