<figure><img src="image_1.png" /></figure>

JCMT BISTRO Survey: Magnetic Fields within the Hub-filament Structure in IC 5146

Jia-Wei Wang<fnref n="1" />©, Shih-Ping Lai<fnref n="1" />©, Chakali Eswaraiah<fnref n="1" />,<fnref n="2" />©, Kate Pattle<fnref n="1" />©, James Di Francesco<fnref n="3" />,<fnref n="4" />©, Doug Johnstone<fnref n="3" />,<fnref n="4" />©, Patrick M. Koch<fnref n="5" />©, Tie Liu<fnref n="6" />,<fnref n="7" />©, Motohide Tamura<fnref n="8" />,<fnref n="9" />,<fnref n="10" />©, Ray S. Furuya<fnref n="11" />,<fnref n="12" />©, Takashi Onaka<fnref n="9" />©, Derek Ward-Thompson<fnref n="13" />©, Archana Soam<fnref n="6" />,<fnref n="14" />©, Kee-Tae Kim<fnref n="6" />©, Chang Won Lee<fnref n="6" />,<fnref n="15" />©, Chin-Fei Lee<fnref n="5" />©, Steve Mairs<fnref n="7" />©, Doris Arzoumanian<fnref n="16" />©, Gwanjeong Kim<fnref n="17" />©, Thiem Hoang<fnref n="6" />,<fnref n="15" />©, Jihye Hwang<fnref n="6" />,<fnref n="15" />©, Sheng-Yuan Liu<fnref n="5" />©, David Berry<fnref n="7" />©, Pierre Bastien<fnref n="18" />©, Tetsuo Hasegawa<fnref n="8" />©, Woojin Kwon<fnref n="6" />,<fnref n="15" />©, Keping Qiu<fnref n="19" />,<fnref n="20" />©, Philippe André<fnref n="21" />©, Yusuke Aso<fnref n="9" />©, Do-Young Byun<fnref n="6" />,<fnref n="15" />©, Huei-Ru Chen<fnref n="1" />,<fnref n="5" />©, Michael C. Chen<fnref n="4" />©, Wen Ping Chen<fnref n="22" />©, Tao-Chung Ching<fnref n="2" />,<fnref n="23" />©, Jungyeon Cho<fnref n="24" />©, Minho Choi<fnref n="6" />©, Antonio Chrysostomou<fnref n="25" />©, Eun Jung Chung<fnref n="6" />©, Simon Coudé<fnref n="26" />©, Yasuo Doi<fnref n="27" />©, C. Darren Dowell<fnref n="28" />©, Emily Drabek-Maunder<fnref n="29" />©, Hao-Yuan Duan<fnref n="1" />©, Stewart P. S. Eyres<fnref n="13" />©, Sam Falle<fnref n="30" />©, Lapo Fanciullo<fnref n="5" />©, Jason Fiege<fnref n="31" />©, Erica Franzmann<fnref n="31" />©, Per Friberg<fnref n="7" />©, Rachel K. Friesen<fnref n="32" />©, Gary Fuller<fnref n="33" />©, Tim Gledhill<fnref n="25" />©, Sarah F. Graves<fnref n="7" />©, Jane S. Greaves<fnref n="29" />©, Matt J. Griffin<fnref n="29" />©, Qilao Gu<fnref n="34" />©, Ilseung Han<fnref n="6" />,<fnref n="15" />©, Jennifer Hatchell<fnref n="35" />©, Saeko S. Hayashi<fnref n="36" />©, Wayne Holland<fnref n="37" />,<fnref n="38" />©, Martin Houde<fnref n="39" />©, Tsuyoshi Inoue<fnref n="16" />©, Shu-ichiro Inutsuka<fnref n="16" />©, Kazunari Iwasaki<fnref n="40" />©, Il-Gyo Jeong<fnref n="6" />©, Yoshihiro Kanamori<fnref n="27" />©, Ji-hyun Kang<fnref n="6" />©, Miju Kang<fnref n="6" />©, Sung-ju Kang<fnref n="6" />©, Akimasa Kataoka⁴¹©, Koji S. Kawabata⁴²,⁴³,⁴⁴©, Francisca Kemper<fnref n="5" />©, Jongsoo Kim<fnref n="6" />,<fnref n="15" />©, Kyoung Hee Kim⁴⁵©, Mi-Ryang Kim<fnref n="6" />©, Shinyoung Kim<fnref n="6" />,<fnref n="15" />©, Jason M. Kirk<fnref n="13" />©, Masato I. N. Kobayashi⁴⁶©, Vera Konyves<fnref n="13" />©, Jungmi Kwon<fnref n="9" />©, Kevin M. Lacaille⁴⁷,⁴⁸©, Hyeseung Lee<fnref n="24" />©, Jeong-Eun Lee⁴⁹©, Sang-Sung Lee<fnref n="6" />,<fnref n="15" />©, Yong-Hee Lee⁴⁹©, Dalei Li⁵⁰©, Di Li<fnref n="2" />,⁵¹©, Hua-bai Li<fnref n="34" />©, Hong-Li Liu⁵²©, Junhao Liu<fnref n="19" />,<fnref n="20" />©, A-Ran Lyo<fnref n="6" />©, Masafumi Matsumura⁵³©, Brenda C. Matthews<fnref n="3" />,<fnref n="4" />©, Gerald H. Moriarty-Schieven<fnref n="3" />©, Tetsuya Nagata⁵⁴©, Fumitaka Nakamura⁴¹,⁵⁵©, Hiroyuki Nakanishi⁵⁶©, Nagayoshi Ohashi<fnref n="36" />©, Geumsook Park<fnref n="6" />©, Harriet Parsons<fnref n="7" />©, Enzo Pascale<fnref n="29" />©, Nicolas Peretto<fnref n="29" />©, Andy Pon<fnref n="39" />©, Tae-Soo Pyo<fnref n="36" />,⁵⁵©, Lei Qian<fnref n="2" />©, Ramprasad Rao<fnref n="5" />©, Mark G. Rawlings<fnref n="7" />©, Brendan Retter<fnref n="29" />©, John Richer⁵⁷,⁵⁸©, Andrew Rigby<fnref n="29" />©, Jean-François Robitaille⁵⁹©, Sarah Sadavoy⁶⁰©, Hiro Saito⁶¹©, Giorgio Savini⁶²©, Anna M. M. Scaife<fnref n="33" />©, Masumichi Seta⁶³©, Hiroko Shinnaga⁵⁶©, Ya-Wen Tang<fnref n="5" />©, Kohji Tomisaka⁴¹,⁵⁵©, Yusuke Tsukamoto⁵⁶©, Sven van Loo⁶⁴©, Hongchi Wang⁶⁵©, Anthony P. Whitworth<fnref n="29" />©, Hsi-Wei Yen<fnref n="5" />©, Hyunju Yoo<fnref n="6" />©, Jinghua Yuan<fnref n="23" />©, Hyeong-Sik Yun⁴⁹©, Tetsuya Zenko⁵⁴©, Chuan-Peng Zhang<fnref n="23" />©, Guoyin Zhang<fnref n="2" />©, Ya-Peng Zhang<fnref n="34" />©, Jianjun Zhou⁵⁰©, and Lei Zhu<fnref n="21" />

<footnote n="1">¹Institute of Astronomy and Department of Physics, National Tsing Hua University, Hsinchu 30013, Taiwan</footnote>

<footnote n="2">²CAS Key Laboratory of FAST, National Astronomical Observatories, Chinese Academy of Sciences, Beijing, People's Republic of China</footnote>

<footnote n="3">³NRC Herzberg Astronomy and Astrophysics, 5071 West Saanich Road, Victoria, BC V9E 2E7, Canada</footnote>

<footnote n="4">⁴Department of Physics and Astronomy, University of Victoria, Victoria, BC V8P 5C2, Canada</footnote>

<footnote n="5">⁵Academia Sinica Institute of Astronomy and Astrophysics, P.O. Box 23-141, Taipei 10617, Taiwan</footnote>

<footnote n="6">⁶Korea Astronomy and Space Science Institute, 776 Daedeokdae-ro, Yuseong-gu, Daejeon 34055, Republic of Korea</footnote>

<footnote n="7">⁷East Asian Observatory, 660 N. A'ohōkū Place, University Park, Hilo, HI 96720, USA</footnote>

<footnote n="8">⁸National Astronomical Observatory of Japan, National Institutes of Natural Sciences, Osawa, Mitaka, Tokyo 181-8588, Japan</footnote>

<footnote n="9">⁹Department of Astronomy, Graduate School of Science, The University of Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo 113-0033, Japan</footnote>

<footnote n="10">¹⁰Astrobiology Center, National Institutes of Natural Sciences, 2-21-1 Osawa, Mitaka, Tokyo 181-8588, Japan</footnote>

<footnote n="11">¹¹Tokushima University, Minami Jousanajima-machi 1-1, Tokushima 770-8502, Japan</footnote>

<footnote n="12">¹²Institute of Liberal Arts and Sciences Tokushima University, Minami Jousanajima-machi 1-1, Tokushima 770-8502, Japan</footnote>

<footnote n="13">¹³Jeremiah Horrocks Institute, University of Central Lancashire, Preston PR1 2HE, UK</footnote>

<footnote n="14">¹⁴SOFIA Science Centre, USRA, NASA Ames Research Centre, MS N232, Moffett Field, CA 94035, USA</footnote>

<footnote n="15">¹⁵University of Science and Technology, Korea, 217 Gajeong-ro, Yuseong-gu, Daejeon 34113, Republic of Korea</footnote>

<footnote n="16">¹⁶Department of Physics, Graduate School of Science, Nagoya University, Furo-cho, Chikusa-ku, Nagoya 464-8602, Japan</footnote>

<footnote n="17">¹⁷Nobeyama Radio Observatory, National Astronomical Observatory of Japan, National Institutes of Natural Sciences, Nobeyama, Minamimaki, Minamisaku, Nagano 384-1305, Japan</footnote>

<footnote n="18">¹⁸Centre de recherche en astrophysique du Québec and Département de physique, Université de Montréal, C.P. 6128, Succ. Centre-ville, Montréal, QC H3C 3J7, Canada</footnote>

<footnote n="19">¹⁹School of Astronomy and Space Science, Nanjing University, 163 Xianlin Avenue, Nanjing 210023, People's Republic of China</footnote>

<footnote n="20">²⁰Key Laboratory of Modern Astronomy and Astrophysics (Nanjing University), Ministry of Education, Nanjing 210023, People's Republic of China</footnote>

<footnote n="21">²¹Laboratoire AIM CEA/DSM-CNRS-Université Paris Diderot, IRFU/Service d'Astrophysique, CEA Saclay, F-91191 Gif-sur-Yvette, France</footnote>

<footnote n="22">²²Institute of Astronomy, National Central University, Chung-Li 32054, Taiwan</footnote>

<footnote n="23">²³National Astronomical Observatories, Chinese Academy of Sciences, A20 Datun Road, Chaoyang District, Beijing 100012, People's Republic of China</footnote>

<footnote n="24">²⁴Department of Astronomy and Space Science, Chungnam National University, 99 Daehak-ro, Yuseong-gu, Daejeon 34134, Republic of Korea</footnote>

<footnote n="25">²⁵School of Physics, Astronomy & Mathematics, University of Hertfordshire, College Lane, Hatfield, Hertfordshire AL10 9AB, UK</footnote>

<footnote n="26">²⁶SOFIA Science Center, Universities Space Research Association, NASA Ames Research Center, M.S. N232-12, Moffett Field, CA 94035, USA</footnote>

<footnote n="27">²⁷Department of Earth Science and Astronomy, Graduate School of Arts and Sciences, The University of Tokyo, 3-8-1 Komaba, Meguro, Tokyo 153-8902, Japan</footnote>

<footnote n="28">²⁸Jet Propulsion Laboratory, M/S 169-506, 4800 Oak Grove Drive, Pasadena, CA 91109, USA</footnote>

<footnote n="29">²⁹School of Physics and Astronomy, Cardiff University, The Parade, Cardiff CF24 3AA, UK</footnote>

<footnote n="30">³⁰Department of Applied Mathematics, University of Leeds, Woodhouse Lane, Leeds LS2 9JT, UK</footnote>

<footnote n="31">³¹Department of Physics and Astronomy, The University of Manitoba, Winnipeg, MB R3T 2N2, Canada</footnote>

<footnote n="32">³²National Radio Astronomy Observatory, 520 Edgemont Road, Charlottesville, VA 22903, USA</footnote>

<footnote n="33">³³Jodrell Bank Centre for Astrophysics, School of Physics and Astronomy, University of Manchester, Oxford Road, Manchester, M13 9PL, UK</footnote>

<footnote n="34">³⁴Department of Physics, The Chinese University of Hong Kong, Shatin, N.T., Hong Kong</footnote>

<footnote n="35">³⁵Physics and Astronomy, University of Exeter, Stocker Road, Exeter EX4 4QL, UK</footnote>

<footnote n="36">³⁶Subaru Telescope, National Astronomical Observatory of Japan, 650 N. A'ohōkū Place, Hilo, HI 96720, USA</footnote>

<footnote n="37">³⁷UK Astronomy Technology Centre, Royal Observatory, Blackford Hill, Edinburgh EH9 3HJ, UK</footnote>

<footnote n="38">³⁸Institute for Astronomy, University of Edinburgh, Royal Observatory, Blackford Hill, Edinburgh EH9 3HJ, UK</footnote>

<footnote n="39">³⁹Department of Physics and Astronomy, The University of Western Ontario, 1151 Richmond Street, London, ON N6A 3K7, Canada</footnote>

<footnote n="40">⁴⁰Department of Environmental Systems Science, Doshisha University, Tatara, Miyakodani 1-3, Kyotanabe, Kyoto 610-0394, Japan</footnote>
41 Division of Theoretical Astronomy, National Astronomical Observatory of Japan, Mitaka, Tokyo 181-8588, Japan

42 Hiroshima Astrophysical Science Center, Hiroshima University, Kagamiyama 1-3-1, Higashi-Hiroshima, Hiroshima 739-8526, Japan

43 Department of Physics, Hiroshima University, Kagamiyama 1-3-1, Higashi-Hiroshima, Hiroshima 739-8526, Japan

44 Core Research for Energetic Universe (CORE-U), Hiroshima University, Kagamiyama 1-3-1, Higashi-Hiroshima, Hiroshima 739-8526, Japan

