JCMT POL-2 and ALMA polarimetric observations of 6000-100 au scales in the protostar B335: linking magnetic field and gas kinematics in observations and MHD simulations
Hsi-Wei Yen, Bo Zhao, I-Ta Hsieh, Patrick Koch, Ruben Krasnopolsky, Chin-Fei Lee, Zhi-Yun Li, Sheng-Yuan Liu, Nagayoshi Ohashi, Shigehisa Takakuwa, Ya-Wen Tang
DD raft version J anuary
3, 2019Typeset using L A TEX twocolumn style in AASTeX62
JCMT POL-2 and ALMA polarimetric observations of 6000–100 au scales in the protostar B335: linking magnetic field and gaskinematics in observations and MHD simulations H si -W ei Y en , B o Z hao , I-T a H sieh , P atrick K och , R uben K rasnopolsky , C hin -F ei L ee , Z hi -Y un L i , S heng -Y uan L iu , N agayoshi O hashi , S higehisa T akakuwa , and Y a -W en T ang Academia Sinica Institute of Astronomy and Astrophysics, P.O. Box 23-141, Taipei 10617, Taiwan European Southern Observatory (ESO), Karl-Schwarzschild-Str. 2, D-85748 Garching, Germany Max-Planck-Institut f¨ur extraterrestrische Physik (MPE), Garching, Germany, 85748 Astronomy Department, University of Virginia, Charlottesville, VA 22904, USA Subaru Telescope, National Astronomical Observatory of Japan, 650 North A’ohoku Place, Hilo, HI, 96720, USA Department of Physics and Astronomy, Graduate School of Science and Engineering, Kagoshima University, 1-21-35 Korimoto, Kagoshima, Kagoshima890-0065, Japan
ABSTRACTWe present our analysis of the magnetic field structures from 6000 au to 100 au scales in the Class 0 protostarB335 inferred from our JCMT POL-2 observations and the ALMA archival polarimetric data. To interpret theobservational results, we perform a series of (non-)ideal MHD simulations of the collapse of a rotating non-turbulent dense core, whose initial conditions are adopted to be the same as observed in B335, and generatesynthetic polarization maps. The comparison of our JCMT and simulation results suggests that the magneticfield on a 6000 au scale in B335 is pinched and well aligned with the bipolar outflow along the east–westdirection. Among all our simulations, the ALMA polarimetric results are best explained with weak magneticfield models having an initial mass-to-flux ratio of 9.6. However, we find that with the weak magnetic field,the rotational velocity on a 100 au scale and the disk size in our simulations are larger than the observationalestimates by a factor of several. An independent comparison of our simulations and the gas kinematics inB335 observed with the SMA and ALMA favors strong magnetic field models with an initial mass-to-flux ratiosmaller than 4.8. We discuss two possibilities resulting in the di ff erent magnetic field strengths inferred fromthe polarimetric and molecular-line observations, (1) overestimated rotational-to-gravitational energy in B335and (2) additional contributions in the polarized intensity due to scattering on a 100 au scale. Keywords:
Stars: formation - ISM: kinematics and dynamics - ISM: individual objects (B335) - ISM: magneticfields INTRODUCTIONStars form via gravitational collapse of dense cores (Shu etal. 1987), which are magnetized (Crutcher 2012). During thecollapse, the magnetic field lines are expected to be draggedinward by collapsing material, and the magnetic flux in theinner envelope around a central protostar increases (e.g., Li& Shu 1996; Galli et al. 2006). As the magnetic field strengthincreases in the inner envelope, the magnetic field can slowdown the infalling and rotational motions of the collapsingmaterial more e ffi ciently (e.g., Allen et al. 2003), if the fieldand matter remain well coupled. As a result, the collapsingmaterial is expected to infall toward the center at a velocity Corresponding author: Hsi-Wei [email protected] slower than the free-fall velocity, and its angular momentumis transferred outward, leading to suppression of formationand growth of a rotationally-supported disk (e.g., Mellon &Li 2008). Signs of infalling motion slower than free fall andremoval of angular momentum of collapsing material havebeen observed in several Class 0 and I protostars, such asL1527 (Ohashi et al. 2014), B335 (Yen et al. 2015b), TMC-1A (Aso et al. 2015), L1551 IRS 5 (Chou et al. 2014), andHH 111 (Lee et al. 2016). These results suggest that the mag-netic field could play an importance role in the dynamics incollapsing dense cores. Nevertheless, it remains unclear asto how e ffi ciently the magnetic field a ff ects the star-formingprocess.Polarized thermal dust continuum emission at (sub-)millimeter wavelengths can be adopted to trace magneticfield structures on scales of hundreds to thousands of au indense cores (Crutcher 2012), where dust grains are expected a r X i v : . [ a s t r o - ph . S R ] J a n to preferentially align their long axis perpendicular to themagnetic field (Lazarian & Hoang 2007). Thus, the magneticfield orientation can be inferred by rotating the polarizationorientation by 90 ◦ . In addition, in young protostellar sources,the sizes of dust grains on these scales are likely smaller than100 µ m (e.g., Li et al. 2017), so that dust scattering unlikelyinduces any significant polarized intensity (Kataoka et al.2015; Yang et al. 2016, 2017). With polarimetric observa-tions at (sub-)millimeter wavelengths, magnetic field linesbeing dragged to form an hour-glass morphology by col-lapse or being wrapped by rotational motion have been seenin protostellar envelopes around several protostars (Girartet al. 2006; Attard et al. 2009; Rao et al. 2009; Hull et al.2014; Davidson et al. 2014; Cox et al. 2018; Sadavoy et al.2018; Lee et al. 2018; Maury et al. 2018; Kwon et al. 2018)as well as in high-mass star-forming regions (Girart et al.2009; Qiu et al. 2013, 2014; Ward-Thompson et al. 2017;Pattle et al. 2017). These results suggest an interplay be-tween the magnetic field and the gas motions. Therefore,linking observational results of magnetic field structures andgas kinematics could shed light on the role of the magneticfield in the dynamics of collapsing cores.B335 is an isolated Bok globule with an embedded Class0 protostar at a distance of 100 pc (Keene et al. 1980, 1983;Stutz et al. 2008; Olofsson & Olofsson 2009). The size ofthe dense core in B335 observed at millimeter wavelengths is ∼ / submillimeter Array (ALMA) at an angular resolution of0 . (cid:48)(cid:48) ff ects of the magnetic fieldon the gas kinematics in B335. In addition, ALMA obser-vations in the C O and H CO + lines show no detectabledi ff erence in the infalling velocities of neutral and ionizedgas on a 100 au scale with a constraint on the upper limit ofthe ambipolar drift velocity of 0.3 km / s, suggesting that themagnetic field likely remains well coupled with the matter inthe inner envelope in B335 (Yen et al. 2018). The magneticfield structures on a 1000 au scale in B335 also show signsof being dragged toward the center and become pinched, asinferred from the ALMA polarimetric observations (Maury et al. 2018). Therefore, B335 is an excellent target to investi-gate the interplay between the magnetic field and gas motionsand the e ff ects of the magnetic field on the dynamics in col-lapsing dense cores.Theoretical simulations show that the importance of themagnetic field on the dynamics in collapsing dense cores isclosely related to the magnetic field strength, coupling be-tween the magnetic field and matter, and alignment betweenthe magnetic field and rotational axis (Mellon & Li 2008,2009; Li et al. 2011, 2013; Dapp et al. 2012; Joos et al.2012; Padovani et al. 2014; Tsukamoto et al. 2015; Zhao etal. 2016, 2018a; Masson et al. 2016). In the present work,we study these three physical parameters by linking the in-formation from the polarimetric and molecular-line observa-tions of B335 and by comparing them with theoretical sim-ulations. We conducted polarimetric observations at sub-millimeter wavelength with the James Clerk Maxwell Tele-scope (JCMT) to map the magnetic field structures on a scaleof thousands of au in B335. We additionally obtained theALMA polarimetric data at millimeter wavelength from thepublic archive in order to study the magnetic field structureson scales from the dense core to the inner envelope in B335.The details of the observations and the observational resultsare presented in Section 2. To analyze the observed magneticfield structures, we performed a series of (non-)ideal mag-netohydrodynamics (MHD) simulations and generated syn-thetic polarization maps to compare with the observations.The simulation setup and results are described in Section 3.In Section 4, we compare the magnetic field structures fromthe observations and the simulations. With these polarimet-ric data and our MHD simulations as well as the observa-tional results of the gas kinematics from our previous studieswith the SMT, SMA, and ALMA, we discuss the magneticfield strength and the coupling between the field and matterin B335 in Section 5. POLARIMETRIC OBSERVATIONS2.1.
JCMT POL-2 observations
The observations with the JCMT were conducted on April18 and May 5 and 12 in 2017. During the observations, the225 GHz opacity ranged from 0.03 to 0.08. Polarized contin-uum emission was observed at 850 µ m and 450 µ m simulta-neously with the polarimeter POL-2 (Friberg et al. 2016) andthe continuum receiver SCUBA-2 (Holland et al. 2013). ThePOL-2 observing mode at 450 µ m is not fully commissioned.In this paper, we present the results at 850 µ m. The angularresolution of JCMT at 850 µ m is 14 (cid:48)(cid:48) . The total on-sourceobserving time was 6.8 hours. The Daisy observing modewas adopted. With this observing mode, the exposure timeof the central region within a radius of 3 (cid:48) is above 80% of thetotal observing time. Table 1.
JCMT POL-2 detectionsRA Dec PP ∆ PP PA ∆ PA I p ∆ I p (%) (%) ( ◦ ) ( ◦ ) (mJy / beam) (mJy / beam)19:37:02.28 + + −
19 8 4.8 1.519:37:00.66 + + + − + + ote —PP, PA, and I p present the polarization percentage, the position angle of the polarization orientation, and the polarized intensity,respectively. ∆ PP, ∆ PA, and ∆ I p are their uncertainties. The data were reduced with the software
Starlink (Currieet al. 2014) and the task pol2map . The data were first re-duced with the default pixel size of 4 (cid:48)(cid:48) and the procedure makemap to obtain initial Stokes I maps. These initial Stokes I maps were adopted as the reference to correct the instru-mental polarization and to generate masks in the subsequentdata reduction process. Then, the final Stokes IQU mapswere generated with the procedure skyloop , and the polar-ized intensity was debiased in this process. We binned upthe final Stokes
IQU maps from the default pixel size of 4 (cid:48)(cid:48) to 12 (cid:48)(cid:48) to improve the signal-to-noise ratio (S / N). The finalpixel size is slightly smaller than the angular resolution of14 (cid:48)(cid:48) . The achieved noise level is 1.6 mJy Beam − in Stokes I and 1.5 mJy Beam − in the polarized intensity ( I p ) withthe pixel size of 12 (cid:48)(cid:48) . Our detection criteria of the polarizedemission are I p over its uncertainty ( ∆ I p ), I p / ∆ I p >
3, and theS / N in Stokes I larger than 5. This leads to seven detections(Table 1). The uncertainties in the polarization orientations,which are ∆ I p / I p (e.g., Hull et al. 2014), range from 6 ◦ to8 ◦ . 2.2. ALMA observations
The ALMA polarimetric data at 1.3 mm analyzed herewere retrieved from the archive (project code: 2013.1.01380.S).The details of the ALMA observations were described inMaury et al. (2018). We reduced the data with the calibrationscript provided in the archive using the software CommonAstronomy Software Applications (CASA; McMullin et al.2007) of version 4.7.0, and we further performed self cali-bration on the phase. The images were generated with Briggsweighting with a robust parameter of 0.5. The achieved angu-lar resolution is 0 . (cid:48)(cid:48) × . (cid:48)(cid:48)
47 ( ∼
70 au ×
50 au). The achievednoise level is 50 µ Jy in Stokes I and 15 µ Jy in both Stokes Q and U . We debiased I p with I p = (cid:112) Q + U − σ Q , U ,where σ Q , U is the noise level in Stokes Q and U (Wardle& Kronberg 1974; Simmons & Stewart 1985). We binnedup the Stokes Q and U maps to have a pixel size of 0 . (cid:48)(cid:48) / N higher than 3 in Stokes I and I p / ∆ I p > (cid:46) ◦ . 2.3. Observational results
Figure 1a presents the maps of the Stokes I and polarizedintensity at 850 µ m of B335 observed with the JCMT. Thediameter of the observed dense core is ∼ (cid:48)(cid:48) re-gion. All the other detections are at outer radii larger than20 (cid:48)(cid:48) . Figure 1b shows the polarization percentage as a func-tion of Stokes I intensity. The polarization percentage clearlydecreases with increasing Stokes I intensity, as reported inseveral dense cores (e.g., Wolf et al. 2003). The polarizationprecentage close to the center is 0.7%. This is much lowerthan the median polarization percentage of 10% on a scaleof 1000 au (10 (cid:48)(cid:48) ), which is comparable to the JCMT beamsize of 14 (cid:48)(cid:48) , detected with the SMA at an angular resolutionof 5 (cid:48)(cid:48) (Galametz et al. 2018). As the ALMA observations atan angular resolution of 0 . (cid:48)(cid:48) ff erent polarizationorientations within the larger JCMT beam.Figure 2a presents the orientations of the magnetic fieldsegments inferred by rotating the polarization orientationsby 90 ◦ . The inferred magnetic field orientations are orga-nized and are along the east–west direction, except for onesegment in the northwest with an orientation clearly deviatedfrom the others. Figure 2c presents the number distributionof the inferred magnetic field orientations from the JCMT po-larization detections. The mean orientation of the magneticfield segments is 99 ◦ with a standard deviation of 20 ◦ . Ex-cluding the segment in the northwest with a position angle (a) (b) Polarization segments P o l a r i z e d i n t e n s i t y ( m J y / B e a m ) Stokes I (mJy/Beam) P o l a r i z a t i o n p e r c e n t a g e ( % ) Figure 1.
JCMT POL-2 results of B335 at 850 µ m. (a) Stokes I map(contours) overlaid on polarized intensity map (grey scale). Seg-ments show the polarization orientations detected at a level above3 σ , and their lengths are proportional to the polarization percentage.The center of the map is 19 h m . s + d m . s
8. (b) Polarizationpercentage as a function of Stokes I intensity. (PA) of 137 ◦ , the mean orientation becomes 92 ◦ with a stan-dard deviation of 12 ◦ , comparable to the 2 σ uncertainty inthe polarization angle. The mean orientation of the magneticfield is consistent with the direction of the outflow having aPA of 90 ◦ –99 ◦ observed in the CO lines in B335 (Hirano etal. 1988; Yen et al. 2010; Hull et al. 2014).We note that the inferred magnetic field orientations fromour POL-2 observations are di ff erent from the observationalresults obtained with SCUPOL (Wolf et al. 2003; Matthewset al. 2009). In the SCUPOL results, the magnetic field orien-tations have a mean PA of 3 ◦ and a standard deviation of 36 ◦ ,and are perpendicular to the direction of the outflow (Wolfet al. 2003). Only two SCUPOL detections in the north-west show orientations along the east–west direction, simi-lar to our POL-2 detections (Wolf et al. 2003; Matthews etal. 2009). In addition, the polarization percentage observedwith SCUPOL is a factor of two to three higher than that in our POL-2 observations (Wolf et al. 2003). A detailed com-parison between the performance of POL-2 and SCUPOLhas been presented in Ward-Thompson et al. (2017) with theobservations of OMC 1. The comparison shows that the ob-served polarization orientations and percentages with POL-2and SCUPOL are consistent in the central bright regions ofOMC 1. On the other hand, in the outer faint regions, theSCUPOL detections tend to show higher polarization per-centages, and the di ff erence in the polarization orientationscan be as large as 30 ◦ to 50 ◦ (Fig. 7 and 8 in Ward-Thompsonet al. 2017). This discrepancy is possibly due to the low S / Nof the SCUPOL data in the faint region (Ward-Thompsonet al. 2017). As the noise level of our POL-2 observationsis a factor of three lower than that of the SCUPOL data(Matthews et al. 2009), the di ff erence between our POL-2and the SCUPOL results can be attributed to the higher noiselevel in the SCUPOL data.Figure 2b presents the Stokes I map of the 1.3 mm con-tinuum emission on a 1000 au scale in B335 overlaid withthe magnetic field orientations from the ALMA polarimet-ric observations. The details of the ALMA results have beendescribed by Maury et al. (2018). The polarization detec-tions obtained with ALMA are primarily along the wall ofthe outflow cavity, and there is an additional patch of detec-tions in the north. The area of the polarized emission de-tected with ALMA is within the JCMT beam at the inten-sity peak position. The number distribution of the magneticfield orientations inferred from the ALMA polarization de-tections is shown in Fig. 2d. The distribution has a doublepeak at a PA of 70 ◦ and 160 ◦ , as shown in Maury et al. (2018).This number distribution with the double peaks is consistentwith the expectation from a pinched magnetic field along theeast–west direction (Li & Shu 1996; Galli et al. 2006), asdiscussed and demonstrated in Maury et al. (2018) with theresults from their non-ideal MHD simulations.In summary, our JCMT results and the ALMA archivaldata show that the magnetic field on scales from 1000 au to6000 au in B335 is most likely along the east–west direction,which is aligned with the outflow direction within 10 ◦ on theplane of the sky, and the magnetic field is likely pinched ona 1000 au scale in B335. We note that the SMA polarimetricobservations at an angular resolution of ∼ (cid:48)(cid:48) show orangizedmagnetic field orientations tilted from the outflow axis by35 ◦ on a 1000 au scale in B335 (Galametz et al. 2018). Thiscould be due to the limited angular resolution and sensitivityof the SMA that likely only detected a part of the pinchedmagnetic field, leading to the mean orientation misalignedwith the outflow axis. MHD SIMULATIONSTo study the strength of the magnetic field and its e ff ect onthe gas kinematics in B335, we carried out three-dimensional (a) (b)(c) (d) JCMT B segmentsJCMT B segments ALMA B segmentsALMA B segments
Figure 2. (a) & (b) Magnetic field orientations (orange segments) inferred by rotating the polarization orientations by 90 ◦ obtained with theJCMT POL-2 and ALMA observations, respectively. Contours shows the Stokes I maps. Blue and red arrows present the directions of theblue- and redshifted outflows. Blue ellipses at the bottom right corners show the angular resolutions. (c) and (d) Number distributions of themagnetic field orientations. Dashed vertical lines denote the PA of the outflow axis. (3D) non-ideal MHD simulations using ZeusTW code(Krasnopolsky et al. 2010) to model the collapse of a mag-netized non-turbulent dense core, and generated syntheticpolarization maps from the simulations to compare with theobservations. 3.1. Simulation setup
In our simulations, among the three non-ideal MHD ef-fects, we only consider ambipolar di ff usion (AD) because itis the most e ffi cient mechanism in magnetic field decouplingon a scale from 100 au to thousands of au in a collapsingdense core (Zhao et al. 2016, 2018a). The magnetic di ff u-sivity of ambipolar di ff usion is self-consistently computed atrun-time using the equilibrium chemical network from Zhaoet al. (2016), and is closely related to size of dust grains andcosmic-ray ionization (Zhao et al. 2016, 2018b). For com- parison, we also performed ideal MHD simulations, whichare independent on grain size and cosmic-ray ionization.We prescribe the initial conditions of the dense core in oursimulations to approximately match the observed propertiesof B335. In B335, the mass of the dense core within a radiusof 6000 au is estimated to be 0.5–1.8 M (cid:12) based on the single-dish observational results in the millimeter continuum emis-sion and in the C O and H CO + lines (Saito et al. 1999;Harvey et al. 2003; Evans et al. 2005; Kurono et al. 2013).Thus, in our simulations, the initial core was set to have a to-tal mass of 1.0 M (cid:12) and an outer radius of 10 cm (6684 au),and its density distribution is uniform and spherical. We note We corrected the values in the literature for the distance from 250 pc to100 pc (Stutz et al. 2008; Olofsson & Olofsson 2009). that the outer radius of the initial core in our simulations issmaller than the core radius of ∼
10 000–14 000 au in B335(Saito et al. 1999). Nevertheless, the adopted outer radiusof the initial core is a factor of two larger than the estimatedradius of the dynamically infalling region in B335 (Choi etal. 1995; Harvey et al. 2003; Evans et al. 2005, 2015), wherethe interplay between the magnetic field and the gas kinemat-ics becomes prominent. The initial core is rotating as a solidbody in our simulations. The angular speed of its solid-bodyrotation is adopted to be 4 × − s − , which is the same asthe core rotation in B335 measured from the velocity gradientin the C O and H CO + emission observed with single-dishtelescopes (Saito et al. 1999; Yen et al. 2011; Kurono et al.2013). The corresponding ratio of rotational-to-gravitationalenergy β rot is around 0.4%.In our simulations, the magnetic field uniformly threadsthe initial core along the rotation axis. Thus, our simula-tions are initially axisymmetric. We focus on axisymmetricmodels because the misalignment between the magnetic fieldand the rotational axis is observed to be small on the planeof the sky ( < ◦ ; Section 2.3). Three di ff erent magneticfield strengths B of 42.5 µ G, 21 . µ G, and 10 . µ G wereadopted for strong, intermediate, and weak field cases. Thisgives a dimensionless mass-to-flux ratio λ ( ≡ M c π R B π √ G ) of2.4, 4.8, and 9.6, respectively. The adopted magnetic fieldstrength is the same as the observational estimate of 10–40 µ G on a 0.1 pc scale in B335 from the infrared polarimet-ric observations using the Chandrasekhar & Fermi (1953)method by Bertrang et al. (2014). After having narroweddown the main physical quantities (see Section 4), we alsocarried out two additional simulations with a slight misalign-ment with an angle of 15 ◦ between the magnetic field and therotation axis.In addition, we varied the minimum grain size a min andthe cosmic-ray ionization rate ζ because they play an im-portant role in the magnetic field decoupling in the collaps-ing envelope, which in turn determines the field strength andthe degree of pinch during the evolution as well as the disksize formed eventually. In our simulations, the slope of thegrain size distribution and the maximum grain size a max arefixed at − . µ m, respectively, and we adopted twodi ff erent values of a min , 0.005 µ m and 0.1 µ m. Those mod-els with a min = µ m and 0.1 µ m are denoted as MRNand tr-MRN (Mathis-Rumpl-Nordsieck; Mathis et al. 1977),respectively. As shown by Zhao et al. (2016), the tr-MRNcase strongly promotes disk formation due to the enhancedAD in the collapsing envelope. The simulations without thenon-ideal MHD e ff ect of AD are denoted as ideal. The mod-els with a typical ζ of 10 − s − are denoted as CR1, andthose with a higher ζ of 5 × − s − as CR5. The num-bers in the names of our models present the initial mass-to-flux ratios λ in those models. All our simulation models are Table 2.
Summary of MHD simulationsModel λ ζ M (cid:63) + disk Disk radius(10 s − ) ( M (cid:12) ) (AU)2.4Ideal 2.4 1 0.02 –4.8Ideal 4.8 1 0.02 –9.6Ideal 9.6 1 0.02 –2.4MRN-CR1 2.4 1 0.04 –4.8MRN-CR1 4.8 1 0.04 –9.6MRN-CR1 9.6 1 0.04 (cid:46)
10 au → < (cid:46) <
15 au4.8tr-MRN-CR5 4.8 5 0.04 < (cid:46)
10 au9.6 5 0.1 (cid:46)
10 au9.6tr-MRN-CR10 9.6 10 0.1 < → ote —Ideal MHD simulations and non-ideal MHD simulationswith the minimum grain sizes of 0.005 µ m and 0.1 µ m are denotedas ideal, MRN, and tr-MRN, respectively. λ is the initialmass-to-flux ratio. ζ is the cosmic-ray ionziation rate. M (cid:63) + disk isthe total mass of the central star + disk system when we extractedthe simulation data, and the disk radius is the radius of therotationally supported disk (if forms) at that time. Arrows indicatetransient disks. Model mis-tr-MRN-CR5 is the simulation with themisaligned magnetic field and rotational axis by 15 ◦ . In the modelswith the misalignment, the central disks are additionallysurrounded by spiral structures which extend to tens of au (e.g., Liet al. 2013). summarized in Table 2. The results of our non-ideal MHDsimulations were extracted when the total mass of the centralprotostar and (if any) rotationally supported disk ( M (cid:63) + disk )reached 0.04 M (cid:12) and 0.1 M (cid:12) , which is in the range of the ob-servational estimate of 0.04–0.2 M (cid:12) (Yen et al. 2015b; Evanset al. 2015). The results of the ideal MHD simulations wereextracted when M (cid:63) + disk was 0.02 M (cid:12) because these simula-tions fail to evolve longer due to the build-up and growth ofMHD instabilities close to the center (Zhao et al. 2011).3.2. Simulation results
Disk formation
We found that in the adopted range of λ , persistent disksonly form in the tr-MRN ( a min = . µ m) cases, as discussedin Zhao et al. (2016, 2018a). The ideal MHD models allfail to evolve to later stages after the formation of the pro- − −
50 0 50 100AU − − A U − −
50 0 50 100AU − − -15.5 -14.9 -14.3 -13.7 -13.1 -12.5 -11.8 -11.2 -10.6 -10.0 log 10 -16.4 -15.7 -15.0 -14.3 -13.6 -12.9 -12.1 -11.4 -10.7 -10.0 log 10 M ⋆ = 0.041 M ⊙ M d = 0.008 M ⊙ − −
50 0 50 100AU − − A U − −
50 0 50 100AU − − A U -15.6 -15.0 -14.4 -13.9 -13.3 -12.8 -12.2 -11.6 -11.1 -10.5 log10 ρ -16.5 -15.8 -15.2 -14.5 -13.8 -13.2 -12.5 -11.8 -11.2 -10.5 log10 ρ M ⋆ = 0.039 M ⊙ M d = 0.000 M ⊙ Figure 3.
Comparison of the disk and outflow morphologies between model 4.8MRN-CR1 (upper panels) and 4.8tr-MRN-CR1 (lower panels),along the equatorial (left) and meridian (right) planes, respectively. Color and arrows present the density and velocity in the simulations,respectively. The masses of the central protostar ( M (cid:63) ) and the disk ( M d ) are labeled below the panels. The disk and outflow structures arewell-defined in 4.8tr-MRN-CR1 but not in 4.8MRN-CR1. tostar, due to the build-up and growth of MHD instabilitiesnear the central stellar sink (Zhao et al. 2011). Thus, no diskis formed. The disk formation is also hindered in the MRNmodels. Among all the MRN models, only the one with theweak field and the low ionization, 9.6MRN-CR1, forms atransient ∼
10 au disk that shrinks over time due to influx ofgas with a low specific angular momentum (see also Zhaoet al. 2018a). Rotationally-supported disks generally form inour tr-MRN models. In these models, small dust grains withsizes of tens to hundreds of Å, which are well coupled to themagnetic field and also exert strong drag to neutral gas, areremoved. As a result, the ambipolar di ff usivity is enhancedby one to two orders of magnitude (Zhao et al. 2016). Withthe enhanced AD, magnetic braking becomes less e ffi cient,leading to the formation of a persistent disk. The radius ofthe disk typically remains below ∼
10 au for most tr-MRNmodels due to the slow core rotation, except for model 9.6tr-MRN-CR1 where both magnetic field strength and cosmic-ray ionization rate are low. In Fig. 3, we compare the typi-cal MRN and tr-MRN models in terms of disk and outflowmorphologies. In the 4.8tr-MRN-CR1 case (bottom panels),a small rotationally supported disk with a radius of few auforms and drives a clear bipolar outflow. In contrast, in the4.8MRN-CR1 case, no disk is formed, and there is also nowell-defined outflow.3.2.2.
Pinched magnetic field
In our MHD simulations, the magnetic field morphologygenerally follows an hour-glass configuration. We find thatthe degree of the pinch of the magnetic field lines is relatedto evolutionary time, the magnetic field strength in the initialcore, and the level of AD. Figure 4 presents the magnetic fieldlines in the tr-MRN models as an example to demonstrate thedependance of the degree of pinch on the initial magneticfield strength and evolutionary time.As shown in Fig. 4, the magnetic field lines are morepinched at the later time with M (cid:63) + disk of 0.1 M (cid:12) than at theearly time with M (cid:63) + disk of 0.04 M (cid:12) . This is because the mag-netic field lines are dragged more severely as more material iscollapsing into the center. The trend is more clearly revealedin the radial profiles of the pinch angles shown in Fig. 5. Thepinch angle is defined as the angle between the vertical axisand the direction of the magnetic field just above and belowthe midplane, and it is azimuthally averaged. Thus, the pinchangle here presents the degree of the magnetic field beingbent from the original orientation, which is perpendicular tothe midplane. A pinch angle closer to 90 ◦ means that themagnetic field is more pinched and become more parallel tothe midplane. Fig. 5 shows that the magnetic field becomespinched at the inner radii of a few hundred au because of thecollapse, and the pinch angle at these inner radii increases by 10 ◦ –20 ◦ from the earlier to later times, meaning that themagnetic field is more pinched at the later time.In addition, the comparison of the pinch angles in the threemodels in Fig. 5 shows that the magnetic field lines pinchmore severely in the models with the weaker initial fieldthan the stronger initial field. We find that the pinch angleis closely regulated by the Lorentz force. The radial pro-files of the Lorentz force are also presented in Fig. 5. TheLorentz force along the midplane (pseudo-disk) and pointingaway from the center is proportional to B z δ B r δ z in the cylindri-cal coordinate, where B z is the poloidal component and B r isthe radial component. In our simulations, the Lorentz forceincreases with decreasing radii, as the magnetic field linesare dragged inward by the collapse. Eventually, the Lorentzforce becomes approximately balanced with the gravity ona 100 au scale such that the AD drift velocity is high andthe velocity of the magnetic field is low. The simulation re-sults shown in Fig. 5 were extracted at the same evolution-ary time, when M (cid:63) + disk are 0.04 M (cid:12) (upper panels) or 0.1 M (cid:12) (lower panels). In these models at the same evolutionarytime, their gravitational forces on a 100 au scale are compa-rable. Consequently, the Lorentz forces on a 100 au scale inthese models with di ff erent magnetic field strengths all reachapproximately the same magnitude (Fig. 5). Given the sameLoreztz force, in the models with the weaker initial magneticfield and thus lower B z , the magnetic field lines pinch moreseverely to have a larger δ B r /δ z along the midplane (pseudo-disk), resulting in a larger pinch angle.In Fig. 6, we compare the 2.4MRN-CR1 model with lesse ffi cient AD and the 2.4trMRN-CR1 model with enhancedAD to study the dependence of the pinch angle on magneticdi ff usion, which is in the form of AD in our simulations. Themagnetic field lines in the MRN model show a sharp kinkat a radius of ∼
200 au, where the pinch angle is as large as50 ◦ –60 ◦ . Compared to the MRN model, the tr-MRN modelshows less pinched field lines with a smaller pinch angle of20 ◦ –30 ◦ at the same radius. Because of the enhanced AD,magnetic di ff usion is more e ffi cient in the tr-MRN modelthan in the MRN model. Thus, the magnetic fields in the tr-MRN model are dragged toward the center less e ffi ciently bythe collapsing material, which quenches the increase of theLorentz force, compared to the MRN model. As a result, theMRN model has a stronger Lorentz force on an inner scaleof 100–300 au than the tr-MRN model (Fig. 6). For a givenstrength of the poloidal component B z , the magnetic field ismore pinched when the Lorentz force is stronger because of alarger δ B r /δ z . In these two models, the field strengths of theinitial cores are both λ = .
4. Therefore, the magnetic fieldlines are more pinched in the MRN model with less e ffi cientmagnetic di ff usion than in the tr-MRN model.In summary, our MHD simulations suggest that the mag-netic field lines pinch more severely in a collapsing dense Figure 4.
Comparison of the structures of the magnetic field lines along the meridian plane (white lines) between cases with di ff erent fieldstrengths and at di ff erent evolutionary times. The model name is labeled at the upper left corner in each panel. Upper and lower panels presentthe models at the earlier and later evolutionary times when the star + disk masses are 0.04 M (cid:12) and 0.1 M (cid:12) , respectively. Color represents thedensity distributions. core with weaker magnetic field, longer time evolution, andlower magnetic di ff usion (or AD in our simulations). COMPARISON BETWEEN SIMULATIONS ANDOBSERVATIONSTo compare the magnetic field structures inferred from thepolarimetric observations with the simulations, we generatedsynthetic Stokes
IQU maps from our MHD simultions usingthe radiative transfer code, Simulation Package for Astro-nomical Radiative Xfer (SPARX; https: // sparx.tiara.sinica.edu.tw / ). SPARX adopts the properties of dust grains and for-mulization in Lee & Draine (1985) to compute Stokes param-eters and generate synthetic images. In the present paper, weonly compare the polarization orientations from the observa-tions and the synthetic images, and we do not compare the polarization percentage. To properly compute polarizationpercentage in our models requires more sophisticated cal-culations of mechanisms and e ffi ciency of grain alignments(e.g., Lazarian & Hoang 2007), which is beyond the scope ofthe present paper. Thus, in our calculations of radiative trans-fer, the polarization e ffi ciency is assumed to be constant, andwe do not consider polarization due to dust scattering (e.g.,Kataoka et al. 2015; Yang et al. 2016, 2017).For each model, we generated Stokes IQU maps at twowavelengths, 850 µ m and 1.3 mm, the same as the obser-vations using SPARX. We convolved the model Stokes IQU maps at 850 µ m and 1.3 mm with the beams of the JCMT andALMA observations, respectively. We have also performedALMA imaging simulations on our models with the most andleast pinched magnetic field, and found that there is no sig-0 Radius (AU)0102030405060708090 P i n c h A ng l e ( D eg ) Radius (AU)10 -29 -28 -27 -26 -25 -24 -23 -22 -21 -20 -19 -18 R ad i a l Lo r en t z F o r c e g ·( c m · s ) − Radius (AU)0102030405060708090 P i n c h A ng l e ( D eg ) Radius (AU)10 -29 -28 -27 -26 -25 -24 -23 -22 -21 -20 -19 -18 R ad i a l Lo r en t z F o r c e g ·( c m · s ) − Figure 5.
Profiles of the pinch angle (left) and radial Lorentz force (right) across the equatorial plane in the models with di ff erent magneticfield strengths and at di ff erent evolutionary times. A larger pinch angle denotes that the magnetic field lines pinch more severely. Black solid,blue dashed, and red dashed-dotted lines present model 2.4tr-MRN-CR1, 4.8tr-MRN-CR1, and 9.6tr-MRN-CR1, respectively. Upper and lowerpanels present the models at the earlier and later evolutionary times when the star + disk masses are 0.04 M (cid:12) and 0.1 M (cid:12) , respectively. Radius (AU)010203040506070 P i n c h A ng l e ( D eg ) Radius (AU)10 -29 -28 -27 -26 -25 -24 -23 -22 -21 -20 -19 -18 R ad i a l Lo r en t z F o r c e g ·( c m · s ) − Figure 6.
Dependence of the pinch of the magnetic field on ambipolar di ff usion. Upper panels: comparsion of the structures of the magneticfield lines along the meridian plane (white lines) between the MRN (left) and tr-MRN (right) models with the same initial λ of 2.4 when thecentral star + disk mass is 0.1 M (cid:12) . Color scales present the density distributions. Lower panels: profiles of the pinch angle (left) and radialLorentz force (right) across the equatorial plane. A larger pinch angle denotes that the magnetic field lines pinch more severely. ff erence between the polarization orientations ex-tracted from the simulated and convolved model maps in theregions where there are ALMA polarization detections (Ap-pendix A). Therefore, we adopt the convolved model mapsfor the comparison between the models and the observations.We re-gridded the convolved model maps to have the sameimage and pixel sizes as our observational maps. When wecompared the number distributions of the magnetic field ori-entations between the observations and models, the modelfield orientations were extracted from the model maps at thesame positions as the detections in the observational maps.Thus, the total numbers of the segments are the same in thecomparison.4.1. Magnetic field on a 6000 au scale
Figure 7 presents two examples of the comparison betweenthe magnetic field orientations from the JCMT observationsand from the models. Two extreme cases are shown here,2.4Ideal having the strongest magnetic field and coupling be-tween the field and the matter, and 9.6tr-MRN-CR1 havingthe weakest field and coupling, among our models. The PAof the magnetic field axis in the synthetic maps is adopted tobe 90 ◦ . Figure 7 shows that the magnetic field structures on a6000 au scale in B335 observed with the JCMT can be welldescribed with the pinched magnetic field in both simula-tions. Although the degrees of the pinch are di ff erent amongour MHD simulations, there is no detectable di ff erence on a6000 au scale with the angular resolution of the JCMT obser-vations. Thus, all our MHD simulations can well explain themagnetic field structures observed with the JCMT, as demon-strated in Fig. 7. We note that the magnetic field orientationdetected at the o ff set of ( − (cid:48)(cid:48) , 16 (cid:48)(cid:48) ) has a relatively large de-viation of 20 ◦ –25 ◦ from all the models. For the other detec-tions, the mean di ff erence between the observed orientationsand the models is 5 ◦ –9 ◦ , comparable to the observational un-certainty of 6 ◦ –8 ◦ , when the PA of the magnetic field axis inour synthetic maps is adopted to be 87 ◦ –94 ◦ . Therefore, thecomparison between the JCMT observations and our MHDsimulations suggest that in B335, the magnetic field on a6000 au scale is pinched and well aligned with the directionof the bipolar outflow within a few degrees on the plane ofthe sky. 4.2. Magnetic field on a 1000 au scale
Figure 8a and b present the comparison between the mag-netic field orientations on a 1000 au scale from the ALMAobservations and from our models 9.6tr-MRN-CR5 with M (cid:63) + disk of 0.1 M (cid:12) and 2.4Ideal with M (cid:63) of 0.02 M (cid:12) , whichare the models showing the most and least pinched fieldamong our models, respectively. The overall morphologyof the observed magnetic field orientations is similar tothe pinched magnetic field in our MHD simulations, except (a) 2.4ideal(b) 9.6tr-MRN-CR1 Figure 7.
Comparison of the magnetic field orientations from thePOL-2 observations (orange segments) and from the synthetic maps(green segments) of our models, (a) 2.4Ideal and (b) 9.6tr-MRN-CR1. There is no detectable di ff erence between the two syntheticmaps with the JCMT resolution. for the detections in the northern and northwestern patches(shown as blue segments). The northern patch of the mag-netic field segments is oriented along the north–south direc-tion (PA of 0 ◦ ), which is the direction of the midplane. Al-though our model 9.6tr-MRN-CR5 also shows the magneticfield segments along the direction of the midplane, these seg-ments are located close to the midplane and within a radius3 (a) 2.4ideal (M * = 0.02 Msun) (b) 9.6tr-MRN-CR5 (M *+d = 0.1 Msun)(c) 2.4ideal (M * = 0.02 Msun)(e) 2.4ideal (M * = 0.02 Msun) (d) 9.6tr-MRN-CR5 (M *+d = 0.1 Msun)(f) 9.6tr-MRN-CR5 (M *+d = 0.1 Msun) northern patch northwestern patchB parallel to midplane Figure 8. (a) & (b) Comparison of the magnetic field orientations from the ALMA observations (blue and orange segments) and from thesynthetic maps (black segments) of our models, 2.4Ideal with the stellar mass of 0.02 M (cid:12) (left) and 9.6tr-MRN-CR5 with the star + disk mass of0.1 M (cid:12) (right), respectively. The center of the map, (0, 0), is the protostellar position. (c) and (d) Number fraction distributions of the magneticfield orientations from the ALMA observations (orange histogram) and the models (black dotted histograms). Only the orange segments in (a)and (b) and the model segments at the same positions are included in these plots. (e) & (f) Number fraction distributions of the PA di ff erencein the magnetic field orientations from the ALMA observations and the models. The northern and northwestern patches discussed in Section 4and 5 are labeled in (a). A red circle in (a) denotes the region showing the magnetic field orientation almost parallel to the midplane. Figure 9.
Same as Fig. 8 (a) and (b) but for model mis-tr-MRN-CR5, where the magnetic field and rotation axis are misaligned by15 ◦ . of 2 (cid:48)(cid:48) (200 au). All the segments at a radius larger than 200au are tilted from the direction of the midplane in our models.On the other hand, the orientations of the northwestern mag-netic field segments change from the northwest–southeastdirection with a PA of ∼ ◦ to the northeast–southwestdirection with a PA of ∼ ◦ as the distance to the centerincreases. At the outer radii larger than 250 au (RA o ff sets < − . (cid:48)(cid:48) ◦ in thenumber distribution. These segments can be well explainedwith the pinched magnetic field in our models. The segmentsobserved around the center are mostly oriented close to thedirection of the midplane. Such an orientation of the mag-netic field segments around the center can also be explainedwith our weak field models with initial λ of 9.6 at the laterevolutionary time when M (cid:63) + disk is 0.1 M (cid:12) . We found that allof our models with the stronger field of initial λ of 2.4 and4.8 do not show such magnetic field segments along the mid-plane direction even at the later evolutionary time. In thesemodels, the magnetic field segments around the center aremore perpendicular to the midplane, as the case of 2.4Idealshown in Fig. 8c. The observed magnetic field segmentsin the southeast are oriented along the northwest–southeastdirection, corresponding to the peak at PA close to 150 ◦ in the number distribution. Although the overall orientation ofthese segments is similar to those in our models, the seg-ments in our models are more tilted away from the midplaneand have a PA of 100 ◦ , as seen in the number distribution(Fig. 8d). The magnetic field segments with PA of ∼ ◦ are also present in the southeast in our model map of theweak field case with initial λ of 9.6 (Fig. 8b), but they arelocated closer to the midplane and do not extended to 1 (cid:48)(cid:48) –2 (cid:48)(cid:48) from the midplane as the observations. Additionally, wealso performed simulations with the magnetic field and rota-tion axis misaligned by 15 ◦ . We rotated the magnetic fielddirection to have a PA of 90 ◦ , when we generated the syn-thetic maps of these simulations with the misalignment. Al-though the model maps from the simulations with the mis-alignment show di ff erent PAs of the magnetic field orienta-tions in the northwest and the southeast, compared to thosewith the aligned magnetic field and rotational axis (Fig. 9),incorporating the misalignment still cannot explain the ob-served magnetic field orientations in the northwest and thesoutheast.In summary, there are two patches of the magnetic field ori-entations observed in the north and the northwest that cannotbe explained with our MHD simulations. The possible rea-sons are discussed in Section 5. The other observed magneticfield orientations can be explained with the pinched magneticfield in our simulations, and we find that the simulationswith the weak field of initial λ of 9.6 explain the observa-tions better than our other simulations with the stronger field(e.g., Fig. 8e and f). Especially, the presence of the magneticfield orientations parallel to the midplane observed aroundthe center can only be explained with the models with initial λ of 9.6 at the later evolutionary time when M (cid:63) + disk is 0.1 M (cid:12) among all of our models. This is similar to the results byMaury et al. (2018). They found that the observed magneticfield orientations are similar to their models with the largestmass-to-flux ratio µ of 6, corresponding to λ of 6.3. Nev-ertheless, the observed magnetic field segments in the south-east are oriented closer to the midplane by ∼ ◦ compared tothose in our models. In addition, we note that the observedmagnetic field structures are not symmetric with respect tothe magnetic field direction of PA of ∼ ◦ . The number dis-tribution of the observed magnetic field orientations has twopeaks at PA of 70 ◦ and 150 ◦ (Fig. 8), and thus, the symmetricaxis is at PA of 110 ◦ . This distribution is di ff erent from thatin our models having the peaks at PA of 70 ◦ and 105 ◦ withthe magnetic field direction of PA of 90 ◦ . DISCUSSIONS Because of the di ff erent definitions of the mass-to-flux ratio, the mass-to-flux ratio µ in Maury et al. (2018) is √ λ . Figure 10.
Three-dimensional structures of the magnetic field lines in our model 9.6tr-MRN-CR5 (green lines). Color scale presents thedensity.
Discrepancy between observations and simulations
Although the overall morphology of the observed magneticfield in B335 is similar to the pinched field in our MHD sim-ulations (Fig. 7), the symmetric axis of the magnetic field ori-entations on a 1000 au scale observed with ALMA is likelyat PA of ∼ ◦ , di ff erent from the magnetic field directionof 87 ◦ –94 ◦ on a 6000 au scale inferred from the JCMT ob-servations. In addition, a part of the magnetic field segmentsin the northwest observed with ALMA are almost perpendic-ular to those in our simulations (Fig. 8). We found that inour MHD simulations, the magnetic field is highly wrappedaround the bipolar outflow, and there are significant toroidalcomponents of the magnetic field in the wall of the outflowcavity (Fig. 10). Because the outflow in our simulations aremore or less axisymmetric, these toroidal components arecancelled out after the integration along the line of sight inour radiative transfer calculation, and only the poloidal com-ponents that are aligned with or moderately tilted from theaxis of the bipolar outflow are seen in the region of the out-flow in our synthetic maps (Fig. 8). In the ALMA observa-tions, the polarized emission is primarily detected in the wallof the outflow cavity. The structures of the bipolar outflow in B335 are clearly asymmetric as observed in the CO line (Yenet al. 2010) and in the continuum (Fig. 2b Maury et al. 2018).Thus, the inhomogeneous density structures (if any) in thewall of the outflow cavity could a ff ect the observed polariza-tion orientations in B335 because the magnetic field directioncould change a lot along the line of sight passing through theoutflow, as seen in our MHD simulations. As presented inFig. 10, in our simulations, there are indeed magnetic fieldlines perpendicular to the wall of the outflow cavity. Thiscould explain the observed magnetic field segments in thenorthwestern region in B335 with ALMA. Therefore, the dis-crepancy in the magnetic field orientations in the northwestin the observations and the simulations could be due to the ef-fect of the integration along the line of sight passing throughthe asymmetric outflow.The northern patch of the magnetic field orientations ob-served with ALMA, which are almost parallel to the mid-plane, cannot be reproduced in our MHD simulations regard-less of the field strengths considered. As suggested by ourJCMT results, the magnetic field direction on the large scaleis most likely perpendicular to the midplace. To form mag-netic field lines severely pinched and parallel to the midplane6on a 1000 au scale as observed with ALMA, one possibil-ity is that there is an accretion flow with a high mass in-falling rate, a weak field strength, and low magnetic di ff u-sion from the northern region in B335, as discussed in Sec-tion 3.2.2. The severely pinched magnetic field formed by astrong accretion flow in B335 was also suggested by Mauryet al. (2018). Nevertheless, such magnetic field orientationsparallel to the midplane are not detected in the southern re-gion. If the magnetic field orientations parallel to the mid-plane in the north are indeed formed by a strong accretionflow, the results could hint that the distribution of the collaps-ing material is not symmetric. An asymmetric accretion flowdragging the magnetic field line to form a severely pinchedmagnetic field at a radius of 200–500 au is not present inour initially axisymmetric MHD simulations or our simula-tions with the misaligned magnetic field and rotational axis.Asymmetric accretion flows have been seen in other MHDsimulations with turbulence, such as those in Li et al. (2014),where the pseudo disk and thus the accretion flow are warpedbecause of the turbulence. On the other hand, we note thatthere is no polarized emission detected in the northeasternand southern regions with our JCMT observations. This non-detection could hint at the presence of complex magneticfield structures in those regions. Consequencely, the polar-ized emission with di ff erent polarization orientations couldbe cancelled out. ALMA polarimetric observations with awider field of view and a high resolution are needed to imagecomplete field structures in B335.5.2. Magnetic field strength and gas kinematics in B335
The magnetic field structures on a 6000 au scale observedwith the JCMT can be well explained with the pinched fieldin our MHD simulations, but the angular resolution of theJCMT observations is insu ffi cient to distinguish di ff erent de-grees of the pinched field in our simulations with di ff erentfield strengths (Fig. 7). The comparison between our MHDsimulations and the ALMA observations shows that the ob-served magnetic field structures on a 1000 au scale are betterexplained with our simulations with the weak magnetic fieldof initial λ of 9.6 than with the strong magnetic field of initial λ of 2.4 and 4.8 (Fig. 8). Especially, the observed magneticfield orientations close to the center within a radius of 200 auare almost parallel to the midplane. Such field orientationsare only seen in our simulations with the weak field at thelater evolutionary time when M (cid:63) + disk is 0.1 M (cid:12) . Thus, theseresults favor a weak magnetic field of initial λ of 9.6 in B335,as also suggested by Maury et al. (2018).However, we also found that with the weak magnetic field,a Keplerian disk with a radius of 10 au quickly forms andgrows to have a radius of 50 au when M (cid:63) + disk is 0.1 M (cid:12) in ourtr-MRN simulations with a typical cosmic-ray ionzation rateof 10 − s − , while no disk forms in our ideal MHD or MRN simulations. The ALMA observations in the C O and SOlines at an angular resolution of 0 . (cid:48)(cid:48) ffi ciency of magnetic braking.Our tr-MRN simulation with the weak magnetic field and ahigher cosmic-ray ionzation rate of 5 × − s − forms asmall Keplerian disk with a radius smaller than 10 au, andcan explain the observed magnetic field parallel to the mid-plane within a radius of 200 au (Fig. 7b).In Fig. 11, we compare the radial profiles of the infallingand rotational velocities along the midplane in the simula-tions and the observational estimates. The infalling veloci-ties on a scale of 100–250 au in the simulations with di ff erentmagnetic field strengths are all consistent with the observa-tions within the uncertainties because in these simulations, M (cid:63) + disk is in the range of the observational estimates. In con-trast to that, the rotational velocities in the simulations withthe weak magnetic field and the high cosmic-ray ionizationrate are higher than those estimated from the observations bya factor of four to ten, suggesting that with the weak magneticfield of initial λ of 9.6, magnetic braking is not e ffi cient toslow down the rotational motion. As also discussed in Yen etal. (2015b), the observed specific angular momentum of thecollapsing material on a 100 au scale in B335 is a factor oftwo lower than the expectation from the inside-out collapsemodel where the angular momentum is conserved, sugges-tive of e ffi cient magnetic braking in the inner envelope. Wenote that even in our ideal MHD simulations, the rotationalvelocity in the simulations is still higher than in the obser-vations by a factor of three, when the initial magnetic fieldstrength is λ of 9.6. Thus, the higher rotational velocity andine ffi cient magnetic braking in these simulations are due tothe magnetic field strength but not the non-ideal MHD ef-fect. Only when the initial magnetic field strength is strongerthan λ of 4.8, the rotational velocities in the simulations andthe observations become comparable within a factor of two.7 (a)(b) Figure 11.
Radial profiles of infalling (solid lines) and rotational(dashed lines) velocities along the midplane from the models with(a) the weak magnetic field of initial λ of 9.6 and (b) the strongermagnetic field of initial λ of 2.4 and 4.8 in comparison with theobservational measurements (data points) obtained with the SMAand ALMA. Di ff erent models are denoted with di ff erent colors andare labeled at the bottom right corner in each panel. Because B335 isalmost edge on, the uncertainty of the estimated rotational velocity,which is 10%–20%, is smaller than that of the estimated infallingvelocity (Yen et al. 2015b). Thus, the error bars of the rotationalvelocity cannot be seen clearly in the figures. Therefore, the comparison of the rotational velocities in thesimulations and the observations could favor the strong mag-netic field of initial λ < . ff erent fromthe inferred magnetic field strength of λ of 9.6 based on thepolarimetric observations described above.There are two possibilities resulting in the di ff erentmagnetic field strengths inferred from the polarimetricand molecular-line observations: (1) the rotational-to-gravitational energy β rot is overestimated and (2) additionalcontributions in the polarized intensity from other mecha-nisms, such as dust scattering. In our MHD simulations, β rot is adopted to be 0.4% based on the observational esti-mates of the core mass of ∼ M (cid:12) and the angular speed ofthe core rotation of 4 × − s − . The angular speed wasestimated based on the global velocity gradient along themajor axis of the dense core observed with the single-dishtelescopes (Saito et al. 1999; Yen et al. 2011; Kurono et al.2013). Numerical simulations of dense cores including syn- thetic observations show that the specific angular momentumderived from the synthetic images of the dense cores canbe a factor of 8–10 higher than their actual specific angularmomentum computed by the sum of the angular momentacontributed by the individual gas parcels in the dense cores(Dib et al. 2010). In addition, if there are filamentary struc-tures in the dense core in B335, which could not be resolvedwith the single-dish observations, infalling motions along thefilamentary structures could also contribute to the observedvelocity gradient, leading to an overestimated angular speedof the core rotation (Tobin et al. 2012). We have also per-formed our simulations with a lower β rot , and we find thatthe rotational velocity on a 100 au scale in the simulationsdecreases with decreasing β rot . Thus, the discrepancy in themagnetic field strengths inferred from the field structures andthe gas kinematics can be reconciled, if the core rotation inB335 is overestimated by a factor of a few in the observa-tions, and these results would suggest a weak magnetic fieldof initial λ of 9.6 in B335. Further observations combiningsingle dishes and interferometers to have a high spatial dy-namical range and to map the velocity structures of the entiredense core in B335 at a high angular resolution are needed tostudy coherent velocity features and provide a better estimateof the core rotation.The other possibility is that the polarized emission ob-served close to the center could have contributions from dustscattering. If large dust grains with a size of the order of100 µ m are present in the inner dense region on a 100 auscale, where the density is higher than 10 cm − (Harvey etal. 2003; Evans et al. 2015), scattering of anisotropic con-tinuum emission by these dust grains could induce polarizedemission with a polarization percentage of ∼
1% (Kataoka etal. 2015; Yang et al. 2016, 2017), similar to that observedaround the center with ALMA (Maury et al. 2018). In ad-dition, the polarization orientations close to the midplanearound the center observed with ALMA are almost perpen-dicular to the midplane, and thus, the inferred magnetic fieldorientations are parallel to the midplane (Maury et al. 2018),as shown in Fig. 2. The presence of the polarization orien-tations perpendicular to the midplane are consistent with theexpectation of scattering-induced polarization in an edge-ondisk (Yang et al. 2016, 2017), and B335 is indeed an edge-onsource (Hirano et al. 1988; Stutz et al. 2008). Signatures ofscattering-induced polarized emission have been observed ina few embedded young protostars (Stephens et al. 2017; Coxet al. 2018; Lee et al. 2018). Thus, because of the contri-bution of the scattering, the actual magnetic field structuresclose to center may not be as severely pinched as discussedin Section 4.2. In this case, these results would favor thestrong magnetic field of initial λ < . >
200 au in the north(Section 5.1), where the grains are unlikely large enough toinduce any significant polarzation through scattering. FutureALMA polarimetric observations at higher angular resolu-tions and at di ff erent wavelengths are required to resolve thepolarization emission on a 100 au scale in B335 and to studyits nature. SUMMARYWe present the results and analysis of our JCMT POL-2observations at 850 µ m and the ALMA archival polarimetricdata at 1.3 mm of B335. In addition, we carried out a seriesof (non-)ideal MHD simulations of the collapse of a rotatingnon-turbulent dense core, whose mass and angular speed ofthe rotation are adopted to be the same as those observed inB335. We generated synthetic polarization maps from thesesimulations to compare with the polarimetric data. Our mainresults are summarized below.1. We carried out MHD simulations with di ff erent mag-netic field strengths of λ of 2.4, 4.8, and 9.6, cosmic-ray ionization rates of 1, 5, and 10 × − s − , anddistributions of dust grain sizes with minimum sizes of0.005 µ m and 0.1 µ m. The latter two parameters de-termine the magnetic di ff usivity. With the initial corerotation and mass set to be the same as in B335, wefind that a persistent disk only forms in the simulationswith the minimum grain size of 0.1 µ m, where ambipo-lar di ff usion is enhanced because of the removal of thesmall grains, regardless of other parameters. In oursimulations, the magnetic field lines generally show anhour-glass morphology. We find that the degree of thepinch of the magnetic field is regulated by the Lorentzforce, and that the magnetic field in the simulationswith weaker field strength, lower magnetic di ff usivity,and longer evolutionary time is pinched more severely.2. The magnetic field orientations on a 6000 au scale inB335 inferred from the JCMT POL-2 observations arealong the east–west direction, well aligned with the di-rection of the outflow within 10 ◦ on the plane of thesky. The observed magnetic field structures with theJCMT can be well explained with all our simulationswith di ff erent magnetic field strengths. The compari-son between our JCMT and simulation results suggestthat the magnetic field on a 6000 au scale in B335 ispinched. The JCMT resolution is not su ffi cient to dis-tinguish di ff erent magnetic field structures in our sim-ulations with di ff erent field strengths.3. The ALMA polarization results show the signatures ofa pinched magnetic field on a 1000 au scale in B335,which are similar to our MHD simulations. In addi-tion, the magnetic field orientations close to the center and at radii of 200–500 au in the north are almost par-allel to the midplane. The observed field orientationsin the north cannot be explained with our simulations,where the magnetic field orientations are tilted fromthe midplane. These magnetic field structures couldbe caused by a strong accretion flow from the north.Among all our MHD simulations, the observed mag-netic field orientations almost parallel to the midplanein the central region can only be explained with theweak field models having an initial mass-to-flux ratioof 9.6 at the late evolutionary time when the mass ofthe central star + disk system is 0.1 M (cid:12) .4. The comparison between the ALMA and simulationresults favor weak field models with an initial mass-to-flux ratio of 9.6. Nevertheless, with such a weakmagnetic field, the rotational velocity on a 100 au scaleand the size of the Keplerian disk formed in our sim-ulations are a factor of several higher than the obser-vational estimates. We find that when the cosmic-rayionization rate is increased by a factor of five or moreto enhance the field–matter coupling and thus the e ffi -ciency of magnetic braking, the disk size in our weakfield simulations is reduced and consistent with the ob-servational estimate, while the rotational velocity on a100 scale in the simulations is still higher than in theobservations. On the other hand, the rotational veloc-ity on a 100 au scale and the disk size in the simula-tions both become comparable to the observations inour stronger field models with an initial mass-to-fluxratio smaller than 4.8, regardless of magnetic di ff usiv-ity. Thus, the comparison between the observed gaskinematics in B335 and our MHD simulations favorsthe stronger field models.5. There are two possibilities resulting in the di ff erentmagnetic field strengths inferred from the magneticfield structures and the gas kinematics: (1) an overes-timated rotational-to-gravitational energy β rot and (2)additional contributions in the polarized intensity fromdust scattering. The rotational velocity on a 100 auscale in the simulations is proportional to β rot of theinitial core. The presence of incoherent gas motions orinfalling motion along filamentary structures (if any)in the dense core in B335 could lead to an overesti-mated β rot from the global velocity gradient observedwith single-dish telescopes with limited resolutions. Inthis case, our results favor the weak field models forB335. On the other hand, dust scattering in an edge-on source like B335 tends to induce a polarization ori-entation perpendicular to the midplane. This e ff ectcould make the inferred magnetic field orientations inthe central dense region more parallel to the midplane,9and the actual magnetic field structures might not beseverely pinched. In this case, based on the observed gas kinematics, our results favor the strong field mod-els for B335.REFERENCES Allen, A., Li, Z.-Y., & Shu, F. H. 2003, ApJ, 599, 363Aso, Y., Ohashi, N., Saigo, K., et al. 2015, ApJ, 812, 27Attard, M., Houde, M., Novak, G., et al. 2009, ApJ, 702, 1584Bertrang, G., Wolf, S., & Das, H. S. 2014, A&A, 565, A94Blandford, R. D., & Payne, D. G. 1982, MNRAS, 199, 883Caselli, P., Walmsley, C. M., Terzieva, R., Herbst, E. 1998, ApJ,499, 234Chandrasekhar, S., & Fermi, E. 1953, ApJ, 118, 116Choi, M., Evans, N. J., II, Gregersen, E. M., & Wang, Y. 1995,ApJ, 448, 742Chou, T.-L., Takakuwa, S., Yen, H.-W., Ohashi, N., & Ho, P. T. P.2014, ApJ, 796, 70Cox, E. G., Harris, R. J., Looney, L. W., et al. 2018, ApJ, 855, 92Crutcher, R. M. 2012, ARA&A, 50, 29Currie, M. J., Berry, D. S., Jenness, T., et al. 2014, AstronomicalData Analysis Software and Systems XXIII, 485, 391Dapp, W. B., Basu, S., & Kunz, M. W. 2012, A&A, 541, A35Davidson, J. A., Li, Z.-Y., Hull, C. L. H., et al. 2014, ApJ, 797, 74Dib, S., Hennebelle, P., Pineda, J. E., et al. 2010, ApJ, 723, 425Evans, N. J., II, Di Francesco, J., Lee, J.-E., et al. 2015, ApJ, 814,22Evans, N. J., II, Lee, J.-E., Rawlings, J. M. C., & Choi, M. 2005,ApJ, 626, 919Friberg, P., Bastien, P., Berry, D., et al. 2016, Proc. SPIE, 9914,991403Galametz, M., Maury, A., Girart, J. M., et al. 2018,arXiv:1804.05801Galli, D., Lizano, S., Shu, F. H., & Allen, A. 2006, ApJ, 647, 374Girart, J. M., Beltr´an, M. T., Zhang, Q., Rao, R., & Estalella, R.2009, Science, 324, 1408Girart, J. M., Rao, R., & Marrone, D. P. 2006, Science, 313, 812Harvey, D. W. A., Wilner, D. J., Myers, P. C., & Tafalla, M. 2003,ApJ, 596, 383Hirano, N., Kameya, O., Nakayama, M., & Takakubo, K. 1988,ApJL, 327, L69Holland, W. S., Bintley, D., Chapin, E. L., et al. 2013, MNRAS,430, 2513Hull, C. L. H., Plambeck, R. L., Kwon, W., et al. 2014, ApJS, 213,13Joos, M., Hennebelle, P., & Ciardi, A. 2012, A&A, 543, A128Kataoka, A., Muto, T., Momose, M., et al. 2015, ApJ, 809, 78Keene, J., Davidson, J. A., Harper, D. A., et al. 1983, ApJL, 274,L43 Keene, J., Hildebrand, R. H., Whitcomb, S. E., & Harper, D. A.1980, ApJL, 240, L43Krasnopolsky, R., Li, Z.-Y., & Shang, H. 2010, ApJ, 716, 1541Kurono, Y., Saito, M., Kamazaki, T., Morita, K.-I., & Kawabe, R.2013, ApJ, 765, 85Kwon, W., Stephens, I., Tobin, J., et al. 2018, arXiv:1805.07348Lazarian, A., & Hoang, T. 2007, MNRAS, 378, 910Lee, H. M., & Draine, B. T. 1985, ApJ, 290, 211Lee, C.-F., Hwang, H.-C., & Li, Z.-Y. 2016, ApJ, 826, 213Lee, C.-F., Li, Z.-Y., Ching, T.-C., Lai, S.-P., & Yang, H. 2018,ApJ, 854, 56Li, J. I., Liu, H. B., Hasegawa, Y., & Hirano, N. 2017, ApJ, 840, 72Li, Z.-Y., Krasnopolsky, R., & Shang, H. 2011, ApJ, 738, 180Li, Z.-Y., Krasnopolsky, R., & Shang, H. 2013, ApJ, 774, 82Li, Z.-Y., Krasnopolsky, R., Shang, H., & Zhao, B. 2014, ApJ, 793,130Li, Z.-Y., & Shu, F. H. 1996, ApJ, 472, 211Masson, J., Chabrier, G., Hennebelle, P., Vaytet, N., &Commerc¸on, B. 2016, A&A, 587, A32Mathis, J. S., Rumpl, W., & Nordsieck, K. H. 1977, ApJ, 217, 425Matthews, B. C., McPhee, C. A., Fissel, L. M., & Curran, R. L.2009, ApJS, 182, 143Maury, A. J., Girart, J. M., Zhang, Q., et al. 2018, MNRAS, 477,2760McMullin, J. P., Waters, B., Schiebel, D., et al. 2007, AstronomicalData Analysis Software and Systems XVI (ASP Conf. Ser. 376),ed. R. A. Shaw, F. Hill, & D. J. Bell (San Francisco, CA: ASP),127Mellon, R. R., & Li, Z.-Y. 2008, ApJ, 681, 1356Mellon, R. R., & Li, Z.-Y. 2009, ApJ, 698, 922Motte, F., & Andr´e, P. 2001, A&A, 365, 440Ohashi, N., Saigo, K., Aso, Y., et al. 2014, ApJ, 796, 131Olofsson, S., & Olofsson, G. 2009, A&A, 498, 455Padovani, M., Galli, D., Hennebelle, P., Commerc¸on, B., & Joos,M. 2014, A&A, 571, A33Pattle, K., Ward-Thompson, D., Berry, D., et al. 2017, ApJ, 846,122Qiu, K., Zhang, Q., Menten, K. M., et al. 2014, ApJL, 794, L18Qiu, K., Zhang, Q., Menten, K. M., Liu, H. B., & Tang, Y.-W.2013, ApJ, 779, 182Rao, R., Girart, J. M., Marrone, D. P., Lai, S.-P., & Schnee, S.2009, ApJ, 707, 921Sadavoy, S. I., Myers, P. C., Stephens, I. W., et al. 2018, ApJ, 859,165 Saito, M., Sunada, K., Kawabe, R., Kitamura, Y., & Hirano, N.1999, ApJ, 518, 334Shirley, Y. L., Evans, N. J., II, & Rawlings, J. M. C. 2002, ApJ,575, 337Shu, F. H., Adams, F. C., & Lizano, S. 1987, ARA&A, 25, 23Simmons, J. F. L., & Stewart, B. G. 1985, A&A, 142, 100Stephens, I. W., Yang, H., Li, Z.-Y., et al. 2017, ApJ, 851, 55Stutz, A. M., Rubin, M., Werner, M. W., et al. 2008, ApJ, 687,389-405Tobin, J. J., Hartmann, L., Bergin, E., et al. 2012, ApJ, 748, 16Tsukamoto, Y., Iwasaki, K., Okuzumi, S., Machida, M. N., &Inutsuka, S. 2015, MNRAS, 452, 278Ward-Thompson, D., Pattle, K., Bastien, P., et al. 2017, ApJ, 842,66Wardle, J. F. C., & Kronberg, P. P. 1974, ApJ, 194, 249Wolf, S., Launhardt, R., & Henning, T. 2003, ApJ, 592, 233 Yang, H., Li, Z.-Y., Looney, L. W., Girart, J. M., & Stephens, I. W.2017, MNRAS, 472, 373Yang, H., Li, Z.-Y., Looney, L., & Stephens, I. 2016, MNRAS,456, 2794Yen, H.-W., Koch, P. M., Takakuwa, S., et al. 2015, ApJ, 799, 193Yen, H.-W., Takakuwa, S., Koch, P. M., et al. 2015b, ApJ, 812, 129Yen, H.-W., Takakuwa, S., & Ohashi, N. 2010, ApJ, 710, 1786Yen, H.-W., Takakuwa, S., & Ohashi, N. 2011, ApJ, 742, 57Yen, H.-W., Zhao, B., Koch, P. M., et al. 2018, A&A, 615, A58Zhao, B., Caselli, P., Li, Z.-Y., et al. 2016, MNRAS, 460, 2050Zhao, B., Caselli, P., Li, Z.-Y., & Krasnopolsky, R. 2018, MNRAS,473, 4868Zhao, B., Caselli, P., & Li, Z.-Y. 2018, MNRAS, 478, 2723Zhao, B., Li, Z.-Y., Nakamura, F., Krasnopolsky, R., & Shang, H.2011, ApJ, 742, 10Zhou, S. 1995, ApJ, 442, 685Zhou, S., Evans, N. J., II, Koempe, C., & Walmsley, C. M. 1993,ApJ, 404, 232 / JAO.ALMA / NRAOand NAOJ. The JCMT data were obtained under program ID M17AP067. JCMT is operated by the East Asian Observatoryon behalf of The National Astronomical Observatory of Japan; Academia Sinica Institute of Astronomy and Astrophysics; theKorea Astronomy and Space Science Institute; the Operation, Maintenance and Upgrading Fund for Astronomical Telescopesand Facility Instruments, budgeted from the Ministry of Finance (MOF) of China and administrated by the Chinese Academy ofSciences (CAS), as well as the National Key R&D Program of China (No. 2017YFA0402700). Additional funding support isprovided by the Science and Technology Facilities Council of the United Kingdom and participating universities in the UnitedKingdom and Canada. We thank all the ALMA and JCMT sta ff supporting this work. P.M.K. acknowledges support from MOST107-2119-M-001-023 and from an Academia Sinica Career Development Award. ZYL is supported in part by NSF AST-1716259and 1815784 and NASA 80NSSC18K1095 and NNX14AB38G. S.T. acknowledges a grant from JSPS KAKENHI Grant Num-bers JP16H07086 and JP18K03703 in support of this work. This work was supported by NAOJ ALMA Scientific Research GrantNumbers 2017-04A. APPENDIX A. SIMULATED AND CONVOLVED MODEL MAPSWe performed imaging simulations of the ALMA observations on the synthetic Stokes
IQU maps from two models, 2.4Idealand 9.6tr-MRN-CR5, using the CASA simulator. These two models show the least and most pinched magnetic field among ourmodels, respectively. The same array configuration as that in the ALMA observations was adopted in the imaging simulations.The synthesized beam in our simulated maps is 0 . (cid:48)(cid:48) × . (cid:48)(cid:48)
47 comparable to that in the observations. We extracted the polarizationorientations from the simulated Stokes Q and U maps and rotated them by 90 ◦ to compare with those from the convolved modelmaps. Figure 12 compares the magnetic field orientations extracted from the simulated and convolved model maps. The numbersof the magnetic field orientations extracted from the simulated model maps are smaller than those from the convolved modelmaps because the extended emission in the outer region along the outflow axis is filtered out in the simulated observations. Formodel 2.4Ideal, there is no significant di ff erence in the magnetic field orientations extracted from the simulated and convolvedmodel maps. For model 9.6tr-MRN-CR5, the magnetic field orientations extracted from the simulated model map are generallyconsistent with those from the convolved model maps. Few segments at o ff sets close to (1 (cid:48)(cid:48) , 0 (cid:48)(cid:48) ) extracted from the simulatedmodel maps are orientated more along the direction of the mid-plane, while those from the convolved model maps are moreperpendicular to the mid-plane. These segments are located near the boundary between the outflow and the inner flatten envelope,where the directions of the magnetic field lines change significantly (e.g., Fig. 4). Nevertheless, in the region where thereare ALMA polarization detections, there is no significant di ff erence between the magnetic field orientations extracted fromthe simulated and convolved model maps for model 9.6tr-MRN-CR5. Therefore, our comparison between the models and theobservations is not a ff ected by the filtering e ff ect of the ALMA observations.2 (a) 2.4ideal (b) 9.6tr-MRN-CR5 simulatedconvolved Figure 12.
Comparison between the magnetic field orientations extracted from the convolved (black segments) and simulated (red segments)model maps for the two cases with the least and most pinched magnetic field, (a) 2.4Ideal and (b) 9.6tr-MRN-CR5. Blue contours delineate theregion where there are ALMA polarization detections. Because of the filtering e ffff