# Filamentary Network and Magnetic Field Structures Revealed with BISTRO in the High-mass Star-forming Region NGC 2264: Global Properties and Local Magnetogravitational Configurations

Jia-Wei Wang, Patrick M. Koch, Seamus D. Clarke, Gary Fuller, Nicolas Peretto, Ya-Wen Tang, Hsi-Wei Yen, Shih-Ping Lai, Nagayoshi Ohashi, Doris Arzoumanian, et al.

► To cite this version:

Jia-Wei Wang, Patrick M. Koch, Seamus D. Clarke, Gary Fuller, Nicolas Peretto, et al. Filamentary Network and Magnetic Field Structures Revealed with BISTRO in the High-mass Star-forming Region NGC 2264: Global Properties and Local Magnetogravitational Configurations. The Astrophysical Journal, 2024, 962, 10.3847/1538-4357/ad165b. insu-04836900

HAL Id: insu-04836900

https://insu.hal.science/insu-04836900v1

Submitted on 13 Dec 2024

**HAL** is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire **HAL**, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

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

Distributed under a Creative Commons Attribution 4.0 International License
Filamentary Network and Magnetic Field Structures Revealed with BISTRO in the High-
mass Star-forming Region NGC 2264: Global Properties and Local Magnetogravitational
Configurations

Jia-Wei Wang<sup>1,97</sup>, Patrick M. Koch<sup>1</sup>, Seamus D. Clarke<sup>1</sup>, Gary Fuller<sup>2,3</sup>, Nicolas Peretto<sup>4</sup>, Ya-Wen Tang<sup>1</sup>,
Hsi-Wei Yen<sup>1</sup>, Shih-Ping Lai<sup>1,5</sup>, Nagayoshi Ohashi<sup>1</sup>, Doris Arzoumanian<sup>6</sup>, Doug Johnstone<sup>7,8</sup>, Ray Furuya<sup>9</sup>,
Shu-ichiro Inutsuka<sup>10</sup>, Chang Won Lee<sup>11,12</sup>, Derek Ward-Thompson<sup>13</sup>, Valentin J. M. Le Gouellec<sup>14</sup>, Hong-Li Liu<sup>15</sup>,
Lapo Fanciullo<sup>16</sup>, Jihye Hwang<sup>11</sup>, Kate Pattle<sup>17</sup>, Frédérick Poidevin<sup>18,19</sup>, Mehrnoosh Tahani<sup>20</sup>, Takashi Onaka<sup>21</sup>,
Mark G. Rawlings<sup>22</sup>, Eun Jung Chung<sup>23</sup>, Junhao Liu<sup>24</sup>, A-Ran Lyo<sup>11</sup>, Felix Priestley<sup>4</sup>, Thiem Hoang<sup>11,12</sup>,
Motohide Tamura<sup>21,25,26</sup>, David Berry<sup>24</sup>, Pierre Bastien<sup>27</sup>, Tao-Chung Ching<sup>28</sup>, Simon Coudé<sup>29,30</sup>,
Woojin Kwon<sup>31,32</sup>, Mike Chen<sup>33</sup>, Chakali Eswaraiah<sup>34</sup>, Archana Soam<sup>35</sup>, Tetsuo Hasegawa<sup>25</sup>, Keping Qiu<sup>36,37</sup>,
Tyler L. Bourke<sup>2,38</sup>, Do-Young Byun<sup>11,12</sup>, Zhiwei Chen<sup>39</sup>, Huei-Ru Vivien Chen<sup>1,5</sup>, Wen Ping Chen<sup>40</sup>,
Jungyeon Cho<sup>23</sup>, Minho Choi<sup>11</sup>, Yunhee Choi<sup>11</sup>, Youngwoo Choi<sup>41</sup>, Antonio Chrysostomou<sup>38</sup>, Sophia Dai<sup>42</sup>,
James Di Francesco<sup>7,8</sup>, Pham Ngoc Diep<sup>43</sup>, Yasuo Doi<sup>44</sup>, Yan Duan<sup>42</sup>, Hao-Yuan Duan<sup>40</sup>, David Eden<sup>45</sup>,
Jason Fiege<sup>46</sup>, Laura M. Fissel<sup>33</sup>, Erica Franzmann<sup>46</sup>, Per Friberg<sup>24</sup>, Rachel Friesen<sup>47</sup>, Tim Gledhill<sup>48</sup>,
Sarah Graves<sup>24</sup>, Jane Greaves<sup>4</sup>, Matt Griffin<sup>4</sup>, Qilao Gu<sup>49</sup>, Ilseung Han<sup>11,12</sup>, Saeko Hayashi<sup>50</sup>, Martin Houde<sup>51</sup>,
Tsuyoshi Inoue<sup>52</sup>, Kazunari Iwasaki<sup>53</sup>, Il-Gyo Jeong<sup>11,54</sup>, Vera Könyves<sup>13</sup>, Ji-hyun Kang<sup>11</sup>, Miju Kang<sup>11</sup>,
Janik Karoly<sup>13</sup>, Akimasa Kataoka<sup>55</sup>, Koji Kawabata<sup>56,57,58</sup>, Zacariyya Khan<sup>17</sup>, Mi-Ryang Kim<sup>59</sup>, Kee-Tae Kim<sup>11,12</sup>,
Kyoung Hee Kim<sup>60</sup>, Shinyoung Kim<sup>11</sup>, Jongsoo Kim<sup>11,12</sup>, Hyosung Kim<sup>31</sup>, Gwanjeong Kim<sup>61</sup>,
Florian Kirschlager<sup>62</sup>, Jason Kirk<sup>13</sup>, Masato I. N. Kobayashi<sup>6</sup>, Takayoshi Kusune<sup>63</sup>, Jungmi Kwon<sup>21</sup>,
Kevin Lacaille<sup>64,65</sup>, Chi-Yan Law<sup>66,67</sup>, Sang-Sung Lee<sup>11,12</sup>, Hyeseung Lee<sup>68</sup>, Jeong-Eun Lee<sup>69</sup>, Chin-Fei Lee<sup>1</sup>,
Dalei Li<sup>70</sup>, Hua-bai Li<sup>66</sup>, Guangxing Li<sup>71</sup>, Di Li<sup>72,73</sup>, Sheng-Jun Lin<sup>1</sup>, Tie Liu<sup>74</sup>, Sheng-Yuan Liu<sup>1</sup>, Xing Lu<sup>49</sup>,
Steve Mairs<sup>24</sup>, Masafumi Matsumura<sup>75</sup>, Brenda Matthews<sup>7,8</sup>, Gerald Moriarty-Schieven<sup>7</sup>, Tetsuya Nagata<sup>76</sup>,
Fumitaka Nakamura<sup>55,77</sup>, Hiroyuki Nakanishi<sup>78</sup>, Nguyen Bich Ngoc<sup>43,79</sup>, Geumsook Park<sup>11</sup>, Harriet Parsons<sup>24</sup>,
Tae-Soo Pyo<sup>50,77</sup>, Lei Qian<sup>72</sup>, Ramprasad Rao<sup>1</sup>, Jonathan Rawlings<sup>17</sup>, Brendan Retter<sup>4</sup>, John Richer<sup>80,81</sup>,
Andrew Rigby<sup>4</sup>, Sarah Sadavoy<sup>33</sup>, Hiro Saito<sup>82</sup>, Giorgio Savini<sup>83</sup>, Masumichi Seta<sup>84</sup>, Ekta Sharma<sup>72</sup>,
Yoshito Shimajiri<sup>85</sup>, Hiroko Shinnaga<sup>78</sup>, Xindi Tang<sup>86</sup>, Hoang Duc Thuong<sup>87</sup>, Kohji Tomisaka<sup>55</sup>, Le Ngoc Tram<sup>88</sup>,
Yusuke Tsukamoto<sup>78</sup>, Serena Viti<sup>89</sup>, Hongchi Wang<sup>39</sup>, Anthony Whitworth<sup>4</sup>, Jintai Wu<sup>36</sup>, Jinjin Xie<sup>42</sup>,
Meng-Zhe Yang<sup>5</sup>, Hyunju Yoo<sup>23</sup>, Jinghua Yuan<sup>42</sup>, Hyeong-Sik Yun<sup>11</sup>, Tetsuya Zenko<sup>76</sup>, Chuan-Peng Zhang<sup>42,72</sup>,
Yapeng Zhang<sup>90</sup>, Guoyin Zhang<sup>72</sup>, Jianjun Zhou<sup>70</sup>, Lei Zhu<sup>72</sup>, Ilse de Looze<sup>62</sup>, Philippe André<sup>91</sup>, C. Darren Dowell<sup>92</sup>,
Stewart Eyres<sup>93</sup>, Sam Falle<sup>94</sup>, Jean-François Robitaille<sup>95</sup>, and Sven van Loo<sup>96</sup>

<footnote n="1">¹ Academia Sinica Institute of Astronomy and Astrophysics, No. 1, Sec. 4, Roosevelt Road, Taipei 10617, Taiwan; jwwang@asiaa.sinica.edu.tw</footnote>

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

<footnote n="3">³ Physikalisches Institut, University of Cologne, Zülpicher Str. 77, D-50937 Köln, Germany</footnote>

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

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

<footnote n="6">⁶ Division of Science, National Astronomical Observatory of Japan, 2-21-1 Osawa, Mitaka, Tokyo 181-8588, Japan</footnote>

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

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

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

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

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

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

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

<footnote n="14">¹⁴ NASA Ames Research Center, Space Science and Astrobiology Division, MS 245-6, Moffett Field, CA 94035, USA</footnote>

<footnote n="15">¹⁵ School of Physics and Astronomy, Yunnan University, Kunming, 650091, People's Republic of China</footnote>

<footnote n="16">¹⁶ National Chung Hsing University, 145 Xingda Road, South District, Taichung City 402, Taiwan</footnote>

<footnote n="17">¹⁷ Department of Physics and Astronomy, University College London, WC1E 6BT London, UK</footnote>

<footnote n="18">¹⁸ Instituto de Astrofísica de Canarias, E-38200 La Laguna, Tenerife, Canary Islands, Spain</footnote>

<footnote n="19">¹⁹ Departamento de Astrofísica, Universidad de La Laguna (ULL), E-38206 La Laguna, Tenerife, Spain</footnote>

<footnote n="20">²⁰ Banting and KIPAC Fellow: Kavli Institute for Particle Astrophysics and Cosmology (KIPAC), Stanford University, Stanford, CA, USA</footnote>

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

<footnote n="22">²² Gemini Observatory/NSF's NOIRLab, 670 North A'ohōkū Place, Hilo, HI 96720, USA</footnote>

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

<footnote n="24">²⁴ East Asian Observatory, 660 North A'ohōkū Place, University Park, Hilo, HI 96720, USA</footnote>

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

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

<footnote n="27">²⁷ Centre de recherche en astrophysique du Québec & département de physique, Université de Montréal, 1375, Avenue Thérèse-Lavoie-Roux, Montréal, QC H2V 0B3, Canada</footnote>

<footnote n="28">²⁸ National Radio Astronomy Observatory, 1003 Lopezville Road, Socorro, NM 87801, USA</footnote>

<footnote n="29">²⁹ Department of Earth, Environment, and Physics, Worcester State University, Worcester, MA 01602, USA</footnote>

<footnote n="30">³⁰ Center for Astrophysics | Harvard & Smithsonian, 60 Garden Street, Cambridge, MA 02138, USA</footnote>

<footnote n="31">³¹ Department of Earth Science Education, Seoul National University, 1 Gwanak-ro, Gwanak-gu, Seoul 08826, Republic of Korea</footnote>
32 SNU Astronomy Research Center, Seoul National University, 1 Gwanak-ro, Gwanak-gu, Seoul 08826, Republic of Korea

33 Department for Physics, Engineering Physics and Astrophysics, Queen's University, Kingston, ON K7L 3N6, Canada

34 Indian Institute of Science Education and Research (IISER) Tirupati, Rami Reddy Nagar, Karakambadi Road, Mangalam (P.O.), Tirupati 517 507, India

35 Indian Institute of Astrophysics, II Block, Koramangala, Bengaluru 560034, India

36 School of Astronomy and Space Science, Nanjing University, 163 Xianlan Avenue, Nanjing 210023, People's Republic of China

37 Key Laboratory of Modern Astronomy and Astrophysics (Nanjing University), Ministry of Education, Nanjing 210023, People's Republic of China

38 SKA Observatory, Jodrell Bank, Lower Withington, Macclesfield SK11 9FT, UK

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

40 Institute of Astronomy, National Central University, Zhongli 32001, Taiwan

41 Department of Physics and Astronomy, Seoul National University, 1 Gwanak-ro, Gwanak-gu, Seoul 08826, Republic of Korea

42 National Astronomical Observatories, Chinese Academy of Sciences, A20 Datun Road, Chaoyang District, Beijing 100012, People's Republic of China

43 Vietnam National Space Center, Vietnam Academy of Science and Technology, 18 Hoang Quoc Viet, Hanoi, Vietnam

44 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

45 Armagh Observatory and Planetarium, College Hill, Armagh BT61 9DG, UK

46 Department of Physics and Astronomy, The University of Manitoba, Winnipeg, MB R3T2N2, Canada

47 National Radio Astronomy Observatory, 520 Edgemont Road, Charlottesville, VA 22903, USA

48 School of Physics, Astronomy & Mathematics, University of Hertfordshire, College Lane, Hatfield, Hertfordshire AL10 9AB, UK

49 Shanghai Astronomical Observatory, Chinese Academy of Sciences, 80 Nandan Road, Shanghai 200030, People's Republic of China

50 Subaru Telescope, National Astronomical Observatory of Japan, 650 North A'ohōku Place, Hilo, HI 96720, USA

51 Department of Physics and Astronomy, The University of Western Ontario, 1151 Richmond Street, London, ON N6A 3K7, Canada

52 Department of Physics, Konan University, Okamoto 8-9-1, Higashinada-ku, Kobe 658-8501, Japan

53 Center for Computational Astrophysics, National Astronomical Observatory of Japan, Mitaka, Tokyo 181-8588, Japan

54 Department of Astronomy and Atmospheric Sciences, Kyungpook National University, Daegu 41566, Republic of Korea

55 Division of Theoretical Astronomy, National Astronomical Observatory of Japan, Mitaka, Tokyo 181-8588, Japan

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

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

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

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

60 Ulsan National Institute of Science and Technology (UNIST), UNIST-gil 50, Eonyang-eup, Ulju-gun, Ulsan 44919, Republic of Korea

61 Nobeyama Radio Observatory, National Astronomical Observatory of Japan, National Institutes of Natural Sciences, Nobeyama, Minamimaki, Minamisaku, Nagano 384-1305, Japan

62 Sterrenkundig Observatorium, Ghent University, Krijgslaan 281-S9, B-9000 Ghent, Belgium

63 Astronomical Institute, Graduate School of Science, Tohoku University, Aoba-ku, Sendai, Miyagi 980-8578, Japan

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

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

66 Department of Physics, The Chinese University of Hong Kong, Shatin, N.T., Hong Kong

67 Department of Space, Earth & Environment, Chalmers University of Technology, SE-412 96 Gothenburg, Sweden

68 Ulsan National Institute of Science and Technology (UNIST), 50 UNIST-gil, Ulsan 44919, Republic of Korea

69 Astronomy Program, Department of Physics and Astronomy, Seoul National University, 1 Gwanak-ro, Gwanak-gu, Seoul 08826, Republic of Korea

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

71 Department of Astronomy, Yunnan University, Kunming, 650091, People's Republic of China

72 CAS Key Laboratory of FAST, National Astronomical Observatories, Chinese Academy of Sciences, People's Republic of China

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

74 Key Laboratory for Research in Galaxies and Cosmology, Shanghai Astronomical Observatory, Chinese Academy of Sciences, 80 Nandan Road, Shanghai 200030, People's Republic of China

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

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

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

78 Department of Physics and Astronomy, Graduate School of Science and Engineering, Kagoshima University, 1-21-35 Korimoto, Kagoshima, Kagoshima 890-0065, Japan

79 Graduate University of Science and Technology, Vietnam Academy of Science and Technology, 18 Hoang Quoc Viet, Cau Giay, Hanoi, Vietnam

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

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

82 Faculty of Pure and Applied Sciences, University of Tsukuba, 1-1-1 Tennodai, Tsukuba, Ibaraki 305-8577, Japan

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

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

85 Kyushu Kyoritsu University, 1-8, Jiyugaoka, Yahatanishi-ku, Kitakyushu-shi, Fukuoka 807-8585, Japan

86 Xinjiang Astronomical Observatory, Chinese Academy of Sciences, 830011 Urumqi, People's Republic of China

87 Kavli Institute for the Physics and Mathematics of the Universe (WPI), UTIAS, The University of Tokyo, Kashiwa, Chiba 277-8583, Japan

88 University of Science and Technology of Hanoi, Vietnam Academy of Science and Technology, 18 Hoang Quoc Viet, Hanoi, Vietnam

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

90 Department of Astronomy, Beijing Normal University, Beijing 100875, People's Republic of China

91 Laboratoire AIM CEA/DSM-CNRS-Université Paris Diderot, IRFU/Service d'Astrophysique, CEA Saclay, F-91191 Gif-sur-Yvette, France

92 Jet Propulsion Laboratory, M/S 169-506, 4800 Oak Grove Drive, Pasadena, CA 91109, USA

93 University of South Wales, Pontypridd, CF37 1DL, UK

94 Department of Applied Mathematics, University of Leeds, Woodhouse Lane, Leeds LS2 9JT, UK

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

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

Received 2023 August 13; revised 2023 December 14; accepted 2023 December 15; published 2024 February 14

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

<footnote n="97">⁹⁷ Corresponding author.</footnote>
Abstract