45 Department of Earth Science Education, Korea National University of Education, 250 Taeseongtabyeon-ro, Grangnae-myeon, Heungdeok-gu, Cheongju-si, Chungbuk 28173, Republic of Korea

46 Department of Earth and Space Science, Graduate School of Science, Osaka University, 1-1 Machikaneyama-cho, Toyonaka, Osaka, Japan

47 Department of Physics and Astronomy, McMaster University, Hamilton, ON L8S 4M1, Canada

48 Department of Physics and Atmospheric Science, Dalhousie University, Halifax, NS B3H 4R2, Canada

49 School of Space Research, Kyung Hee University, 1732 Deogyeong-daero, Giheung-gu, Yongin-si, Gyeonggi-do 17104, Republic of Korea

50 Xinjiang Astronomical Observatory, Chinese Academy of Sciences, 150 Science 1-Street, Urumqi 830011, Xinjiang, People's Republic of China

51 University of Chinese Academy of Sciences, Beijing 100049, People's Republic of China

52 Chinese Academy of Sciences, South America Center for Astrophysics, Camino El Observatorio #1515, Las Condes, Santiago, Chile

53 Faculty of Education and Center for Educational Development and Support, Kagawa University, Saiwai-cho 1-1, Takamatsu, Kagawa, 760-8522, Japan

54 Department of Astronomy, Graduate School of Science, Kyoto University, Sakyo-ku, Kyoto 606-8502, Japan

55 SOKENDAI (The Graduate University for Advanced Studies), Hayama, Kanagawa 240-0193, Japan

56 Amanoga Galaxy Astronomy Research Center (AGARC), Kagoshima University, 1-21-35 Korimoto, Kagoshima 890-0065, Japan

57 Astrophysics Group, Cavendish Laboratory, J J Thomson Avenue, Cambridge CB3 0HE, UK

58 Kavli Institute for Cosmology, Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge CB3 0HA, UK

59 Univ. Grenoble Alpes, CNRS, IPAG, F-38000 Grenoble, France

60 Harvard-Smithsonian Center for Astrophysics, 60 Garden Street, Cambridge, MA 02138, USA

61 Department of Astronomy and Earth Sciences, Tokyo Gakugei University, Koganei, Tokyo 184-8501, Japan

62 OSL, Physics & Astronomy Department, University College London, London WC1E 6BT, UK

63 Department of Physics, School of Science and Technology, Kwansei Gakuin University, 2-1 Gakuen, Sanda, Hyogo 669-1337, Japan

64 School of Physics and Astronomy, University of Leeds, Woodhouse Lane, Leeds LS2 9JT, UK

65 Purple Mountain Observatory, Chinese Academy of Sciences, 2 West Beijing Road, 210008 Nanjing, People's Republic of China

Received 2018 December 13; revised 2019 March 22; accepted 2019 March 25; published 2019 May 1

## Abstract

We present the 850 µm polarization observations toward the IC 5146 filamentary cloud taken using the Submillimetre Common-User Bolometer Array 2 (SCUBA-2) and its associated polarimeter (POL-2), mounted on the James Clerk Maxwell Telescope, as part of the B-fields In STar forming Regions Observations. This work is aimed at revealing the magnetic field morphology within a core-scale ($\lesssim$1.0 pc) hub-filament structure (HFS) located at the end of a parsec-scale filament. To investigate whether the observed polarization traces the magnetic field in the HFS, we analyze the dependence between the observed polarization fraction and total intensity using a Bayesian approach with the polarization fraction described by the Rice likelihood function, which can correctly describe the probability density function of the observed polarization fraction for low signal-to-noise ratio data. We find a power-law dependence between the polarization fraction and total intensity with an index of 0.56 in $A_V \sim 20-300$ mag regions, suggesting that the dust grains in these dense regions can still be aligned with magnetic fields in the IC 5146 regions. Our polarization maps reveal a curved magnetic field, possibly dragged by the contraction along the parsec-scale filament. We further obtain a magnetic field strength of $0.5 \pm 0.2$ mG toward the central hub using the Davis-Chandrasekhar-Fermi method, corresponding to a mass-to-flux criticality of $\sim 1.3 \pm 0.4$ and an Alfvénic Mach number of $<0.6$. These results suggest that gravity and magnetic field are currently of comparable importance in the HFS and that turbulence is less important.

**Key words:** ISM: individual objects (IC 5146) – ISM: magnetic fields – ISM: structure – polarization – radio continuum: ISM – stars: formation

## 1. Introduction

Observations over the last few decades have revealed that stars predominantly form within magnetized and turbulent molecular clouds (Crutcher 2012). How magnetic fields regulate star formation, however, is still poorly understood. Theoretical works have suggested that magnetic fields could be important in supporting molecular clouds, suppressing the star formation rate (e.g., Nakamura & Li 2008; Price & Bate 2008), and removing angular momentum (e.g., Mestel 1985; Mouschovias & Paleologou 1986). Nonetheless, measurements of magnetic field morphologies and strength are still too rare to test these theories (Tamura et al. 1987; Kwon et al. 2015; Tamura & Kwon 2015, and references therein).

Recently, much attention has been drawn to filamentary molecular clouds, which are suggested as the key progenitors of star formation (André et al. 2010). Li et al. (2013) found that the intercloud media magnetic fields, traced by optical polarimetry, were often oriented either parallel or perpendicular

to the filamentary clouds. Based on the recent Planck data, the Planck Collaboration et al. (2016) showed that the relative orientations between magnetic fields and filamentary density structure changed systematically from parallel to filaments in low column density areas to perpendicular in high column density areas, with a switch point at $N_H \sim 10^{21.7}$ cm$^{-2}$. The observed alignment between magnetic fields and filaments is consistent with theoretical works suggesting that magnetic fields play an important role in guiding the gravitational or turbulence-driven contraction and also in supporting filaments from contraction along the filament major axis (e.g., Nakamura & Li 2008; Busquet et al. 2013; Inutsuka et al. 2015).

Within dense filamentary clouds, morphological configurations named “hub-filament structures” (HFSs) are commonly seen. Such structures consist of a central dense hub ($N_H > 10^{22}$ cm$^{-2}$) with several converging filaments surrounding the hub (Myers 2009; Li et al. 2014). The central hubs of HFSs often host most of the star formation in a filamentary cloud and
hence are the potential sites for cluster formation. Observations
found that converging filaments connecting HFSs often have
similar orientations and spacings (Myers 2009). Polarization
observations toward the HFS G14.225-0.506 found that its
converging filaments are perpendicular to the large-scale
magnetic fields (Busquet et al. 2013; Santos et al. 2016). To
explain these features, theoretical works have suggested that
HFSs are formed via layer fragmentation of clouds threaded by
magnetic fields. In this model, the local densest regions
collapse quickly and form dense hubs, and then the surround-
ing material tends to fragment along magnetic field lines and
become parallel layers, because the gravitational instability
grows faster along the magnetic fields (Nagai et al. 1998;
Myers 2009; Van Loo et al. 2014).

Kinematic analyses have shown that the surrounding filaments within HFSs might indeed be infalling material (e.g., Liu et al. 2015; Juárez et al. 2017; Yuan et al. 2018), attributed either to accretion flows attracted by the dense hubs (e.g., Friesen et al. 2013; Kirk et al. 2013) or to gravitationally collapsing filaments (e.g., Pon et al. 2011, 2012). As an example, spiral arm-like converging filaments with significant velocity gradients were revealed in G33.92+0.11 using ALMA, supporting the idea that these filaments were eccentric accretion flows preserving high angular momentum and were previously fragmented from a rotating clump (Liu et al. 2012, 2015). From this point of view, the formation of HFSs is dominated by gravity, and magnetic fields are merely dragged by the accretion flows. Therefore, the magnetic field morphologies are parallel to the converging filaments and different from the larger scale magnetic fields, which have been seen in NGC 6334 V via Submillimeter Array (SMA) polarimetry (Juárez et al. 2017). Because most of the current HFS formation models were based on the large-scale magnetic field morphologies, it is still challenging to explain how the observed infalling features evolve at a smaller scale. Hence, more observations on the scale of HFS are essential to complete evolutionary models.

Submillimeter dust polarization is commonly used to
measure the magnetic field morphology. Whether or not
submillimeter polarization can really trace the dust grains in
dense clouds, however, is still being debated. Current radiative
torque dust alignment (RATs) theory (Lazarian et al. 1997;
Lazarian & Hoang 2007a; Hoang & Lazarian 2009) suggests
that dust grains in high-extinction regions cannot be efficiently
aligned with magnetic fields due to the lack of radiation fields.
Observationally, polarization efficiency (PE), defined as a ratio
of absorption polarization percentage to visual extinction A<sub>v</sub>, is
commonly used to evaluate whether dust grains within clouds
could be aligned. Past observations found that the PE in high-
extinction regions decreases with A<sub>v</sub> by a power-law index of
-1, indicating that the polarization only comes from the
surface of the clouds (e.g., Andersson et al. 2015; Jones et al.
2015), and also support the prediction of the RATs theory.
Nevertheless, some observations do show flatter PE–A<sub>v</sub>
relations (e.g., Jones et al. 2016; Wang et al. 2017), suggesting
that the dust grains in dense regions do contribute to
polarization. It is still unclear what mechanism can efficiently
align dust grains in high-extinction regions and what environ-
ment the mechanism requires. More measurements of PE in
different environments are needed to settle this debate.

The IC 5146 system is a nearby star-forming region in
Cygnus, consisting of an H II region, known as the Cocoon

Nebula, and a long dark cloud extending from the H II region. The distance of the IC 5146 cloud is ambiguous. Harvey et al. (2008) estimated a distance of 950 pc based on the comparison of zero-age main sequence among the Orion Nebula Cluster and the B-type members of IC 5146. In contrast, Lada et al. (1999) derived a distance of $460_{-60}^{+40}$ pc by comparing the number of foreground stars, identified by near-infrared extinction measurements, to those expected from galactic models. Dzib et al. (2018) estimated a distance of $813 \pm 106$ pc based on the Gaia second data release (Gaia Collaboration et al. 2018) parallax measurements toward the embedded young stellar objects (YSOs) within the Cocoon Nebula. In this paper, we assume a default distance of $813 \pm 106$ pc for consistency.

The Herschel Gould Belt Survey (André et al. 2010) revealed a complex network of filaments within the IC 5146 dark clouds (Arzoumanian et al. 2011): several diffuse subfilaments extend from its main filamentary structures, and two HFSs are located at the ends of the main filaments. The main filament is a known active star-forming region, where more than 200 YSOs have been identified with Spitzer (Harvey et al. 2008; Dunham et al. 2015). The variety of filamentary features in the IC 5146 system suggests it as an ideal target for investigating the formation and evolution of these filaments (Johnstone et al. 2017). Wang et al. (2017) (hereafter WLE17) measured the optical and near-infrared starlight polarization across the whole IC 5146 cloud and showed that the large-scale magnetic fields are uniform and perpendicular to the main filaments, suggesting that the large-scale filaments were formed under strong magnetic field conditions. Since the large-scale magnetic fields have been well probed, the IC 5146 system is an excellent target to perform further submillimeter polarimetry to reveal the role of the magnetic field to smaller scales.

In this paper, we report the 850 µm polarization observations toward the brightest HFS in the IC 5146 system taken with Submillimetre Common-User Bolometer Array 2 (SCUBA-2, Holland et al. 2013) and its associated polarimeter (POL-2, Friberg et al. 2016; P. Bastien et al. 2019, in preparation), mounted on the James Clerk Maxwell Telescope (JCMT), as part of the B-fields In STar forming Regions Observations (BISTRO) (Ward-Thompson et al. 2017; Kwon et al. 2018; Pattle et al. 2018; Soam et al. 2018). The target HFS has a physical size less than ~1.0 pc and a total mass of ~100 $M_{\odot}$ (Harvey et al. 2008), lower than the commonly seen parsec-scale HFS such as NGC 1333 or IC 348, and thus hereafter we named our target “core-scale HFS” to distinguish it from other HFSs with much larger physical scale. In Section 2, we address the details of our observations and data reduction. In Section 3, we discuss the magnetic field morphology revealed by the polarization map, and we present an analysis of the dependency between PE and $A_V$, the magnetic field morphology, and magnetic strength. Our interpretations of the observed polarization data are discussed in Section 4. A summary of our conclusions is given in Section 5.

2. Observations

2.1. *Data Acquisition and Reduction Techniques*

Our polarimetric continuum observations toward the
IC 5146 dark cloud system were carried out between 2016
May and 2017 April. The observed field targeted the brightest
HFS located at the eastern end of the IC 5146 main filament, as
shown in Figure 1. We performed 20 sets of 40-minute
<figure><img src="image_2.png" /><figcaption>Figure 1. The IC 5146 field observed in BISTRO overlaid on the Herschel 250 µm image. The blue solid circle represents the field of view of the POL-2 850 µm polarimetry observations, and the dashed circle indicates the inner 3' region with the best sensitivity. The green vectors show the optical and infrared polarization measurements (Wang et al. 2017). The white circle at the bottom right corner shows the Herschel 250 µm FWHM beam size. We note that the Cocoon Nebula is about 1° to the east of this field.</figcaption></figure>

observations toward the IC 5146 region with $\tau_{225 GHz}$ ranging from 0.04 to 0.07.

The POL-2 observations were made using POL-2 DAISY scan mode (Friberg et al. 2016; P. Bastien et al. 2019, in preparation), producing a fully sampled circular region of 11 arcmin diameter. Within the DAISY map, the noise is uniform and lowest in the central 3 arcmin diameter region and increases to the edge of the map. The POL-2 data were simultaneously taken at 450 µm with a resolution of 9.6 arcsec and at 850 µm with a resolution of 14.1 arcsec. The 450 µm data are not reported in this paper since the 450 µm instrumental polarization (IP) model has been only recently commissioned.

The IC 5146 data were reduced in a three-stage process using pol2map,<fnref n="66" /> a script recently added to the SCUBA-2 mapmaking routine SMURF (Berry et al. 2005; Chapin et al. 2013).

In the first stage, the raw bolometer time streams for each observation are converted into separate Stokes Q, U, and I time streams using the process calcqu.

