Distortion of Magnetic Fields in Barnard 335

In this study, the detailed magnetic field structure of the dense protostellar core Barnard 335 (B335) was revealed, based on near-infrared polarimetric observations of background stars to measure dichroically polarized light produced by magnetically aligned dust grains in the core. Magnetic fields pervading B335 were mapped using 24 stars after subtracting unrelated ambient polarization components, revealing that they have an axisymmetrically distorted hourglass-shaped structure toward the protostellar core. On the basis of simple two- and three-dimensional magnetic field modeling, magnetic inclination angles in the plane-of-sky and line-of-sight directions were determined to be 90° ± 7° and 50° ± 10°, respectively. The total magnetic field strength of B335 was determined to be 30.2 ± 17.7 μG. The critical mass of B335, evaluated using both magnetic and thermal/turbulent support against collapse, was determined to be Mcr = 3.37 ± 0.94 M⊙, which is identical to the observed core mass of Mcore = 3.67 M⊙. We thus concluded that B335 started its contraction from a condition near equilibrium. We found a linear relationship in the polarization versus extinction diagram, up to AV ∼ 15 mag toward the stars with the greatest obscuration, which verified that our observations and analysis provide an accurate depiction of the core.


Introduction
Barnard 335 (B335) is one of the most well-defined isolated dense cores. The core harbors far-infrared (FIR) source IRAS 19347+0727, which is associated with a dense molecular gas envelope (e.g., Frerking et al. 1987;Menten et al. 1989;Saito et al. 1999;Kurono et al. 2013). The FIR source shows strong submillimeter (submm) emissions (Chandler et al. 1990) and is categorized as a Class 0 protostar based on its spectral energy distribution (Barsony 1994). The protostar is thought to be the driver of the bipolar outflow streaming in an east-west direction (e.g., Frerking & Langer 1982;Hirano et al. 1988Hirano et al. , 1992Yen et al. 2010). Molecular line spectra display an inverse P-Cygni profile toward the central region of the core (e.g., Zhou et al. 1993;Choi et al. 1995;Evans et al. 2005Evans et al. , 2015, which is evidence of gas moving inwardly toward the core center, as expected for a protostar. The core is seen in optical images as a dark, mostly opaque globule that obscures background starlight and is elongated toward the north-south direction (see Figure3 of Gålfalk & Olofsson 2007, see also Stutz et al. 2008 for the detection of flattened core).
Owing to its isolated geometry, simple shape, and rich stellar backdrop, the density structure of B335 was investigated with near-infrared (NIR) extinction measurements (Harvey et al. 2001).
The Bonnor-Ebert sphere model (Ebert 1955;Bonnor 1956), which describes a pressure-confined, self-gravitating isothermal gas sphere in hydrostatic equilibrium, was well suited for describing the density structure of the core using a dimensionless radius of ξ max =12.5±2.6 (Harvey et al. 2001). This ξ max value, which is a measure of the gas sphere's stability against gravitational collapse, clearly exceeds the stability criterion of 6.5 (Bonnor 1956), resulting in an unstable equilibrium solution. As shown by Kandori et al. (2005), the density profile of a collapsing gas sphere mimics a series of unstable Bonnor-Ebert solutions. The solution ξ max =12.5 is thus consistent with the existence of a protostar-forming inside B335. Subsequent density structure observations at radio wavelengths suggested a density profile slope of p∼−1.5 for the r µ r p relationship in the inner region and p∼−2 in the outer region (Harvey et al. 2003;Kurono et al. 2013). This is consistent with the shape of the density profile of an inside-out collapse solution (Shu 1977), as well as an unstable Bonnor-Ebert sphere with a rather flat inner region (plateau) and p∼−2 profile in the outer region. Since these profiles cannot be discerned using observational data (Harvey et al. 2001), we treat the shape of density structure of B335 as a Bonnor-Ebert sphere with ξ max =12.5 in this study.
On the basis of the NIR extinction data and the distance obtained by Tomita et al. (1979), the mass of B335 inside its radius (R=125″, obtained by the Bonnor-Ebert fitting, Harvey et al. 2001) was determined to be M core =14 M e (d 2 /250 pc; Harvey et al. 2001). Using the newly revised distance of d=105±15 pc (Olofsson & Olofsson 2009) determined using the photometric distance of foreground and background stars, the mass and radius of B335 are M core =3.67 M e and R=13,100 au (0.0656 pc), respectively. Zhou et al. (1990) reported the gas kinematic temperature T k =13 K and effective sound speed = C 0.23 s,eff km s −1 (i.e., turbulent velocity dispersion σ turb =0.085 km s −1 ) for B335, through the microturbulent radiative transfer modeling of H 2 CO 6 cm molecular line observations. Using the same parameter in theoretical modeling, gas infall spectra of various molecules were well explained (Zhou et al. 1993;Choi et al. 1995;Evans et al. 2005Evans et al. , 2015. Note that NH 3 observations produced temperatures that were slightly higher (T k =15 K toward the center of the core, from Kurono et al. 2013) and lower (T k =10-12 K for the bulk of gas, from Menten et al. 1984).
The magnetic field is the last physical property of B335 that needs to be explored, and it has been the subject of several studies. Wolf et al. (2003) conducted submm dust emission polarimetry toward the core to obtain a magnetic field strength of -+ 134 39 46 μG and mean magnetic field direction of θ mag =3°. This magnetic field direction is almost perpendicular to the outflow axis. Bertrang et al. (2014) conducted NIR (J s band) dust extinction polarimetry toward the core to obtain a magnetic field strength of 12-40 μG and mean magnetic field direction of θ mag =115°±6°. Optical polarimetry of the B335 region showed a similar mean magnetic field direction of θ mag =111°±4° (Vrba et al. 1986). Comparing these observations, mean magnetic field directions obtained at submm wavelengths are almost perpendicular to the direction obtained at optical or NIR wavelengths (note that recently Yen et al. 2019 reported that the core-scale magnetic field of B335 is in the east-west direction, based on the submm polarimetry using JCMT/POL-2). Furthermore, the estimates of magnetic field strength vary, although the uncertainty in both observations are large.
There are several possible reasons for these differences. First, the field of view of the submm dust emission polarimetry was small, less than 1′, so the submm results could be affected by the magnetic fields associated with the disk-like structure around the protostar. Furthermore, magnetic field strength was estimated using the data covering only the central region of the core, resulting in a large value, not an estimate of the mean magnetic field strength for the entire core. Second, the submm polarization map clearly contains a "polarization hole." The obtained submm polarization degree is small toward the center and gets larger toward the outer region. Therefore, it is anticorrelated with column density and does not reflect the overall magnetic field structure of the core. This may be due to the effect of optical depth and/or grain growth (Brauer et al. 2016). The same effect is seen in the subsequent dust emission polarimetry of B335 (Hull et al. 2014) and other dense cores (Matthews et al. 2009). Finally, though the NIR extinction polarimetry covered the entire core extent of ∼4′, the field of view was insufficient to map the outer ambient polarization distribution. Thus, the NIR data toward B335 are a superposition of the B335 polarization and the polarization arising from the ambient medium. This effect makes the observed polarization distribution of the core deviate from the true values. Note that optical polarimetry only traced ambient fields or the periphery of the core region. For these reasons, the magnetic field structure of B335 was not fully revealed based on the previous observations. This study conducted wide-field background star linear polarimetry at NIR wavelengths for B335. The plane-of-sky magnetic field structure was revealed using several tens of stars in and around core radius, and the ambient field component was subtracted. The total magnetic field strength was estimated based on the Davis-Chandrasekhar-Fermi method (Davis 1951;Chandrasekhar & Fermi 1953) and three-dimensional (3D) magnetic field modeling of the core. The magnetic field information was then used to determine the stability of B335 and describe the geometry among magnetic fields, outflow axis, and core elongation. The polarization versus extinction (P-A) relationship in B335 was constructed using corrections for (1) ambient polarization, (2) depolarization effect caused by a distorted magnetic field shape, and (3) line-of-sight magnetic inclination angle. The obtained linear P-A relationship verifies that the polarizations reported here reflect the overall magnetic field structure in the core.

Observations and Data Reduction
We observed B335 using the JHK s -simultaneous imaging camera SIRIUS (Nagayama et al. 2003) and its polarimetry mode SIRPOL (Kandori et al. 2006) on the Infrared Survey Facility (IRSF) 1.4-m telescope at the South African Astronomical Observatory (SAAO). IRSF/SIRPOL is one of the most useful instruments for NIR polarization surveys, providing deep-and wide-( ¢´¢ 7.7 7.7 with a scale of  0. 45 pixel −1 ) field polarization images. The uncertainty due to sky variation during exposures is typically 0.3% in polarization degree. The uncertainty of polarization angle due to the uncertainty in the determination of the polarization angle origin of the polarimeter is less than 3° (Kandori et al. 2006 and updates for Kusune et al. 2015).
Observations were conducted on the nights of 2017 June 30 and July 2 using the linear polarimetry mode of SIRPOL. Tensecond exposures at four half-waveplate angles (in the sequence 0°, 45°, 22°. 5, and 67°.5) were performed at 10 dithered positions (1 set). The total integration time was 2000 s (20 sets) per waveplate angle. The typical seeing during the observations was ∼  1. 35 (3.0 pixels) in the H band. We reduced the observed data in the same manner as described in Kandori et al. (2007) using the Interactive Data Language (IDL) software (flat-field correction with twilight flat frames, median sky subtraction, and frame combination after registration). Note that the accuracy of calibration based on twilight flat is discussed in the Appendix of Kandori et al. (2020b). Software aperture polarimetry was carried out for a number of sources in the field of view. The point sources having a peak intensity greater than 10σ above local sky background were detected on the Stokes I image. The number of detected source is 2089, 2181, and 1385 in the J, H, and K s bands, respectively. The limiting magnitudes are 18.5, 18.0, and 17.2 mag in the J, H, and K s bands, respectively. The local background was subtracted using the mean of an annulus around the source on the original image. This process was carried out for each position angle image (I 0 , I 45 , I 22.5 , and I 67.5 ). The aperture radius was the same as the full width at half maximum of stars (3.0, 3.0, and 2.8 pixels for the J, H, and K s bands, respectively), and the sky annulus was set to 10 pixels with a 5 pixel width. This relatively small aperture radius was used to suppress flux contamination from neighboring stars in the crowded field. Another method to avoid stellar contaminations is the use of the psf-fitting photometry on each waveplate angle image. In this case, however, the goodness of psf-fitting for each star varies on each waveplate angle image, and this can be the source of systematic errors in measured polarization signals. Thus, we do not use the psf-fitting method for the polarization measurements of B335. Note that though we use relatively small radius for the aperture polarimetry, our important conclusions in this paper do not change even if we change the size of the aperture.
The sources with a photometric error greater than 0. and θ=0.5 atan(U/Q). Since the polarization degree, P, is a positive quantity, the derived P values tend to be overestimated, especially for low S/N sources. We corrected the bias using d = -P P P db 2 2 (Wardle & Kronberg 1974). The H band observations provided polarizations for a total of 108 stars with P/δP5 in the field of view.
In this study, we used the results only for the H band, where NIR extinction by dust is less severe than in the J band and polarization efficiency is greater than that for the K s band. Note that the polarization vectors detected in the J band and K s band are roughly coincident with those in the H band. The correlation coefficients of P/δP>10 sources for polarization angle are 0.62 (J versus H) and 0.82 (K s versus H), and the coefficients for polarization degree are 0.79 (J versus H) and 0.50 (K s versus H). Figure 1 shows the observed polarization vector map of B335 in the H band. B335 is visible near the center of the image as a dark region that obscures the stars behind it. The white circle marks the core boundary (R=125″, Harvey et al. 2001). Inside the core radius, the bent structure of the polarization vectors can be seen, particularly in the northern part of the core. Outside of the core radius (i.e., the off-core region), polarization vectors generally flow from northwest to southeast. These polarization components can be regarded as off-core polarization, located along the same line of sight but unrelated to the B335 core.

Distortion of Magnetic Fields
Considering the direction of galactic longitude (PA∼30°in the equatorial coordinates at the position of B335), the off-core vectors deviate from the polarization direction of the galactic plane. As reported by Frerking et al. (1987), the B335 core is accompanied by diffuse envelope extending more than 20′ in east-west direction (see Figure1 of Frerking et al. 1987). The position angle of off-core vectors is roughly perpendicular to the elongation axis of the diffuse envelope, and may be associated with this structure.
The off-core vectors have , where x and y are the pixel coordinates, and A, B, and C are the parameters to be fitted. The estimated off-core vectors are shown in Figure 4. The regression vectors of off-core components were subtracted from the original vectors to isolate the polarization vectors associated with B335 ( Figure 5). Comparing Figures 1 and 5, polarization vectors toward B335 are slightly rotated by the effect of ambient subtraction. After subtraction, polarization degrees of off-core vectors were successfully suppressed toward 0% (Figure 2), and approximately random distribution of polarization angles was realized ( Figure 3). Some on-core polarization vectors increase in polarization degree after the subtraction of off-core vectors. This occurs for the on-core vectors with angles perpendicular to the off-core vectors, and is not an artificial effect. The increase of polarization degree for some on-core vectors is the result of the suppression of the depolarization effect by off-core vectors.
In Figure 5, the north part of off-core polarizations is not fully subtracted out and seems to be associated with the polarization pattern of the core region. As described above, these vectors may be associated with the diffuse envelope structure surrounding the B335 core.
In Figure 5, polarization vectors toward B335 are clearly oriented east-west, and a structure reflecting the bending of magnetic fields can be seen in the northern and (possibly) southern part of the core. The number of polarization vectors within the core radius of R125″ and with P H 1% is 24. This number is small, but sufficient to draft the magnetic field lines.

Parabolic Model
The most probable configuration of the magnetic field lines, estimated using a parabolic function and its rotation, is shown in Figure 6 (solid white lines). The field of view matches the diameter of B335 (250″), and the 24 polarization vectors having P H 1% are shown in the figure. The threshold of P H 1% was determined to avoid the systematic error caused by the subtraction of off-core polarizations. Note that most of the residual polarization vectors in the off-core region after the subtraction analysis are less than 1%, as shown in Figure 2. The coordinate origin of the parabolic function is fixed to the central position of the core/ protostar system measured on the Atacama Large Millimeter Array (ALMA) dust emission map (R.A.=19 h 37 m 00 89, decl.=+7°3 4′10 0, J2000, Evans et al. 2015). The best-fit parameters are θ mag =90°±7°and C=1.43 (±0.49)×10 −5 pixel −2 (=7.06×10 −5 arcsec −2 ) for the parabolic function = + y g gCx 2 , where g specifies magnetic field lines, θ mag denotes the position angle of the magnetic field direction (from north through east), and C determines the degree of curvature in the parabolic function. In the fitting procedure, the observational error for each star was taken into account in the calculations of ( 1 obs,i model 2 2 , where n is the number of stars, x and y are the star coordinates, θ obs and θ model denote the polarization angle from observations and the model, respectively, and δθ i is the observational error. The parabolic fitting seems reasonable, since the standard deviation of residual angles, θ res =θ obs −θ fit , is smaller for the parabolic function (dq =    21 . 54 0 . 99 res , Figure 7) than for the uniform field case (dq =    24 . 94 1 . 00 uni ). When dq res and dq uni are compared, the difference is statistically significant at the 3σ level, and the existence of the distorted field is most likely real.
B335 is the first protostellar core to be associated with axisymmetrically distorted hourglass-shaped magnetic fields, and is the third dense core to be associated with such fields following the results of FeSt 1-457 (Paper I) and B68 (Kandori et al. 2020c), which lack protostars. The obtained magnetic curvature value C=7.06×10 −5 arcsec −2 =6.4×10 −9 au −2 is similar to that obtained for FeSt 1-457 (C=5.14×10 −5 arcsec −2 =3.0× 10 −9 au −2 , Paper I) and B68 (C=1.09×10 −4 arcsec −2 = 7.0×10 −9 au −2 , Kandori et al. 2020c), suggesting the existence of similar mechanisms that create hourglass-shaped field distortions in dense cores/globules. The magnetic curvature of dense  cores can be produced with the drag of magnetic fields by materials concentrating toward center, and the degree of curvature can be determined by the moving distance of mass that drags magnetic fields. The similar curvature values among the dense cores indicate that the spatial scale of core formation is similar in these cores. Core formation scale, i.e., initial contraction radius R 0 , may be similar in dense cores/globules. Further observational studies to measure magnetic curvature of dense cores are needed.
The intrinsic dispersion, ( ) dq dq dq = int res 2 err 2 1 2 , estimated using the parabolic fitting, is 21°.34±1°.00 (0.3725 radian), where δθ err is the standard deviation of the observational error in polarization measurements. Note that the choice of a function form better than the function used here (i.e., parabolic form) can reduce the dispersion at residual angles (e.g., a mathematical model for hourglass fields by Ewertowski & Basu 2013;Myers et al. 2018).
If the magnetic field is assumed to be frozen in the medium, the intrinsic dispersion of the magnetic field direction, δθ int , can be attributed to an Alfvén wave perturbed by turbulence. The strength of the plane-of-sky magnetic field (B pos ) can be estimated from the relation ( ) pr s dq = B C 4 pos corr 1 2 turb int , where ρ and σ turb are the mean density of the core and turbulent velocity dispersion, respectively (Davis 1951;Chandrasekhar & Fermi 1953) and C corr is a correction factor suggested by theoretical studies. In the original formulation, C corr =1 (Davis 1951;Chandrasekhar & Fermi 1953), whereas in this study, a value of C corr =0.5 (Ostriker et al. 2001, see also, Heitsch et al. 2001Padoan et al. 2001;Heitsch 2005;Matsumoto et al. 2006) was adopted. Using the data for mean density (ρ=2.30 (±1.74)×10 −19 g cm −3 ) calculated from Harvey et al. (2001), turbulent velocity dispersion (σ turb =0.085 km s −1 ) from Zhou et al. (1990), and dq int derived here, a relatively weak magnetic field was obtained as a lower limit of total field strength (| | B ): B pos =19.4±11.4 μG. We compared the obtained plane-of-sky magnetic field strength with previous studies ( -+ 134 39 46 μG, Wolf et al. 2003;12-40 μG, Bertrang et al. 2014). Note that these previous studies used no theoretical correction factor in the Davis-Chandrasekhar-Fermi calculation. If we apply C corr =0.5 to their data, the values of the magnetic field strength are reduced to -+ 67 20 23 μG (Wolf et al. 2003) and 6-20 μG (Bertrang et al. 2014), respectively. The value by Bertrang et al. (2014) is roughly consistent with our result, and the value by Wolf et al. (2003) is larger than our result. Note that the submm polarimetry by Wolf et al. (2003) covered only the central region of the core, and the obtained magnetic field strength is not the mean value for the entire core.

3D Magnetic Field
For the 3D magnetic field modeling, we followed the procedure described in our previous papers (Kandori et al. 2017b, hereafter Paper II, see also Kandori et al. 2020a, hereafter Paper VI). The 3D version of the simple parabolic function employed in Paper I, ( ) j = + z r g g gCr , , 2 in the cylindrical coordinate ( ) j r z , , , was used for modeling the core magnetic fields, where g specifies the magnetic field line, C is the curvature of lines, and j is the azimuth angle (measured in the plane perpendicular to r). In this function, the shape of Figure 4. Estimated off-core polarization vectors superimposed on an H band intensity image. Note that these vectors were not obtained directly from observations. Off-core vectors estimated by fitting are plotted at the position of each star. The white circle marks the core boundary (radius of 125″, Harvey et al. 2001). The scale of the 5% polarization degree is shown above the image. magnetic field lines is axially symmetric around the r axis. Thus, the function ( ) j z r g , , has no dependence on the parameter j.
The 3D model was virtually observed after rotating in the lineof-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 parabolic model for various viewing angles ( g  -90 mag ) are shown in Figure S1 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 =12.5 (Harvey et al. 2001). Since we obtained a linear polarization-extinction relationship for B335 (see Section 3.5), it is reasonable to assume that the polarization arisen in a cell in the model core is proportional to the density in the cell. The model map can be compared with observations for each parameter set (see Figure  S1 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. Note that the resolution of the final model 3D map is enough for comparison with observations. The nearest neighbor distance of stars in observations is 44 pixels on the average, whereas the sampling grid width of the model is 9 pixels (≈4″).
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 by 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 . In Figure 8, the minimization point is γ mag =45°. It is clear that the high-γ mag region (especially γ mag  80°) is unlikely. Note that the large γ mag model approaches the pole-on magnetic field geometry and shows a radial polarization pattern that does not match observations.
To check consistency of the χ 2 analysis based on polarization angle, we made χ 2 calculations using both polarization angles and degrees. We determined the best model parameters for each γ mag by minimizing the difference in polarization angles, and we calculated where P obs and P model show the polarization degree from observations and the model, and dP i is the observational error in polarization degrees (Figure 9). In the procedure, the relationship between the model core column density and polarization degree was scaled to be consistent with observations. In Figure 9, the minimization point is γ mag =55°. It is clear that the low-γ mag region (especially γ mag  30°) is unlikely. Note that the small γ mag model approaches the edge-on magnetic field geometry and shows a hourglass-like field structure with weak depolarization pattern.
The distributions of c q 2 and c P 2 show relatively large scatter, especially for small or large γ mag . This can be due to the complicated patterns of hourglass-like fields projected onto the sky plane, and the relatively small number of data points (N = 24) of the B335 background stars. The 1σ error estimated at the minimum c q 2 and c P 2 points are 12°.5 and 4°. 7, respectively. These values are similar to the difference between χ 2 minimization points based on c q 2 and c P 2 (10°). We conclude that the most likely value for γ mag is 50°w ith uncertainty of 10°. The magnetic curvature obtained at γ mag =50°is C=1.73×10 −4 arcsec −2 . Note that the χ 2 values in Figures 8 and 9 are relatively large. This seems to come from the existence of polarization angle scattering mainly caused by Alfvén waves, which cannot be included in the observational error term in the calculation of χ 2 . Figure 6. Polarization vectors after subtraction of the off-core component. The field of view is 250″ or 0.131 pc at a distance of 105 pc, which is equal to the diameter of the core. The background image is the optical (H α ) image from Gålfalk & Olofsson (2007), and the polarization vectors are from Figure 5. The white lines indicate the direction of the magnetic field inferred from parabolic fitting. The scale of the 5% polarization degree is shown above the image.  . The best magnetic curvature parameter (C) is determined at each γ mag . γ mag =0°and 90°c orrespond to the edge-on and pole-on geometry in the magnetic axis. Figure 10 shows the best-fit 3D parabolic model with the observed polarization vectors. The direction of polarization vectors in the model generally agrees with observations at the same level compared with the results of 2D fitting, although there is some deviations in between model and observations particularly in the north part of the core. The standard deviation of the angular difference in plane-of-sky polarization angles between the 3D model and observations is 15°. 12, which is less than the fitting result in 2D, as well as in the uniform field case of δθ uni =24°.94. Figure 11 is the same observational data as Figure 10, but with the background image processed using the line integral convolution (LIC) technique (Cabral & Leedom 1993). We used the publicly available IDL code developed by Diego Falceta-Gonçalves. The direction of the LIC "texture" is parallel to the direction of magnetic fields, and the background image is based on the polarization degree of model core.

Magnetic Properties of the Core
Using the magnetic inclination angle θ inc of 50°±10°obtained as described in the previous section, we determined that the mean magnetic field strength of B335 is 30.2 17.7 μG. The ability of the magnetic field to support the core against gravity was investigated using the parameter obs critical , which represents the ratio of the observed mass-to-magnetic flux ratio to a critical value, ( ) p -2 G 1 2 1 , suggested by theory (Mestel & Spitzer 1956;Nakano & Nakamura 1978). We found that λ=3.24±1.52 (magnetically supercritical) and that the magnetic critical mass of the core is 1.13±0.69 M e , which is lower than the observed core mass of M core =3.67 M e .
The critical mass of B335, evaluated using both magnetic and thermal/turbulent support against collapse, M cr ; M mag +M BE (Mouschovias & Spitzer 1976;Tomisaka et al. 1988;McKee 1989), is 1.13+2.24=3.37±0.94 M e , where 2.24±0.64 M e is the Bonnor-Ebert mass calculated using the kinematic temperature of the core of 13 K, turbulent velocity dispersion of 0.085 km s −1 (Zhou et al. 1990), and external pressure P ext of 1.35(±0.77)×10 5 K cm −3 calculated from Harvey et al. (2001). Though the obtained M cr is slightly less than the core mass M core , we concluded that the stability of B335 is near the critical state, with M cr ∼M core .
Thus, the protostellar core B335 appears to have started its contraction from a condition near equilibrium between gravity and magnetic and thermal/turbulent support. We speculate that, in general, the (spontaneous) low-mass star formation in globules is initiated from the state close to the critical condition. This is consistent with the observations of subsonic infall motion in dense cores (Lee et al. 1999(Lee et al. , 2001. If the formation of stars were to start from a condition far out of equilibrium, the dense core would collapse suddenly with strong acceleration, and supersonic infall motion would be widely observed toward dense cores (see, e.g., Aikawa et al. 2005). The infalling gas motion in B335 was successfully modeled (e.g., Zhou et al. 1993;Choi et al. 1995;Evans et al. 2005Evans et al. , 2015 using the spherical inside-out collapse model (Shu 1977), in which the collapse starts from the end of quasistatic evolution. This is consistent with our conclusion that B335 is in a nearly equilibrium state.
We evaluated the stability parameter l º = M M cr core cr core mag BE , and found that λ cr is 1.09 for B335, which can be compared to λ cr of 0.94 for FeSt 1-457 (Papers I, II, VI) and 0.91 for B68 (Kandori et al. 2020c). Two starless cores and a protostellar core all show λ cr values close to unity, again supporting our conclusion that the initial condition of star formation in isolated globules may be in the condition very close to the critical state. If this is true, a slight decrease in turbulence (Nakano 1998) and/or ambipolar diffusion can initiate the onset of star formation. However, the formation mechanism of nearly critical dense cores/ globules remains unknown and is an open problem.
The relative importance of magnetic fields in supporting the core against gravity was investigated using the ratio of the thermal and turbulent energy to the magnetic energy, b º C V 3 s 2 A 2 and b s s where C s , σ turb , and V A denote the isothermal sound speed at 13 K, turbulent velocity dispersion, and Alfvén velocity, respectively. These ratios were found to be β≈4.39 and β turb ≈0.69, respectively. Thus, B335 is dominated by thermal support with a smaller contribution from static magnetic fields. Turbulence seems dissipated and makes a contribution comparable with or smaller than that of magnetic fields. These values were previously derived for FeSt 1-457 (β≈1.27 and β turb ≈0.12 with the application of the core's inclination angle; Paper I; Paper II) and B68 (β≈3.15 and β turb ≈0.70, Kandori et al. 2020c). Comparison of the results for all three dense cores indicates that they are mainly thermally supported with comparably smaller contributions from magnetic fields. Though the support from turbulence seems to be minor, this does not necessarily mean that turbulence is unimportant. Considering the nearly critical state of dense cores, even slight further dissipation of turbulence can initiate core contraction (e.g., Nakano 1998). Figure 12 shows the relationship among the core's elongation axis (θ elon ∼0°), outflow (θ out ∼90°), and magnetic field (θ mag =90°) superimposed on the map obtained by submm dust emission polarimetry (Wolf et al. 2003). Like FeSt 1-457 and B68, the orientation of the elongation of the optical obscuration of B335, as well as in the H 13 CO + (J=1−0) line map (Kurono et al. 2013), is clearly perpendicular to the Figure 9. χ 2 distribution for the polarization degree (c P 2 ). γ mag =0°and 90°c orrespond to the edge-on and pole-on geometry in the magnetic axis. Calculations of χ 2 in polarization degree were performed after determining the best magnetic curvature parameter (C) which minimizes χ 2 in the polarization angle. This calculation was carried out at each γ mag . Figure 11. Same as Figure 10, but the background image was made using the line integral convolution (LIC) technique (Cabral & Leedom 1993). The direction of the LIC "texture" is parallel to the direction of magnetic fields, and the background image is based on the polarization degree of model core. magnetic axis. This geometrical relationship is consistent with the picture of mass accretion along magnetic field lines suggested by theories (e.g., Galli & Shu 1993a, 1993b. This geometry is also consistent with the magneto-hydrostatic configuration (e.g., Tomisaka et al. 1988) as the starting condition for core contraction. The outflow axis is toward eastwest (Hirano et al. 1988) on the plane of sky and is nearly perpendicular to the line of sight, with an inclination angle of 87° (Stutz et al. 2008). Though the outflow axis and magnetic axis are aligned perfectly on the plane of sky, the magnetic inclination angle of γ mag =50°±10°means that a gap probably exists between the two axes in the line of sight.

Polarization-Extinction Relationship
As described in the previous paper (Kandori et al. 2018a, hereafter Paper III, see also Paper VI), the relationship between dust dichroic polarization (P) and extinction (A) in molecular clouds and their cores is important for (1) interpreting the obtained angle of interstellar polarizations, which is closely related to the direction of magnetic fields pervading the observed region, and (2) investigating dust alignment mechanisms in the dense environment. The former point is essential to ensure that our polarization observations can trace the dust alignment and magnetic fields deep inside B335. The latter point is important for the comparison of observations with dust alignment theories.
The observed polarizations toward B335 are a superposition of the polarization from the core and the ambient medium surrounding the core. The polarizations arising from the ambient off-core medium must be subtracted from the observed polarizations in order to extract the polarization associated with the core (Section 3.1). As described in Section 3.3, distorted magnetic fields surrounding the core can produce depolarization, Figure 12. Relationship among the axis of elongation of the core (θ elon ∼0°), the outflow axis (θ rot ∼90°, Hirano et al. 1988), and the magnetic field axis (θ mag =90°) superimposed on the submm dust emission polarimetry map from Wolf et al. (2003). particularly at the equatorial plane of the core. Furthermore, the line-of-sight inclination angle in magnetic axis weakens the observed plane-of-sky polarizations. These effects must be corrected in order to determine the true relationship between polarization and extinction in dense cores. Figure 13(a) shows the observed P H versus -H K s relationship with no correction. The stars with R125″ and P H /δP H 6 were used. We did not employ the threshold of P H 1% as used in the parabolic fit analysis in Section 3.2. Though several stars with P H 1% were included in the plot, this does not change the important conclusions. Note that the strongest polarization vector in Figure 1 was dropped in Figure 13(a) due to the large uncertainty in P H . In the figure, the polarization generally increases with increasing extinction up to H−K s ∼1.2 mag. The slope of the relationship is 3.30±0.14% mag −1 . The observed polarization is the superposition of the polarization from the core and ambient medium. Stars located outside of the core radius, R>125″, were used to estimate the off-core polarization vectors. After subtraction, a relatively linear P H versus -H K s relationship was obtained, as shown in Figure 13(b). The slope of the relationship is 4.67±0.17% mag −1 , which is close to the average interstellar polarization slope (Jones 1989). Note that, since the ambient vectors are roughly north-south and the B335 vectors are roughly east-west, the slope of the P-A relationship after subtraction is greater than the original value.
B335 is associated with inclined distorted magnetic fields, as discussed in Sections 3.1-3.3, which provides depolarization effects inside the core. Based on the known 3D magnetic field structure, a depolarization correction factor was estimated, as shown in Figure 14. This figure was created simply by dividing the polarization degree map of the distorted field model with edge-on geometry (γ mag =0°) by the inclined distorted field model (γ mag =50°). The correction factor map can thus simultaneously correct for both the depolarization and line-ofsight inclination angles. In Figure 14, the factors distributed around the equatorial plane (north-south direction) 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). Figure 13(c) shows the depolarization and inclination corrected P-A relationship obtained by dividing the Figure 13(b) relationship by the correction factor map ( Figure 14). The slope of the relationship is 12.73±0.43% mag −1 . The final corrected slope is larger than those obtained for FeSt 1-457 (6.60 ± 0.41% mag −1 ,  Paper III and Paper IV) and B68 (4.57 ± 0.11% mag −1 , Kandori et al. 2020c), and it is comparable to the statistically estimated upper limit of interstellar polarization efficiency (Jones 1989, -

H H K s
). The correlation coefficients obtained in these analysis are 0.80 (ambient polarization subtraction, Figure 13(b)) and 0.87 (ambient polarization, depolarization, and inclination correction, Figure 13(c)). These are comparable to the value of 0.85 for the original relationship. For B335, our analysis did not greatly improve the tightness of the P-A correlation, but it changed the slope value.
Comparison of Figures 13(a) and (c) shows that the corrections produced a dramatic change, particularly in slope. The final result is a relatively linear relationship between polarization and extinction in the range up to A V ∼15 mag. This linear relationship verifies that our NIR polarimetric observations trace the overall polarization (magnetic field) structure of B335. Figure 15 shows the P H /A V versus A V diagram. Reflecting the relatively linear slope in Figures 13(c), distributions of data points seems flat especially for A V  5 mag. To reject too high P H /A V data for very low A V , the data points with > P A 4 H V % mag −1 were not plotted. The dotted line shows the fitting of the data using the power law µ a P A A H V V , resulted in the α index of −0.22±0.22. The relatively shallow α index indicates that the polarization (and magnetic field) structure toward B335 is traced in our observations. The dotted-dashed line shows an observational upper limit by Jones (1989). The relation was calculated based on the equation , and the parameter η is set to 0.875 (Jones 1989). τ K denotes the optical depth in the K band. The polarization efficiency of B335 is quite high, as high as the interstellar upper limit value.

Summary and Conclusion
This study revealed the detailed magnetic field structure of the dense protostellar core B335 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 B335 were mapped using 24 stars, and axisymmetrically distorted hourglass-shaped magnetic fields were identified in a protostellar core. On the basis of simple 2D and 3D magnetic field modeling, magnetic inclination angles in the plane-of-sky and line-of-sight directions were determined as    90 7 and 50°±10°, respectively. We used the obtained magnetic inclination angle and the Davis-Chandrasekhar-Fermi  μG. The magnetic critical mass of the core, M mag =1.13± 0.69 M e , is less than the observed core mass, M core = 3.67 M e , suggesting a magnetically supercritical state with a ratio of observed mass-to-flux ratio to the critical value λ=3.24±1.52. The critical mass of B335, evaluated using both magnetic and thermal/turbulent support,  + M M cr mag M BE , is 1.13+2.24=3.37±0.94 M e , where 2.24± 0.64 M e is the Bonnor-Ebert mass calculated using the effective sound speed and external pressure. Though the obtained M cr is slightly less than the core mass M core , we concluded that the stability of B335 is in a condition near the critical state, with M cr ∼M core . Thus, B335 is considered to have started its contraction from the condition near the equilibrium. Further, we speculate that the (spontaneous) low-mass star formation in globules is generally initiated in the state close to the critical condition. The direction of planeof-sky magnetic fields of B335 (∼90°) is perpendicular to the core's elongation axis (∼0°), and is parallel to the outflow axis (∼90°). This geometrical relationship is consistent with the picture of mass accretion along magnetic field lines suggested by theories of isolated star formation. We found a linear relationship in the polarization versus extinction diagram, up to A V ∼15 mag toward the stars with greatest obscuration. The linear relationship indicates that the observed polarizations reflect the actual overall magnetic field structure of the core. However, this linear nature also raises questions about the alignment of dust grains in the cold, dense, and low radiation environment of B335. Further theoretical and observational studies are needed to explain the dust alignment.