We report 850 µm continuum polarization observations toward the filamentary high-mass star-forming region NGC 2264, taken as part of the B-fields In STar forming Regions Observations large program on the James Clerk Maxwell Telescope. These data reveal a well-structured nonuniform magnetic field in the NGC 2264C and 2264D regions with a prevailing orientation around 30° from north to east. Field strength estimates and a virial analysis of the major clumps indicate that NGC 2264C is globally dominated by gravity, while in 2264D, magnetic, gravitational, and kinetic energies are roughly balanced. We present an analysis scheme that utilizes the locally resolved magnetic field structures, together with the locally measured gravitational vector field and the extracted filamentary network. From this, we infer statistical trends showing that this network consists of two main groups of filaments oriented approximately perpendicular to one another. Additionally, gravity shows one dominating converging direction that is roughly perpendicular to one of the filament orientations, which is suggestive of mass accretion along this direction. Beyond these statistical trends, we identify two types of filaments. The type I filament is perpendicular to the magnetic field with local gravity transitioning from parallel to perpendicular to the magnetic field from the outside to the filament ridge. The type II filament is parallel to the magnetic field and local gravity. We interpret these two types of filaments as originating from the competition between radial collapsing, driven by filament self-gravity, and longitudinal collapsing, driven by the region's global gravity.

*Unified Astronomy Thesaurus concepts: Collapsing clouds (267); Interstellar filaments (842); Molecular clouds (1072); Interstellar magnetic fields (845); Interstellar medium (847); Polarimetry (1278); Interstellar dynamics (839); Star forming regions (1565); Submillimeter astronomy (1647); Young massive clusters (2049)*

1. Introduction

Observations over the last decade have shown that interstellar filaments are ubiquitous within molecular clouds and that they can be major sites of star formation (Schneider & Elmegreen 1979; Myers 2009; André et al. 2010, 2014; Arzoumanian et al. 2011; Chung et al. 2019; Kumar et al. 2020). Understanding how these filaments form is therefore crucial to determining the initial stage of star formation. Hub-filament structures (HFSs), consisting of a massive hub connecting with several converging filaments, are special filamentary configurations that have recently drawn much attention, because they are now commonly identified in nearby massive star-forming regions (e.g., Myers 2009; Schneider et al. 2012; Motte et al. 2018; Liu et al. 2023). Kumar et al. (2020) show that all massive clumps within 2 kpc detected in the Herschel Hi-GAL survey are located in the hubs of HFSs. Recent Atacama Large Millimeter/submillimeter Array observations suggest that self-similar hierarchical HFSs, where small-scale HFSs are embedded in the hub of large-scale HFSs, are possibly common in massive protoclusters (Zhou et al. 2022). All these observational results point to HFSs being an essential step during the formation of massive stars and clusters.

Theoretically, the formation mechanism of filaments and their assembly into a hub remains a topic of debate, with several proposed processes including layer fragmentation threaded by magnetic fields (Myers 2009; Van Loo et al. 2014), magnetic field-channeled gravitational collapse (e.g., Nakamura & Li 2008), shock/turbulence compression (e.g., Padoan et al. 2001; Federrath 2016), hierarchical collapse in molecular clouds (e.g., Gómez & Vázquez-Semadeni 2014; Vázquez-Semadeni et al. 2017; Gómez et al. 2018), and filament-filament collisions (e.g., Dobashi et al. 2014; Nakamura et al. 2014; Frau et al. 2015; Kumar et al. 2020). The magnetic field can be one of the key factors determining the dominating mechanism in these processes. A strong large-scale magnetic field might be important in guiding large-scale accretion flows (Palmeirim et al. 2013; Li et al. 2014), gravitational collapse driven by magnetohydrodynamic (MHD)

turbulence (Nakamura & Li 2008), or layer fragmentation (Myers 2009; Van Loo et al. 2014).

Observationally, magnetic fields are commonly found to correlate with the orientation of interstellar filaments (Heyer et al. 2008; Li & Houde 2008; Busquet et al. 2013; Palmeirim et al. 2013; Planck Collaboration et al. 2016; Hwang et al. 2022), suggesting that they play an important role in the filament formation. The correlation between filaments and magnetic fields likely varies with column densities along or across filaments (Planck Collaboration et al. 2016; Pillai et al. 2020; Arzoumanian et al. 2021; Kwon et al. 2022), implying that the role of the magnetic field possibly evolves locally. These variations could be more significant in HFSs because the local densities within filaments can dramatically change due to the strong gravitational field of a massive hub (Wang et al. 2020a; Arzoumanian et al. 2021; Chung et al. 2022). However, studies addressing the local variations of the magnetic field on subparsec scales are still rare, and hence the exact role of the magnetic field remains unclear.

Most studies on interstellar magnetic fields focus on global, averaged properties of magnetic fields, such as mass-to-magnetic flux ratio, virial parameters, and Alfvénic Mach number. This is due to the difficulties and challenges in both probing the detailed magnetic field morphologies and extracting quantitative information. These analyses on averaged quantities are useful to describe relatively simple systems but become insufficient for complex systems, such as hierarchical filamentary networks within HFSs. Recently, new analytical methods have been proposed to extract local properties of magnetic fields. These include the polarization-intensity gradient method (Koch et al. 2012a, 2012b, 2013), the histogram of relative orientation (HRO) analysis (Soler et al. 2013; Planck Collaboration et al. 2016; Soler 2019; Kwon et al. 2022), spatial distributions of magnetic field strengths (Wang et al. 2020b; Hwang et al. 2021, 2022), and relative orientations among gravity, magnetic fields, and density structures (Koch et al. 2018, 2022; Wang et al. 2020a, 2022; Liu et al. 2023). In this paper, we aim at building an analysis scheme highlighting the local physical conditions within HFSs in order to reveal the variations in the role of the magnetic field.
NGC 2264 is an active cluster-forming region populated by hundreds of young stellar objects (YSOs; Sung et al. 2009; Rapson et al. 2014) embedded in the Mon OB1 molecular cloud complex. It is located at a Gaia DR2-determined distance of $715_{-42}^{+81}$ pc (Zucker et al. 2020). Buckle et al. (2012) identified six parsec-scale filaments converging toward NGC 2264 from $^{12}$CO (3–2) and H$_2$ 1–0 S(1) wide-field images, revealing that NGC 2264 is the hub center of a 10 pc-sized HFS. Kumar et al. (2020) further identified filamentary structures of young stars joining at NGC 2264. Since the radial velocities of these stars in different filaments fall into different velocity subgroups (Tobin et al. 2015), they proposed that NGC 2264 originated from a collision of stellar filaments. Within the hub center, high-resolution millimeter and sub-millimeter continuum observations reveal an HFS-like morphology, where filamentary structures with lengths around 0.5–2 pc converge toward dense cores (Peretto et al. 2006; Burkhart et al. 2015). Hence, this likely points at self-similar hierarchical HFSs. Based on this, NGC 2264 is an ideal source to investigate the formation of hierarchical HFSs.

NGC 2264 primarily consists of two massive, supervirial cluster-forming clumps (Peretto et al. 2006) embedded in the hub of the 10 pc-sized HFS. The southern 1650 M$_{\odot}$ clump (Peretto et al. 2006) hosts the bright 9.5 M$_{\odot}$ B2 zero-age main-sequence star IRAS 06384+0932 (IRS1), also known as Allen's source (Allen 1972), which is associated with the molecular outflow 2264C (Margulis et al. 1988). The northern 1310 M$_{\odot}$ clump (Peretto et al. 2006) contains a class I YSO, IRAS 06382+0939 (IRS2), which is associated with the molecular outflow NGC 2264D (Margulis et al. 1988). These two clumps are commonly named after their two dominating sources (IRS1 and IRS2) or their associated outflows (NGC 2264C and NGC 2264D). For clarity, we use NGC 2264C and NGC 2264D (or 2264C and 2264D) when referring to the two massive clumps and IRS1 and IRS2 when referring to the two sources.

In this paper, we report 850 µm continuum polarization observations toward NGC 2264 using the James Clerk Maxwell Telescope (JCMT) POL-2 polarimeter as part of the B-fields In STar forming Regions Observations (BISTRO) large program (Ward-Thompson et al. 2017). These observations, with a resolution of 14", allow us to probe the detailed magnetic field morphology within 2264C and 2264D. The goal of the paper is to study how a hierarchical HFS can be formed by investigating how gravity and magnetic fields interact and link to hubs and filaments. In Section 2, we present the observations and data reduction. Section 3 reports the observed intensity and magnetic field morphology. In Section 4, we present how to extract spatial parameters from the observed data. Statistical analyses are performed to identify possible trends and correlations of these parameters using both maps probing the local spatial properties and histograms highlighting the statistical properties. The roles of gravity and magnetic field in the formation of NGC 2264 are discussed in Section 5. Our conclusions are summarized in Section 6.

## 2. Observations

### 2.1. JCMT POL-2 Observations

We carried out polarization continuum observations toward NGC 2264 with SCUBA-2 and POL-2 mounted on the JCMT (project code M17BL011) between 2017 November and 2019