In the second stage, an initial Stokes I map is created from the I time stream from each observation using the iterative map-making routine makemap. For each reduction, areas of astrophysical emission are defined using a signal-to-noise ratio (S/N) based mask determined iteratively by makemap. Areas outside this masked region are set to zero after each iteration until the final iteration of makemap (see Mairs et al. 2015 for a detailed description of the role of masking in SCUBA-2 data reduction). Convergence is reached when successive iterations of the mapmaker produce pixel values that differ by <5% on average. Each map is compared to the first map in the sequence to determine a set of relative pointing corrections. The individual I maps are then coadded to produce an initial I map of the region.

In the third stage, the final Stokes I, Q, and U maps are created. The initial I map described above is used to generate a fixed S/N-based mask for all further iterations of makemap. The pointing corrections determined in the second stage are applied during the map-making process. In this stage, skyloop,

a variant mode of makemap, is invoked. In this mode, rather than each observation being reduced consecutively as is the standard method, one iteration of the mapmaker is performed on each of the observations. At the end of each iteration, all the maps created are coadded. The coadded maps created after successive iterations are compared, and when these coadded maps on average vary by <5% between successive iterations, convergence is reached. Using skyloop typically improves the mapmaker's ability to recover faint extended structure, at the expense of additional memory usage and processing time. The mapmaker was run three times successively to produce the output I, Q, and U maps from their respective time streams. The Q and U data were corrected for IP using the final output I map and the latest IP model (2018 January) (Friberg et al. 2016, 2018).

In all pol2map instances of makemap, the polarized sky background is estimated by doing a principal component analysis of the I, Q, and U time streams to identify components that are common to multiple bolometers. In the first run of makemap (stage 2), the 50 most correlated components are removed at each iteration. In the second run (stage 3), 150 components are removed at each iteration, resulting in smaller changes in the map between iterations and lower noise in the final map.

The output I, Q, and U maps were calibrated in units of Jy beam⁻¹, using a flux conversion factor of 725 Jy pW⁻¹—the standard 850 µm SCUBA-2 flux conversion factor multiplied by 1.35 to account for additional losses due to POL-2 (Dempsey et al. 2013; Friberg et al. 2016).

Finally, a polarization vector catalog was created from the coadded Stokes I, Q, and U maps. To improve the sensitivity, we binned the coadded Stokes I, Q, and U maps into 12" pixels, and the binned data reached rms noise levels of 1.1 mJy beam⁻¹ for Stokes Q and U.

We calculated the polarization fractions and orientations in the 12" pixel map. We debiased the former with the asymptotic estimator (Wardle & Kronberg 1974) as

$$P = \frac{1}{I} \sqrt{(U^2 + Q^2) - \frac{1}{2}(\delta Q^2 + \delta U^2)}, \quad (1)$$

<footnote n="66">⁶⁶ http://starlink.eao.hawaii.edu/docs/sc22.pdf</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 2.** B-field orientation map sampled on a 12" grid shown on the 850 µm dust continuum map, sampled on a 4" grid, of IC 5146 region. The vectors are selected by $I/\delta I > 10$ and $P/\delta P > 2$ and rotated by 90° to represent magnetic field orientations. The yellow and black vectors show the greater than 3σ and 2–3σ polarization detections. The green vectors represent the H-band starlight polarization. The white circle at the top right corner shows the POL-2 850 µm beam size of 14 arcsec. The zoom-in to the red box is shown in Figure 3.</figcaption></figure>

where P is the debiased polarization percentage, and I, Q, U, δI, δQ, and δU are the Stokes I, Q, U, and their uncertainties. The uncertainty of polarization fraction was estimated using

$$ \delta P = \sqrt{\frac{(Q^2\delta Q^2 + U^2\delta U^2)}{I^2(Q^2 + U^2)}} + \frac{\delta I^2(Q^2 + U^2)}{I^4}. \quad (2) $$

The polarization position angle (PA) was calculated as

$$ \mathrm{PA} = \frac{1}{2} \tan^{-1} \left( \frac{U}{Q} \right), \quad (3) $$

and its corresponding uncertainties were estimated using

$$ \delta \mathrm{PA} = \frac{1}{2} \sqrt{\frac{(Q^2 \delta U^2 + U^2 \delta Q^2)}{(Q^2 + U^2)^2}}. \quad (4) $$

The magnetic field orientations used in this paper were assumed to be PA + 90°.

## 2.2. CO Contamination

The SCUBA-2 850 µm waveband covers the wavelength of the CO (J = 3–2) rotational line, and thus our measured continuum flux could be affected by CO line emission (e.g., Drabek et al. 2012; Coudé et al. 2016). Furthermore, the CO (J = 3–2) rotational line is known to be polarized via the Goldreich–Kylafis effect (Goldreich & Kylafis 1981, 1982). For example, the typical polarization fraction of CO (J = 3–2) could be ≲3% in dense clouds and outflows (Ching et al. 2016), calculated using the formulation in Deguchi & Watson (1984)

and Cortes et al. (2005). The polarization angle of CO line is either parallel or perpendicular to the magnetic fields depending on optical depth and the relative angle between magnetic fields and gas velocity fields (Cortes et al. 2005). If a typical polarization fraction of 2% for CO (J = 3–2) is assumed, the polarized intensity from the line would be only 0.02%–0.14% of the total 850 µm flux, which is insignificant compared with the uncertainties of polarization, ≳0.2%–0.5%, in the central hub.

The CO contamination in total intensity might also decrease the observed polarization fraction. Johnstone et al. (2017) calculated the fraction of CO (J = 3–2) line emission to the total flux in the JCMT 850 µm waveband toward several clumps in the IC 5146 system. The fraction of CO (J = 3–2) to total flux in our target region is mostly ~1%–3%, but a higher fraction of ~7% was found in the central hub. Hence, the CO contamination would reduce the measured polarization fraction by a factor of 1%–7%. Nevertheless, this effect is insignificant to our analysis because the S/Ns of our polarization detections are typically only ~2–4.

# 3. Results and Analysis

## 3.1. Magnetic Field Morphology

We show the observed magnetic field orientations traced by POL-2 850 µm polarization, with pixel size of 12", overlaid on the Stokes I map, with pixel size of 4", in Figure 2. We selected the 139 vectors with $I/\delta I > 10$ to ensure that the selected data are associated with the target core-scale HFS. Montier et al. (2015a) suggest that the uncertainty in Stokes I may enhance
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Same as Figure 2, but zoom-in to the HFS region. The vectors are plotted with constant lengths to emphasize their orientations.</figcaption></figure>

the bias in polarization fraction for data with $I/\delta I < 10$, and thus our $I/\delta I > 10$ selection criterion could exclude these biased data. Among these $I/\delta I > 10$ data, 30 of them have $2 < P/\delta P < 3$ and 42 of them have $P/\delta P > 3$. In order to better probe the magnetic fields, we further added the $P/\delta P > 2$ criterion to exclude the samples with higher uncertainties in PA, and the final selected samples have a maximum $\delta$PA of 12.7° and a mean $\delta$PA of 8.5°. Figure 3 shows the zoom-in polarization map toward the HFS and our final selected samples. We note that the CO contamination in Stokes *I* only has insignificant effect on our sample selection. Assuming the CO contamination in total intensity is 7% everywhere, as the worst case, the number of $P/\delta P > 2$ vectors would only decrease to 68 from 72.

The Stokes *I* map shows a central massive clump in which three filaments intersect. The observed morphology is consistent with the typical HFS. The central massive clump hosts ~80% of the total intensity within the system, and thus we recognize the central massive clump as the hub of the HFS. Three filaments are identified extending from the central hub to the north, east, and south. The magnetic field revealed by our polarization map seems to have small angular dispersion but also shows a change of orientations from the north to the south.