February. The two-pointing mosaic observations targeted 2264C (south) and 2264D (north) with the reference positions (R.A., decl.) = (6$^{h}$41$^{m}$11$^{s}$97, 9$^{\circ}$29'44"5) and (6$^{h}$41$^{m}$01$^{s}$5, 9$^{\circ}$35'39"5). We performed 21 sets of 40 minute observations for each pointing under a weather condition with $\tau_{225 GHz}$ ranging from 0.02 to 0.06. The POL-2 DAISY scan mode (Friberg et al. 2016) was adopted, producing a fully sampled circular region with a diameter of 11' for each pointing and a resolution of 14"1 (corresponding to 0.05 pc at a source distance of 715 pc). Polarization was simultaneously observed in both the 450 and 850 µm continuum bands. This paper focuses on the 850 µm data.

The POL-2 850 µm polarization data were reduced with pol2map$^{98}$ in the SMURF package$^{99}$ (Berry et al. 2005; Chapin et al. 2013). The skyloop mode was invoked, which is a script that runs mapmaking on the full set of observations in order to find a solution that minimizes residuals across the full set of maps. The MAPVARS mode was activated to estimate the uncertainties from the standard deviation among the individual observations, accounting for the instrumental and mapmaking uncertainties. The details of the data reduction steps and procedure are described in previous BISTRO papers (e.g., Wang et al. 2019; Pattle et al. 2021). The POL-2 data reduction was done with a 4" pixel size because larger pixel sizes might increase the uncertainty during the mapmaking process.

The output Stokes I, Q, and U images were calibrated in units of Jy arcsec$^{-2}$ using an aperture flux conversion factor (FCF) of 2.79 Jy pW$^{-1}$ arcsec$^{-2}$ (including a factor of 1.35 for POL-2, equivalent to 630 Jy pW$^{-1}$ beam$^{-1}$) for extended sources (Mairs et al. 2021) and binned to a pixel size of 12" to improve the sensitivity. The typical rms noise of the final Stokes I, Q, and U maps is ~0.01 mJy arcsec$^{-2}$ at the map center and gradually increases to ~0.04 mJy arcsec$^{-2}$ near the edge of the mapped region. The Stokes I image has a higher rms noise of up to ~0.07 mJy arcsec$^{-2}$ for pixels with large intensities ($I > 6$ mJy arcsec$^{-2}$), where the rms noise is dominated by the mapmaking uncertainties. The calculated polarization fraction P was debiased using the asymptotic estimator (Wardle & Kronberg 1974) as

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

with a polarization uncertainty $\sigma_P$ calculated as

$$\sigma_P = \sqrt{\frac{Q^2\sigma_Q^2 + U^2\sigma_U^2}{(Q^2 + U^2)I^2} + \frac{\sigma_I^2(Q^2 + U^2)}{I^4}}, \quad (2)$$

where $\sigma_I$, $\sigma_Q$, and $\sigma_U$ are the uncertainties in the I, Q, and U Stokes parameters. The polarization position angle (PA) is

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

and the corresponding uncertainty is estimated as

$$\delta 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 in this paper are inferred to be PA + 90°.

<footnote>$^{98}$ http://starlink.eao.hawaii.edu/docs/sc22.pdf</footnote>

<footnote>$^{99}$ Version 2021 December 2.</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Magnetic field map (yellow segments) toward NGC 2264 based on POL-2 850 µm continuum polarization observations, overlaid on Stokes I 850 µm continuum in color. Segments are rotated by 90° with respect to the originally detected polarization orientations (Figure 17). White contours mark 10σ (0.1 mJy arcsec⁻²) in Stokes I. Green segments are the larger-scale magnetic field from Planck at 353 GHz (850 µm). A scale bar and the JCMT POL-2 beam are shown in the lower right corner.</figcaption></figure>

## 2.2. JCMT HARP Observations

¹³CO (3–2) and C¹⁸O (3–2) molecular line observations were carried out simultaneously using the Heterodyne Array Receiver Program (HARP) instrument (Buckle et al. 2009) on the JCMT toward 2264D in 2022 February (project code: M22AP045; PI: Jia-Wei Wang). The raster scan produced an 8.5 × 5.5 map centered on (R.A., decl.) = (6h41m00s, 9°34′43″2). The observations with a half-power beamwidth of 14.1 were carried out with a spectral resolution of 50 kHz, yielding 0.066 km s⁻¹, and a bandwidth of 625 MHz. The data reduction was performed using the Starlink package using the default ORAC-DR pipeline (Cavanagh 2008). To enhance the sensitivity, we binned 2 × 2 pixels of the reduced data, originally Nyquist sampled with a pixel size of 7", resulting in a pixel size of 14". Additionally, we smoothed the spectral resolution with a Gaussian kernel to 0.11 km s⁻¹.

In order to obtain a complete ¹³CO (3–2) and C¹⁸O (3–2) map over the entire NGC 2664 system, we additionally included archival JCMT HARP data toward 2264C (project code: M08BU18) for these two lines, with a spectral resolution of 0.11 km s⁻¹. These data cover an 8.5 × 5.5 area centered on (R.A., decl.) = (6h41m10s, 9°29′40″2). The archival data were rebinned to the same pixel grid as our above 2264D data using the Starlink package *wcsalign*. The two data sets were then combined to obtain a full map over NGC 2264. The final rms noise of these data is ~0.5 and 0.2 K for NGC 2264C and D, respectively. This paper focuses on the velocity dispersion from these two lines in order to estimate a magnetic field strength. A more detailed analysis of these velocity data sets will be reported in a later paper.

# 3. Results

## 3.1. POL-2 850 µm Continuum Map

Figure 1 shows the observed 850 µm continuum map in a logarithmic color scale, selected to display all reliable continuum detection (>10σ). The structures within NGC 2264 span a wide dynamical range, from 0.1 mJy arcsec⁻² (signal-to-noise ratio, S/N = 10) in the faint outskirts to 16 mJy arcsec⁻² (S/N = 1600) in the brightest clumps. To better visualize the filamentary and clumpy structures, Figure 2 shows the 850 µm continuum map with different color scales, highlighting the bright compact clumps in the top panels and the filamentary extended structures in the bottom panels. In 2264C, at least four dense clumps are found lined up along an east–west filamentary structure. A number of fainter streamer-like structures are revealed extending from this east–west filament to the northern and southern outskirts. In contrast to this, the bright clumps in 2264D are more scattered, and they are likely connected with a fainter filamentary network. A number of fainter streamer-like structures can also be seen in 2264D extending from the filamentary network to the fainter outskirts. These filament/streamer–clump configurations are self-similar to the parsec-scale hub–filament configuration around NGC 2264 (Kumar et al. 2020), appearing as typical hierarchical structures. We will further extract these filament-like structures in Section 4.1.1.

## 3.2. Magnetic Field from 850 µm Continuum Polarization

In order to ensure the robustness of polarization measurements, we select polarization detections based on the criteria $I/\sigma_I > 20$ and $P/\sigma_P > 3$. We further exclude detections with large uncertainties in polarization fraction ($\sigma_P > 5\%$). As a result, a total of 622 polarization measurements are selected across NGC 2264. For these selected data, the maximum and median uncertainties in the polarization PA are 9°:54 and 4°:5. In order to examine whether or not the observed polarization originates from magnetically aligned dust grains, Appendix A shows the polarization fraction map and discusses how the polarization fraction correlates with the total intensity. Based on a Bayesian analysis, we find that, indeed, dust grains are likely magnetically aligned in NGC 2264.

Figure 1 shows the observed magnetic field. Its morphology appears organized over the entire region. In the northern 2264D, the overall magnetic field is oriented with a PA around 30° (measured counterclockwise from north to east). In some of the dense regions, the magnetic field locally shows curved patterns, which seem to be associated with the elongated dense structures. In the southern 2264C, an overall magnetic field orientation around a similar PA ~30° is observed. Additionally, an hourglass-like pattern can be seen along the brightest filamentary structure (along a southeast–northwest direction, roughly perpendicular to the prevailing 30° field orientation). As a result, the prevailing magnetic field orientation has a wider range than in 2264D. At the northeastern end of 2264C, the magnetic field is flipped by almost 90° to an orientation around ~120°, perpendicular to the overall prevailing orientation, becoming parallel to the central filamentary structure. This flipped magnetic field morphology is possibly influenced by large-scale accreting filaments, or it might be part of the hourglass-like morphology but twisted by projection effects.

The larger-scale magnetic field, traced by the Planck 353 GHz data (beam size of 5'), reveals an organized, nearly uniform field with a mean orientation of 133° across the entire
<figure><img src="image_4.png" /><figcaption>**Figure 2.** NGC 2264 continuum maps on different intensity scales. The cyan contours mark 0.1, 0.2, 0.5, 1.0, 2.0, 5.0, and 10 mJy arcsec⁻². (Top panels) High-intensity compact clumps are distributed along a filamentary structure in 2264C. The bright clumps in 2264D are more scattered. (Bottom panels) A number of streamer-like structures are found in 2264C extending from the massive filament with the dense clumps. In 2264D, a more complex filamentary network is revealed connecting the bright clumps, with fainter streamers extending to the low-intensity outskirts.</figcaption></figure>

NGC 2264 region (overlaid in Figure 1). The clear difference in orientations between large- and small-scale magnetic fields suggests that the subparsec-scale magnetic field traced by POL-2 has been decoupled from the larger-scale field traced by Planck.

### 3.3. Velocity Dispersion Traced by ¹³CO and C¹⁸O Lines

We extracted the velocity dispersion from our molecular line data adopting the scousePy package (Henshaw et al. 2016a, 2019).

This package is a Python implementation of the spectral line-fitting IDL code SCOUSE (Henshaw et al. 2016b), which can efficiently identify multiple velocity components in a large spectral cube. We applied scousePy on both the ¹³CO (3–2) and C¹⁸O (3–2) data pixel by pixel, where the pixel size is 14" and the channel width is 0.11 km s⁻¹. The pixel rms noise was calculated using the line-free channels, and the identified velocity components were selected using the criteria that the peak brightness temperature be greater than 5σ and the velocity dispersion less than 2 km s⁻¹.
<figure><img src="image_5.png" /><figcaption>**Figure 3.** (a) C¹⁸O (3–2) integrated intensity map and (b) velocity dispersion map. The black and white contours mark the major clumps in 2264C and 2264D, extracted using a dendrogram algorithm. The 0.1 mJy arcsec⁻² contour is shown in yellow and red, as in Figure 1. Examples of spectra are extracted at positions (a), (b), (c), and (d).</figcaption></figure>

The C¹⁸O (3–2) line has a maximum brightness temperature of 6 K and is mostly below 4 K. This is much lower than the possible gas temperature in molecular clouds (~10–20 K), and it is thus likely optically thin. We further found that the intensity ratio of ¹³CO (3–2) and C¹⁸O (3–2) is around a constant of 3.7 over the entire cloud (see details in Appendix B, Figure 22), except for the few brightest pixels. This indicates that the ¹³CO (3–2) line is also optically thin except for the brightest pixels.

In the ¹³CO data cube, we commonly detected two to three velocity components across the source (for more information, see Appendix B). In contrast to this, only one major component was detected in the C¹⁸O cube. Since both ¹³CO and C¹⁸O are optically thin, this difference possibly results from velocity components in the low-density outskirts, seen only in ¹³CO but not traced by C¹⁸O. Since we aim at estimating the magnetic field strength from combining polarization and molecular line data, we mainly adopt the C¹⁸O data for the following analyses. This is because they better match the POL-2 polarization data, in which the large-scale, low-intensity structures are filtered out. Nevertheless, we also perform the analysis using the ¹³CO data, with velocity dispersions higher by 30%–50% (Appendix F) to evaluate the possible impact of the choice of gas tracers.

Figure 3 displays the C¹⁸O integrated intensity and extracted velocity dispersion maps. The moment maps are made using the moment-masking technique included in the BTS package (Clarke et al. 2018) with the masking threshold $T_C = 8\sigma$ and $T_L = 3\sigma$. Full details of the technique may be found in the code documentation.<fnref n="100" /> The integrated intensity map shows

structures visually similar to those in the continuum map. In order to identify the major components of 2264C and 2264D from the 850 µm continuum image, we performed a dendrogram analysis using the astrodendro package (Goodman et al. 2009) with a minimum significance and threshold of 0.05 mJy arcsec⁻² (5σ). The two branches matching the major clumps in 2264C and 2264D are selected and plotted in Figure 3 with intensity levels of 2.0 and 1.3 mJy arcsec⁻², respectively. The estimated velocity dispersion seems to increase toward the intensity peaks, and thus a range of velocity dispersion values is present in each of the major clumps. The amplitude-weighted mean velocity dispersion in the 2264C and 2264D major clump is 0.94 ± 0.24 and 1.04 ± 0.23 km s⁻¹, where the uncertainties are the standard deviations of velocity dispersion within the two clumps. These values are much more turbulent than other low-mass filamentary clouds observed in the BISTRO survey (e.g., IC 5146, 0.1–0.3 km s⁻¹, Wang et al. 2019; Ophiuchus C, 0.13 km s⁻¹, Liu et al. 2019; Serpens Main, 0.2–0.3 km s⁻¹, Kwon et al. 2022) but similar to other massive regions (e.g., Orion A, 1.33 ± 0.31 km s⁻¹, Pattle et al. 2017; Mon R2, 0.45–0.73 km s⁻¹, Hwang et al. 2022).

## 4. Analysis

One fundamental question in the formation of an HFS with its filamentary network is how the density structures are shaped by gravity and magnetic field. In order to answer this question, we aim to extract spatial features (Section 4.1) from both the dust emission continuum and the polarization map and investigate how these features possibly correlate. Figure 4 shows a flowchart of our analysis scheme. In Section 4.1.1, we

<footnote n="100">¹⁰⁰ https://github.com/SeamusClarke/BTS</footnote>
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Analysis flowchart. Our analysis consists of two aspects: local measures revealing the spatially localized information on maps and statistical measures showing the overall trends from histograms. Our analysis focuses on individual parameters in Section 4.2, extends to the correlation between two parameters in Section 4.3, and further includes information on densities in Section 4.4.</figcaption></figure>

identify filaments (F) and estimate their orientations. Additionally, we estimate the projected gravitational force (G) at every pixel to probe the gravitational vector field (Section 4.1.2), and we use the magnetic field (B) orientations as traced by the polarization data. With these spatially resolved parameters (F, G, B), we examine their spatial distributions on the maps (hereafter local measures) and their overall tendencies using histograms (hereafter statistical measures). Section 4.2 focuses on the properties of these individual parameters, and Section 4.3 investigates trends and possible correlations between these parameters. Section 4.4 further addresses how these correlations evolve with the local densities. Finally, we conduct an analysis of the stability of 2264C and D in Section 4.5, focusing on their global properties.

## 4.1. Introducing Local Measures

### 4.1.1. Filament Identification

In order to identify the filamentary structures within NGC 2264, we apply the DisPerSE algorithm (Sousbie 2011) on the 4" pixel JCMT 850 µm Stokes I image. We use a persistent threshold of 0.8 mJy arcsec⁻² (8σ), which evaluates the intensity difference between two connecting critical points. The identified filaments are smoothed by 28" (2× beam size). Filaments with lengths shorter than 3× beam size are excluded to ensure a minimum length. We further examine the reliability and robustness of the extracted filaments by performing Gaussian fits on the radial intensity profile at each pixel along the filaments. This shows that at least 70% of the radial profiles within each identified filament can be well described by a Gaussian, where the major peak can be identified and the fitted filament width is above 3σ. We note that those profiles that are not well fit are usually asymmetric profiles due to overlapping clumps or filaments. The identified filaments are plotted in Figure 5. The detailed filament properties are given in Appendix C. The identified filaments have typical lengths of around 0.2–1.0 pc and widths of around 0.1 pc. We note that this is a much smaller aspect ratio than seen for the typical

interstellar filaments identified by Herschel (e.g., André et al. 2010; Arzoumanian et al. 2011).

In 2264C, a bright 0.5 pc long filament (Fil 1) with an orientation around 120° connects most of the central bright structures. Two long but fainter filaments (Fil 2 and 3) extend to the south from the two ends of the brightest filament (Fil 1). In 2264D, the longest (0.4–0.7 pc) filaments (Fil 4, 5, 6, and 7) seem nearly parallel to each other with orientations around 120°, with numerous shorter (0.2–0.3 pc) filaments extending roughly orthogonally from the long ones.

To estimate the local orientations of the identified filament ridges, we extract the orientation of the tangent in each pixel along the filaments. For the ith pixel along a filament, we fit the positions of the (i − 2)th to the (i + 2)th consecutive pixels along the filament with a straight line. With this, we find a median fitting error in the local orientation of a filament of ~2°.

### 4.1.2. Local Gravitational Vector Field

In order to evaluate whether the density structures in NGC 2264 are primarily driven by gravity, we estimate the projected gravitational vector field from the 4" pixel JCMT 850 µm continuum data, following the polarization-intensity gradient technique in Koch et al. (2012a, 2012b). This is done by calculating the vector sum (F<sub>G,i</sub>) of all gravitational forces from masses at all pixels j over a map<fnref n="101" /> (Wang et al. 2020a), acting at a pixel location i as

$$F_{G,i} = kI_i \sum_{j=1}^{n} \frac{I_j}{r_{i,j}^2} \hat{r}, \quad (5)$$

where k is a factor accounting for the gravitational constant and conversion from emission to total column density. I<sub>i</sub> and I<sub>j</sub> are the intensity at the pixel positions i and j, and n is the total number of pixels within the area of relevant gravitational influence. r<sub>i,j</sub> is the plane-of-sky projected distance between

<footnote n="101">¹⁰¹ We use the 4" pixel continuum map to better sample the gravitational field. We note that while the oversampled pixels may exhibit partial correlations, this does not impact the determined orientations of the gravitational field due to the isotropic nature of the beam pattern.</footnote>
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Filaments (green lines) overlaid on 850 µm continuum toward (a) 2264C and (b) 2264D. The cyan contours mark 0.1, 0.2, 0.5, 1.0, 2.0, 5.0, and 10 mJy arcsec⁻². The red and yellow stars label IRS1 (zero-age main-sequence) and IRS2 (class I), the dominating sources in 2264C and 2264D, respectively. We note that IRS2 is offset from any filament or intensity peak, while IRS1 is at the converging point of multiple filaments. The main filaments along the prevailing directions are labeled.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 6.** Projected gravitational vector field (red) and filaments (green lines) overlaid on 850 µm continuum toward (a) 2264C and (b) 2264D. The gravity vectors are calculated per pixel but only shown per 3 × 3 pixels for simplicity and with uniform length. The major gravitational converging centers are labeled C1–C5.</figcaption></figure>

pixels *i* and *j*, and $\hat{r}$ is its unit vector. We only focus on the directions of the local gravitational forces and not their absolute magnitudes, and hence a constant factor *k* = 1 is adopted. We apply a mask with a threshold of 0.1 mJy arcsec⁻² (10σ) on the continuum map, because the gravitational force originating from the diffuse and extended surrounding material is both small and rather symmetrical and thus assumed to be mostly

canceled out. Assuming that the intensity distribution is a fair approximation for the distribution of the total mass, and that these mass components in NGC 2264 are roughly at the same distance, the calculated $F_{G,i}$ can be used as a proxy for the projected gravitational vector field.

Figure 6 displays the local gravitational vector field in NGC 2264. In 2264C, the gravitational field is largely dominated by
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Histograms of filament, magnetic field, local gravity, and local intensity gradient orientations. The green and red histograms represent 2264C and 2264D, respectively. The blue filled histograms are for the two combined. The vertical black dashed line labels the large-scale mean magnetic field orientation of 144°, as traced by Planck across the entire NGC 2264 region (Figure 1). All four parameters show common peaks around 20°–50° for both 2264C and 2264D separately as well as for the combined data. Only the histograms of filament orientations reveal additional peaks around ~150°. The intensity gradient distribution in the rightmost panel is further discussed in Appendix D.</figcaption></figure>

the densest elongated clump. Thus, the gravitational field visually tends to be either parallel or perpendicular to it. One clear gravitational converging point (C1) can be found at the massive major filament Fil 1. The connection between the gravitational vector field and filaments Fil 2 and 3 extending from the major filament appears more irregular. In 2264D, the gravitational field reveals four major converging points (C2–C5) within four east-west filaments (Fil 4, 5, 6, and 7). This suggests that 2264D is structured into four major filaments parallel to each other.

## 4.2. One-parameter Analysis: Globally Prevailing Orientations

We plot the histograms of the orientations of the filament (F), B-field (B), and gravity (G) in Figure 7, both individually for 2264C and 2264D and for them combined. These histograms reveal the global tendencies of these parameters. None of the histograms is random. They show either a single-peaked or a bimodal distribution.

A common prevailing orientation, between about 20° and 50°, is found for all parameters. This possibly implies that there is one underlying mechanism that is regulating and driving these distributions. Different from the single-peaked distributions found for the B-field and gravity, the filament orientations appear bimodal, peaking around 20°–50° and additionally also around 150°. The two peaks are separated by almost 90° and hence indicate two prevailing orientations that are nearly perpendicular. The peak observed around 150° is specifically linked to the brightest filaments that are nearly parallel to each other (Fil 1, 4, 5, 6, and 7) as depicted in Figure 5. Therefore, we will henceforth refer to these filaments as the “main filaments.” We note the close similarity in the B-field and gravity distributions, with the B-field distribution being a little wider. Its wider peak is visible from the 0°/180° wrapping, i.e., an increase around 180° that is connecting to and continuing at 0°.

In summary, the distributions of the orientations of the filament, B-field, and gravity in 2264C and 2264D show similarities and a prominent approximate 90° spacing in the filament orientations. This might suggest an evolution and connection between these parameters driven by the same process, as we further elaborate in Sections 5.1 and 5.2.

<figure><img src="image_10.png" /><figcaption>**Figure 8.** Histograms of relative orientations between the spatial parameters F, G, and B for all associated pairs across the entire NGC 2264 region. The data are separated into low- and high-intensity regions with an intensity threshold of 2.0 mJy arcsec⁻². The *p*-value from Rayleigh’s test indicates the similarity between an observed and a uniform distribution, and *p* < 0.05 favors a nonuniform distribution. PRS $Z_x$ indicates whether the relative orientations tend to be parallel ($Z_x > 0$) or perpendicular ($Z_x < 0$).</figcaption></figure>

## 4.3. Two-parameter Analysis: Relative Orientations among Local Parameters

In this section, we examine the relative orientations among the parameters F, G, and B. In other words, we assess how closely aligned or clearly misaligned any two parameters are. Section 4.3.1 will first search for statistical trends, while Section 4.3.2 will look at the local alignment or misalignment on the maps, with a focus on B and G.

### 4.3.1. Statistical Measure

In order to define relative orientations, the different parameters are matched with each other as follows. For each estimated parameter, we select a nearest estimate of another parameter within a distance of 9'' (~√2/2 × beam size). These two estimates are then defined as one associated pair. Figure 8 shows the relative orientations (ΔPA) for all associated pairs for all three combinations among F, G, and B for the entire NGC 2264 region. All combinations show significant nonuniform distributions,
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Histograms of relative orientations between magnetic field and local gravity in 2264C and 2264D. The data are separated into low- and high-intensity regions with an intensity threshold of 2.0 mJy arcsec⁻². The *p*-value from the Rayleigh's test indicates the similarity between an observed and a uniform distribution, and *p* < 0.05 favors a nonuniform distribution. PRS *Z*<sub>x</sub> indicates whether the relative orientations tend to be parallel (*Z*<sub>x</sub> > 0) or perpendicular (*Z*<sub>x</sub> < 0).</figcaption></figure>

indicated by Rayleigh's test (*p* < 0.05).<fnref n="102" /> The projected Rayleigh statistic (PRS; Jow et al. 2018) is a modification of Rayleigh's test to further examine the alternative hypothesis whether a distribution is rather parallel- (*Z*<sub>x</sub> > 0) or perpendicular-like (*Z*<sub>x</sub> < 0), with *Z*<sub>x</sub> defined as

$$Z_x = \frac{\sum_{i}^{n} \cos \theta_i}{\sqrt{n/2}}, \qquad (6)$$

and the related uncertainties determined by

$$\sigma_{Z_x}^2 = \frac{2\sum_{i}^{n}(\cos \theta_i)^2 - (Z_x)^2}{n}, \qquad (7)$$

where $\theta_i = 2\Delta PA_i$, leading to $\theta_i \in [0, 180^\circ]$ for our measured $PA_i \in [0, 90^\circ]$. This change of variable is sufficient to match the domain of PRS (0°–360°), as cosine is an even function and allows us to examine the parallel ($\theta_i = 0^\circ$) and perpendicular cases ($\theta_i = 180^\circ$) simultaneously.

The PRS results suggest that most of the distributions are parallel-like, except for *B* versus *F* in high-intensity regions. This is consistent with the one-parameter analysis (Section 4.2), which finds very similar prevailing orientations for *F*, *G*, and *B*, implying that, at least statistically, one should also expect closely aligned orientations among these parameters. However, the two-parameter analysis reveals differences once low- and high-density regions are separated. This is very obvious for *B* versus *F*, with a clear change from rather aligned in the diffuse to rather orthogonal in the dense regions, and possibly more subtle for *G* versus *B*, with a possible hint of becoming more aligned in the dense regions (Figure 8). While these differences start to show statistically in these histograms, this also raises the question of whether their origin can be localized and understood in a comprehensive scenario. This will be further explored in the next section and Section 5.2.

If NGC 2264 is further separated into its northern and southern complexes, further nuances can be seen (Figure 9). We focus on *G* versus *B*, as this forms the basis for our analysis in the next sections and Section 5.2. While the diffuse regions in NGC 2264C and *D* show rather flat distributions (similar to the

entire region as seen in Figure 8), the dense regions evolve differently. ΔPA peaks around 0°–20° in 2264C, indicating a statistically growing alignment between magnetic field and gravity, and around 50°–70° in 2264D, indicating a statistically growing misalignment. This difference might be driven by the environment in the two regions, which will be examined with a global stability analysis in Section 4.5. The additional combinations of relative orientations are presented in Appendix E.

### 4.3.2. Local Alignment Maps

Figure 10 displays maps of relative orientations for the three combinations among *F*, *G*, and *B*. We first note that the areas of alignment or misalignment appear systematic and not random. They seem organized, sometimes confined in localized regions, sometimes transitioning toward or along filaments.

Motivated by the question of whether magnetic fields are sufficiently strong to support a system against gravitational collapse, we turn our attention to the ΔPA maps for *G* versus *B* (first row of panels in Figure 10). In 2264D, several stretches along filaments show relative orientations of nearly 90° (red). They appear to have transitioned from surrounding regions where ΔPA is around 0°–50° (blue to green). These regions coincide with those filaments along which the main gravitational converging points are located. A similar transition occurs in the western part of 2264C around the dominating source IRS1. The largest regions of close alignment between *G* and *B* are the northern end of 2264C and the nearby southern end of 2264D. Along the outer boundaries of the combined 2264C and 2264D complex tend to be larger regions of misalignment.

The third row of panels in Figure 10 is motivated by the question of how the observed filamentary structures are possibly shaped by gravity or magnetic field. In both clouds, the ΔPA maps for *G* versus *F* reveal that gravity is often parallel to filaments, especially toward the gravitational converging centers. The ΔPA maps for *B* versus *F* show certain filaments with *B* parallel to *F*, while others show *B* perpendicular to *F*.

This peculiar feature will be discussed more in Section 5.2, where we will argue that the alignment of *G* versus *F* and *B* versus *F* reveals whether gravity and magnetic field tend to radially compress a filament (ΔPA ~ 90°) or pull/guide the gas flow along a filament (ΔPA ~ 0°). With this, the detailed inspection of the relative orientation, i.e., alignment or misalignment between *F*, *G*, and *B* in the maps, can provide insight that is not accessible from the statistical measures presented in the previous section. We finally note that many of the alignment trends seen in the northern end of 2264C appear to connect and continue in the southern part of 2264D.

### 4.4. Three-parameter Analysis: Dependence between Relative Orientations and Local Densities

As noted in Section 4.3.1, the histograms reveal that the trends between *F*, *G*, and *B* appear to evolve with local intensity. To examine whether these trends are statistically significant, we adopt the HRO technique (Soler et al. 2013; Planck Collaboration et al. 2016). This technique probes how the shape of an HRO’s ΔPA changes with local density. The shape is described by the parameter

$$\xi = \frac{A_c - A_e}{A_c + A_e}, \qquad (8)$$

<footnote n="102">¹⁰² Rayleigh’s test is defined in the domain of 0°–360°. Since the range for our *PA* is 0°–90°, we use 4ΔPA as the input of this test. Note that the uniformity of 4ΔPA remains unchanged compared to ΔPA, and as a result, this change of variable preserves the validity of Rayleigh’s test. We note that we additionally also ran Kolmogorov–Smirnov tests on all distributions, and we found identical results, all pointing at nonuniform distributions.</footnote>
<figure><img src="image_12.png" /><figcaption>**Figure 10.** Alignment maps between local gravity (G), magnetic field (B), and filament orientations (F) toward NGC 2264.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>Figure 11. Histogram shape parameter $\xi$ vs. local intensity for all parameter pairs toward 2264C (red) and 2264D (blue). All pairs, except for B vs. F, show a clear tendency of increasing alignment as the local intensity increases until certain turnover points. An opposite trend is seen in B vs. F, which is more aligned in low-intensity regions and becoming more perpendicular in high-intensity regions.</figcaption></figure>

where $A_c$ and $A_e$ are the areas of a histogram around its center ($0^\circ < |\Delta PA| < 22^\circ.5$) and at its extreme ($67^\circ.5 < |\Delta PA| < 90^\circ$). A positive $\xi$ indicates relative orientations close to 0, while a negative $\xi$ indicates perpendicular orientations. The uncertainty $\sigma_\xi$ is obtained as

$$\sigma_{\xi} = \frac{4(A_e^2 \sigma_{A_c}^2 - A_c^2 \sigma_{A_e}^2)}{(A_c + A_e)^4}, \qquad (9)$$

with the variances $\sigma_{A_c}^2$ and $\sigma_{A_e}^2$ of $A_c$ and $A_e$ as

$$\sigma_{A_{c,e}}^2 = h_k(1 - h_k/h_{tot}), \qquad (10)$$

where $h_k$ is the number of data points in the central or extreme bins, and $h_{tot}$ is the total number of data points. The HRO analysis was originally used to study the correlation between magnetic field and intensity gradient (Soler et al. 2013). Here, we further include local gravity and filament orientation in this analysis in order to obtain a more complete picture of the physical parameters.

Figure 11 displays $\xi$ as a function of local intensity. In both 2264C and 2264D, the $\xi$ for these pairs becomes more positive with local intensity in low-intensity regions but transitions to decrease in high-intensity regions. This suggests that these orientation pairs are getting more aligned with growing intensity, but this alignment breaks down in the densest areas. These breaking points vary, with a range of $1-5$ mJy arcsec$^{-2}$ ($N_{H_2} = 3 \times 10^{22}-10^{23}$ cm$^{-2}$) for different pairs, and they occur at higher intensities in 2264C than in 2264D. To examine whether the turnover of $\xi$ in these pairs is occurring for the same filaments, we plot the 2D kernel density distributions (KDEs) of $\Delta PA$ for the G-B and B-F pairs in Figure 12, with a bandwidth of $2^\circ$ determined by Scott's rule. These pairs are selected from filaments with an intensity threshold of $2.0$ mJy arcsec$^{-2}$ to exclude diffuse outskirts regions. The $\Delta PAs$ for G-B and B-F show a positive correlation, seen with a Pearson correlation coefficient of $0.31 \pm 0.02$, estimated using the bootstrap method to account for the observational

uncertainty, and a $p$-value $<0.01$ against the null hypothesis of no correlation. This distribution further reveals two clustering regions: one in the upper right corner ($B \perp F$ and $G \perp B$, which will later be introduced as the type I configuration in Section 5.2) and one in the lower left corner ($B \parallel F$ and $G \parallel B$, hereafter the type II configuration). The possible implications and origins of these configurations are discussed in Section 5.2, which will show how local measures can reveal important features and differences that are unnoticed in global statistics.

### 4.5. Stability of Major Clumps in NGC 2264

In Section 4.3, we have shown that the magnetic field orientations can become perpendicular to local gravity within massive clumps and filaments. In such a configuration, magnetic fields might be important to support clumps and filaments against gravitational collapse, if the magnetic fields are sufficiently strong. To evaluate this further, we estimate the magnetic field strength using the Davis-Chandrasekhar-Fermi (DCF) method (Davis 1951; Chandrasekhar & Fermi 1953). Based on a comparison between the observed polarization angular dispersion $\delta\phi$ and the line-of-sight nonthermal velocity dispersion $\sigma_{v,NT}$, and assuming equipartition of kinetic and magnetic energy, the DCF method yields a plane-of-sky magnetic field strength $B_{pos}$ as

$$B_{pos} = Q \sqrt{4\pi\rho \frac{\sigma_{v,NT}}{\delta\phi}}, \qquad (11)$$

where $\rho$ is the gas volume density, and $Q=0.5$ is a factor accounting for complex magnetic field and inhomogeneous density structures (Ostriker et al. 2001).

We select the polarization segments within the two major clumps in NGC 2264C and D (Figure 3) and calculate the angular dispersion using pairs from only adjacent segments as

$$\delta\phi = \frac{\Sigma(PA_i - PA_j)}{N_{pairs}}, \qquad (12)$$

where $PA_i - PA_j$ is the smaller angle between the two adjacent segments i and j, and $N_{pairs}$ is the total number of adjacent pairs in a clump. This is a simplified version of the angular dispersion function in Hildebrand et al. (2009), aimed at removing the angular dispersion contributed by large-scale magnetic field structures but only focused on the shortest resolved scale. We estimate the total mass of the two major clumps, assuming a constant dust temperature of 15 K (Ward-Thompson et al. 2000) and a dust opacity $\kappa$ of $0.0125 \text{ cm}^2 \text{ g}^{-1}$ at $850 \text{ µm}$ (Hildebrand 1983). The total mass of a clump is then obtained by integrating the column density over the selected region, and the clump volume is estimated using the spherical volume with the effective radius $\sqrt{A/\pi}$, where A is the clump area. From the observed $C^{18}O$ velocity dispersion (Section 3.3), we remove the thermal velocity component and derive the nonthermal velocity dispersion as

$$\sigma_{v,NT}^2 = \sigma_{obs}^2 - \frac{k_B T_{kin}}{m_{C^{18}O}}, \qquad (13)$$

where the kinematic energy $T_{kin}$ is assumed to be 15 K (Peretto et al. 2006).

With the above estimates, the derived magnetic field strengths are $0.24 \pm 0.06$ and $0.22 \pm 0.05$ mG for 2264C and 2264D (Table 1). The mass-to-flux criticality $\lambda_{obs}$, commonly
<table><caption>**Table 1**
Magnetic Field Strengths and Mass-to-flux Ratios in NGC 2264</caption><thead><tr><th>Region</th><th>n<sub>H2</sub><br>(cm<sup>-3</sup>)</th><th>σ<sub>v,NT</sub><br>(km s<sup>-1</sup>)</th><th>δφ<br>(deg)</th><th>B<sub>pos</sub><br>(mG)</th><th>λ</th><th>T/W</th><th>M/W</th></tr></thead><tbody><tr><td>NGC 2264C</td><td>(1.5 ± 0.1) × 10<sup>5</sup></td><td>0.94 ± 0.24</td><td>30.2 ± 0.3</td><td>0.24 ± 0.06</td><td>1.4<sup>+0.4</sup><sub>-0.4</sub></td><td>0.3<sup>+0.2</sup><sub>-0.2</sub></td><td>0.2<sup>+0.1</sup><sub>-0.1</sub></td></tr><tr><td>NGC 2264D</td><td>(5.1 ± 0.4) × 10<sup>4</sup></td><td>1.04 ± 0.23</td><td>21.5 ± 0.4</td><td>0.22 ± 0.05</td><td>0.6<sup>+0.1</sup><sub>-0.2</sub></td><td>0.9<sup>+0.4</sup><sub>-0.4</sub></td><td>0.8<sup>+0.4</sup><sub>-0.4</sub></td></tr></tbody></table>

**Note.** Magnetic field strengths and mass-to-flux ratios are derived from the DCF method. The uncertainties listed here are obtained from propagating the observational uncertainties through the corresponding equations using a Monte Carlo approach. Possible additional systematic uncertainties due to, e.g., the unknown dust opacity are not included.

<figure><img src="image_14.png" /><figcaption>**Figure 12.** 2D KDE of relative orientations between gravity (G)–magnetic field (B) and magnetic field (F)–filament (F). The red dots are the samples estimated from individual pixels. The colored contours show the sample density (blue: largest; white: lowest). A Pearson correlation coefficient of 0.31 ± 0.02 and a p-value < 0.01 (against the null hypothesis of no correlation) indicate a positive correlation. The additional clustering in peak densities localizes the type I (upper right corner) and type II (lower left corner) filaments, as discussed in Section 5.2.</figcaption></figure>

used to evaluate the relative importance between magnetic field and gravity (Nakano & Nakamura 1978), is estimated as

$$ \lambda_{\text{obs}} = 2\pi\sqrt{G} \frac{\mu m_{\text{H}} N_{\text{H}_2}}{B_{\text{pos}}}, \qquad (14) $$

where $\mu = 2.8$ is the mean molecular weight per H$_2$ molecule, G is the gravitational constant, and $N_{H_2}$ is the molecular hydrogen column density. We adopt a statistical average factor of 1/3 (Crutcher et al. 2004) to account for the projection effect for oblate structures flattened perpendicular to the magnetic field. The corrected mass-to-flux ratio $\lambda$ becomes

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

Since the large uncertainties in $\sigma_{v,NT}$ can cause unequal upper and lower uncertainties in $\lambda$, we use a Monte Carlo approach with 10,000 samplings to represent the distribution of the input uncertainties and handle the error propagation. An uncertainty of 7% in the intensity measurements, due to the uncertainty in the FCF (Mairs et al. 2021), is included. The resulting upper and lower uncertainties (Table 1) are defined as the difference between the mean value and the 16th and 84th percentiles. The estimated $\lambda$ are $1.4_{-0.4}^{+0.4}$ and $0.6_{-0.2}^{+0.1}$ for clouds C and D.

In order to evaluate the stability of the major clumps, we perform a virial analysis using the above derived quantities. The virial theorem can be written as

$$ \frac{1}{2} \ddot{I} = 2(T - T_s) + M + W \qquad (16) $$

(e.g., Mestel & Spitzer 1956; McKee & Ostriker 2007), where $I$ is a quantity proportional to the trace of the inertia tensor of a cloud. The sign of $\ddot{I}$ determines the acceleration of the expansion or contraction of the spherical clump. The term

$$ T = \frac{3}{2} M \sigma_{\text{tot}}^2 \qquad (17) $$

is the total kinetic energy, where M is the total mass and $\sigma_{\text{tot}}$ is the total velocity dispersion in a clump, estimated as $\sigma_{\text{tot}}^2 = \sigma_{v,NT}^2 + \sigma_{\text{thermal}}^2$. The thermal velocity dispersion $\sigma_{\text{thermal}}^2$ for mean free particles is calculated assuming a temperature of 2 K and a mass of mean free particles = 2.33. We neglect the surface kinetic term $T_s$ because it cannot be determined from the current observations. We note that the presence of substantial external pressure might still influence a cloud’s instability. The magnetic energy term is

$$ M = \frac{1}{2} M v_A^2, \qquad (18) $$

where $v_A = B/\sqrt{4\pi\rho}$ is the Alfvén velocity and $\rho$ is the mean density as in Equation (11). Since we estimate the magnetic field strengths using the DCF method (Equation (11)), the Alfvén velocity is determined by the polarization angular dispersion $\delta\phi$ and the line-of-sight nonthermal velocity dispersion $\sigma_{v,NT}$. As the DCF method only constrains the plane-of-sky magnetic field component, we use the statistical average to correct and estimate the total magnetic field strength as $B = (4/\pi)B_{\text{pos}}$ (Crutcher et al. 2004). The term

$$ W = -\frac{3}{5} \frac{GM^2}{R} \qquad (19) $$

is the gravitational potential of a sphere with a uniform density $\rho$ and a radius $R$. The derived $T/W$ and $M/W$ are listed in Table 1.

Our mass-to-flux criticality and virial analysis both show that the magnetic energy is smaller than the gravitational and kinematic energy in 2264C ($T: M: W = 0.3^{+0.2}_{-0.1}: 0.2^{+0.1}_{-0.1}$: 1, normalized to $W$). In contrast to this, the constituents are comparably strong within the uncertainties in 2264D ($T: M: W = 0.9^{+0.4}_{-0.4}: 0.8^{+0.6}_{-0.4}$: 1). This could explain the difference in distributions of $\Delta$PA for G–B in 2264C and 2264D (Figure 9); since 2264C is dominated by gravity, the magnetic field is expected to be more aligned (parallel) with gravity toward
higher densities. On the other hand, since none of the constituents dominates in 2264D, the magnetic field and gravity do not simply appear parallel or perpendicular. As the stability analysis here is based on global and averaged properties (Table 1), a different analysis (Section 5.2) is still necessary to shed light on the more detailed and local role of gravity and magnetic field.

We note that a major source of uncertainty in our estimates is the possible complexity of velocity structures. Multiple velocity components, possibly tracing the low-density outskirts, are identified in $^{13}$CO. On the other hand, $C^{18}$O appears to trace the major high-density components, which are more likely associated with the regions traced by POL-2. In Appendix B, we show that the velocity dispersion from $^{13}$CO is systematically larger than the one from $C^{18}$O by 30%–50%. Additionally, the velocity dispersion estimated from $N_2H^+$ (1–0) lines in Peretto et al. (2006), tracing higher densities, is systematically lower than our estimates by ~30%. This suggests that the lower a density is traced, the larger the velocity dispersion can be. This is possible if the diffuse gas is distributed over a larger volume and thus mixed with more complex velocity structures. Given that the polarization efficiency, determined by unknown dust properties and radiation field together with the SCUBA-2/POL-2 large-scale filtering effect (Sadavoy et al. 2013), is also rather unclear, it is not obvious which gas tracer is the best to use together with continuum polarization in the DCF method. To illustrate these caveats, we present tables of the estimated magnetic field strengths, mass-to-flux ratios, and energy ratios for all structures identified by the dendrogram and for the three different gas tracers ($C^{18}$O, $^{13}$CO, $N_2H^+$) in Appendix F. Nevertheless, regardless of whether the $^{13}$CO or $C^{18}$O line is chosen, our conclusion regarding the relative importance of energy sources remains valid.

Finally, we note that outflow activities were extensively observed in NGC 2264 (e.g., Cunningham et al. 2016). Maury et al. (2009) conducted an analysis focusing on the equilibrium between the force exerted by outflows in 2264C and the gravitational force. Their findings indicate that these outflows fall short by a factor of ~25 in terms of providing significant support against the gravitational collapse of the 2264C protocluster. Consequently, it is unlikely that the outflows will have a substantial influence on the stability analysis presented here.

# 5. Discussion

## 5.1. Statistical Measures of the Filamentary Network

The filamentary network in NGC 2264 is embedded in the hub of a 10 pc–scale hub–filament system (Kumar et al. 2020). This network is therefore denser ($N_{H_2}$ = $10^{22}$–$10^{24}$ cm$^{-2}$) and more compact (length = 0.2–1.0 pc) than the typical filamentary networks discovered by Herschel, i.e., $N_{H_2}$ = $10^{21}$–$10^{22}$ cm$^{-2}$ and length = 1–10 pc in IC 5146 (Arzoumanian et al. 2011). A major scientific question that we are addressing in the following is how the network in NGC 2264 is formed and shaped.

### 5.1.1. One-parameter Statistics: Global Trends from Prevailing Orientations

The one-parameter statistics reveal that filaments, gravity, and magnetic fields share a common prevailing orientation

(Section 4.2, Figure 7). It appears that the filamentary network primarily consists of filaments with orientations approximately perpendicular to each other (~30° and ~150°). This arrangement is similar to the so-called main-filament and subfilament (or striation) structures found in other filamentary clouds, such as Serpens South (Sugitani et al. 2011) and Taurus B211/213 (Palmeirim et al. 2013). The one prevailing orientation in gravity (around ~40°–50°) indicates one dominating converging direction for gravity along a ~40°–50° orientation. The finding that one of the two prevailing filament orientations, namely, the one around ~150°, is approximately perpendicular to the converging gravity direction is suggestive of a picture where gravity (on a global scale across both 2264C and 2264D) is predominantly driving mass accretion along this direction (Figure 13). The second prevailing filament orientation around ~30° aligns roughly with the prevailing gravity orientation, though gravity shows a broader distribution in orientations. This broader distribution (also seen in the intensity gradient) likely is the result of small (local) deviations of gravity from its prevailing global direction toward the filaments.

The B-field orientations are again around a prevailing orientation of ~40°–50°, similar to gravity yet with an even broader distribution. This broadening with respect to gravity is likely capturing a lag in alignment between B-field and gravity. The broader wings in the distribution likely result from the more diffuse regions where the B-field orientation is less aligned with gravity. This is also seen from the two-parameter statistics discussed in the next section. Since the polarization measurements are independent of the intensity maps (which are used for the filament extraction and the derivation of gravity), these similar orientations around ~50° imply that the magnetic field morphology is strongly connected to the density structures. This observed filament–magnetic field configuration can initially be generated by either a cloud collapsing/fragmenting along the magnetic field (Nakamura & Li 2008; Li et al. 2013; Van Loo et al. 2014) or a bow-shaped magnetic field associated with filaments formed via the compression of ISM bubbles (Inutsuka et al. 2015; Inoue et al. 2018; Tahani et al. 2019, 2022a, 2022b). The cartoon in Figure 13 summarizes the suggested scenario based on these global trends.

### 5.1.2. Two-parameter Statistics: Global Alignment Statistics

The relative orientation analysis (Figure 8) suggests that filament orientation and local gravity are both correlated with the magnetic field orientation; i.e., their relative alignment is not random but systematic. Moreover, low- and high-intensity regions appear different, with a trend that denser regions display a closer alignment, i.e., a shift in the distributions toward smaller $\Delta$PA. The only exception is the alignment between the magnetic field and filaments, which shows an opposite trend, the magnetic field orientation becoming more orthogonal to the filaments in denser regions. Finally, all these trends among all parameters also remain when binned to finer intensity steps, except for the turnover trends only seen in regions with the highest densities (Figure 11). This indicates that these alignment correlations are a function of the local density.

Previous statistical analyses (e.g., Planck Collaboration et al. 2016; Kwon et al. 2022) based on polarization observations toward a number of filamentary systems have also shown that magnetic fields tend to be parallel to less-dense filamentary
<figure><img src="image_15.png" /><figcaption>**Figure 13.** Cartoon illustrating the correlation between filaments, gravity, and magnetic field in NGC 2264. The green lines, red arrows, and yellow segments represent the filament ridges, direction of gravity, and magnetic field orientation. The bottom right corner shows a zoom-in to the filament ridge in the black rectangle. The right panels are the expected histograms of the one- and two-parameter statistics for this illustrated picture. All three parameters share a common prevailing orientation around ~50°, while the filaments show an additional prevailing orientation around ~140°. Gravity displays a broader distribution due to the feature in the zoom-in. The magnetic field distribution is even wider due to the 0°/180° wrapping.</figcaption></figure>

structures but become perpendicular toward denser filaments. As the filaments identified in our work are distinct from filaments identified using Planck and Herschel data, our results suggest that the alignment in NGC 2264 transitions in 0.1 pc–scale regions with much higher densities (~$10^4–10^5$ cm$^{-3}$ or $10^{22}–10^{23}$ cm$^{-2}$) than the transition density of $10^{21}–10^{21.5}$ cm$^{-2}$ reported in Planck Collaboration et al. (2016).

In order to explain the observed transition in alignment, recent numerical simulations have proposed a scenario where such a transition originates from changing physical processes in the filament formation, namely, from the phase dominated by magnetically regulated MHD turbulence in low-density regions to the phase dominated by gravitational collapse or converging flow compression in high-density regions (e.g., Soler & Hennebelle 2017; Seifried et al. 2020; Ibáñez-Mejía et al. 2022; Ganguly et al. 2023). However, these simulations predict a wide range of transitional density thresholds, from $10^2–10^3$ cm$^{-3}$ (Seifried et al. 2020; Ganguly et al. 2023), to $10^3–10^4$ cm$^{-3}$ (Ibáñez-Mejía et al. 2022), to $10^4–10^5$ cm$^{-3}$ (Soler & Hennebelle 2017). Our results favor a higher threshold of ~$10^4–10^5$ cm$^{-3}$.

## 5.2. Local Measures of the Filamentary Network: Two Different Types of Filaments?

The analysis presented in Section 5.1 has revealed statistical trends globally across the entire NGC 2264 region. It is evident from the alignment maps in Figure 10 that the relative orientations $\Delta$PA between any two parameters are not random. They seem to appear in certain spatial locations, possibly in organized structures and possibly even with characteristic length scales. The detailed local connections among the magnetic field, gravity, and filament are largely unexplored. In this section, we focus on a particular aspect, namely, the combined magnetic field–gravity connection in shaping

filaments. This is unseen in global statistics but apparent in local structures.

Motivated by the two types of configurations revealed in the $\Delta$PA distribution of B–F and B–G (Figure 12), we extract two representative regions in NGC 2264 associated with these two types (Figure 14). Configuration I (left panels) shows local gravity (G) converging from the lower-density surroundings toward the filament (F), and then once close enough, G aligns with the filament ridge. On the resolved scale in these data, the magnetic field (B) remains perpendicular to the filament. The magnetogravitational configuration shows G parallel to B at larger distances from F but then transitioning to G becoming increasingly perpendicular to B closer to F. This configuration is seen in both lower- and higher-density areas. Configuration II (right panels) shows local gravity being mostly along the filament and not yet converging toward the filament ridge (at the resolved scale in these data). Additionally, the magnetic field is also predominantly aligned with both local gravity and the filament. Here, the magnetogravitational configuration shows G parallel to B both at larger distances and close to F. This configuration is also seen in both lower- and higher-density areas.

In the following, we discuss the possible implications and origins of these two different magnetogravitational configurations. These two systematically different configurations might be suggestive of two different types of filament. A quantitative measure to assess differences in these constellations is the $\sin\omega$ measure (Koch et al. 2018, 2022). With the angle $\omega$ between B and G, $\sin\omega$ quantifies a fraction (between 0 and 1) of the magnetic field tension force that can work against a local gravitational pull. Small values in $\omega$, i.e., close alignment between gravity and magnetic field, allow gravity to maximally accelerate inflowing material. Large values in $\omega$, i.e., gravity and magnetic field close to orthogonal, indicate maximal
<figure><img src="image_16.png" /><figcaption>**Figure 14.** Illustration of the role of gravity and magnetic field in two types of filaments (type I, left panels; type II, right panels). The two maps in the top row are two representative regions as observed in NGC 2264C and D, where the green lines, red arrows, and yellow segments are the filament, gravity, and magnetic field orientations. Middle row: for type I filaments, the radial component of the gravitational field ($\partial U/\partial r$) is greater than or comparable to the longitudinal component ($\partial U/\partial \ell$) in the outer regions, but the longitudinal component starts to dominate at the filament ridges. The magnetic field is therefore parallel to gravity in the outer regions but becomes more perpendicular toward the ridges. In contrast to this, the gravitational field in type II filaments is always dominated by the longitudinal component. Bottom row: $\sin \omega$ maps for the two selected regions in the top row.</figcaption></figure>

obstruction by the magnetic field and hence a reduced acceleration. Consequently, the type I filament (left panels in Figure 14) will initially experience fast accretion from the outside (small $\omega$; middle panels in Figure 14) and hence short radial timescales toward the filament but then have a longer longitudinal timescale along the filament ridge (larger $\omega$). The type II filament (right panels) will favor a short longitudinal timescale for any motion or collapse along the filament (small $\omega$). On a more speculative note, these findings might further imply that the type I filament is able to accrete more material.

The two different types of filament and magnetogravitational configurations can also be explained in terms of their gravitational potentials (middle panels in Figure 14). The gravitational potential $U$ of a filament within an HFS is expected to be a valley. Radially, within a filament

($r < \text{filament width}$), $U$ becomes shallower toward its bottom at the ridge (e.g., a filament with a Plummer-like density profile; Arzoumanian et al. 2011). Longitudinally, driven by the global gravitational field influenced by the nearest massive hub, $U$ smoothly decreases toward the major potential well. Whether material is accreted onto a filament is determined by the competition between the radial and longitudinal components of the gradient of the gravitational potential. The type I filament (left panels in Figure 14) will have a radial gravitational potential profile ($\partial U/\partial r$) that is steeper (as a result of a higher density or a steeper density profile) than its longitudinal gravitational potential profile ($\partial U/\partial \ell$; as a result of mass and distance to the major potential well). Thus, this configuration favors accretion from the ambient material. Once
<figure><img src="image_17.png" /><figcaption>**Figure 15.** YSO KDE maps for class 0/I YSOs (left panel) and class II/TD sources (right panel). The cyan lines are the identified filaments. The yellow and green squares label the positions of class 0/I and class II/TD sources. The obvious YSO clusters are labeled as clusters 1–10.</figcaption></figure>

approaching the filament ridge, where the radial gravitational potential approaches its bottom and $\partial U/\partial r$ becomes 0, the longitudinal motion along the filament will start to dominate. This then causes a transition from small to larger $\omega$ at the observed resolution. In contrast to this, type II filaments (right panels in Figure 14) have a radial gravitational potential profile that is shallower than the longitudinal profile. Therefore, the pull from the global gravitational field is always more efficient. This picture has also been a subject of discussion in theoretical studies of hierarchical collapse within molecular clouds (e.g., Gómez & Vázquez-Semadeni 2014; Vázquez-Semadeni et al. 2019), where both longitudinal and radial gravitational accretion modes of filaments, linked to the global and local cloud collapse, can occur simultaneously and impact the filaments’ stationary, growing, or dissipating states (Vázquez-Semadeni et al. 2019; Naranjo-Romero et al. 2022).

The role of magnetic fields can easily be incorporated into the above picture. Since the magnetic field tension force is perpendicular to field lines, a magnetic field component perpendicular to a filament will suppress the longitudinal potential, while a component parallel to a filament will suppress the radial potential. In other words, a magnetic field more perpendicular to a filament will favor the type I filament, whereas a field more parallel to a filament will favor the type II filament. This reveals itself as the positive correlation we found in the ΔPA of B–F and G–B pairs in Figure 12, which identifies the type I and type II filaments.

In short, the two different magnetogravitational configurations have different effective ways of driving motions, either more effectively perpendicular to a filament (type I) or more

effectively parallel to a filament (type II). Based on all of the above, we propose two different types of filament resulting from two different magnetogravitational configurations.

### 5.3. Star Formation in NGC 2264

YSOs are observed to form within filaments (André et al. 2014; Pineda et al. 2023). It is thus a question of which type of filament tends to form stars. In order to trace the star formation activity in NGC 2264, we adopt the YSO catalog in Rapson et al. (2014), which is based on the Spitzer and Two Micron All Sky Survey data. We select the YSOs classified as class 0/I, class II, and transitional disks (TD) from this catalog and plot the YSO KDE for the young YSOs (class 0/I) and evolved YSOs (class II and TD) in Figure 15. We do not include the class III YSOs, because those diskless YSOs are difficult to distinguish from field stars using only infrared data (Rapson et al. 2014). Since this YSO catalog covers an area greater than NGC 2264, we only select YSOs within R.A. = 6h40m36s–6h41m37s and decl. = 9°22′23″–9°42′11″.

The KDE uses a Gaussian kernel with an FWHM width of 200". This is chosen because it is comparable to the median proper-motion velocity of 1.1 mas yr⁻¹ (Buckner et al. 2020) times the freefall dynamical timescale of 1.7 × 10⁵ yr for the NGC 2264 clumps (Peretto et al. 2006). In this way, this represents the spatial probability distribution of each YSO, accounting for its possible migration. We note that this FWHM width is merely an order-of-magnitude estimate. A more accurate YSO migration would need to be described with an N-body model also considering the interaction between YSOs.
<figure><img src="image_18.png" /><figcaption>**Figure 16.** Cumulative histogram of distances between YSOs and the nearest filament. The blue and orange histograms are for class 0/I and class II/TD, respectively. The 90th percentile of the nearest distance is 0.07 pc for class 0/I and 0.45 pc for class II/TD YSOs.</figcaption></figure>

The derived YSO KDE maps show that the distributions of young and evolved YSOs are clearly different (Figure 15). The young YSOs are more tightly clustered. In both 2264C and 2264D, the class 0/I YSOs are clustered closely either on filaments or on filament converging points. In contrast to this, the overall distribution of the evolved YSOs is more scattered. The few evolved YSO clusters that can still be identified have their density peaks offset from filaments. In 2264C, these clusters are located to the north of the major filaments (cluster 8), and only the two clusters (clusters 9 and 10) in the southern 2264D are likely associated with filaments. In 2264D, evolved YSOs (cluster 7) are concentrated in a filament-surrounded cavity, where no filaments are identified at the center.

In order to further confirm that the evolved YSOs are more distant from filaments, we show cumulative histograms of the nearest distance between YSOs and filaments in Figure 16. The 90th percentile of the nearest distance is 0.07 pc for class 0/I YSOs and 0.45 pc for class II/TD YSOs. This suggests that nearly all class 0/I YSOs are located very close to filaments, while evolved YSOs tend to be more distant from filaments, up to 0.5–1 pc in some cases.

In the following, we discuss two possibilities that can cause offsets between evolved YSOs and filaments: (1) YSO migration and (2) change of filamentary structures.

(1) Due to their older age, evolved YSOs can migrate further from their primordial birth sites than younger YSOs, leading to a spatially more scattered distribution. This possibility has been used to explain the different spatial distribution between protostars and disk sources in Southern Orion A (Mairs et al. 2016). Buckner et al. (2020) found a correlation between the evolutionary stage and source concentration in NGC 2264: 69.4% of class 0/I, 27.9% of class II, and 7.7% of class III objects are found to be clustered. Based on this, migration of evolved YSOs was used to explain the differences in concentration in NGC 2264 (Buckner et al. 2020). Nevertheless, no information on cloud structure was included in their study. In addition, the observed offset direction is not fully consistent with the direction of proper motions of the evolved YSOs. For example, the proper-motion direction of the evolved

YSOs is ~240°–270° in cluster 8 and ~270°–290° in cluster 7 (Buckner et al. 2020), which likely does not trace back to existing filaments.

(2) The evolved YSOs might have formed from filaments that have been dissipated. Inutsuka et al. (2016) argue that present YSOs can have partially formed within filaments that existed in the past. This is suggested because the total mass of the observed YSOs often exceeds the product of the total filament mass and the star formation efficiency. Since the freefall dynamical timescale (~2 × 10⁵ yr: Peretto et al. 2006) of the clumps in NGC 2264 protoclusters is shorter than the typical age of class II YSOs (1–2 Myr; Evans et al. 2009), those filaments or clumps forming the evolved YSOs might have dissipated in between. In addition, the star formation feedback, such as outflows or expanding shells, might also clear the primordial local density structures. This possibility is supported by the presence of cavities associated with evolved YSO clusters (e.g., clusters 7 and 10), and it also aligns with the scenario proposed by Kumar et al. (2020), which suggests that NGC 2264 originated from a collision between filaments on parsec scales. In this scenario, the first-generation massive stars are formed at the filament collision center (cluster 8). The feedback from the massive star subsequently clears the material in the gap between 2264C and 2264D, while 2264C and 2264D are the remnants of the filaments that collided. Finally, it is worth mentioning that according to Nony et al. (2021), the weakly embedded class II sources (such as cluster 10) have the potential to migrate ~6 pc if they move freely without being bound to any structures. Thus, those class II sources that are significantly separated from filaments could also originate from distant regions.

## 6. Conclusions

This study utilizes JCMT SCUBA-2/POL-2 850 µm continuum polarization observations to investigate the hub-filament system NGC 2264. An organized magnetic field morphology is uncovered with the following main outcome.
1. Orientations of filaments (F), local gravity (G), and magnetic field (B) are estimated. These spatial parameters are found to share common prevailing orientations around 20°–50°, with only the filaments showing a second prevailing orientation around 150°. This suggests that the evolution of these four parameters is possibly regulated by the same physical processes.

2. The pairwise analysis of the relative orientations between these parameters indicates that they are locally correlated. Most of the parameter pairs tend to become more tightly aligned toward higher intensities, suggesting that these parameters are likely aligned by the increased gravity. A transition of the relative orientation becoming more perpendicular at very high densities is found in the HRO diagram. This is likely associated with the transition between two types of filaments.

3. A significant number of the filaments in NGC 2264 can be categorized into two distinct groups. This is based on their alignment properties as revealed in the 2D KDE map for B×F versus B×G: B⊥F and B⊥G (type I filament) and B||F and B||G (type II filament).

4. We examine the spatial arrangements of the magneto-gravitational configuration locally on the maps associated with the two types of filaments. These maps reveal that the magnetogravitational configuration in type I filaments shows *G* parallel to *B* at larger distances from *F* but then transitions to *G* becoming increasingly perpendicular to *B* closer to *F*. In contrast to that, the type II filaments display *G* parallel to *B* both at larger distances and close to *F*.

5. We propose a picture to explain the observed features of the two types of filaments. This involves the competition between the radial and longitudinal collapsing timescale. We interpret type I filaments as potentially being able to accumulate more nearby material, while type II filaments are instead shaped by a dominating longitudinal pull from the larger-scale global gravitational field.

6. Magnetic field strengths are estimated for the major clumps in 2264C and 2264D, together with a virial analysis. The estimated ratios of kinematic energy (T): magnetic energy (M):gravitational energy (W) are 0.3$^{\pm0.2}_{-0.1}$; 0.2$^{\pm0.1}_{-0.1}$; 1 in 2264C and 0.9$^{\pm0.4}_{-0.4}$; 0.8$^{\pm0.4}_{-0.4}$; 1 in 2264D. This indicates that the global evolution of 2264C is dominated by gravity, while the three factors appear comparably important in 2264D.

7. The class 0/I YSOs are found to be closer in distance to the identified filaments or filament converging points. Most of the class II and TD YSOs are located at systematically larger distances from filaments.

## Acknowledgments

The James Clerk Maxwell Telescope is operated by the East Asian Observatory on behalf of the National Astronomical Observatory of Japan, the Academia Sinica Institute of Astronomy and Astrophysics, the Korea Astronomy and Space Science Institute, the National Astronomical Research Institute of Thailand, and the Center for Astronomical Mega-Science (as well as the National Key R&D Program of China with No. 2017YFA0402700). Additional funding support is provided by the Science and Technology Facilities Council of the United Kingdom and participating universities and organizations 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 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. and P.M.K. acknowledge support from the National Science and Technology Council (NSTC) in Taiwan through grants NSTC 112-2112-M-001-049 and NSTC 111-2112-M-001-039. K.Q. is partially supported by the National Key R&D Program of China Nos. 2022YFA1603100 and 2017YFA0402604. K.Q. acknowledges National Natural Science Foundation of China (NSFC) grant U1731237 and the science research grant from the China Manned Space Project with No. CMS-CSST-2021-B06. F.P. acknowledges support from the Agencia Canaria de Investigación, Innovación y Sociedad de la Información (ACIISI) under the European FEDER (FONDO EUROPEO DE DESARROLLO REGIONAL) de Canarias 2014-2020 grant No. PROID2021010078. The work of M.G.R. is supported by NOIRLab, which is managed by the Association of Universities for Research in Astronomy (AURA) under a cooperative agreement with the National Science Foundation. L.F. acknowledges support from the Ministry of Science and Technology of Taiwan under grant Nos. 111-2811-M-005-007 and 109-2112-M-005-003-MY3. D.J. is supported by NRC Canada and an NSERC Discovery Grant. K.P. is a Royal Society University Research Fellow, supported by grant No. URF\R1\211322. C.E. acknowledges the financial support from grant RJF/2020/000071 as a part of the Ramanujan Fellowship awarded by the Science and Engineering Research Board (SERB), Department of Science and Technology (DST), Government of India. C.W.L. was supported by the Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Education, Science and Technology (NRF-2019R1A2C1010851) and by the Korea Astronomy and Space Science Institute grant funded by the Korean government (MSIT; project No. 2023-1-84000). W.K. was supported by the National Research Foundation of Korea (NRF) grant funded by the Korean government (MSIT; NRF-2021R1F1A1061794). M.T. is supported by JSPS KAKENHI grant No. 18H05442.

Facility: JCMT.

Software: Aplpy (Robitaille & Bressert 2012; Robitaille 2019), Astrodendro (Goodman et al. 2009), Astropy (Astropy Collaboration et al. 2013, 2018), BTS (Clarke et al. 2018), DisPerSE (Sousbie 2011), Filchap (Suri & Sánchez-Monge 2019), NumPy (Oliphant 2006), SciPy (Virtanen et al. 2020), scousePy (Henshaw et al. 2016a, 2019), Smurf (Berry et al. 2005; Chapin et al. 2013), Starlink (Currie et al. 2014).

## Appendix A
### Polarization Properties

Figure 17 shows the 850 µm POL-2 polarization map, where the lengths of the segments are scaled with the polarization fraction. The resulting histogram is reported in Figure 18. The polarization fraction, with a median value of 2.4%, is mostly below 10%. Values of 10%–20% are rare, typically occuring near the edge of the clouds.

The correlation between Stokes *I* and polarization fraction *P* is commonly used to examine whether the observed
<figure><img src="image_19.png" /><figcaption>**Figure 17.** POL-2 polarization segments overlaid on 850 µm dust continuum toward NGC 2264. The yellow segments show polarization detections selected by $I/\sigma_I > 20$ and $P/\sigma_P > 3$. The lengths of the segments are scaled with the polarization fractions. The yellow scale bar at the top left shows a 20% polarization fraction. The red and yellow stars label IRS1 (zero-age main-sequence) and IRS2 (class I), the dominating sources in 2264C and 2264D.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>**Figure 18.** Histogram of the debiased polarization fraction of the selected samples. Most of the segments show a polarization fraction smaller than 10% with a median value of 2.4%.</figcaption></figure>

polarization originated from aligned dust grains. This is essential to confirm that the polarization measurements can trace the magnetic fields within molecular clouds (e.g., Andersson et al. 2015; Jones et al. 2015, 2016; Planck Collaboration et al. 2015, 2020; Pattle et al. 2019; Wang et al. 2019). The $I-P$ relation can be described by a power law, $P = PI/I \propto I^{-\alpha}$, where $PI$ is the polarized intensity. If the dust grains are not aligned with the magnetic field, $PI$ is independent of Stokes $I$, and $\alpha$ is 1. An $\alpha$ smaller than 1 indicates that the polarization fraction results from magnetically aligned dust grains.

<figure><img src="image_21.png" /><figcaption>**Figure 19.** 850 µm total intensity $I$ vs. polarization fraction $P$. The green points are the nondebiased POL-2 polarization measurements, selected by $I/\sigma_I > 10$ and $\sigma_I < 0.02$ mJy arcsec$^{-2}$. The colored regions indicate the 68%, 95%, and 98% CRs predicted by the Bayesian model. The black line indicates the posterior mean. Most of the data points are within the 98% CR of our prediction.</figcaption></figure>

The observed $I-P$ relation for NGC 2264 is shown in Figure 19. We follow the Bayesian analysis in Wang et al. (2019) to determine the power-law index $\alpha$ using the model

$$P = \beta I^{-\alpha}. \qquad (A1)$$

The Bayesian model allows us to directly model the Ricean distribution (Wardle & Kronberg 1974) as the error function of the observed $P$,

$$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 (A2)$$

where $P$ is the observed polarization fraction, $P_0$ is the real polarization fraction, $\sigma_P$ is the Ricean dispersion in the polarization fraction, and $I_0$ is the zeroth-order modified Bessel function. We further assume that the uncertainty in the polarization fraction is dominated by the uncertainties in Stokes $Q$ and $U$,

$$\sigma_P = \sigma_Q/I, \qquad (A3)$$

where $\sigma_Q$ is the dispersion in Stokes $Q$, assuming the same noise in Stokes $Q$ and $U$.

We reselect the polarization data based on the criteria $\sigma_I < 0.02$ mJy arcsec$^{-2}$ and $I/\sigma_I > 10$ in order to have sufficient data points to model the noise-dominated regime. Additionally, the nondebiased polarization fractions are used because the Ricean noise is well accounted for in this model. The $I-P$ relation predicted by our Bayesian model is shown in Figure 19. Most samples are well covered by the predicted 98% confidence region (CR). Nevertheless, some data points with polarization fractions higher than predicted by the model can be seen in the high-intensity regions ($I > 2.0$ mJy arcsec$^{-2}$). This can be due to different dust properties or additional radiative alignment torques caused by the embedded sources. The computed posterior distributions of the Bayesian model are reported in Figure 20. $\alpha$ is constrained to $0.29 \pm 0.01$, which is significantly smaller than 1. Consequently, this suggests that the observed polarization likely originates from magnetically aligned dust grains and hence traces the magnetic field morphology within NGC 2264.
<figure><img src="image_22.png" /><figcaption>**Figure 20.** Posterior distributions of our Bayesian model for *I-P* distribution. The black lines delineate the 68% (1σ) highest density interval (HDI). The most probable α of 0.29 ± 0.01 is significantly less than unity, suggesting that the observed polarization likely originates from aligned dust grains in NGC 2264.</figcaption></figure>

## Appendix B
### $^{13}$CO (3–2) and C$^{18}$O (3–2) Molecular Lines

In this section, we discuss the possible influence of the choice of the $^{13}$CO (3–2) or C$^{18}$O (3–2) line to estimate the velocity dispersion. Figure 21 shows example spectra of $^{13}$CO (3–2) in NGC 2264 at different positions. In most of the locations in NGC 2264, we identified multiple velocity components or asymmetric line profiles in the spectra. The separations of those multiple components are usually less than or comparable to their line widths.

The fitted amplitude of the velocity components in $^{13}$CO (3–2) and C$^{18}$O (3–2) are reported in Figure 22. These components are matched if their centroid velocity is close within twice the fitted uncertainties in the velocity component identification. The black dashed line represents the best-fit

<figure><img src="image_23.png" /><figcaption>**Figure 21.** $^{13}$CO (3–2) integrated intensity map with spectra extracted at different pixels. Double or triple velocity components separated by only a few km s⁻¹ are seen.</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>**Figure 22.** Comparison of the fitted peak amplitudes of the matched velocity components in $^{13}$CO (3–2) and C$^{18}$O (3–2). C$^{18}$O (3–2) is likely optically thin because the brightness temperature is much lower than the possible gas temperature in molecular clouds. The dashed line marks the best-fit intensity ratio of 3.7. Most points follow the fitted intensity ratio, except for a few pixels where the $^{13}$CO amplitude is above 10 K. This suggests that $^{13}$CO (3–2) is likely optically thin in most of the regions.</figcaption></figure>

$T(^{13}\text{CO})/T(\text{C}^{18}\text{O})$ ratio of 3.7. Most points follow the fitted line well, suggesting that $^{13}$CO (3–2) and C$^{18}$O (3–2) are both optically thin.

Figure 23 shows a comparison between the velocity dispersion estimated from the two tracers $^{13}$CO (3–2) and C$^{18}$O (3–2). These velocity components are matched as in Figure 22 but with an additional criterion of $2.5 < T(^{13}\text{CO})/T(\text{C}^{18}\text{O}) < 5.5$ to ensure that the components are really associated. The $^{13}$CO velocity dispersion is systematically larger by 30%–50%. We speculate that $^{13}$CO might trace more of the relatively low-density outskirts of NGC 2264, which could be affected by the large-scale accreting filaments, hence being more turbulent and leading to larger velocity dispersion. Therefore, we conclude that C$^{18}$O is more reliable to trace the velocity dispersion within the NGC 2264 clumps. The five outliers with $\sigma_V(\text{C}^{18}\text{O})/\sigma_V(^{13}\text{CO}) > 1.2$ are from the faint pixels, where the velocity components are not well-constrained by the C$^{18}$O data, and thus we flagged the five pixels from calculating the mean velocity dispersion in the DCF analysis.
<figure><img src="image_25.png" /><figcaption>Figure 23. Comparison of the fitted velocity dispersion of the matched velocity components in $^{13}$CO (3–2) and C$^{18}$O (3–2). The dashed line marks where the two dispersion values are equal. The velocity dispersion values extracted from $^{13}$CO are systematically wider than the ones from C$^{18}$O by 30%–50%.</figcaption></figure>

Appendix C

Filament Properties: Width, Line Mass, and Filament-to-cloud Fraction

We derive the properties of the identified filaments by fitting the radial intensity profile at each pixel along the filaments using the bootstrap method. A fitting result is selected if (1) one major peak is found, and (2) the fitted width is above $3\sigma$. We further derive the central column densities based on the fitted amplitude (assuming a constant dust temperature of 15 K, Ward-Thompson et al. 2000, and a dust opacity $\kappa$ of 0.0125 cm$^2$ g$^{-1}$ at 850 µm, Hildebrand 1983) and a line mass ($M_{\text{line}}$; also known as mass per unit length) by multiplying the central column density with the filament width. Figure 24 shows a histogram of the FWHM filament width. The mean and median filament widths are both 0.13 pc, which is very similar to the characteristic width of 0.1 pc found in Arzoumanian et al. (2011). The shape of the filament width distribution with a rapid decline to larger widths is consistent with the filament formation via supersonic turbulence (Priestley & Whitworth 2022). Figure 25 plots the line mass along each filament. A gradual increase is seen, ranging from 10 to 50 $M_{\odot}$ pc$^{-1}$ at the outskirts to over 100 $M_{\odot}$ pc$^{-1}$ at the converging points. Notably, the main filament within 2264C has the highest line mass, reaching

approximately 600 $M_{\odot}$ pc$^{-1}$. This suggests that these filaments are actively accumulating ambient material and possibly transport the accumulated mass toward the converging points.

The critical filament line mass (Fiege & Pudritz 2000) can be estimated by

$$M_{\text{line,critical}} = 2(c_s^2 + \sigma_{v,NT}^2)/G, \quad (C1)$$

where $c_s$ is the sound speed and $\sigma_{v,NT}$ is the nonthermal velocity dispersion. Assuming a nonthermal velocity dispersion of 1.0 km s$^{-1}$ and a gas temperature of 15 K, the $M_{\text{line,critical}}$ is about 400 $M_{\odot}$ pc$^{-1}$. This is higher than the derived line mass in most of the regions. However, it is important to note that the critical line mass is determined for filaments assumed to be in hydrostatic balance. This might not be the case for most of the observed filaments with accretion onto or along them as indicated by velocity gradients (e.g., Peretto et al. 2013; Storm et al. 2016; Williams et al. 2018; Zhou et al. 2022). These motions could then change the energy balance. In addition, HFSs, with the presence of dense hubs, behave as major gravitational potential wells. Thus, the gravitational potential in a filament is a superposition of both the filament’s self-gravity and the global mass distribution. This is much more complex than the case of an isolated, static filament. Consequently, a dynamic filament model is desirable to more accurately describe the stability of this system.

Assessing the proportion of mass within a cloud that resides in filaments is essential to evaluate the final star formation efficiency (Inoue & Inutsuka 2012; Inutsuka et al. 2016). To accomplish this, we employ two different methods to calculate the total mass of the identified filaments. First, we integrate the products of $M_{\text{line}}$ and $\Delta L$, where $\Delta L$ represents the length of each filament segment. According to this calculation, the total filament mass amounts to 757 $M_{\odot}$, which corresponds to 26% of the total cloud mass of 2960 $M_{\odot}$ (Peretto et al. 2006). Second, we calculate the total mass associated with the pixels belonging to filaments, selected if pixels are within the fitted $2\sigma$ width. To estimate the large-scale background, we use a Gaussian kernel with a $\sigma$ of 5 pixels to interpolate the intensity map of the unselected pixels. After subtracting the contribution of the large-scale background, the total filament mass is determined to be 567 $M_{\odot}$. This is equivalent to 19% of the total cloud mass. While the first method might be an upper limit (because filaments partially overlap near the converging points), the second method

<figure><img src="image_26.png" /><figcaption>Figure 24. Histogram of FWHM filament widths estimated from individual intensity radial profiles. The red dashed line is the median width of 0.13 pc.</figcaption></figure>
<figure><img src="image_27.png" /><figcaption>Figure 25. Maps of filament line mass. The color indicates the line mass at each position, estimated from the radial intensity profile. The filament line mass gradually increases from a range of [10, 50] $M_{\odot}$ pc$^{-1}$ in the outskirts to above 100 $M_{\odot}$ pc$^{-1}$ in the converging points. The main filament in 2264C has the highest line mass of ~600 $M_{\odot}$ pc$^{-1}$.</figcaption></figure>

might be a lower limit (due to a possible oversubtraction of the background). We note that the major uncertainty in this estimation is the large-scale filtering effect. This is because the total mass recovered by our POL-2 continuum map is 1150 $M_{\odot}$, which is only 32% of the total mass estimated from the unfiltered 1.2 mm continuum data in Peretto et al. (2006).

## Appendix D: Intensity Gradient

In addition to filament extraction algorithms (Section 4.1.1), the intensity gradient can also characterize features in density structures. Koch et al. (2012a, 2012b) show that the local intensity gradient is a measure for the resulting direction of motion in the MHD force equation, and hence it can be used to evaluate the relative importance of the magnetic field together with other constituents (such as local gravity; Section 4.1.2). One outcome of their polarization-intensity gradient technique is the magnetic field significance $\Sigma_B = \frac{\sin\Phi}{\cos\delta}$, which is a ratio that describes the relative importance between the magnetic field and gravity solely based on measurable angles. $\Phi$ is the angle between the intensity gradient and local gravity (Section 4.1.2), and $\delta$ is the angle between the intensity gradient and magnetic field. This method has been applied to a 50-source sample of star-forming regions to quantify the importance of the magnetic field in these sources (Koch et al. 2014). An additional method to characterize the magnetic field and density morphologies is histograms of relative orientations (Soler et al. 2013), which have been used for a variety of systems, from larger-scale to filamentary systems (e.g., Planck Collaboration et al. 2016; Kwon et al. 2022), finding that the relative orientation (which is identical to the angle $\delta$) can be either parallel, perpendicular, or transitioning, depending on the local density.

The local intensity gradients in our 4" pixel 850 µm continuum image are derived. We note that even though the oversampled pixels may be partially correlated, the derived orientations of the intensity gradients remain unaffected. This is

because of the isotropic beam pattern. The gradients are calculated by fitting at each grid point over 3 × 3 pixels in the continuum image (Goodman et al. 1993) with

$$I(x, y) = c_x x + c_y y + c_0, \qquad (D1)$$

where $x$ and $y$ are the coordinates of each pixel, and $c_0$, $c_x$, and $c_y$ are free parameters representing the first-order expansion of the intensity field.

The resulting intensity gradients are shown as arrows in Figure 26, with a median orientation uncertainty of 6°.3. Overall, the estimated intensity gradients are pointing toward both the massive clumps and filament ridges. This suggests that the local mass is pulled not only by the filaments themselves but also by global gravity toward the cloud's main mass centers. As a result, the intensity gradients often show offset angles of ~30°–60° to the filament ridges, instead of simply being parallel or perpendicular. The histogram of intensity gradients (right panel in Figure 7) displays a less pronounced but still visible broad peak around 30°–60°. This is a similar prevailing orientation as seen for the other parameters (filament, magnetic field, and gravity) but flatter with a very broad shoulder.

This is possibly because our resolved physical scale captures both filamentary structures and clumps, which are both delineated by closed intensity contours. These closed contours typically lead to a broad distribution in intensity gradient orientations with a broad peak that is reflecting a prevailing shape. Figure 27 shows histograms of the relative orientations between the intensity gradient $\nabla I$ and other spatial parameters, with the corresponding visualization maps in Figure 28. All the plots reveal a trend that $\nabla I$ is more aligned with filaments, magnetic field, and local gravity toward high-intensity regions. This trend is most obvious with local gravity, suggesting that the density structures are shaped by an increasingly dominating gravitational force. Figure 29 presents how the histogram shape changes with local intensity for the relative alignments between the intensity gradient and the other spatial parameters. $\nabla I$ is preferentially aligned with all parameters.
<figure><img src="image_28.png" /><figcaption>**Figure 26.** Filaments (green lines) and intensity gradients (magenta arrows) overlaid on 850 µm continuum toward (a) 2264C and (b) 2264D. The intensity gradients are calculated per pixel but only shown per 3 × 3 pixels and with uniform length for simplicity.</figcaption></figure>

<figure><img src="image_29.png" /><figcaption>**Figure 27.** Histograms of relative orientations between $\nabla I$ and other spatial parameters over all samples. The samples are separated into low- and high-intensity regions with an intensity threshold of 2.0 mJy arcsec$^{-2}$. The *p*-value from Rayleigh's test indicates the similarity between an observed and a uniform distribution, and *p* < 0.05 favors a nonuniform distribution. PRS $Z_x$ indicates whether the relative orientations tend to be parallel ($Z_x > 0$) or perpendicular ($Z_x < 0$).</figcaption></figure>
<figure><img src="image_30.png" /><figcaption>**Figure 28.** Maps of relative orientations between intensity gradient ($\nabla I$), filament ($F$), local magnetic field ($B$), and local gravity ($G$) orientations toward NGC 2264.</figcaption></figure>
<figure><img src="image_31.png" /><figcaption>**Figure 29.** Histogram shape parameter $\xi$ vs. local intensity for $\nabla I$ and other parameters toward 2264C (red) and D (blue). All pairs tend to be aligned and show a tendency of increasing alignment as the local intensity increases, except for $\nabla I$ vs. B in NGC 2264D.</figcaption></figure>
## Appendix E
### Histograms of Relative Orientations in 2264C and 2264D

Figures 30 and 31 present the histograms of relative orientations between all spatial parameters selected in 2264C

and 2264D separately. For most of the pairs, the overall trends in 2264C and 2264D are not significantly different, and the only pair showing an obvious difference is *G* versus *B*, which is further discussed in Sections 4.3.1 and 5.1.2.

<figure><img src="image_32.png" /><figcaption>**Figure 30.** Histograms of relative orientations between all spatial parameters selected in 2264C. The samples are separated into low-density and high-intensity regions with an intensity threshold of 2.0 mJy arcsec⁻². The *p*-value from Rayleigh’s test indicates the similarity between an observed and a uniform distribution, and *p* < 0.05 favors a nonuniform distribution. PRS *Z*<sub>x</sub> indicates whether the relative orientations tend to be parallel (*Z*<sub>x</sub> > 0) or perpendicular (*Z*<sub>x</sub> < 0).</figcaption></figure>
<figure><img src="image_33.png" /><figcaption>Figure 31. Identical to Figure 30 but for 2264D.</figcaption></figure>
## Appendix F

### Magnetic Field Strength Estimates Using Different Molecular Lines

We estimate magnetic field strengths, mass-to-flux ratios, and energy ratios for all structures identified in NGC 2264 with different gas tracers. The structures identified from the dendrogram algorithm are labeled in Figure 32. They are indicated by three symbols: trunk (N, S<sub>1</sub>, or S<sub>2</sub>) + branch (a

or b) + leaf (1, 2, or 3). If a level is not identified, it is omitted.

In addition to the ¹³CO and C¹⁸O data, we further adopt the N₂H⁺ data in Peretto et al. (2006) toward 2264D. The estimates using ¹³CO, C¹⁸O, and N₂H⁺ are listed in Tables 2–4, respectively. In general, the estimated magnetic field strengths increase by 10%–50% if ¹³CO is used instead of C¹⁸O, and they decrease by 10%–70% if N₂H⁺ is used instead of C¹⁸O.

<figure><img src="image_34.png" /><figcaption>**Figure 32.** Identified clumps in NGC 2264. Contours are the boundaries of the clumps identified using the dendrogram. The green, cyan, and red labels indicate whether the clumps correspond to trunks, branches, or leaves.</figcaption></figure>
<table><caption>**Table 2**
Magnetic Field Strengths and Mass-to-flux Ratios in NGC 2264 Using <sup>13</sup>CO</caption><thead><tr><th>Leaves</th><th>n<sub>H2</sub><br>(cm<sup>-3</sup>)</th><th>σ<sub>v,NT</sub><br>(km s<sup>-1</sup>)</th><th>δφ<br>(deg)</th><th>B<sub>pos</sub><br>(mG)</th><th>λ</th><th>T/W</th><th>M/W</th></tr></thead><tbody><tr><td>S<sub>1</sub>a1</td><td>(3.7 ± 0.2) × 10<sup>5</sup></td><td>1.21 ± 0.43</td><td>34.0 ± 0.3</td><td>0.43 ± 0.16</td><td>1.5<sup>+0.3</sup><sub>-0.5</sub></td><td>0.6<sup>+0.4</sup><sub>-0.4</sub></td><td>0.2<sup>+0.1</sup><sub>-0.1</sub></td></tr><tr><td>S<sub>1</sub>a2</td><td>(6.2 ± 0.4) × 10<sup>5</sup></td><td>0.92 ± 0.28</td><td>9.5 ± 1.4</td><td>1.56 ± 0.54</td><td>0.2<sup>+0.7</sup><sub>-0.1</sub></td><td>1.4<sup>+0.7</sup><sub>-0.7</sub></td><td>7.0<sup>+4.3</sup><sub>-4.3</sub></td></tr><tr><td>S<sub>2</sub>1</td><td>(5.6 ± 0.3) × 10<sup>4</sup></td><td>0.95 ± 0.26</td><td>14.3 ± 2.1</td><td>0.32 ± 0.10</td><td>0.2<sup>+0.1</sup><sub>-0.1</sub></td><td>3.9<sup>+1.8</sup><sub>-1.8</sub></td><td>8.3<sup>+4.7</sup><sub>-4.6</sub></td></tr><tr><td>N1</td><td>(1.5 ± 0.1) × 10<sup>5</sup></td><td>0.90 ± 0.20</td><td>8.5 ± 0.4</td><td>0.82 ± 0.18</td><td>0.2<sup>+0.1</sup><sub>-0.1</sub></td><td>1.3<sup>+0.5</sup><sub>-0.5</sub></td><td>7.3<sup>+3.1</sup><sub>-3.0</sub></td></tr><tr><td>Na1</td><td>(2.2 ± 0.2) × 10<sup>5</sup></td><td>0.93 ± 0.46</td><td>16.4 ± 2.0</td><td>0.53 ± 0.27</td><td>0.5<sup>+0.1</sup><sub>-0.3</sub></td><td>1.8<sup>+1.3</sup><sub>-1.3</sub></td><td>2.8<sup>+2.3</sup><sub>-2.3</sub></td></tr><tr><td>Na2</td><td>(1.4 ± 0.1) × 10<sup>5</sup></td><td>1.01 ± 0.21</td><td>11.4 ± 0.7</td><td>0.67 ± 0.15</td><td>0.3<sup>+0.1</sup><sub>-0.1</sub></td><td>1.4<sup>+0.6</sup><sub>-0.6</sub></td><td>4.5<sup>+1.9</sup><sub>-1.9</sub></td></tr><tr><td>Na3</td><td>(1.6 ± 0.1) × 10<sup>5</sup></td><td>1.10 ± 0.36</td><td>9.3 ± 0.9</td><td>0.93 ± 0.32</td><td>0.2<sup>+0.1</sup><sub>-0.1</sub></td><td>3.4<sup>+1.9</sup><sub>-1.9</sub></td><td>17.1<sup>+10.3</sup><sub>-10.2</sub></td></tr><tr><td>Nb1</td><td>(8.7 ± 0.6) × 10<sup>4</sup></td><td>1.04 ± 0.08</td><td>28.0 ± 2.7</td><td>0.22 ± 0.03</td><td>0.4<sup>+0.1</sup><sub>-0.1</sub></td><td>3.4<sup>+0.5</sup><sub>-0.5</sub></td><td>1.9<sup>+0.5</sup><sub>-0.5</sub></td></tr><tr><td>S<sub>1</sub>a</td><td>(1.5 ± 0.1) × 10<sup>5</sup></td><td>1.09 ± 0.39</td><td>30.2 ± 0.3</td><td>0.28 ± 0.10</td><td>1.3<sup>+0.5</sup><sub>-0.5</sub></td><td>0.50<sup>+0.3</sup><sub>-0.3</sub></td><td>0.2<sup>+0.1</sup><sub>-0.1</sub></td></tr><tr><td>Na</td><td>(5.1 ± 0.4) × 10<sup>4</sup></td><td>1.01 ± 0.32</td><td>21.5 ± 0.4</td><td>0.21 ± 0.07</td><td>0.7<sup>+0.2</sup><sub>-0.2</sub></td><td>0.9<sup>+0.5</sup><sub>-0.5</sub></td><td>0.8<sup>+0.5</sup><sub>-0.5</sub></td></tr><tr><td>Nb</td><td>(2.5 ± 0.2) × 10<sup>4</sup></td><td>1.09 ± 0.12</td><td>23.0 ± 1.8</td><td>0.15 ± 0.02</td><td>0.3<sup>+0.1</sup><sub>-0.1</sub></td><td>3.9<sup>+0.8</sup><sub>-0.8</sub></td><td>3.1<sup>+0.8</sup><sub>-0.8</sub></td></tr><tr><td>S<sub>1</sub></td><td>(4.7 ± 0.3) × 10<sup>4</sup></td><td>1.06 ± 0.37</td><td>31.1 ± 0.3</td><td>0.14 ± 0.05</td><td>1.4<sup>+0.3</sup><sub>-0.6</sub></td><td>0.5<sup>+0.3</sup><sub>-0.3</sub></td><td>0.2<sup>+0.1</sup><sub>-0.1</sub></td></tr><tr><td>S<sub>2</sub></td><td>(2.4 ± 0.2) × 10<sup>4</sup></td><td>0.74 ± 0.29</td><td>31.4 ± 1.3</td><td>0.07 ± 0.03</td><td>0.7<sup>+0.3</sup><sub>-0.3</sub></td><td>1.8<sup>+1.1</sup><sub>-1.2</sub></td><td>0.8<sup>+0.5</sup><sub>-0.5</sub></td></tr><tr><td>N</td><td>(3.4 ± 0.2) × 10<sup>4</sup></td><td>1.00 ± 0.31</td><td>22.6 ± 0.4</td><td>0.16 ± 0.05</td><td>0.8<sup>+0.2</sup><sub>-0.3</sub></td><td>0.8<sup>+0.4</sup><sub>-0.4</sub></td><td>0.6<sup>+0.4</sup><sub>-0.4</sub></td></tr></tbody></table>

**Note.** Magnetic field strengths and mass-to-flux ratios derived from the DCF method. The uncertainties listed here are obtained from propagating the observational uncertainties through the corresponding equations using a Monte Carlo approach. Possible additional systematic uncertainties due to, e.g., the unknown dust opacity are not included.

<table><caption>**Table 3**
Magnetic Field Strengths and Mass-to-flux Ratios in NGC 2264 Using C<sup>18</sup>O</caption><thead><tr><th>Leaves</th><th>n<sub>H2</sub><br>(cm<sup>-3</sup>)</th><th>σ<sub>v,NT</sub><br>(km s<sup>-1</sup>)</th><th>δφ<br>(deg)</th><th>B<sub>pos</sub><br>(mG)</th><th>λ</th><th>T/W</th><th>M/W</th></tr></thead><tbody><tr><td>S<sub>1</sub>a1</td><td>(3.7 ± 0.3) × 10<sup>5</sup></td><td>0.97 ± 0.25</td><td>34.0 ± 0.4</td><td>0.35 ± 0.09</td><td>1.6<sup>+0.4</sup><sub>-0.5</sub></td><td>0.4<sup>+0.2</sup><sub>-0.2</sub></td><td>0.1<sup>+0.1</sup><sub>-0.1</sub></td></tr><tr><td>S<sub>1</sub>a2</td><td>(6.2 ± 0.4) × 10<sup>5</sup></td><td>0.69 ± 0.18</td><td>9.5 ± 1.4</td><td>1.16 ± 0.35</td><td>0.3<sup>+0.1</sup><sub>-0.1</sub></td><td>0.8<sup>+0.4</sup><sub>-0.3</sub></td><td>3.8<sup>+2.1</sup><sub>-2.1</sub></td></tr><tr><td>S<sub>2</sub>1</td><td>(5.6 ± 0.4) × 10<sup>4</sup></td><td>0.58 ± 0.09</td><td>14.3 ± 2.1</td><td>0.19 ± 0.04</td><td>0.3<sup>+0.1</sup><sub>-0.1</sub></td><td>1.5<sup>+0.4</sup><sub>-0.4</sub></td><td>2.9<sup>+1.2</sup><sub>-1.2</sub></td></tr><tr><td>N1</td><td>(1.5 ± 0.1) × 10<sup>5</sup></td><td>0.90 ± 0.07</td><td>8.5 ± 0.4</td><td>0.82 ± 0.08</td><td>0.2<sup>+0.1</sup><sub>-0.1</sub></td><td>1.2<sup>+0.2</sup><sub>-0.2</sub></td><td>7.0<sup>+1.3</sup><sub>-1.3</sub></td></tr><tr><td>Na1</td><td>(2.2 ± 0.2) × 10<sup>5</sup></td><td>0.92 ± 0.11</td><td>16.4 ± 2.0</td><td>0.55 ± 0.09</td><td>0.4<sup>+0.1</sup><sub>-0.1</sub></td><td>1.4<sup>+0.3</sup><sub>-0.3</sub></td><td>2.3<sup>+0.7</sup><sub>-0.7</sub></td></tr><tr><td>Na2</td><td>(1.4 ± 0.1) × 10<sup>5</sup></td><td>1.22 ± 0.15</td><td>11.4 ± 0.7</td><td>0.80 ± 0.11</td><td>0.2<sup>+0.1</sup><sub>-0.1</sub></td><td>2.0<sup>+0.5</sup><sub>-0.5</sub></td><td>6.4<sup>+1.7</sup><sub>-1.7</sub></td></tr><tr><td>Na3</td><td>(1.6 ± 0.1) × 10<sup>5</sup></td><td>0.78 ± 0.07</td><td>9.3 ± 0.9</td><td>0.67 ± 0.09</td><td>0.2<sup>+0.1</sup><sub>-0.1</sub></td><td>1.7<sup>+0.3</sup><sub>-0.3</sub></td><td>8.0<sup>+2.0</sup><sub>-2.0</sub></td></tr><tr><td>Nb1</td><td>(8.7 ± 0.6) × 10<sup>4</sup></td><td>0.70 ± 0.08</td><td>28.0 ± 2.7</td><td>0.15 ± 0.02</td><td>0.6<sup>+0.1</sup><sub>-0.1</sub></td><td>1.6<sup>+0.3</sup><sub>-0.3</sub></td><td>0.8<sup>+0.2</sup><sub>-0.2</sub></td></tr><tr><td>S<sub>1</sub>a</td><td>(1.5 ± 0.1) × 10<sup>5</sup></td><td>0.94 ± 0.24</td><td>30.2 ± 0.3</td><td>0.24 ± 0.06</td><td>1.4<sup>+0.4</sup><sub>-0.4</sub></td><td>0.3<sup>+0.2</sup><sub>-0.2</sub></td><td>0.2<sup>+0.1</sup><sub>-0.1</sub></td></tr><tr><td>Na</td><td>(5.1 ± 0.4) × 10<sup>4</sup></td><td>1.04 ± 0.23</td><td>21.5 ± 0.4</td><td>0.22 ± 0.05</td><td>0.6<sup>+0.1</sup><sub>-0.2</sub></td><td>0.9<sup>+0.4</sup><sub>-0.4</sub></td><td>0.8<sup>+0.4</sup><sub>-0.4</sub></td></tr><tr><td>Nb</td><td>(2.5 ± 0.2) × 10<sup>4</sup></td><td>0.72 ± 0.09</td><td>23.0 ± 1.8</td><td>0.10 ± 0.01</td><td>0.5<sup>+0.1</sup><sub>-0.1</sub></td><td>1.8<sup>+0.4</sup><sub>-0.4</sub></td><td>1.3<sup>+0.4</sup><sub>-0.4</sub></td></tr><tr><td>S<sub>1</sub></td><td>(4.7 ± 0.3) × 10<sup>4</sup></td><td>0.91 ± 0.25</td><td>31.1 ± 0.3</td><td>0.13 ± 0.04</td><td>1.3<sup>+0.5</sup><sub>-0.3</sub></td><td>0.4<sup>+0.2</sup><sub>-0.2</sub></td><td>0.2<sup>+0.1</sup><sub>-0.1</sub></td></tr><tr><td>S<sub>2</sub></td><td>(2.4 ± 0.2) × 10<sup>4</sup></td><td>0.57 ± 0.09</td><td>31.4 ± 1.3</td><td>0.06 ± 0.01</td><td>0.8<sup>+0.1</sup><sub>-0.2</sub></td><td>1.0<sup>+0.3</sup><sub>-0.3</sub></td><td>0.4<sup>+0.1</sup><sub>-0.1</sub></td></tr><tr><td>N</td><td>(3.4 ± 0.2) × 10<sup>4</sup></td><td>1.01 ± 0.23</td><td>22.6 ± 0.4</td><td>0.16 ± 0.04</td><td>0.7<sup>+0.2</sup><sub>-0.2</sub></td><td>0.7<sup>+0.3</sup><sub>-0.3</sub></td><td>0.6<sup>+0.2</sup><sub>-0.2</sub></td></tr></tbody></table>

**Note.** Identical to Table 2.
<table><caption>**Table 4**

Magnetic Field Strengths and Mass-to-flux Ratios in NGC 2264 Using N₂H⁺</caption><thead><tr><th>Leaves</th><th>n<sub>H2</sub> (cm<sup>-3</sup>)</th><th>σ<sub>v,NT</sub> (km s<sup>-1</sup>)</th><th>δφ (deg)</th><th>B<sub>pos</sub> (mG)</th><th>λ</th><th>T/W</th><th>M/W</th></tr></thead><tbody><tr><td>N1</td><td>(1.5 ± 0.2) × 10<sup>5</sup></td><td>0.76 ± 0.04</td><td>8.5 ± 0.4</td><td>0.68 ± 0.05</td><td>0.2<sup>+0.1</sup><sub>-0.1</sub></td><td>0.9<sup>+0.1</sup><sub>-0.1</sub></td><td>4.9<sup>+0.7</sup><sub>-0.7</sub></td></tr><tr><td>Na1</td><td>(2.2 ± 0.2) × 10<sup>5</sup></td><td>0.61 ± 0.12</td><td>16.4 ± 2.0</td><td>0.34 ± 0.08</td><td>0.6<sup>+0.1</sup><sub>-0.1</sub></td><td>0.7<sup>+0.2</sup><sub>-0.2</sub></td><td>1.0<sup>+0.4</sup><sub>-0.4</sub></td></tr><tr><td>Na2</td><td>(1.4 ± 0.1) × 10<sup>5</sup></td><td>1.06 ± 0.18</td><td>11.4 ± 0.7</td><td>0.70 ± 0.13</td><td>0.2<sup>+0.1</sup><sub>-0.1</sub></td><td>1.5<sup>+0.5</sup><sub>-0.5</sub></td><td>4.9<sup>+1.7</sup><sub>-1.7</sub></td></tr><tr><td>Na3</td><td>(1.6 ± 0.1) × 10<sup>5</sup></td><td>0.59 ± 0.04</td><td>9.3 ± 0.9</td><td>0.48 ± 0.06</td><td>0.3<sup>+0.1</sup><sub>-0.1</sub></td><td>1.0<sup>+0.1</sup><sub>-0.1</sub></td><td>4.2<sup>+0.9</sup><sub>-0.9</sub></td></tr><tr><td>Nb1</td><td>(8.7 ± 0.6) × 10<sup>4</sup></td><td>0.41 ± 0.07</td><td>28.0 ± 2.7</td><td>0.09 ± 0.02</td><td>1.0<sup>+0.2</sup><sub>-0.2</sub></td><td>0.7<sup>+0.2</sup><sub>-0.2</sub></td><td>0.3<sup>+0.1</sup><sub>-0.1</sub></td></tr><tr><td>Na</td><td>(5.1 ± 0.4) × 10<sup>4</sup></td><td>0.86 ± 0.13</td><td>21.5 ± 0.4</td><td>0.18 ± 0.03</td><td>0.7<sup>+0.1</sup><sub>-0.1</sub></td><td>0.6<sup>+0.2</sup><sub>-0.2</sub></td><td>0.6<sup>+0.2</sup><sub>-0.2</sub></td></tr><tr><td>Nb</td><td>(2.5 ± 0.2) × 10<sup>4</sup></td><td>0.53 ± 0.15</td><td>23.0 ± 1.8</td><td>0.07 ± 0.02</td><td>0.7<sup>+0.2</sup><sub>-0.2</sub></td><td>1.1<sup>+0.5</sup><sub>-0.5</sub></td><td>0.8<sup>+0.4</sup><sub>-0.4</sub></td></tr><tr><td>N</td><td>(3.4 ± 0.2) × 10<sup>4</sup></td><td>0.87 ± 0.13</td><td>22.6 ± 0.4</td><td>0.14 ± 0.02</td><td>0.8<sup>+0.1</sup><sub>-0.1</sub></td><td>0.5<sup>+0.2</sup><sub>-0.2</sub></td><td>0.4<sup>+0.1</sup><sub>-0.1</sub></td></tr></tbody></table>

Note. Identical to Table 2.

## ORCID iDs

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

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

Seamus D. Clarke https://orcid.org/0000-0001-9751-4603

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

Nicolas Peretto https://orcid.org/0000-0002-6893-602X

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

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

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

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

Doris Arzoumanian https://orcid.org/0000-0002-1959-7201

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

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

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

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

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

Valentin J. M. Le Gouellec https://orcid.org/0000-0002-5714-799X

Hong-Li Liu https://orcid.org/0000-0003-3343-9645

Lapo Fanciullo https://orcid.org/0000-0001-9930-9240

Jihye Hwang https://orcid.org/0000-0001-7866-2686

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

Frédérick Poidevin https://orcid.org/0000-0002-5391-5568

Mehrnoosh Tahani https://orcid.org/0000-0001-8749-1436

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

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

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

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

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

Felix Priestley https://orcid.org/0000-0002-5858-6265

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

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

David Berry https://orcid.org/0000-0001-6524-2447

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

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

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

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

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

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

Tetsuo Hasegawa https://orcid.org/0000-0003-1853-0184

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

Tyler L. Bourke https://orcid.org/0000-0001-7491-0048

Do-Young Byun https://orcid.org/0000-0003-1157-4109

Zhiwei Chen https://orcid.org/0000-0003-0849-0692

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

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

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

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

Sophia Dai https://orcid.org/0000-0002-7928-416X

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

Pham Ngoc Diep https://orcid.org/0000-0002-2808-0888

Yasuo Doi https://orcid.org/0000-0001-8746-6548

Yan Duan https://orcid.org/0000-0003-3758-7426

Hao-Yuan Duan https://orcid.org/0000-0002-7022-4742

David Eden https://orcid.org/0000-0002-5881-3229

Laura M. Fissel https://orcid.org/0000-0002-4666-609X

Erica Franzmann https://orcid.org/0000-0003-2142-0357

Per Friberg https://orcid.org/0000-0002-8010-8454

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

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

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

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

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

Ilseung Han https://orcid.org/0000-0002-9143-1433

Saeko Hayashi https://orcid.org/0000-0001-5026-490X

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

Tsuyoshi Inoue https://orcid.org/0000-0002-7935-8771

Kazunari Iwasaki https://orcid.org/0000-0002-9892-1881

Il-Gyo Jeong https://orcid.org/0000-0002-5492-6832

Vera Könyves https://orcid.org/0000-0002-3746-1498

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

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

Janik Karoly https://orcid.org/0000-0001-5996-3600

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

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

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

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

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

Shinyoung Kim https://orcid.org/0000-0001-9333-5608

Jongsoo Kim https://orcid.org/0000-0002-1229-0426

Gwanjeong Kim https://orcid.org/0000-0003-2011-8172
Florian Kirschslager https://orcid.org/0000-0002-3036-0184

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

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

Takayoshi Kusune https://orcid.org/0000-0002-9218-9319

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

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

Chi-Yan Law https://orcid.org/0000-0003-1964-970X

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

Hyeseung Lee https://orcid.org/0000-0003-3465-3213

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

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

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

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

Sheng-Jun Lin https://orcid.org/0000-0002-6868-4483

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

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

Xing Lu https://orcid.org/0000-0003-2619-9305

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

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

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

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

Tetsuya Nagata https://orcid.org/0000-0001-9264-9015

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

Nguyen Bich Ngoc https://orcid.org/0000-0002-5913-5554

Geumsook Park https://orcid.org/0000-0001-8467-3736

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

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

Jonathan Rawlings https://orcid.org/0000-0001-5560-1303

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

Andrew Rigby https://orcid.org/0000-0002-3351-2200

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

Giorgio Savini https://orcid.org/0000-0003-4449-9416

Ekta Sharma https://orcid.org/0000-0002-4541-0607

Yoshito Shimajiri https://orcid.org/0000-0001-9368-3143

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

Xindi Tang https://orcid.org/0000-0002-4154-4309

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

Le Ngoc Tram https://orcid.org/0000-0002-6488-8227

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

Serena Viti https://orcid.org/0000-0001-8504-8844

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

Anthony Whitworth https://orcid.org/0000-0002-1178-5486

Jintai Wu https://orcid.org/0000-0001-7276-3590

Jinjin Xie https://orcid.org/0000-0002-2738-146X

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

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

Hyeong-Sik Yun https://orcid.org/0000-0001-6842-1555

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

Yapeng Zhang https://orcid.org/0000-0002-5102-2096

Jianjun Zhou https://orcid.org/0000-0002-3536-818X

Ilse de Looze https://orcid.org/0000-0001-9419-6355

Philippe André https://orcid.org/0000-0002-3413-2293

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

Sam Falle https://orcid.org/0000-0002-9829-0426

Jean-François Robitaille https://orcid.org/0000-0001-5079-8573

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

## References

* Allen, D. A. 1972, ApJ, 172, L55

* 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. (Tuscon, 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

* Arzoumanian, D., Furuya, R. S., Hasegawa, T., et al. 2021, A&A, 647, A78

* Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123

* 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. A. Adamson et al. (San Francisco, CA: ASP), 71

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

* Buckle, J. V., Richer, J. S., & Davis, C. J. 2012, MNRAS, 423, 1127

* Buckner, A. S. M., Khorrami, Z., González, M., et al. 2020, A&A, 636, A80

* Burkhart, B., Collins, D. C., & Lazarian, A. 2015, ApJ, 808, 48

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

* Cavanagh, B. 2008, StaUN, 260

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

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

* Chung, E. J., Lee, C. W., Kim, S., et al. 2019, ApJ, 877, 114

* Chung, E. J., Lee, C. W., Kwon, W., et al. 2022, AJ, 164, 175

* Clarke, S. D., Whitworth, A. P., Spowage, R. L., et al. 2018, MNRAS, 479, 1722

* Crutcher, R. M., Nutter, D. J., Ward-Thompson, D., & Kirk, J. M. 2004, ApJ, 600, 279

* Cunningham, N., Lumsden, S. L., Cyganowski, C. J., Maud, L. T., & Purcell, C. 2016, MNRAS, 458, 1742

* 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

* Dobashi, K., Matsumoto, T., Shimoikura, T., et al. 2014, ApJ, 797, 58

* Evans, N. J., II, Dunham, M. M., Jørgensen, J. K., et al. 2009, ApJS, 181, 321

* Federrath, C. 2016, MNRAS, 457, 375

* Fiege, J. D., & Pudritz, R. E. 2000, MNRAS, 311, 105

* Frau, P., Girart, J. M., Alves, F. O., et al. 2015, A&A, 574, L6

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

* Ganguly, S., Walch, S., Seifried, D., Clarke, S. D., & Weis, M. 2023, MNRAS, 525, 721

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

* Gómez, G. C., Vázquez-Semadeni, E., & Zamora-Avilés, M. 2018, MNRAS, 480, 2939

* Goodman, A. A., Benson, P. J., Fuller, G. A., & Myers, P. C. 1993, ApJ, 406, 528

* Goodman, A. A., Rosolowsky, E. W., Borkin, M. A., et al. 2009, Natur, 457, 63

* Henshaw, J. D., Ginsburg, A., Haworth, T. J., et al. 2019, MNRAS, 485, 2457

* Henshaw, J. D., Longmore, S. N., Kruijssen, J. M. D., et al. 2016a, MNRAS, 457, 2675

* Henshaw, J. D., Longmore, S. N., Kruijssen, J. M. D., et al., 2016b, SCOUSE: Semi-automated multi-COmponent Universal Spectral-line fitting Engine, Astrophysics Source Code Library, ascl:1601.003

* Heyer, M., Gong, H., Ostriker, E., & Brunt, C. 2008, ApJ, 680, 420

* Hildebrand, R. H. 1983, QJRAS, 24, 267

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

* Hwang, J., Kim, J., Pattle, K., et al. 2021, ApJ, 913, 85

* Hwang, J., Kim, J., Pattle, K., et al. 2022, ApJ, 941, 51

* Ibáñez-Méjía, J. C., Mac Low, M.-M., & Klessen, R. S. 2022, ApJ, 925, 196

* Inoue, T., Hennebelle, P., Fukui, Y., et al. 2018, PASJ, 70, S53

* Inoue, T., & Inutsuka, S.-i. 2012, ApJ, 759, 35

* Inutsuka, S.-i., Inoue, T., Iwasaki, K., & Hosokawa, T. 2015, A&A, 580, A49
Inutsuka, S.-I., Inoue, T., Iwasaki, K., Hosokawa, T., & Kobayashi, M. I. N. 2016, in IAU Symp. 315, From Interstellar Clouds to Star-forming Galaxies: Universal Processes?, ed. P. Jablonka, P. André, & F. van der Tak (Cambridge: Cambridge Univ. Press), 61

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

Jow, D. L., Hill, R., Scott, D., et al. 2018, MNRAS, 474, 1018

Koch, P. M., Tang, Y.-W., & Ho, P. T. P. 2012a, ApJ, 747, 79

Koch, P. M., Tang, Y.-W., & Ho, P. T. P. 2012b, ApJ, 747, 80

Koch, P. M., Tang, Y.-W., & Ho, P. T. P. 2013, ApJ, 775, 77

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

Koch, P. M., Tang, Y.-W., Ho, P. T. P., et al. 2022, ApJ, 940, 89

Kumar, M. S. N., Palmeirim, P., Arzoumanian, D., & Inutsuka, S. I. 2020, A&A, 642, A87

Kwon, W., Pattle, K., Sadavoy, S., et al. 2022, ApJ, 926, 163

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., & Houde, M. 2008, ApJ, 677, 1151

Liu, H.-L., Tej, A., Liu, T., et al. 2023, MNRAS, 522, 3719

Liu, J., Qiu, K., Berry, D., et al. 2019, ApJ, 877, 43

Liu, J., Zhang, Q., Koch, P. M., et al. 2023, ApJ, 945, 160

Mairs, S., Dempsey, J. T., Bell, G. S., et al. 2021, AJ, 162, 191

Mairs, S., Johnstone, D., Kirk, H., et al. 2016, MNRAS, 461, 4022

Margulis, M., Lada, C. J., & Snell, R. L. 1988, ApJ, 333, 316

Maury, A. J., André, P., & Li, Z. Y. 2009, A&A, 499, 175

McKee, C. F., & Ostriker, E. C. 2007, ARA&A, 45, 565

Mestel, L., & Spitzer, L. 1956, MNRAS, 116, 503

Motte, F., Bontemps, S., & Louvet, F. 2018, ARA&A, 56, 41

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

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

Nakamura, F., Sugitani, K., Tanaka, T., et al. 2014, ApJL, 791, L23

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

Naranjo-Romero, R., Vázquez-Semadeni, E., & Loughnane, R. M. 2022, MNRAS, 512, 4715

Nony, T., Robitaille, J. F., Motte, F., et al. 2021, A&A, 645, A94

Oliphant, T. 2006, NumPy: A Guide to NumPy (USA: Trelgol Publishing), https://web.mit.edu/dvp/Public/numpybook.pdf

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

Padoan, P., Juvela, M., & Goodman, A. A. 2001, ApJ, 553, 227

Palmeirim, P., André, P., Kirk, J., et al. 2013, A&A, 550, A38

Pattle, K., Lai, S.-P., Di Francesco, J., et al. 2021, ApJ, 907, 88

Pattle, K., Lai, S.-P., Hasegawa, T., et al. 2019, ApJ, 880, 27

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

Peretto, N., André, P., & Belloche, A. 2006, A&A, 445, 979

Peretto, N., Fuller, G. A., Duarte-Cabral, A., et al. 2013, A&A, 555, A112

Pillai, T. G. S., Clemens, D. P., Reissl, S., et al. 2020, NatAs, 4, 1195

Pineda, J. E., Arzoumanian, D., André, P., et al. 2023, in ASP Conf. Ser. 534, Protostars and Planets VII, ed. S.-i. Inutsuka (San Francisco, CA: ASP), 233

Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2015, A&A, 576, A105

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

Planck Collaboration, Aghanim, N., Akrami, N., et al. 2020, A&A, 641, A12

Priestley, F. D., & Whitworth, A. P. 2022, MNRAS, 509, 1494

Rapson, V. A., Pipher, J. L., Gutermuth, R. A., et al. 2014, ApJ, 794, 124

Robitaille, T. 2019, APLpy v2.0: The Astronomical Plotting Library in Python, v2.0, Zenodo, doi:10.5281/zenodo.2567476

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

Sadavoy, S. I., Di Francesco, J., Johnstone, D., et al. 2013, ApJ, 767, 126

Schneider, N., Csengeri, T., Hennemann, M., et al. 2012, A&A, 540, L11

Schneider, S., & Elmegreen, B. G. 1979, ApJS, 41, 87

Seifried, D., Walch, S., Weis, M., et al. 2020, MNRAS, 497, 4196

Soler, J. D. 2019, A&A, 629, A96

Soler, J. D., & Hennebelle, P. 2017, A&A, 607, A2

Soler, J. D., Hennebelle, P., Martin, P. G., et al. 2013, ApJ, 774, 128

Sousbie, T. 2011, MNRAS, 414, 350

Storm, S., Mundy, L. G., Lee, K. I., et al. 2016, ApJ, 830, 127

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

Sung, H., Stauffer, J. R., & Bessell, M. S. 2009, AJ, 138, 1116

Suri, S., Sánchez-Monge, Á., Schilke, P., et al. 2019, A&A, 623, A142

Tahani, M., Glover, J., Lupypciw, W., et al. 2022a, A&A, 660, L7

Tahani, M., Lupypciw, W., Glover, J., et al. 2022b, A&A, 660, A97

Tahani, M., Plume, R., Brown, J. C., Soler, J. D., & Kainulainen, J. 2019, A&A, 632, A68

Tobin, J. J., Hartmann, L., Fűrész, G., Hsu, W.-H., & Mateo, M. 2015, AJ, 149, 119

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

Vázquez-Semadeni, E., González-Samaniego, A., & Colín, P. 2017, MNRAS, 467, 1313

Vázquez-Semadeni, E., Palau, A., Ballesteros-Paredes, J., Gómez, G. C., & Zamora-Avilés, M. 2019, MNRAS, 490, 3061

Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, NatMe, 17, 261

Wang, J.-W., Koch, P. M., Galván-Madrid, R., et al. 2020a, ApJ, 905, 158

Wang, J.-W., Koch, P. M., Tang, Y.-W., et al. 2022, ApJ, 931, 115

Wang, J.-W., Lai, S.-P., Clemens, D. P., et al. 2020b, ApJ, 888, 13

Wang, J.-W., Lai, S.-P., Eswaraiah, C., et al. 2019, ApJ, 876, 42

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

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

Ward-Thompson, D., Zylka, R., Mezger, P. G., & Sievers, A. W. 2000, A&A, 355, 1122

Williams, G. M., Peretto, N., Avison, A., Duarte-Cabral, A., & Fuller, G. A. 2018, A&A, 613, A11

Zhou, J.-W., Liu, T., Evans, N. J., et al. 2022, MNRAS, 514, 6038

Zucker, C., Speagle, J. S., Schlafly, E. F., et al. 2020, A&A, 633, A51