To compare the magnetic fields in the observed HFS with the large-scale magnetic fields shown in Figure 1, we plot a histogram of the PAs of the local magnetic fields from POL-2 and WLE17 data within our field of view (diameter of 11') in Figure 4 with a bin size of 10° that is close to our mean $\delta$PA of 8.5°.

The PA histogram of our data shows two major components separated by a dip at 10°. The PA > 10° component has a mean PA of 37° and a PA dispersion of 15°, which is similar to the large-scale magnetic fields (28° ± 21°). In contrast, the PA < 10°

<figure><img src="image_5.png" /><figcaption>**Figure 4.** Histogram of the magnetic field orientations. The colored histogram shows the 850 µm polarimetry data, and the blue histogram represents the optical/infrared data. The bin size of the histograms is set to 10°, similar to our typical uncertainties in PA. The PA = 0° corresponds to north, and PA = +90° is east. Our data show two major components with mean PA of -27° (red) and 37° (green), separated by a dip at ≈10°, and a minor component peaked at -75° (gray). The black solid and dashed vertical lines label the orientation parallel and perpendicular to the large-scale filament, respectively, and the perpendicular orientation is consistent with the dip between the observed two components.</figcaption></figure>

component has a mean PA of -27° and a PA dispersion of 27°. The PA difference of 64° between the two components is much greater than the PA dispersion for large-scale magnetic fields as
<figure><img src="image_6.png" /></figure>

**Figure 5.** The polarization vectors of the three components shown in Figure 4. The green, red, and gray vectors are associated with the components with $10^\circ < \text{PA} < 90^\circ$, $-50^\circ < \text{PA} < 10^\circ$, and $-90^\circ < \text{PA} < -60^\circ$, respectively. The white contours show the H$_2$ column densities of (0.5, 1, 1.5, and 3) × 10$^{21}$ cm$^{-2}$ (Arzoumanian et al. 2011), indicating the morphology of the large-scale filament. The yellow dashed line shows the direction across the intensity peak of the two clumps along the large-scale filament (PA = −73°), which we used to represent the major axis of the large-scale filament. The vectors from the two major components tend to be distributed in the upper and bottom half of the system and are separated by the dashed line, which favors the possibility that the magnetic field is dragged by the large-scale main filament. The vectors from the minor component seem to be randomly distributed over the area, which are probably small-scale structures that we cannot resolve.

well as for our mean observational uncertainty (8°5), suggesting that the observed magnetic field morphology is significantly different from the large-scale magnetic fields.

Figure 5 shows the locations of these two components. To represent the major axis of the main filament, we plotted the yellow dashed line that shows the direction across the intensity peaks of the two clumps along the parsec-scale filament. This major axis has an orientation of −73° and roughly separates the spatial distribution of the two magnetic field orientation components. Within the HFS, the red and green components tend to be distributed in the northern and southern half of the HFS; the tendency, however, is reversed in the western clump. In addition, the orientation perpendicular to the main filament (17°) is also close to the dip between the two components. These features favor the possibility that the magnetic fields in the HFS are curved along the main filament. In contrast, the WLE17 data only show a major peak (28°) in the PA histogram, which is roughly perpendicular to the large-scale filament but with a ≈10° offset in PA. A minor component peaked at ≈−75° is also shown by our data; however, this component is diffusely distributed over the area, and thus more vectors are needed to reveal these structures.

## 3.2. PE

To investigate whether our polarization data trace the dust grains in high-extinction regions, we plot the 850 µm emission polarization fraction $P_{\text{emit}}$ versus $A_V$ in Figure 6. To reveal the complete $P_{\text{emit}}-A_V$ distribution, this figure includes all the data

<figure><img src="image_7.png" /></figure>

**Figure 6.** PE vs. $A_V$. The green, blue, and black points represent the POL-2 data with S/N > 3, between 2 and 3, and <2, respectively. The optical extinction of POL-2 data is derived from total intensities $I_\mu = \tau B_\mu(T)$ with temperatures derived in Arzoumanian et al. (2011) using the Herschel data. The green and blue dashed lines show the best least-squares fitting to the S/N > 3 and 2 < S/N < 3 data with indices of −1.08 and −1.03, respectively. The magenta points are the mean H-band PE, observed across the whole IC 5146 cloud system (Wang et al. 2017). The $P_{\text{emit}}$ (for POL-2) and PE (for H-band) values are shown in the right and left y-axis, respectively. The magenta line represents the best fit for the H-band data (Wang et al. 2017), and the red line shows the prediction from the mean posterior for our data (Section 3.3). These two fitting results are offset by a factor of 48.3 at $A_V$ = 20 mag, due to the wavelength-dependent optical depth of the aligned dust grains, which we use to scale and match the two data sets.

with $I/\delta I > 10$, and the data points are color coded based on their S/N of $P_{\text{emit}}$.

To estimate the $A_V$, we calculated the $\tau_{850 \mu m}$ from the observed 850 µm intensity using $I_{850 \mu m} = \tau_{850 \mu m} B(T_{\text{dust}})$, assuming that the dust emission is optically thin at 850 µm. We used the dust temperature $B(T_{\text{dust}})$ derived in Arzoumanian et al. (2011) via fitting the Herschel data at five wavelengths with a modified blackbody function assuming a dust emissivity index of 2. The dust temperature map and the Stokes I map were both resampled on a 12" grid to match our polarization catalog. The $\tau_{850 \mu m}$ was converted to $A_V$ using the $R_V = 3.1$ extinction curve in Weingartner & Draine (2001). We note that the extinction curve may vary at dense regions due to grain growth. If $R_V$ changes from 3.1 to 5.5 within the observed regions, we would underestimate the $P_{\text{emit}}$ versus $A_V$ slope by 10%.

The $P_{\text{emit}}$ are equivalent to the extinction polarization percentages divided by the optical depth ($P_{\text{ext}}/\tau_\lambda$) in the optically thin case (Andersson et al. 2015) and so are proportional to the PE (defined as $P_{\text{ext}}/A_V$). Thus, the observed $P_{\text{emit}}$ versus $A_V$ slope is equivalent to the PE versus $A_V$ slope.

We further plotted in Figure 6 the PE versus $A_V$ revealed by WLE17 optical and infrared polarization data to show how PE varies in low $A_V$ regions. The PE at 850 µm are in the form of $P_{\text{ext,850}}/\tau_{850}$, and the PE obtained in the H-band are represented by $P_{\text{ext,H}}/A_V$. Thus, a scaling factor $\frac{P_{\text{ext,850}}}{P_{\text{ext,H}}} \cdot \frac{A_V}{\tau_{850}}$ is required to convert the PE at the two wavelengths, which is determined by the unknown dust properties (Andersson et al. 2015). Via matching the fitting results of PE versus $A_V$ relation in the H-band (WLE17) and in the 850 µm bands (described in Section 3.3) at $A_V$ = 20 mag, we found a scaling factor of 48.3,
which we used to match the two data sets. This scaling factor is not a universal constant, as discussed by Jones et al. (2015), and varies with physical conditions in different clouds.

### 3.3. PE-AV Dependence

To determine the $P_{\text{emit}}$ versus $A_V$ slope, the conventional approach is to apply a least-squares power law fit to data selected by an S/N cut in $P_{\text{emit}}$. Following this approach, we fitted the $P/\delta P > 2$ and $>3$ data with a power-law function. The best-fit functions are shown in Figure 6 by blue and green dashed lines, and the best-fit power-law indices are $-1.02 \pm 0.02$ and $-1.08 \pm 0.02$, respectively. Nevertheless, because Figure 6 shows that the $P_{\text{emit}}-A_V$ distribution is significantly truncated by the S/N cut, and also because the best-fit trends are very similar to the truncated boundary, it raises doubts about whether or not the fitting is biased by the sample selection.

We investigated how the sample selection biased on $P/\delta P$ affects the fitting to the $P_{\text{emit}}$ versus $A_V$ distribution by performing Monte Carlo simulations of data sets with underlying $P_{\text{emit}}$ versus $A_V$ function and randomly generated measurement errors in the Appendix. We found that the fitted power index would be dominated by the S/N cut and approaches $-1$ rather than the true underlying value, if the $P_{\text{emit}}$ versus $A_V$ distribution is significantly truncated by the applied $P/\delta P$ selection criteria. Hence, to obtain an unbiased power index, it is recommended to include the low $P/\delta P$ data, so that a complete probability density function (PDF) of $P_{\text{emit}}$ can be recovered. Nevertheless, the use of low $P/\delta P$ data would break the Gaussian PDF assumption (Wardle & Kronberg 1974; Vaillancourt 2006) required for least-squares fit, and therefore we turn to using a Bayesian approach to fit the observed $P_{\text{emit}}$ versus $A_V$ distribution.

We used a Bayesian approach to apply a non-Gaussian PDF and fit the observed $P_{\text{emit}}$ versus $A_V$ trend with a power-law model. The Bayesian statistical framework provides a model fitting tool based on probability theory (see detailed introduction in Jaynes & Bretthorst 2003). The general form of the Bayesian inference is

$$P(\theta|D) = \frac{P(\theta)P(D|\theta)}{P(D)} \qquad (5)$$

or

$$\text{Posterior} = \frac{\text{Prior} \times \text{Likelihood}}{\text{Evidence}}, \qquad (6)$$

where $D$ is the observed data and $\theta$ represents the model parameters. The posterior $P(\theta|D)$ describes the probability of the model parameters matching the given data, which is what we are interested in. The evidence $P(D)$ is the probability of obtaining the data, which mainly serves as a normalization factor for the posterior. The prior $P(\theta)$ serves as the initial guessed probability of the model parameters based on our prior knowledge. The likelihood $P(D|\theta)$ describes how likely it is for a given model parameter set to match the observed data.

Assuming the measurements in Stokes Q and U have similar and Gaussian distributed noise, the PDF of the observed polarization fraction has been known to follow the Rice

distribution (Rice 1945; Wardle & Kronberg 1974; Simmons & Stewart 1985; Quinn 2012)

$$F(P|P_0) = \frac{P}{\sigma_P^2} \exp \left[ -\frac{P^2 + P_0^2}{2\sigma_P^2} \right] I_0 \left( \frac{PP_0}{\sigma_P^2} \right), \qquad (7)$$

where $P$ is the observed polarization fraction, $P_0$ is the real polarization fraction, $\sigma_P$ is the uncertainty in polarization fraction, and $I_0$ is the zeroth-order modified Bessel function. The likelihood function of polarization measurements is defined as

$$L(P_0) = \prod_{i=1}^{n} F(P_n|P_0), \qquad (8)$$

where $P_n$ represents the nth measurement.

To perform the fit to the $P_{\text{emit}}$ versus $A_V$ trend using a Bayesian approach, we assumed the following power-law model such that

$$P_0 = \beta A_V^{-\alpha}, \qquad (9)$$

where $\alpha$ and $\beta$ are the free model parameters and $A_V$ is the observed visual extinction. The uncertainty in the polarization fraction is

$$\sigma_P = \sigma_{QU}/I. \qquad (10)$$

Here, $I$ is the observed total intensity, and the $\sigma_{QU}$ is a free model parameter describing the dispersion in Stokes Q and U, which has contributions from both the instrumental uncertainty and the intrinsic dispersion due to source properties. We then simply used uniform priors within reasonable limits as

$$\begin{align}
P(\alpha) &= \begin{cases} \text{Uniform} & 0 < \alpha < 2 \\ 0 & \text{otherwise} \end{cases} \nonumber \\
P(\beta) &= \begin{cases} \text{Uniform} & 0 < \beta < 400 \\ 0 & \text{otherwise} \end{cases} \nonumber \\
P(\sigma_{QU}) &= \begin{cases} \text{Uniform} & 0 < \sigma_{QU} < 10 \\ 0 & \text{otherwise.} \end{cases} \tag{11}
\end{align}$$

The Bayesian model fitting was performed with the Python Package PyMC3 (Salvatier et al. 2016) via a Markov Chain Monte Carlo method using the Metropolis-Hastings sampling algorithm. The 12 arcsec pixel data were used for the fitting to ensure that each measurement is independent.

The derived posterior of each model parameter is shown in Figure 7, and the 95% highest posterior density (HPD) interval of each parameter is plotted to represent the uncertainty. The 95%, 68%, and 50% confidence regions (CRs) predicted by the posterior distribution in $P_{\text{emit}}$ versus $A_V$ space are shown in Figure 8, assuming a dust temperature of 13 K. Since the error distribution of $P$ is asymmetric and also varies with $P/\sigma_P$ and $A_V$, the predicted $P$ versus $A_V$ is not simply a straight line on a logarithmic scale, even though the input model is a power law. Almost all the data points are well within the 95% CRs predicted by the posterior.

The derived $\alpha$ has a mean value of 0.56 with a 95% confidence interval from 0.22 to 0.83. The $\alpha$ derived by the Bayesian method is less steep than the $\alpha \approx 1.0$ derived from the conventional approach, confirming that the conventional method is biased (see Figure 6). The $\alpha$ range of 0.22–0.83 includes the index of $0.25 \pm 0.06$ obtained from near-infrared polarization data in $A_V$ of 3–20 mag regions (see Wang et al. 2017), and thus no significant difference in PE was found between the
<figure><img src="image_8.png" /><figcaption>**Figure 7.** PDF of the model parameters derived using Bayesian model fitting to our 12 arcsec data. The 95% HPD intervals are plotted to represent the uncertainties. The derived α has a mean of 0.56 and a 95% confidence interval from 0.22 to 0.83. The α value much lower than 1 suggests that the dust grains in $A_V \sim 20-300$ mag can still be aligned with magnetic fields.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 8.** The comparison between the Bayesian posterior prediction and the observations. The green, blue, and black points represent the POL-2 data with $S/N > 3$, 2–3, and <2, respectively. The black line and colored regions show the mean, 95%, 68%, and 50% CRs, predicted by the posterior shown in Figure 7, assuming a dust temperature of 13 K. Since the polarization error distribution is non-Gaussian and changes with $A_V$, the expected mean polarization is not simply a straight line on the logarithmic scale.</figcaption></figure>

$A_V < 20$ mag and $A_V = 20–300$ mag regions (Figure 6). In addition, the fitted $\sigma_{QU}$ of $1.78$ mJy beam$^{-1}$ is greater than our estimated instrumental noise of $1.1$ mJy beam$^{-1}$, which indicates a significant intrinsic dispersion in PE.

The value of $\alpha$ smaller than unity indicates that the extinction polarization fraction ($P_{ext} = \tau P_{emit}$) increases with column density. Since the extinction polarization fraction is defined as $\tanh(\Delta\tau)$ (Jones 1989), where $\Delta\tau$ is the differential optical depth between two polarization directions, the increase of extinction polarization fraction indicates a higher amount of aligned dust grains. Hence, our results suggest that the dust grains in the IC 5146 dense regions can still be aligned with magnetic fields. The $\alpha$ of ~0.5 is also predicted by the simulations based on RATs theory (e.g., Whittet et al. 2008) in low-density regions, where the radiation field is sufficiently strong to align dust grains.

Three possibilities could explain why the dust grains within these dense regions can still be aligned. First, because our target is an active star-forming region, the embedded young stars could be the sources of radiation needed to align the dust grains in dense regions. Second, WLE17 found that the dust grains in IC 5146 have significantly grown from the diffuse interstellar medium. These large dust grains could be aligned more efficiently by the radiation with longer wavelengths,

<figure><img src="image_10.png" /><figcaption>**Figure 9.** The clumps identified by the JCMT Gould Belt survey (Johnstone et al. 2017) overlaid on our polarization map. The image shows the 850 µm continuum emission. The black contours and circles represent the boundaries and the emission peaks of the identified clumps, respectively. The yellow vectors show the orientation of the major axis for each clump. The green diamonds and yellow boxes label the class 0/I and II/III YSOs identified in Dunham et al. (2015), respectively.</figcaption></figure>

which can penetrate the dense regions (Lazarian & Hoang 2007a; Hoang & Lazarian 2009). Third, the mechanical torques due to infalling gas and outflows in the star-forming regions could possibly align the dust grains in the absence of a radiation field (Lazarian & Hoang 2007b; Hoang et al. 2018). These possibilities will be further investigated in upcoming BISTRO papers probing the PE in different environments.

### 3.4. Orientation of Clumps and Magnetic Fields

To investigate whether magnetic fields influence the clump fragmentation within the IC 5146 HFS, we examined the correlation between the magnetic field orientations and the clump morphologies. Based on the JCMT 850 µm Gould Belt Survey data, Johnstone et al. (2017) identified eight sub-millimeter clumps within the regions where we had polarization detections (Figure 9). To represent the orientation of each clump, we used our total intensity map and performed a 2D Gaussian fit to each clump to find the PA of its FWHM major axis. The 2D Gaussian fit had typical orientation uncertainties of ~10°. The obtained clump orientations are listed in Table 1 and plotted over the polarization map in Figure 9.
<table><caption>**Table 1**

Geometric and Polarization Properties of the Clumps</caption><thead><tr><th>ID<sup>a</sup></th><th>Total Mass<sup>a</sup><br>(M<sub>☉</sub>)</th><th>Major Axis<br>(arcsec)</th><th>Minor Axis<br>(arcsec)</th><th>R<sub>eff</sub><br>(arcsec)</th><th>Clump Orientation<br>(deg)</th><th>B-field PA<sub>peak</sub><sup>b</sup><br>(deg)</th><th>σ<sub>NT</sub><br>(km s<sup>-1</sup>)</th><th>α<sub>vir</sub> = <sup>M<sub>vir</sub></sup>/<sub>M<sub>clump</sub></sub><sup>c</sup></th><th>α<sub>vir,B</sub></th></tr></thead><tbody><tr><td>42</td><td>11 ± 3</td><td>51 ± 2</td><td>33 ± 2</td><td>43 ± 2</td><td>14 ± 10</td><td>134.8 ± 28.5</td><td>0.25 ± 0.01</td><td>0.9 ± 0.2</td><td>...</td></tr><tr><td>43</td><td>2.0 ± 0.5</td><td>29 ± 1</td><td>16 ± 1</td><td>22 ± 1</td><td>4 ± 10</td><td>100.1 ± 6.6</td><td>0.12 ± 0.01</td><td>1.4 ± 0.3</td><td>...</td></tr><tr><td>45</td><td>0.77 ± 0.18</td><td>37 ± 3</td><td>16 ± 1</td><td>24 ± 2</td><td>174 ± 10</td><td>0.7 ± 9.5</td><td>0.18 ± 0.02</td><td>5.1 ± 1.0</td><td>...</td></tr><tr><td>46</td><td>6.4 ± 1.6</td><td>26 ± 2</td><td>23 ± 1</td><td>24 ± 2</td><td>25 ± 10</td><td>51.5 ± 19.7</td><td>0.21 ± 0.03</td><td>0.7 ± 0.1</td><td>...</td></tr><tr><td>47</td><td>85 ± 20</td><td>40 ± 2</td><td>36 ± 2</td><td>38 ± 2</td><td>86 ± 10</td><td>21.0 ± 2.7</td><td>0.36 ± 0.12</td><td>0.2 ± 0.1</td><td>0.2–1.0<sup>d</sup></td></tr><tr><td>48</td><td>6.0 ± 1.5</td><td>35 ± 3</td><td>25 ± 1</td><td>30 ± 3</td><td>164 ± 10</td><td>-4.1 ± 1.6</td><td>0.14 ± 0.01</td><td>0.7 ± 0.1</td><td>...</td></tr><tr><td>50</td><td>0.97 ± 0.24</td><td>29 ± 1</td><td>18 ± 1</td><td>23 ± 1</td><td>9 ± 10</td><td>13.6 ± 11.4</td><td>...</td><td>...</td><td>...</td></tr><tr><td>52</td><td>7.6 ± 1.9</td><td>31 ± 2</td><td>17 ± 1</td><td>23 ± 2</td><td>135 ± 10</td><td>64.7 ± 40.7</td><td>0.29 ± 0.01</td><td>0.9 ± 0.2</td><td>...</td></tr><tr><td>53</td><td>1.7 ± 0.4</td><td>34 ± 3</td><td>18 ± 1</td><td>24 ± 3</td><td>140 ± 10</td><td>-19.7 ± 10.1</td><td>...</td><td>...</td><td>...</td></tr></tbody></table>

**Notes.**

<sup>a</sup> The clumps ID and mass were obtained from Johnstone et al. (2017), but the masses were scaled to a distance of 813 ± 106 pc.

<sup>b</sup> Mean magnetic field orientation averaged using the 3 × 3 pixels at the intensity peaks.

<sup>c</sup> The virial masses of clumps were calculated considering the support from thermal pressure and turbulence.

<sup>d</sup> If the inclination angle of the magnetic field with respect to the line of sight is >15°.

To estimate the local magnetic field orientation within each clump, we calculated the mean magnetic field orientations by averaging the data within 3 × 3 pixels at the clump intensity peaks. The size of 3 × 3 pixels is comparable to the typical radius of these clumps (≈20′′–40′′, see Table 1), so the average represents the mean magnetic field orientation over the dense center of these clumps. In addition, The 3 × 3 pixels also provide an estimation of orientation dispersions, which were used as the uncertainties of the averaged orientations; if only one vector was obtained for a given clump, the instrumental uncertainty would be used. To handle the ±180° ambiguity, the mean and dispersion of the magnetic field PAs were calculated in a new coordinate system, where the PA dispersion was minimized, and the calculated results were converted back to the standard coordinate system.

The derived local magnetic field orientations versus the clump orientations are plotted in Figure 10. The comparison between the clump axis and the magnetic field orientations in the clump is limited by the small number of statistics. Nevertheless, there appears to be a tendency that the observed clumps are likely either parallel or perpendicular to the mean magnetic field orientation within ~20°, as shown in Figure 10. The upcoming BISTRO data would provide a much bigger sample set from various systems to statistically confirm this tendency. If the tendency is real, it would suggest that the magnetic fields are a key factor in the fragmentation of these clumps. We note that clumps 43 and 52 only contain two polarization vectors, which are almost perpendicular to each other, and thus the mean magnetic field orientations are not meaningful for these two cases. Because the orientation of these two clumps is still parallel to one of the vectors and perpendicular to the other, these two clumps are still consistent with the tendency.

We further plot in Figure 10 the mean large-scale magnetic field orientation, 28°, obtained from WLE17. Only clump 47 has a magnetic field orientation similar to the large-scale magnetic fields within 20°. All other clumps are aligned either parallel or perpendicular with the local magnetic field and show no significant correlation with the large-scale magnetic fields. Hence, these clumps are more likely formed after the local magnetic fields were distorted by the process of the clump formation.

<figure><img src="image_11.png" /><figcaption>**Figure 10.** The comparison of orientations between clumps and magnetic fields. The clump orientations are obtained from 2D Gaussian fits to the CO-subtracted intensity, and the magnetic field orientations are averaged from the polarization detection within the 3 × 3 pixels at the clump intensity peak. The gray dashed line labels the mean PA of the large-scale magnetic field from WLE17 data. The magenta region represents where the orientation of clumps and magnetic fields are equal, and the green regions show where the orientation of clumps and magnetic fields are perpendicular. All our clumps are close to either the magenta or green regions.</figcaption></figure>

### 3.5. Magnetic Field Strength in IC 5146

The Davis-Chandrasekhar-Fermi (DCF) method (Davis 1951; Chandrasekhar & Fermi 1953) is commonly used to evaluate the magnetic field strength from dust polarizations. The DCF method assumes that turbulent kinematic energy and the turbulent magnetic energy are in equipartition, and hence the magnetic field strength can be estimated using

$$B_{\text{pos}} = Q \sqrt{4\pi\rho} \frac{\sigma_v}{\delta\phi_{\text{intrinsic}}}, \quad (12)$$

where $\delta\phi$ is the intrinsic angular dispersion of the magnetic fields, $\sigma_v$ is the velocity dispersion, $\rho$ is the gas density, and $Q$ is a factor accounting for the complicated magnetic field and inhomogeneous density structure. Ostriker et al. (2001) suggested that $Q = 0.5$ yields a good estimation of magnetic
<table><caption>**Table 2**
Derived Magnetic Field Strength of Clump 47</caption><thead><tr><th>ID</th><th>σ<sub>NT</sub><br>(km s<sup>-1</sup>)</th><th>δφ<br>(deg)</th><th>n<sub>H<sub>2</sub></sub><br>(cm<sup>-3</sup>)</th><th>B<sub>pos</sub><br>(mG)</th><th>λ</th></tr></thead><tbody><tr><td>47</td><td>0.36 ± 0.12</td><td>17.4 ± 0.6</td><td>(9.8 ± 2.4) × 10<sup>5</sup></td><td>0.5 ± 0.2</td><td>1.3 ± 0.4</td></tr></tbody></table>

**Note.** The magnetic field strength estimated using the DCF method. The variables $σ_{NT}$, $δφ$, $n_{H_2}$, $B_{pos}$, and $λ$ represent the velocity dispersion, magnetic field angular dispersion, H₂ volume density, plane of sky magnetic field strength, and mass-to-flux ratio, respectively.

field strength on the plane of sky if magnetic field angular dispersion is <25°.

### 3.5.1. Magnetic Field Angular Dispersion

We used the 12" pixel polarization data to calculate the magnetic field angular dispersion to ensure that all vectors we used are independent measurements. To avoid small number statistics (fewer than 10 vectors), we only perform the angular dispersion estimation using the polarization vectors (45 vectors) within the central hub (clump 47) (Figure 9).

The DCF method requires an estimation of magnetic field distortion caused by turbulence, and the underlying magnetic field geometry might bias the estimation. Thus, we calculated the magnetic field angular dispersion in a local area to avoid the angular dispersion contribution from the large-scale nonuniform magnetic field geometry. Specifically, we selected each 24" × 24" box (i.e., the width of two independent beams) and calculated the mean and the corresponding sum of squared differences (SSD = $\sum_{i=1}^{n}(\bar{PA} - PA_i)^2$) of the PA using the, at most, nine vectors within each box. The SSD from all boxes were averaged with inverse-variance weighting, and the square root of the mean SSD was taken as the observed angular dispersion. Next, the mean instrumental uncertainties ($δφ_{ins}$) of 8°0 were removed from the observed angular dispersion ($δφ_{obs}$) to obtain the intrinsic dispersion ($δφ_{intrinsic}$) via

$$δφ_{\text{intrinsic}}^2 = δφ_{\text{obs}}^2 - δφ_{\text{ins}}^2 \quad (13)$$

With these corrections, the calculated $δφ_{intrinsic}$ for clump 47 is 17.4 ± 0.6.

### 3.5.2. Velocity Dispersion

To estimate the velocity dispersion, we used the C¹⁸O (J = 3–2) spectral data taken by Graham (2008) with the JCMT HARP receiver (Buckle et al. 2009). CO and its isotopomers are well mixed with H₂ and are commonly used to trace the gas kinematics. The C¹⁸O (J = 3–2) line, in particular, is expected to trace gas with volume density up to ~10⁵ cm⁻³ (e.g., Di Francesco et al. 2007), which is comparable to the densities of our target field. In addition, the C¹⁸O (J = 3–2) line is likely optically thin in this field (Graham 2008), so it traces the kinematics of all the gas in the clump. Therefore, we assumed that the C¹⁸O (J = 3–2) line width can well represent the gas velocity dispersion in our observing regions.

The C¹⁸O data reveal, at least, three velocity components within the central hub, which peaked at 3.7, 4.1, and 4.5 km s⁻¹. Because the three components have very similar velocities, and also multiple YSOs have been identified in the central hub by Harvey et al. (2008), the multiple components are possibly structures within the hub, instead of foreground or background components. We performed a multicomponent

Gaussian fit to estimate the C¹⁸O (J = 3–2) line width, using the python package PySpecKit (Ginsburg & Mirocha 2011). We only accepted the fitted Gaussian components with amplitudes >5σ. To estimate the thermal velocity dispersion, we adopted a gas kinematic temperature ($T_{kin}$) of 10 ± 1 K, which is the same as the excitation temperature estimated from the C¹³CO (J = 3–2) line in Graham (2008), leading to $\sqrt{\frac{k_B T_{kin}}{m_{C^{18}O}}} = 0.05 \pm 0.01 \text{ km s}^{-1}$. The thermal velocity dispersions were then removed from the fitted line widths to obtain the nonthermal velocity dispersions via

$$\sigma_{\text{NT}}^2 = \sigma_{\text{obs}}^2 - \frac{k_B T_{\text{kin}}}{m_{\text{C}^{18}\text{O}}}, \quad (14)$$

where $σ_{NT}$ is the nonthermal velocity dispersion, $σ_{obs}$ is the observed C¹⁸O Gaussian line width, and $m_{C¹⁸O}$ is the molecular weight. The inverse-variance weighted mean of the $σ_{obs}$ and $σ_{NT}$ of all velocity components over the central hub was 0.37 and 0.36 km s⁻¹, respectively, and the dispersion of $σ_{obs}$ of 0.12 km s⁻¹ among the central hub was used as its uncertainty.

### 3.5.3. Volume Density

Johnstone et al. (2017) estimated the total mass of the clumps within the IC 5146 cloud using JCMT 850 µm data assuming a distance of 950 pc. The derived total masses of the central hub were scaled down to a distance of 813 ± 106 pc and are listed in Table 1. We assume that the thickness of the hub is equal to the geometric mean of the observed major and minor axis obtained from the 2D Gaussian fit listed in Table 1, and the uncertainty of the thickness is assumed to be the difference between the observed major and minor axis. The mean volume density of the hub is then estimated using its total mass and ellipsoid volume. The calculated H₂ volume densities ($n_{H_2}$) for clump 47 are (9.8 ± 2.4) × 10⁵ cm⁻³. We note that our estimated radius is underestimated due to the unknown inclination angle (i) of the clump, and thus the volume densities we estimated here are only upper limits.

### 3.5.4. Magnetic Field Strength and Mass-to-flux Ratio

Using Equation (12) and the quantities estimated above (Table 2), the plane-of-sky magnetic field strength ($B_{pos}$) is estimated to be 0.5 ± 0.2 mG. To evaluate the relative importance of magnetic fields and gravity in the central hub, we calculate the mass-to-flux critical ratio via

$$\lambda_{\text{obs}} = \frac{(M/\Phi)_{\text{obs}}}{(M/\Phi)_{\text{cri}}}, \quad (15)$$

where the observed mass-to-flux ratio is

$$(M/\Phi)_{\text{obs}} = \frac{\mu m_{\text{H}} N_{\text{H}_2}}{B_{\text{pos}}}, \quad (16)$$
where $\mu = 2.8$ is the mean molecular weight per H$_2$ molecule and the $(M/\Phi)_{\text{cri}}$ is the critical mass-to-flux ratio defined as

$$ (M/\Phi)_{\text{cri}} = \frac{1}{2\pi\sqrt{G}} \qquad (17) $$

(Nakano & Nakamura 1978). Due to the unknown inclination of the clumps, the observed mass-to-flux ratio $\lambda_{\text{obs}}$ is also only an upper limit. Crutcher (2004) suggested that a statistically average factor of one-third could be used to estimate the real mass-to-flux ratio accounting for the random inclinations for an oblate spheroid core, flattened perpendicular to the orientation of the magnetic field. Since we have shown that the central clump is elongated with its major axis perpendicular to the local magnetic field, we adopt a factor of one-third to estimate the mass-to-flux ratio via

$$ \lambda = \frac{1}{3}\lambda_{\text{obs}}. \qquad (18) $$

The estimated mass-to-flux ratio for the central clump is $1.3 \pm 0.4$. The DCF method often tends to overestimate the magnetic field strength, because the effect of integration over the telescope beam and along the line of sight might smooth out part of the magnetic field structure, resulting in an under-estimated angular dispersion (Heitsch et al. 2001; Ostriker et al. 2001; Crutcher 2012). In addition, our target region has a complicated velocity structure, and therefore the measured velocity dispersion might have contributions from gas accretion or contraction motions instead of only isotropic turbulence, also leading to an overestimated magnetic strength. Hence, our estimation of the mass-to-flux ratio only represents a lower limit. A mass-to-flux ratio $\gtrsim 1.0$ suggests that the central hub is supercritical and that magnetic fields and gravity are comparably important at subparsec scales.

### 3.5.5. Angular Dispersion Function

Hildebrand et al. (2009) developed an alternative method to improve the DCF method using a polarization angular dispersion function (ADF) to accurately extract the turbulent component from the polarization data. Houde et al. (2009) further generalized the ADF method by including the effect of signal integration along the thickness of the clouds and over the telescope beam. In this section, we test whether the magnetic field strength estimated using the ADF method is significantly different from our estimation in Section 3.5.4.

The ADF method assumes that the magnetic fields in clouds are combinations of ordered large-scale component $B_0$ and turbulent component $B_t$, and the ratio of these two components determines the intrinsic polarization angular dispersion that

$$ \delta\phi_{\text{intrinsic}} = \left[ \frac{\langle B_t^2 \rangle}{\langle B_0^2 \rangle} \right]^{\frac{1}{2}}, \qquad (19) $$

where $\langle \dots \rangle$ denotes an average. Hence, the DCF equation (Equation (12)) can be written as

$$ B_{\text{pos}} = \sqrt{4\pi\rho} \sigma_v \left[ \frac{\langle B_t^2 \rangle}{\langle B_0^2 \rangle} \right]^{-\frac{1}{2}}. \qquad (20) $$

The detailed derivation given by Hildebrand et al. (2009) and Houde et al. (2009) shows that the ratio of turbulent to

<figure><img src="image_12.png" /><figcaption>Figure 11. The ADF $1 - \langle \cos[\Delta\Phi(\ell)] \rangle$ as a function of the distance $\ell$. The mean $\Delta\Phi(\ell)$ was calculated in bins of 12 arcsec. The green dashed line shows the best fit of Equation (21) to the data.</figcaption></figure>

magnetic energy can be estimated from the ADF using the following equation:

$$ 1 - \langle \cos[\Delta\Phi(\ell)] \rangle \approx \frac{1}{N} \frac{\langle B_t^2 \rangle}{\langle B_0^2 \rangle} (1 - e^{-\ell^2/2(\delta^2+2W^2)}) + a\ell^2, \qquad (21) $$

where $\Delta\Phi(\ell)$ is the difference in the polarization angle measured at two positions separated by a distance $\ell$. The quantities $\delta$ and $a$ are unknown parameters, representing the turbulent correlation length and the first-order Taylor expansion of the large-scale magnetic field structure. The quantity $W$ is the telescope beam radius, which is 6.2 arcsec at 850 µm. The quantity $N$ is the number of turbulent cells observed along the line of sight and within the telescope beam and can be estimated from

$$ N = \Delta' \frac{\delta^2 + 2W^2}{\sqrt{2\pi}\delta^3}, \qquad (22) $$

where $\Delta'$ is the effective cloud thickness, which is assumed to be the clump effective radius. Via fitting the above equation to the observed $1 - \langle \cos[\Delta\Phi(\ell)] \rangle$ versus $\ell$ distribution, the three unknown parameters $\delta$, $\frac{\langle B_t^2 \rangle}{\langle B_0^2 \rangle}$, and $a$ can be derived.

We applied the ADF method to estimate the magnetic field strength in clump 47 using the same selected polarization vectors as in Section 3.5.1. We calculated the $\cos[\Delta\Phi(\ell)]$ and $\ell$ from each pair of the polarization vectors within clump 47, and the results are averaged in bins of width $\ell = 12$ arcsec to estimate the ADF $1 - \langle \cos[\Delta\Phi(\ell)] \rangle$ versus $\ell$. The calculated ADF is plotted in Figure 11. We fitted the observed ADF using Equations (21) and (22), and the best-fitting parameters are shown in Table 3. The obtained turbulent to magnetic energy ratio $\frac{\langle B_t^2 \rangle}{\langle B_0^2 \rangle}$ is $0.27 \pm 0.03$, suggesting that the turbulent magnetic field component is weaker than the ordered large-scale magnetic field. With the previously derived gas velocity and the volume density (Table 2), the magnetic field strength is estimated to be $0.4 \pm 0.2$ mG, which is consistent with our estimation using the DCF method ($0.5 \pm 0.2$ mG) within the uncertainties.
<table><caption>**Table 3**
Magnetic Field Strength of Clump 47 Using the ADF Method</caption><thead><tr><th colspan="3">Fit Result</th><th colspan="2">Derived Quantities</th></tr><tr><th>δ<br>(arcsec)</th><th>(B<sub>2</sub><sup>2</sup>)/(B<sub>0</sub><sup>2</sup>)</th><th>a<br>(arcsec<sup>-2</sup>)</th><th>N</th><th>B<sub>pos</sub><br>(mG)</th></tr></thead><tbody><tr><td>15.1 ± 3.9</td><td>0.27 ± 0.03</td><td>(1.5 ± 0.3) × 10<sup>-5</sup></td><td>1.73 ± 0.4</td><td>0.4 ± 0.2</td></tr></tbody></table>

### 3.5.6. Alfvénic Mach Number

The turbulent Alfvénic Mach number ($M_A$) describes the relative importance of magnetic fields and turbulence, and hence it is a key factor in most cloud evolution models (e.g., Padoan et al. 2001; Nakamura & Li 2008). In the sub-Alfvénic case ($M_A \le 1$), magnetic fields are strong enough to regulate turbulence and cause an organized magnetic field and cloud structure. In the super-Alfvénic case ($M_A > 1$), the turbulence can efficiently change the magnetic field morphology, and the magnetic field morphology is expected to be random.

The Alfvénic Mach number can be estimated from the angular dispersion of the magnetic field if the same assumptions as for the DCF method are made. In doing so, the definition of the Alfvénic Mach number

$$M_A = \frac{\sigma_{NT}}{V_A} = \frac{\sigma_{NT} \sqrt{4\pi\rho}}{B} \qquad (23)$$

can be combined with the equation of the DCF method (Equation (12)), yielding

$$M_A = \frac{\delta\phi \cdot \sin\theta}{Q}, \qquad (24)$$

where $\theta$ is the inclination of the magnetic fields, with respect to the line of sight, so that $B_{pos} = B \sin\theta$. For $Q = 0.5$, the obtained magnetic field angular dispersion $17^\circ$ corresponds to $M_A$ of $0.6 \sin\theta$, and hence the central hub is likely sub-Alfvénic.

### 3.6. Gravitational Stability of Clumps

In this section, we use the virial analysis to investigate whether thermal pressure, turbulence, and magnetic fields are sufficient to support clumps against gravity. If a clump with uniform density is supported by only thermal pressure and turbulence, the virial mass ($M_{vir}$) is

$$M_{\text{vir}} = \frac{5R_{\text{eff}}}{G}(\sigma_{\text{NT}}^2 + c_s^2) \qquad (25)$$

(Bertoldi & McKee 1992; Pillai et al. 2011; Liu et al. 2018), where $R_{eff}$ is the geometric mean of major and minor radius and $c_s = 0.19 \text{ km s}^{-1}$ is the sound speed for a kinematic temperature of 10 K and mean molecular weight. Virial mass is the maximum mass of a stable clump with the support from kinematic and thermal energy. Hence, a clump mass greater than the virial mass, or a virial parameter ($\alpha_{vir} = M_{vir}/M_{clump}$) less than unity, indicates that the clump is gravitationally unstable. We calculate the virial parameter of each clump and list the results in Table 1. Except for clumps 43 and 45, most of the clumps have $\alpha_{vir}$ less than unity, suggesting that thermal pressure and turbulence are insufficient to support the clump

against gravity. Hence, these clumps require additional support from magnetic fields to stop gravitational collapse.

If support from magnetic fields is taken into account, the virial mass of a clump becomes

$$M_{\text{vir}}^{B} = \frac{5R_{\text{eff}}}{G} \left( \sigma_{\text{NT}}^{2} + c_{s}^{2} + \frac{V_{A}^{2}}{6} \right) \qquad (26)$$

(Bertoldi & McKee 1992; Pillai et al. 2011; Liu et al. 2018), where the additional term $\frac{V_A^2}{6}$ stands for the support from magnetic field pressure. We have estimated an Alfvénic Mach number of $0.6 \sin\theta$ for clump 47 in Section 3.5, which corresponds to an Alfvénic velocity of $0.64/\sin\theta$. With support from magnetic fields, the $\alpha_{vir}$ of clump 47 becomes $0.2-1.0$ for a $\theta$ of $15^\circ-90^\circ$ and greater than unity if $\theta < 15^\circ$. Hence, if the direction of magnetic field is not very close to the line of sight, clump 47 is likely gravitationally unstable, which is consistent with the existence of YSOs in the central hub (Harvey et al. 2008). In addition, the presence of YSOs in clump 47 indicates a density structure more complicated than our simple assumption, which could further reduce the virial mass (Sanhueza et al. 2017), and thus this clump could be even more unstable than our estimation.

## 4. Discussion

### 4.1. The Origin of the Core-scale HFS

In Section 3, we show that the magnetic field orientation around the HFS has two major components, tending to be distributed in the northern and southern part of the system. The two components can be explained by either a curved magnetic field or a foreground/background component overlaid on a uniform magnetic field. Nevertheless, because the $C^{18}\text{O}$ ($J = 3-2$) spectral data taken by Graham (2008) show that all components in the HFS are within a narrow velocities range ($\sim 3-5 \text{ km s}^{-1}$), the first possibility is favored, unless the foreground/background component coincidentally has a velocity very similar to the HFS.

The curved magnetic field could be originated by an uniform large-scale magnetic field dragged by the contraction of the large-scale main filament. The dragging along the major axis of the large-scale filament would cause the single peak in the large-scale magnetic field to broaden and split into two peaks, and thus the center of the splitting shown in the histogram ($\sim 15^\circ$) is consistent with the orientation perpendicular to the large-scale filament. In addition, the spatial distribution tendency of the two components can also be explained, because the contraction along the major axis would lead to an axisymmetric pattern. The feature of parsec-scale magnetic fields being perpendicular to filaments but modified by core collapsing at a smaller scale has been found in other filamentary clouds, such as Serpens South (Sugitani et al. 2011), Orion A (Pattle et al. 2017), and W51 (Koch et al. 2018).

Supercritical main filaments are expected to fragment along their major axis and trigger star formation (e.g., André et al. 2010, 2014; Pon et al. 2011; Miettinen 2012; Clarke et al. 2016), which could be a possible origin of the observed HFS. The main filament connected to our observed HFS has a supercritical mass per unit length ($152 M_{\odot} \text{ pc}^{-1}$, Arzoumanian et al. 2011), and the submillimeter clumps identified in
<figure><img src="image_13.png" /><figcaption>**Figure 12.** Schematic for the magnetic field within an accretion flow, modeled in Gomez et al. (2018). The magnetic field is bent by the ram pressure of the flow and eventually reaches a stationary stage where the ram pressure balances the magnetic field tensor. The relative strength of the two forces determines the curvature radius and angle, $R_c$ and $\alpha$, by Equation (27). This figure is adapted from Figure 5 in Gomez et al. (2018) with permission.</figcaption></figure>

Johnstone et al. (2017) also indicate that some filament fragmentation has already taken place.

Some theoretical work suggests that fragmentation of filaments with aspect ratios greater than 5 tends to first begin at their ends, where the edge-driven collapsing mode is more efficient than the homologous collapse mode over the whole filament (Pon et al. 2011, 2012). In contrast, the centralized collapse mode is more important in shorter filaments with high initial density perturbations or no magnetic support (Seifried & Walch 2015). The edge-driven collapsing mode is consistent with the HFSs found at the end of the main filament in the IC 5146 cloud. In addition, Graham (2008) found the gas velocity within the main filament increases from the center to both ends, based on $^{13}$CO and C$^{18}$O line observations. The velocity gradient suggests that the gas within the main filament is likely fragmented toward the two massive HFSs at the ends.

The center-to-ends filament fragmentation picture might seem inconsistent with the observed magnetic field morphology, which shows a pattern of end-to-center contraction. The curved magnetic field morphology, however, might be shaped at an early evolutionary stage, when the filament was still contracting and accumulating mass until its density was sufficiently high to trigger fragmentation. In addition, the global end-to-center contraction and the local center-to-end fragmentation could be occurring simultaneously but at different scales, as suggested by hierarchical gravitational fragmentation models (Gómez & Vázquez-Semadeni 2014; Gomez et al. 2018).

To explore the magnetic field morphology within collapsing clouds, Gomez et al. (2018) simulated molecular clouds undergoing global, multiscale gravitational collapse. In this simulation, the magnetic fields would be dragged by the gravitational contraction but eventually reach a stationary state in which the ram pressure of the flow balances the magnetic tension. Hence, the model predicts a random magnetic field morphology on parsec scales and a U-shaped magnetic field

within the filaments following the equation

$$ \left( \frac{v_l}{v_A} \right)^2 = 2 \sin(2\alpha), \qquad (27) $$

where $v_l$ is the gas velocity along the filaments, $v_A$ is the Alfvénic velocity and the $\alpha$ is the angle between the magnetic field line and the direction perpendicular to the filament, illustrated in Figure 12. Although the predicted large-scale random magnetic field morphology is inconsistent with the uniform magnetic fields shown by the WLE17 data, a U-shaped magnetic field within the filaments has been observed in our POL-2 data, suggesting that this model might be important when the filaments become dense enough.

The observed $\alpha$ is ~30°, estimated by the two components shown in the PA histogram (Figure 4), so $v_l/v_A$ is expected to be 1.3 by Equation (27). We assume that the $v_l$ is approximately the velocity difference along the filament around the central hub. The line-of-sight C$^{18}$O centroid velocity of the central hub (clump 47) is ~4.1 km s$^{-1}$, and the western clump 42 has a centroid velocity of ~3.8 km s$^{-1}$. Hence, the velocity difference along the filament between clumps 42 and 47 is 0.3/ cos $\phi$ km s$^{-1}$, where $\phi$ is the inclination angle of the filament with respect to the line of sight. With the Alfvénic velocity of 0.64/ sin $\theta$ estimated in Section 3.5, the observed $v_l/v_A$ is ~0.5 sin $\theta$/ cos $\phi$. Due to the unknown inclination angle, we can only speculate that the model expectation would be correct if the filament is nearly perpendicular to the line of sight ($\phi > 67°$).

Based on our observed magnetic field morphologies, we propose a three-stage scenario to explain the origin of the observed HFS, illustrated in Figure 13. In the first stage, the large-scale magnetically subcritical filaments are first formed with dynamically important magnetic fields as described in strong magnetic field filament formation models (e.g., Nakamura & Li 2008; Van Loo et al. 2014), and these filaments appear perpendicular to a uniform large-scale
<figure><img src="image_14.png" /><figcaption>**Figure 13.** A cartoon illustrates the possible formation scenario of the core-scale HFSs. (a) The parsec-scale filaments first form via the contraction and fragmentation along magnetic fields. (b) As the parsec-scale filaments become magnetically and thermally supercritical, the filaments fragment along their major axis, and the most massive components form at the end of filaments. At the same time, the magnetic fields are dragged by the filament contraction. (c) The massive fragment at the end of the parsec-scale filament further fragments along the curved magnetic fields and forms a core-scale HFS with orientation parallel or perpendicular to the local magnetic field instead of the primordial field.</figcaption></figure>

magnetic field, as revealed by the WLE17 data. In the second stage, the large-scale filaments accumulate mass via accretion along magnetic field lines or filament mergers (e.g., Li et al. 2010; André et al. 2014) and eventually become magnetically and thermally supercritical. The contraction of supercritical filaments would bend the uniform primordial magnetic fields, similar to the case in Orion A (Pattle et al. 2017). In the third stage, the dense clumps within filaments, often at the ends of filaments, would tend to fragment along magnetic fields and form second-generation filaments with hub-filament morphologies, because density perturbations parallel to the magnetic fields grow faster than those perpendicular to the fields (e.g., Nagai et al. 1998; Van Loo et al. 2014). The collapse of the cores within the second-generation filaments is also regulated by the bent magnetic fields, and so the cores are oriented either parallel or perpendicular to local magnetic fields, as shown in Figure 10, and are less correlated with the primordial magnetic field.

## 4.2. The Alignment between Local Magnetic Fields and Clumps

Stars form predominantly from high column density filaments (André et al. 2010). Although most filaments are either oriented parallel or perpendicular to the large-scale magnetic fields (Li et al. 2014; Planck Collaboration et al. 2016), only a few young stars have been observed with hourglass magnetic field morphologies, which favor a star formation scenario where the core collapse is regulated by strong magnetic fields (e.g., Girart et al. 2006; Rao et al. 2009; Tang et al. 2009). As a counterexample, ALMA polarization observations toward the embedded source Ser-emb 8 show chaotic magnetic fields (Hull et al. 2017), indicating that this star was formed under weak magnetic field conditions. This difference poses the question of how physical scales and environments generally determine the role of magnetic fields in star formation.

To address the role of magnetic fields in star formation, the SMA polarization survey toward massive cores (Zhang et al. 2014) revealed that magnetic fields on the core scale (0.1–0.01 pc) are mostly either parallel or perpendicular to the magnetic fields on the parsec scales. Li et al. (2015) further analyzed the magnetic field morphologies in NGC 6334 on the 100–0.01 pc scales and found that local magnetic fields on all these scales are either parallel or perpendicular to the local cloud elongation. Both these results suggest that magnetic fields are dynamically important during the collapse and fragmentation of clouds, possibly guiding the contraction of filaments and cores. Koch et al. (2014) further used a large

sample set (50 sources) to examine the bimodal distribution of the relative orientation between the magnetic fields and the density structures and found that the distribution is more scattered than those in previous surveys, although a bimodal distribution cannot be ruled out.

In Section 3.4, we find a tendency toward clumps within the observed HFS having orientations parallel or perpendicular to the local magnetic fields (at the 0.1–0.01 pc scale). The local magnetic fields in many of these clumps, however, have orientations 30°–60° different from the parsec-scale magnetic field, which is inconsistent with the findings of Zhang et al. (2014) and Li et al. (2015). The inconsistent cases are mainly clumps within the extending filaments, which follow the orientation of the curved magnetic fields (see Section 4.1). These clumps are much fainter than those in the central hub, which possibly explains why they were missed in previous surveys. Nevertheless, because we still found the orientations of these clumps to be well coupled with the host filaments and the local magnetic fields, our results support the idea that magnetic fields are important in regulating the core and filament collapse on spatial scales of 0.01–0.1 pc.

## 5. Conclusions

This paper presents the first-look results of SCUBA-2/POL-2 observations at 850 µm toward the IC 5146 filamentary clouds as part of the BISTRO project. Our observations reveal the magnetic field morphology within a core-scale HFS located at the end of a parsec-scale filament. From the analysis of these data, we find the following:

1. The observed polarization fraction is found to vary with Stokes *I* following a power law with an index of ≈0.56, which suggests that dust grains in this $A_V$ ~ 20–300 mag range can be still aligned with magnetic fields.

2. The observed polarization map shows that the magnetic field of the HFS on core scales (~0.05–0.5 pc) is more organized than random. The core-scale magnetic field is likely inherited from a larger scale magnetic field that has been dragged by contraction along the parsec-scale filament.

3. The submillimeter clumps within the observed core-scale HFS tend to be aligned with local magnetic fields, that is, they are oriented within 20° of being either parallel or perpendicular to the local magnetic field direction. This trend may suggest that the core-scale HFS formed after the parsec-scale filament became magnetically supercritical and that the magnetic fields have been
dynamically important during the formation and the following evolution of the core-scale HFS.

4. We propose a scenario to explain the formation of the core-scale HFS: the parsec-scale filaments first form under a strong and uniform magnetic field and start to fragment and locally bend the magnetic field as they becomes magnetically supercritical. The massive clump, formed at the end of the parsec-scale filament, further fragments under the strong magnetic fields and becomes the core-scale HFS.

5. Using the DCF method, the magnetic strength within the central hub is estimated to be 0.5 ± 0.2 mG, and the mass-to-flux ratio is 1.3 ± 0.4 for D = 813 pc. The Alfvénic Mach number estimated using the magnetic field angular dispersion is <0.6. These estimates suggest that gravity and magnetic fields are comparably important in the current core-scale HFS and that turbulence is less important.

The James Clerk Maxwell Telescope is operated by the East Asian Observatory on behalf of the National Astronomical Observatory of Japan; Academia Sinica Institute of Astronomy and Astrophysics; the Korea Astronomy and Space Science Institute; the Operation, Maintenance and Upgrading Fund for Astronomical Telescopes and Facility Instruments, funded by the Ministry of Finance (MOF) of China and administrated by the Chinese Academy of Sciences (CAS), as well as the National Key R&D Program of China (No. 2017YFA0402700). Additional funding support is provided by the Science and Technology Facilities Council of the United Kingdom and participating universities in the United Kingdom and Canada. Additional funds for the construction of SCUBA-2 and POL-2 were provided by the Canada Foundation for Innovation. The Starlink software (Currie et al. 2014) is currently supported by the East Asian Observatory, JCMT Project Code M16AL004. The authors wish to recognize and acknowledge the very significant cultural role and reverence that the summit of Maunakea has always had within the indigenous Hawaiian community. We are most fortunate to have the opportunity to conduct observations from this mountain. J.W.W., S.P.L., K.P., and C.E. are thankful for the support from the Ministry of Science and Technology (MOST) of Taiwan through the grants MOST 105-2119-M-007-021-MY3, 105-2119-M-007-024, and 106-2119-M-007-021-MY3. J.W.W. is a University Consortium of ALMA–Taiwan (UCAT) graduate fellows supported by the MOST of Taiwan through the grant MOST 105-2119-M-007-022-MY3. M.K. was supported by the Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Science, ICT & Future Planning (No. NRF-2015R1C1A1A01052160). C.W.L. is supported by the Basic Science Research Program through the NRF, funded by the Ministry of Education, Science and Technology (NRF-2016R1A2B4012593). W.K. was supported by the Basic Science Research Program through the NRF (NRF-2016R1C1B2013642). D.J. is supported by the National Research Council of Canada and by an NSERC Discovery Grant.

Facility: JCMT.

Software: Aplpy (Robitaille & Bressert 2012), Astropy (Astropy Collaboration et al. 2013), NumPy (van der Walt et al. 2011), PyMC3 (Salvatier et al. 2016), PySpecKit (Ginsburg & Mirocha 2011), SciPy (Jones et al. 2001), Starlink (Currie et al. 2014).

## Appendix

### Bias on Determination of *P* versus *I* Relation

In this section, we use the simulated data to investigate the two possible sources of bias in the estimation of the *P*–*I* relation. The first bias source is the selection criteria *P*/*δP* > 2. Since the uncertainty in the polarization fraction has the following dependence on total intensity of

$$ \delta P \approx \frac{\delta Q}{I} \propto I^{-1}, \qquad (28) $$

assuming $\delta Q \approx \delta U$, the selection criteria would truncate the observed *P*–*I* distribution by a boundary of $P = 2\delta P \propto I^{-1}$. Hence, the selection criteria could cause an artificial trend of $P \propto I^{-1}$ to arise, leading to the erroneous conclusion that the dust cannot be aligned with magnetic fields. This bias source could be avoided simply by including the low S/N data.

The second bias source is the non-Gaussian PDF of the observed polarization fraction. As shown in Section 3.3, the PDF of the observed polarization fraction follows the Rice distribution (Equation (7)). The Rice distribution can be approximated as a Gaussian when $P/\delta P \gtrsim 4$, but it becomes more asymmetric as $P/\delta P$ decreases (Vaillancourt 2006). The misuse of a Gaussian PDF on Rice-distributed data would cause the polarization fraction to be overestimated. Furthermore, the bias is anticorrelated with the S/N of *P*, as well as *I*, and steepens the measured slope. The commonly used debiasing methods, for example, the asymptotic estimator, however, could help remove the bias in the high *P*/*δP* domain, but the PDF of the debiased polarization fraction would be still non-Gaussian (Montier et al. 2015b). To avoid this bias source, using an appropriate PDF, instead of a Gaussian assumption, to analyze the observed polarization fraction is recommended.

We performed Monte Carlo simulations to generate a set of polarization data with a given *P*–*A*<sub>V</sub> relation, and tested how the measured relation is affected by the bias sources discussed above. To simulate the observed polarization data, we first generated a 10,000-element set of *A*<sub>V</sub> values distributed uniformly in logarithm where the polarization fraction of each element was determined by the underlying power law

$$ P \propto A_V^{-0.7}. \qquad (29) $$

For simplicity, we assumed all polarization vectors have PAs of 0°, and we calculated the Stokes *Q* and *U* values for each pair of (*A*<sub>V</sub>, *P*). Here we directly used *A*<sub>V</sub> magnitude as the intensity unit. We also added Gaussian-distributed noise with $\sigma = 1.5 A_V$ mag to both Stokes *Q* and *U* and calculated the debiased polarization fraction from the noise-included Stokes *Q* and *U*. The simulated *P* versus *A*<sub>V</sub> distribution is plotted in Figure 14(a). A least-squares power-law fit to the simulated data returns an index of 0.725 ± 0.002, similar to our input value.

To investigate how the first bias source (selection criteria) affects the *P*–*A*<sub>V</sub> relation, we applied selection criteria *P*/*δP* > 2 and *P*/*δP* > 3 to the simulated data, as shown in Figure 14(b). It is clearly shown that the *P*–*A*<sub>V</sub> distribution is truncated by the applied selection criteria, and the least-squares fits here return power-law indices of 0.796 ± 0.002 and 0.850 ± 0.003 for the data selected by *P*/*δP* > 2 and *P*/*δP* > 3, respectively, suggesting that such selection criteria would significantly steepen the measured slope.
<figure><img src="image_15.png" /><figcaption>**Figure 14.** Monte Carlo simulations to show the possible bias for fitting the P–I relation using $P \propto I^{-\alpha}$. (a) The simulated data assuming $P \propto I^{-0.7}$ with moderate noise of 1.5 mag in Stokes Q and U. The least-squares fitting could recover the $\alpha$ of ~0.7. (b) Sample selection criteria of $P/\delta P > 2$ (blue) and $P/\delta P > 3$ (green) were applied to the simulated data. The derived $\alpha$ of the selected data using least-squares fitting become much higher than those of the input model. (c) The green points represent the simulated data with 5 times higher noise. The derived $\alpha$ for the high noise set using least-squares fitting is almost 1.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Figure 15.** The PDF of model parameters obtained from the Bayesian fitting to the simulated data with high noise (Figure 14 green points). The mean values and 95% HPD regions are labeled in each panel. All the results are consistent with our input model.</figcaption></figure>

To test how the second bias source (non-Gaussianity when $P/\delta P$ is low) affects the P–A<sub>V</sub> relation, we generated another simulated data set with 5 times higher noise in Stokes Q and U than our original data set. Figure 14(c) shows the P–A<sub>V</sub> distribution from the simulated data with both original and higher noise. The simulated data with higher noise show a steeper trend, which mainly results from the positive bias in the low A<sub>V</sub> or low P/δP regimes. A power-law index of 0.950 ± 0.006 for the high noise set was obtained by the least-squares fitting, significantly higher than the 0.725 ± 0.002 derived from the original set. These examples show that the two bias sources both could steepen the measured P–I relation. We will further explore how significant the effects are in a much wider parameter space in our forthcoming paper, K. Pattle et al. (2019, in preparation).

To test whether our Bayesian model, as described in Section 3.3, can avoid the bias due to the non-Gaussianity, we fit the same higher noise simulated data set with our Bayesian model. The derived PDFs of each model parameter are shown in Figure 15. All the derived model parameters are consistent with our input values, and the mean posterior suggests an $\alpha$ of 0.75, which is much more accurate than the $\alpha$ of 0.95 derived from least-squares fitting. In addition, the possibility of an asymmetric PDF is considered in the Bayesian model, and hence it provides a more realistic uncertainty estimation.

## ORCID iDs

Jia-Wei Wang https://orcid.org/0000-0002-6668-974X

Shih-Ping Lai https://orcid.org/0000-0001-5522-486X

Chakali Eswaraiah https://orcid.org/0000-0003-4761-6139

Kate Pattle https://orcid.org/0000-0002-8557-3582

James Di Francesco https://orcid.org/0000-0002-9289-2450

Doug Johnstone https://orcid.org/0000-0002-6773-459X

Patrick M. Koch https://orcid.org/0000-0003-2777-5861

Tie Liu https://orcid.org/0000-0002-5286-2564

Motohide Tamura https://orcid.org/0000-0002-6510-0681

Ray S. Furuya https://orcid.org/0000-0003-0646-8782

Takashi Onaka https://orcid.org/0000-0002-8234-6747

Derek Ward-Thompson https://orcid.org/0000-0003-1140-2761

Archana Soam https://orcid.org/0000-0002-6386-2906

Kee-Tae Kim https://orcid.org/0000-0003-2412-7092

Chang Won Lee https://orcid.org/0000-0002-3179-6334

Chin-Fei Lee https://orcid.org/0000-0002-3024-5864

Steve Mairs https://orcid.org/0000-0002-6956-0730

Gwanjeong Kim https://orcid.org/0000-0003-2011-8172

Thiem Hoang https://orcid.org/0000-0003-2017-0982

Sheng-Yuan Liu https://orcid.org/0000-0003-4603-7119

Pierre Bastien https://orcid.org/0000-0002-0794-3859

Woojin Kwon https://orcid.org/0000-0003-4022-4132

Keping Qiu https://orcid.org/0000-0002-5093-5088

Yusuke Aso https://orcid.org/0000-0002-8238-7709
Do-Young Byun https://orcid.org/0000-0003-1157-4109

Huei-Ru Chen https://orcid.org/0000-0002-9774-1846

Wen Ping Chen https://orcid.org/0000-0003-0262-272X

Tao-Chung Ching https://orcid.org/0000-0001-8516-2532

Jungyeon Cho https://orcid.org/0000-0003-1725-4376

Antonio Chrysostomou https://orcid.org/0000-0002-9583-8644

Eun Jung Chung https://orcid.org/0000-0003-0014-1527

Simon Coudé https://orcid.org/0000-0002-0859-0805

Stewart P. S. Eyres https://orcid.org/0000-0002-6663-7675

Rachel K. Friesen https://orcid.org/0000-0001-7594-8128

Gary Fuller https://orcid.org/0000-0001-8509-1818

Tim Gledhill https://orcid.org/0000-0002-2859-4600

Sarah F. Graves https://orcid.org/0000-0001-9361-5781

Jane S. Greaves https://orcid.org/0000-0002-3133-413X

Qilao Gu https://orcid.org/0000-0002-2826-1902

Jennifer Hatchell https://orcid.org/0000-0002-4870-2760

Martin Houde https://orcid.org/0000-0003-4420-8674

Shu-ichiro Inutsuka https://orcid.org/0000-0003-4366-6518

Kazunari Iwasaki https://orcid.org/0000-0002-2707-7548

Ji-hyun Kang https://orcid.org/0000-0001-7379-6263

Miju Kang https://orcid.org/0000-0002-5016-050X

Sung-ju Kang https://orcid.org/0000-0002-5004-7216

Akimasa Kataoka https://orcid.org/0000-0003-4562-4119

Koji S. Kawabata https://orcid.org/0000-0001-6099-9539

Francisca Kemper https://orcid.org/0000-0003-2743-8240

Kyoung Hee Kim https://orcid.org/0000-0001-9597-7196

Mi-Ryang Kim https://orcid.org/0000-0002-1408-7747

Jason M. Kirk https://orcid.org/0000-0002-4552-7477

Masato I. N. Kobayashi https://orcid.org/0000-0003-3990-1204

Jungmi Kwon https://orcid.org/0000-0003-2815-7774

Kevin M. Lacaille https://orcid.org/0000-0001-9870-5663

Jeong-Eun Lee https://orcid.org/0000-0003-3119-2087

Sang-Sung Lee https://orcid.org/0000-0002-6269-594X

Di Li https://orcid.org/0000-0003-3010-7661

Hua-bai Li https://orcid.org/0000-0003-2641-9240

Junhao Liu https://orcid.org/0000-0002-4774-2998

A-Ran Lyo https://orcid.org/0000-0002-9907-8427

Masafumi Matsumura https://orcid.org/0000-0002-6906-0103

Brenda C. Matthews https://orcid.org/0000-0003-3017-9577

Gerald H. Moriarty-Schieven https://orcid.org/0000-0002-0393-7822

Fumitaka Nakamura https://orcid.org/0000-0001-5431-2294

Nagayoshi Ohashi https://orcid.org/0000-0003-0998-5064

Harriet Parsons https://orcid.org/0000-0002-6327-3423

Andy Pon https://orcid.org/0000-0003-4612-1812

Tae-Soo Pyo https://orcid.org/0000-0002-3273-0804

Lei Qian https://orcid.org/0000-0003-0597-0957

Ramprasad Rao https://orcid.org/0000-0002-1407-7944

Mark G. Rawlings https://orcid.org/0000-0002-6529-202X

John Richer https://orcid.org/0000-0002-9693-6860

Sarah Sadavoy https://orcid.org/0000-0001-7474-6874

Anna M. M. Scaife https://orcid.org/0000-0002-5364-2301

Hiroko Shinnaga https://orcid.org/0000-0001-9407-6775

Ya-Wen Tang https://orcid.org/0000-0002-0675-276X

Kohji Tomisaka https://orcid.org/0000-0003-2726-0892

Yusuke Tsukamoto https://orcid.org/0000-0001-6738-676X

Sven van Loo https://orcid.org/0000-0003-4746-8500

Hongchi Wang https://orcid.org/0000-0003-0746-7968

Hsi-Wei Yen https://orcid.org/0000-0003-1412-893X

Hyunju Yoo https://orcid.org/0000-0002-8578-1728

Jinghua Yuan https://orcid.org/0000-0001-8060-3538

Chuan-Peng Zhang https://orcid.org/0000-0002-4428-3183

## References

Andersson, B.-G., Lazarian, A., & Vaillancourt, J. E. 2015, ARA&A, 53, 501

André, P., Di Francesco, J., Ward-Thompson, D., et al. 2014, in Protostars and Planets VI, ed. H. Beuther et al. (Tucson, AZ: Univ. Arizona Press), 27

André, P., Men'shchikov, A., Bontemps, S., et al. 2010, A&A, 518, L102

Arzoumanian, D., André, P., Didelon, P., et al. 2011, A&A, 529, L6

Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33

Berry, D. S., Gledhill, T. M., Greaves, J. S., & Jenness, T. 2005, in ASP Conf. Ser. 343, Astronomical Polarimetry: Current Status and Future Directions, ed. C. Adamson (San Francisco, CA: ASP), 71

Bertoldi, F., & McKee, C. F. 1992, ApJ, 395, 140

Buckle, J. V., Hills, R. E., Smith, H., et al. 2009, MNRAS, 399, 1026

Busquet, G., Zhang, Q., Palau, A., et al. 2013, ApJL, 764, L26

Chandrasekhar, S., & Fermi, E. 1953, ApJ, 118, 113

Chapin, E. L., Berry, D. S., Gibb, A. G., et al. 2013, MNRAS, 430, 2545

Ching, T.-C., Lai, S.-P., Zhang, Q., et al. 2016, ApJ, 819, 159

Clarke, S. D., Whitworth, A. P., & Hubber, D. A. 2016, MNRAS, 458, 319

Cortes, P. C., Crutcher, R. M., & Watson, W. D. 2005, ApJ, 628, 780

Coudé, S., Bastien, P., Kirk, H., et al. 2016, MNRAS, 457, 2139

Crutcher, R. 2012, ARA&A, 50, 29

Crutcher, R. M. 2004, Ap&SS, 292, 225

Currie, M. J., Berry, D. S., Jenness, T., et al. 2014, in ASP Conf. Ser. 485, Astronomical Data Analysis Software and Systems XXIII, ed. N. Manset & P. Forshay (San Francisco, CA: ASP), 391

Davis, L. 1951, PhRv, 81, 890

Deguchi, S., & Watson, W. D. 1984, ApJ, 285, 126

Dempsey, J. T., Friberg, P., Jenness, T., et al. 2013, MNRAS, 430, 2534

Di Francesco, J., Evans, N. J., II, Caselli, P., et al. 2007, in Protostars and Planets V, ed. B. Reipurth, D. Jewitt, & K. Keil (Tucson, AZ: Univ. Arizona Press), 17

Drabek, E., Hatchell, J., Friberg, P., et al. 2012, MNRAS, 426, 23

Dunham, M. M., Allen, L. E., Evans, N. J., II, et al. 2015, ApJS, 220, 11

Dzib, S. A., Loinard, L., Ortiz-León, G. N., Rodríguez, L. F., & Galli, P. A. B. 2018, ApJ, 867, 151

Friberg, P., Bastien, P., Berry, D., et al. 2016, Proc. SPIE, 9914, 991403

Friberg, P., Berry, D., Savini, G., et al. 2018, Proc. SPIE, 10708, 107083M

Friesen, R. K., Medeiros, L., Schnee, S., et al. 2013, MNRAS, 436, 1513

Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2018, A&A, 616, A1

Ginsburg, A., & Mirocha, J. 2011, PySpecKit: Python Spectroscopic Toolkit, Astrophysics Source Code Library, ascl:1109.001

Girart, J. M., Rao, R., & Marrone, D. P. 2006, Sci, 313, 812

Goldreich, P., & Kylafis, N. D. 1981, ApJL, 243, L75

Goldreich, P., & Kylafis, N. D. 1982, ApJ, 253, 606

Gómez, G. C., & Vázquez-Semadeni, E. 2014, ApJ, 791, 124

Gomez, G. C., Vázquez-Semadeni, E., & Zamora-Aviles, M. 2018, MNRAS, 480, 2939

Graham, S. B. 2008, PhD thesis, Jesus College, Cambridge

Harvey, P. M., Huard, T. L., Jørgensen, J. K., et al. 2008, ApJ, 680, 495

Heitsch, F., Zweibel, E. G., Mac Low, M.-M., Li, P., & Norman, M. L. 2001, ApJ, 561, 800

Hildebrand, R. H., Kirby, L., Dotson, J. L., Houde, M., & Vaillancourt, J. E. 2009, ApJ, 696, 567

Hoang, T., Cho, J., & Lazarian, A. 2018, ApJ, 852, 129

Hoang, T., & Lazarian, A. 2009, ApJ, 695, 1457

Holland, W. S., Bintley, D., Chapin, E. L., et al. 2013, MNRAS, 430, 2513

Houde, M., Vaillancourt, J. E., Hildebrand, R. H., Chitszazzadeh, S., & Kirby, L. 2009, ApJ, 706, 1504

Hull, C. L. H., Mocz, P., Burkhart, B., et al. 2017, ApJL, 842, L9

Inutsuka, S., Inoue, T., Iwasaki, K., & Hosokawa, T., 2015, A&A, 580, A49

Jaynes, E. T., & Bretthorst, G. L. (ed.) 2003, Probability Theory (Cambridge: Cambridge Univ. Press), 758

Johnstone, D., Ciccone, S., Kirk, H., et al. 2017, ApJ, 836, 132

Jones, E., Oliphant, E., Peterson, P., et al. 2001, SciPy: Open Source Scientific Tools for Python, http://www.scipy.org/
Jones, T. J. 1989, ApJ, 346, 728

Jones, T. J., Bagley, M., Krejny, M., Andersson, B.-G., & Bastien, P. 2015, AJ, 149, 31

Jones, T. J., Gordon, M., Shenoy, D., et al. 2016, AJ, 151, 156

Juárez, C., Girart, J. M., Zamora-Avilés, M., et al. 2017, ApJ, 844, 44

Kirk, H., Myers, P. C., Bourke, T. L., et al. 2013, ApJ, 766, 115

Koch, P. M., Tang, Y.-W., Ho, P. T. P., et al. 2014, ApJ, 797, 99

Koch, P. M., Tang, Y.-W., Ho, P. T. P., et al. 2018, ApJ, 855, 39

Kwon, J., Doi, Y., Tamura, M., et al. 2018, ApJ, 859, 4

Kwon, J., Tamura, M., Hough, J. H., et al. 2015, ApJS, 220, 17

Lada, C. J., Alves, J., & Lada, E. A. 1999, ApJ, 512, 250

Lazarian, A., Goodman, A. A., & Myers, P. C. 1997, ApJ, 490, 273

Lazarian, A., & Hoang, T. 2007a, MNRAS, 378, 910

Lazarian, A., & Hoang, T. 2007b, ApJL, 669, L77

Li, H.-b., Fang, M., Henning, T., & Kainulainen, J. 2013, MNRAS, 436, 3707

Li, H.-B., Goodman, A., Sridharan, T. K., et al. 2014, in Protostars and Planets VI, ed. H. Beuther et al. (Tucson, AZ: Univ. Arizona Press), 101

Li, H.-B., Yuen, K. H., Otto, F., et al. 2015, Natur, 520, 518

Li, Z.-Y., Wang, P., Abel, T., & Nakamura, F. 2010, ApJL, 720, L26

Liu, H. B., Galván-Madrid, R., Jiménez-Serra, I., et al. 2015, ApJ, 804, 37

Liu, H. B., Jiménez-Serra, I., Ho, P. T. P., et al. 2012, ApJ, 756, 10

Liu, T., Li, P. S., Juvela, M., et al. 2018, ApJ, 859, 151

Mairs, S., Johnstone, D., Kirk, H., et al. 2015, MNRAS, 454, 2557

Mestel, L. 1985, Protostars and Planets II (Tucson, AZ: Univ. Arizona Press), 320

Miettinen, O. 2012, A&A, 540, A104

Montier, L., Plaszczynski, S., Levrier, F., et al. 2015a, A&A, 574, A135

Montier, L., Plaszczynski, S., Levrier, F., et al. 2015b, A&A, 574, A136

Mouschovias, T. C., & Paleologou, E. V. 1986, ApJ, 308, 781

Myers, P. C. 2009, ApJ, 700, 1609

Nagai, T., Inutsuka, S., & Miyama, S. 1998, ApJ, 506, 306

Nakamura, F., & Li, Z. 2008, ApJ, 687, 354

Nakano, T., & Nakamura, T. 1978, PASJ, 30, 671

Ostriker, E. C., Stone, J. M., & Gammie, C. F. 2001, ApJ, 546, 980

Padoan, P., Goodman, A., Draine, B. T., et al. 2001, ApJ, 559, 1005

Pattle, K., Ward-Thompson, D., Berry, D., et al. 2017, ApJ, 846, 122

Pattle, K., Ward-Thompson, D., Hasegawa, T., et al. 2018, ApJL, 860, L6

Pillai, T., Kauffmann, J., Wyrowski, F., et al. 2011, A&A, 530, A118

Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2016, A&A, 586, A138

Pon, A., Johnstone, D., & Heitsch, F. 2011, ApJ, 740, 88

Pon, A., Toalá, J. A., Johnstone, D., et al. 2012, ApJ, 756, 145

Price, D. J., & Bate, M. R. 2008, MNRAS, 385, 1820

Quinn, J. L. 2012, A&A, 538, A65

Rao, R., Girart, J. M., Marrone, D. P., Lai, S.-P., & Schnee, S. 2009, ApJ, 707, 921

Rice, S. O. 1945, Bell Systems Tech. J., 24, 156

Robitaille, T., & Bressert, E. 2012, APLpy: Astronomical Plotting Library in Python, Astrophysics Source Code Library, ascl:1208.017

Salvatier, J., Wiecki, T. V., & Fonnesbeck, C. 2016, PeerJ Comput. Sci., 2, e55

Sanhueza, P., Jackson, J. M., Zhang, Q., et al. 2017, ApJ, 841, 97

Santos, F. P., Busquet, G., Franco, G. A. P., Girart, J. M., & Zhang, Q. 2016, ApJ, 832, 186

Seifried, D., & Walch, S. 2015, MNRAS, 452, 2410

Simmons, J. F. L., & Stewart, B. G. 1985, A&A, 142, 100

Soam, A., Pattle, K., Ward-Thompson, D., et al. 2018, ApJ, 861, 65

Sugitani, K., Nakamura, F., Watanabe, M., et al. 2011, ApJ, 734, 63

Tamura, M., & Kwon, J. 2015, in Polarimetry of Stars and Planetary Systems, ed. L. Kolokolova, J. Hough, & A. Levasseur-Regourd (Cambridge: Cambridge Univ. Press), 162

Tamura, M., Nagata, T., Sato, S., & Tanaka, M. 1987, MNRAS, 224, 413

Tang, Y.-W., Ho, P. T. P., Koch, P. M., et al. 2009, ApJ, 700, 251

Vaillancourt, J. E. 2006, PASP, 118, 1340

van der Walt, S., Colbert, S. C., & Varoquaux, G. 2011, CSE, 13, 22

Van Loo, S., Keto, E., & Zhang, Q. 2014, ApJ, 789, 37

Wang, J.-W., Lai, S.-P., Eswaraiah, C., et al. 2017, ApJ, 849, 157

Wardle, J. F. C., & Kronberg, P. P. 1974, ApJ, 194, 249

Ward-Thompson, D., Pattle, K., Bastien, P., et al. 2017, ApJ, 842, 66

Weingartner, J. C., & Draine, B. T. 2001, ApJ, 548, 296

Whittet, D. C. B., Hough, J. H., Lazarian, A., & Hoang, T. 2008, ApJ, 674, 304

Yuan, J., Li, J.-Z., Wu, Y., et al. 2018, ApJ, 852, 12

Zhang, Q., Qiu, K., Girart, J. M., et al. 2014, ApJ, 792, 116