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

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.


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;Arzoumanian et al. 2011;André et al. 2014;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. 2023a).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 ALMA observations suggest that self-similar hierarchical HFSs, where smallscale HFSs are embedded in the hub of large-scale HFSs, are possibly common in massive proto-clusters (Zhou et al. 2022).All these observational results point to HFSs being an essential step during the formation of massive stars and clusters.
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;Arzoumanian et al. 2021;Pillai et al. 2020;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 (Arzoumanian et al. 2021;Wang et al. 2020a;Chung et al. 2022).However, studies addressing the local variations of the magnetic field on subparsec-scale 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 polarizationintensity gradient method (Koch et al. 2012a(Koch et al. ,b, 2013)), the histogram of relative orientation 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. 2021Hwang et al. , 2022)), and relative orientations among gravity, magnetic fields, and density structures (Koch et al. 2018;Wang et al. 2020a;Koch et al. 2022;Wang et al. 2022;Liu et al. 2023b).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 +81 −42 pc (Zucker et al. 2020).Buckle et al. (2012) identified six parsec-scale filaments converging towards 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 towards 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.
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, both using 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.

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 November 2017 and February 2019.The twopointing mosaic observations targeted 2264C (south) and 2264D (north) with the reference positions (R.A., Dec.)=(6 h 41 m 11.97 s , 9 • 29 ′ 44.′′ 5) and (6 h 41 m 01.5 s , 9 • 35 ′ 39. ′′ 5).We performed 21 sets of 40-minute observations for each pointing under a weather condition with τ 225GHz 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 µm and 850 µm continuum bands.This paper focuses on the 850 µm data.
The POL-2 850 µm polarization data were reduced with pol2map1 in the smurf package2 (Berry et al. 2005;Chapin et al. 2013).The skyloop mode was invoked, which is a script that runs map-making 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 map-making 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 map-making 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 map-making uncertainties.The calculated polarization fraction P was debiased using the asymptotic estimator (Wardle & Kronberg 1974) as with a polarization uncertainty σ P calculated as where σ I , σ Q , and σ U are the uncertainties in the I, Q, and U Stokes parameters.The polarization position angle (P A) is and the corresponding uncertainty is estimated as The magnetic field orientations in this paper are inferred to be P A + 90 • .
In order to obtain a complete 13 CO (3-2) and C 18 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 −1 .These data cover an 8.5 ′ ×5.5 ′ area centered on (R.A., Dec.)=(6 h 41 m 10 s , 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 K and 0.2 K for NGC 2264 C 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.

POL-2 850 µ Continuum Map
Figure 1 shows the observed 850 µm continuum map in a logarithmic color scale, selected to display all reliable continuum detection (> 10σ I ).The structures within NGC 2264 span over a wide dynamical range, from 0.1 mJy arcsec −2 (S/N = 10) in the faint outskirts to 16 mJy arcsec −2 (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 panel and the filamentary extended structures in the bottom panel.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 hubfilament 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.

Magnetic Field from 850 µ Continuum Polarization
In order to ensure the robustness of polarization measurements, we select polarization detections based on the criteria I/σ I > 20 and P/σ P > 3. We further exclude detections with large uncertainties in polarization fraction (σ P > 5%).As a result, a total of 622 polarization measurements are selected across NGC 2264.For these selected data, the maximum and median uncertainty in the polarization position angle (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 counter-clockwise 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 NGC 2264 region (overlaid in Figure 1).The clear difference in orientations between large-and small-scale magnetic field suggests that the sub-parsec scale magnetic field, traced by POL-2, has been decoupled from the larger-scale field traced by Planck.

Velocity Dispersion traced by 13 O and C 18 lines
We extracted the velocity dispersion from our molecular line data adopting the scousePy package (Henshaw et al. 2016a(Henshaw et al. , 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 13 CO (3-2) and C 18 O (3-2) data pixel-by-pixel, where the pixel size is 14 ′′ and the channel width is 0.11 km s −1 .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 −1 .
The C 18 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 13 CO (3-2) and C 18 O (3-2) is around a constant of 3.7 over the entire cloud (see details in Appendix B in Figure 22), except for the few brightest pixels.This indicates that the 13 CO (3-2) line is also optically thin except for the brightest pixels.
In the 13 CO data cube, we commonly detected 2 to 3 velocity components across the source (for more information see Appendix B).In contrast to this, only one major component was detected in the C 18 O cube.Since both 13 CO and C 18 O are optically thin, this difference possibly results from velocity components in the low-density outskirts, seen only in 13 CO but not traced by C 18 O.Since we aim at estimating the magnetic field strength from combining polarization and molecular line data, we mainly adopt the C 18 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 13 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 18 O integrated intensity and the 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σ and T L = 3σ.Full details of the technique may be found in the code documentation3 .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 −2 (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 −2 , 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.24and 1.04±0.23 km s −1 , 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., IC5146: 0.1-0.3km s −1 (Wang et al. 2019), Ophiuchus C: 0.13 km s −1 (Liu et al. 2019), Serpens Main: 0.2-0.3km s −1 (Kwon et al. 2022)), but similar to other massive regions (e.g.,Orion A: 1.33±0.31km s −1 (Pattle et al. 2017) and Mon R2: 0.45-0.73km s −1 (Hwang et al. 2022)).

ANALYSIS
One fundamental question in the formation of a 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 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.

Filament Identification
In order to identify the filamentary structures within NGC 2264, we apply the DisP erSE algorithm (Sousbie 2011) on the 4 ′′ -pixel JCMT 850 µm Stokes I image.We use a persistent threshold of 0.8 mJy arcsec −2 (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 around 0.2-1.0pc and widths 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 nu-  The cyan contours mark 0.1, 0.2, 0.5, 1.0, 2.0, 5.0, and 10 mJy arcsec −2 .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 labelled.
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 • .

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,b).This is done by calculating the vector sum ( ⃗ F G,i ) of all gravitational forces from masses at all pixels j over a map4 (Wang et al. 2020a), acting at a pixel location i as where k is a factor accounting for the gravitational constant and conversion from emission to total column density.I i and I j are the intensity at the pixel position i and j, and n is the total number of pixels within the area of relevant gravitational influence.r i,j is the plane-ofsky projected distance between pixel i and j, and 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 −2 (10σ) on the continuum map, because the gravitational force originating from the diffuse and extended surrounding material is both small and rather symmetrical, thus assumed to be mostly cancelled 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 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 gravitational vector field and the filaments Fil 2 and 3 extending from the major filament appears more irregular.In 2264D, the gravitational field reveals four major converging points (C2 to 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.

1-Parameter Analysis -Globally Prevailing Orientations
We plot the histograms of the orientations of filament (F), B-field (B), and gravity (G) in Figure 7, both individually for 2264C and 2264D, and also for them combined.These histograms reveal the global tendencies of these parameters.None of the histograms is random.They either show a single-peaked or a bimodal distribution.
A common prevailing orientation, between about 20 • to 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 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 Section 5.1 and Section 5.2.

2-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 how clearly misaligned any two parameters are.Section 4.3.1 will first search for statistical trends, while the later Section 4.3.2 will look  at the local alignment or misalignment on the maps, with a focus on B and G.

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, indicated by the Rayleigh's test (p<0.05)5 .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 x > 0) or perpendicular-like (Z x < 0), with Z x defined as and the related uncertainties determined by where 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 (θ i = 0 • ) and perpendicular cases (θ i = 180 • ) simultaneously.
The PRS results suggest that most of the distributions are parallel-like, except for B vs. F in high-intensity regions.This is consistent with the 1-parameter analysis (Section 4.2) which finds very similar prevailing orientations for F, G, and B, which implies that, at least statistically, one should also expect closely aligned orientations among these parameters.However, the 2-parameter analysis reveals differences once low-and high-density regions are separated.This is very obvious for B vs F -with a clear change from rather aligned in the diffuse to rather orthogonal in the dense regions -and possibly more subtle for G vs 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 whether their origin can be localized and understood in a comprehensive scenario.This will be further explored in the next section and also in Section 5.2.
If NGC 2264 is further separated into its northern and southern complex, further nuances can be seen (Figure 9).We focus on G vs B as this forms the basis for our analysis in the next sections and also Section 5.2.While the diffuse regions in NGC 2264 C 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 gravityand 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.

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 towards 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 vs B (the first row 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 to 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 bottom panels in Figure 10 are 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 vs F reveal that gravity is often parallel to filaments, especially towards the gravitational converging centers.The ∆PA maps for B vs 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 vs F and B vs F reveals whether gravity and magnetic field tend to radially compress a filament (∆P A ∼90 • ) or pull/guide the gas flow along a filament (∆P A ∼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 which 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.

3-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 histogram of relative orientation (HRO) technique (Soler et al. 2013;Planck Collaboration et al. 2016).This technique probes how the shape of a histogram of relative orientations ∆PA  changes with local density.The shape is described by the parameter where A c and A e are the areas of a histogram around its center (0 < |∆P A| <22.• 5) and at its extreme (67.• 5< |∆P A| <90 • ).A positive ξ indicates relative orientations close to 0, while a negative ξ indicates perpendicular orientations.The uncertainty σ ξ is obtained as with the variances σ 2 Ac and σ 2 Ae of A c and A e as 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 ξ as a function of local intensity.In both 2264C and 2264D, ξ 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 H2 = 3 × 10 22 -10 23 cm −2 ) for different pairs, and they occur at higher intensities in 2264C than in 2264D.To examine whether the turn-over of ξ in these pairs is occurring for the same filaments, we plot the 2D kernel density distributions (KDE) of ∆PA for the G-B and B-F pairs in Figure 12, with a bandwidth of 2 • determined by Scott's rule.These pairs are selected from filaments with an intensity threshold of 2.0 mJy arcsec −2 to exclude diffuse, outskirt regions.The ∆PA for G-B and B-F show a positive correlation, seen with a Pearson correlation coefficient of 0.31±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⊥F and G⊥B, which will later be introduced as type I configuration in Section 5.2) and one in the lower left corner (B∥F and G∥B, hereafter 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.

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 δϕ and the line-of-sight non-thermal velocity dispersion σ v,N T , and assuming equipartition of kinetic and magnetic energy, the DCF method yields a planeof-sky magnetic field strength B pos as where ρ 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 2264 C and D (Figure 3), and calculate the angular dispersion using pairs from only adjacent segments as 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 κ of 0.0125 cm 2 /g at 850 µ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 A/π, 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 non-thermal velocity dispersion as 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±0.06and 0.22±0.05mG for 2264C and 2264D (Table 1).The mass-to-flux criticality λ obs , commonly used to evaluate the relative importance between magnetic field and gravity (Nakano & Nakamura 1978), is estimated as where µ=2.8 is the mean molecular weight per H 2 molecule, G is the gravitational constant, and N H2 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 λ becomes Since the large uncertainties in σ v,N T can cause unequal upper and lower uncertainties in λ, 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 flux conversion factor (Mairs et al. 2021), is included.The resulting upper and lower uncertainties (Table 1) are defined as the difference between the mean value, the 16th, and the 84th percentile.The estimated λ are 1.4 +0.4 −0.4 and 0.6 +0.1 −0.2 for cloud 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 (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 Ï determines the acceleration of the expansion or contraction of the spherical clump.The term is the total kinetic energy, where M is the total mass and σ tot is the total velocity dispersion in a clump, estimated as σ 2 tot = σ 2 v,N T +σ 2 thermal .The thermal velocity dispersion σ 2 thermal for mean free particles is calculated assuming a temperature of 2 K and a mass of mean free particle = 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 where v A = B/ √ 4πρ is the Alfvén velocity and ρ 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 δϕ and the line-of-sight nonthermal velocity dispersion σ v,N T .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/π)B pos (Crutcher et al. 2004).The term is the gravitational potential of a sphere with a uniform density ρ 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.2 ±0.1 :1, normalized to W).In contrast to this, the constituents are comparably strong within uncertainties in 2264D (T :M:W = 0.9 ±0.4 :0.8 ±0.4 :1).This could explain the difference in distributions of ∆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, 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.C 18 O on the other hand appears to trace the major highdensity 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 2 H + (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 also the polarization efficiency, determined by unknown dust properties and radiation field together with the SCUBA2/POL-2 large-scale filtering effect (Sadavoy et al. 2013), is 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 2 H + ) 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.

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 H2 = 10 22 -10 24 cm −2 ) and more compact (lengths = 0.2-1.0pc) than the typical filamentary networks discovered by Herschel, i.e., N H2 = 10 21 -10 22 cm −2 and lengths = 1-10 pc in IC5146 (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.

1-Parameter Statistics: Global Trends from Prevailing Orientations
The 1-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 sub-filament (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 • to 50 • ) indicates one dominating converging direction for gravity along a ∼40 • to 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 towards the filaments.
The B-field orientations are again around a prevailing orientation of ∼40 • to 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 2-parameter statistics discussed in the next section.Since the polarization measurements are independent from 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. 2019Tahani et al. , 2022a,b),b).The cartoon in Figure 13 summarizes the suggested scenario based on these global trends.

2-Parameter Statistics: Global Alignment Statistics
The relative orientation analysis (Figure 8) suggests that filament orientation and local gravity are all 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 towards smaller ∆PA.The only exception is the alignment between 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 turn-over 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 structures, but become perpendicular towards 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.1pc-scale regions with much higher densities (∼ 10 4 -10 5 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. 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).

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 ∆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 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 ∆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 panel) shows local gravity (G) converging from the lower-density surrounding 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 magneto-gravitational 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 both in lower-and higherdensity areas.Configuration II (right panel) 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 magneto-gravitational configuration shows G parallel to B both at larger distances and also close to F. This configuration is also seen both in lowerand higher-density areas.
In the following we discuss possible implications and origins of these two different magneto-gravitational 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 ω measure (Koch et al. 2018(Koch et al. , 2022)).With the angle ω between B and G, sin ω quantifies a fraction (between 0 and 1) of the magnetic field tension force that can work against a local gravitational pull.Small values in ω, i.e., close alignment between gravity and magnetic field, allow gravity to maximally accelerate inflowing material.Large values in ω, i.e., gravity and magnetic field close to orthogonal, indicate maximal obstruction by the magnetic field and hence a reduced acceleration.Consequently, the type I filament (left panel in Figure 14) will initially experience fast accretion from the outside (small ω; middle panels in Figure 14), hence short radial timescales towards the filament but then have a longer longitudinal timescale along the filament ridge (larger ω).The type II filament (right panel) will favour a short longitudinal timescale for any motion or collapse along the filament (small ω).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 a HFS is expected to be a valley.Radially, within a filament (r < 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 ( ∂U ∂r ) which is steeper (as a result of a higher density or a steeper density profile) than its longitudinal gravitational potential profile ( ∂U ∂ℓ ) (as a result of mass and distance to the major potential well).Thus, this configuration favors accretion from the ambient material.Once approaching the filament ridge, where the radial gravitational potential approaches its bottom and ∂U ∂r becomes zero, the longitudinal motion along the filament will start to dominate.This then causes a transition from small to larger ω 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 on 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 magneto-gravitational 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.

Star Formation in NGC 2264
YSOs are observed to form within filaments (André et al. 2014;Pineda et al. 2022).It is thus a question 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, 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 disk-less 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 RA (6 h 40 m 36 s ) to (6 h 41 m 37 s ) and Dec (9 The kernel density distribution uses a Gaussian kernel with a 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 free-fall dynamical timescale of 1.7×10 5 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.
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 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 (Cluster 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 of 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 does not likely 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 free-fall dynamical timescale (∼ 2 × 10 5 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., is greater or comparable to the longitudinal component ( ∂U ∂ℓ ) 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 towards the ridges.In contrast to this, the gravitational field in type II filaments is always dominated by the longitudinal component.Bottom row: sin ω maps for the two selected regions in the top row.
Cluster 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 a parsec-scale.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.

CONCLUSIONS
This study utilizes JCMT SCUBA-2/POL-2 850 µm continuum polarization observations to investigate the hubfilament system NGC 2264.An organized magnetic field morphology is uncovered with the following main outcome.
• 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 to 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.
• 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.
• 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 vs B-G: B⊥F and B⊥G (Type I filament) and B∥F and B∥G (Type II filament).
• 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 magneto-gravitational 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 also close to F.
• 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 time scale.We interpret Type I filaments as potentially being able to accumulate more nearby material, while Type II filaments are rather shaped by a dominating longitudinal pull from the larger-scale global gravitational field.
• 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 ±0.2 :0.2 ±0.1 :1 in 2264C and 0.9 ±0.4 :0.8 ±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.
• 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 transition disk YSOs are located at systematically larger distances from filaments.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 polarization is 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;Planck Collaboration et al. 2015;Jones et al. 2016;Wang et al. 2019;Pattle et al. 2019;Planck Collaboration et al. 2020).The I-P relation can be described by a power law, P = P I/I ∝ I −α , where P I is the polarized intensity.If the dust grains are not aligned with the magnetic field, P I is independent of Stokes I and α is 1. α smaller than 1 indicates that the polarization fraction results from magnetically aligned dust grains.
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 α using the model The Bayesian model allows us to directly model the Ricean distribution (Wardle & Kronberg 1974)    where P is the observed polarization fraction, P 0 is the real polarization fraction, σ 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 where σ Q is the dispersion in Stokes Q, assuming the same noise in Stokes Q and U. We re-select the polarization data based on the criteria σ I < 0.02 mJy arcsec −2 and I/σ I > 10, in order to have sufficient data points to model the noise-dominated regime.Additionally, the non-debiased 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.α is constrained to 0.29 ± 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.   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 linewidths.
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 2 times the fitted uncertainties in the velocity component identification.The black dashed line represents the best-fit T( 13 CO)/T(C 18 O) ratio of 3.7.Most points follow well the fitted line, 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 Figure 22, but with an additional criterion 2.5 < T( 13 CO)/T(C 18 O) < 5.5 to ensures that those 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 σ V (C 18 O)/σ V ( 13 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.

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σ.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 κ of 0.0125 cm 2 /g at 850 µm (Hildebrand 1983)) and a line mass (M 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⊙/pc at the outskirts, to over 100 M⊙/pc at the converging points.Notably, the main filament within 2264C has the highest line mass, reaching approximately 600 M ⊙ /pc.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 where c s is the sound speed and σ v,N T is the non-thermal velocity dispersion.Assuming a non-thermal velocity dispersion of 1.0 km/s and a gas temperature of 15 K, the M line,critical is about 400 M ⊙ /pc.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.Firstly, we integrate the products of M line and ∆L, where ∆L represents the length of each filament segment.According to this calculation, the total filament mass amounts to 757 M⊙, which corresponds to 26% of the total cloud mass of 2960 M⊙ (Peretto et al. 2006).Secondly, we calculate the total mass associated with the pixels belonging to filaments, selected if pixels are within the fitted 2σ width.To estimate the large-scale background, we use a Gaussian kernel with a σ 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⊙.This is equivalent to 19% of the total cloud mass.While the first method  (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.might be an upper limit (because filaments partially overlap near the converging points), the second method might be a lower limit (due to a possible over-subtraction 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⊙, which is only 32% of the total mass estimated from the unfiltered 1.2 mm continuum data in Peretto et al. (2006).

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,b) show that the local intensity gradient is a measure for the resulting direction of motion in the magneto-hydrodynamic (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 Σ B = sinΦ cosδ which is a ratio that describes the relative importance between magnetic field and gravity solely based on measurable angles.Φ is the angle between intensity gradient and local gravity (Section 4.1.2) and δ the angle between 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 are 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 δ) can be either parallel or 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 over-sampled 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 , (D5) 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 by not only 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 ∇I and other spatial parameters, with the corresponding visualization maps in Figure 28.All the plots reveal a trend that ∇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 intensity gradient and the other spatial parameters.∇I is preferentially aligned with all parameters.

Figure 1 .
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 −2 ) 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.

Figure 2 .
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 −2 .(Toppanel) High-intensity compact clumps are distributed along a filamentary structure in 2264C.The bright clumps in 2264D are more scattered.(Bottom panel) 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.

Figure 3 .
Figure 3. (a) C 18 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 −2 contour is shown in yellow and red, as in Figure 1.Examples of spectra are extracted at position a, b, c, and d.

Figure 4 .
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 of densities in Section 4.4.

Figure 5 .
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 −2 .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 labelled.

Figure 6 .
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 labelled as C1 to C5.

Figure 7 .
Figure 7. Histograms of filament, magnetic field, local gravity, and local intensity gradient orientations.The red and green 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 • to 50 • , both for 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.

Figure 8 .
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 −2 .The p-value from the Rayleigh's test indicates the similarity between an observed and a uniform distribution, and p< 0.05 favors a non-uniform distribution.PRS Zx indicates whether the relative orientations tend to be parallel (Zx > 0) or perpendicular (Zx < 0).

Figure 9 .
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 −2 .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 Zx indicates whether the relative orientations tend to be parallel (Zx > 0) or perpendicular (Zx < 0).

Figure 11 .
Figure 11.Histogram shape parameter ξ vs. local intensity for all parameter pairs toward 2264C (red) and D (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 and becoming more perpendicular in high-intensity regions.

Figure 12 .
Figure 12. 2D KDE of relative orientations between gravity(G)-magnetic field(B) and magnetic field(B)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.02and 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 filament (upper right corner) and the type II filament (lower left corner), as discussed in Section 5.2.
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 right bottom corner shows a zoom-in to the filament ridge in the black rectangle.The right panels are the expected histograms of the 1-parameter and 2-parameter statistics for this illustrated picture.All three parameters share a common prevailing orientation around ∼50 • , while 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.

Figure 14 .
Figure 14.Illustration of the role of gravity and magnetic field in two types of filament (type I, left panels; type II, right panels).The two maps in the top row are two representative regions as observed in NGC 2264 C 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 ( ∂U ∂r )is greater or comparable to the longitudinal component ( ∂U ∂ℓ ) 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 towards the ridges.In contrast to this, the gravitational field in type II filaments is always dominated by the longitudinal component.Bottom row: sin ω maps for the two selected regions in the top row.

Figure 15 .
Figure 15.YSO kernel density distribution 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 labelled as Cluster 1 to 10.

Figure 16 .
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 90 percentile of the nearest distance is 0.07 pc for Class 0/I and 0.45 pc for Class II/TD YSOs.

Figure 17 .
Figure 17.POL-2 polarization segments overlaid on 850 µm dust continuum toward NGC 2264.The yellow segments show polarization detections selected by I/σI > 20 and P/σ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.

Figure 18 .
Figure 18.Histogram of 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%.
B. 13 CO3-2) AND C 18 O (3-2)OLECULAR 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.Figure21shows example spectra of 13 CO (3-2) in NGC 2264 at different positions.In most of

Figure 19 .
Figure 19.850 µm total intensity I vs polarization fraction P .The green points are the non-debiased POL-2 polarization measurements, selected by I/σI >10 and σI <0.02 mJy arcsec −2 .The colored regions label the 68%, 95%, and 98% confidence regions (CR), predicted by the Bayesian model.The black line indicates the posterior mean.Most of the data points are within the 98% confidence region of our prediction.

Figure 20 .
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.

Figure 21 .
Figure 21. 13CO (3-2) integrated intensity map with spectra extracted at different pixels.Double or triple velocity components separated by only a few km s −1 are seen.

Figure 22 .
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.

Figure 23 .
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%.

Figure 24 .
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.
Figure 30 and Figure 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 vs B, which is further discussed in Section 4.3.1 and Section 5.1.2.

Figure 25 .
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⊙/pc in the outskirts to above 100 M⊙/pc in the converging points.The main filament in 2264C has the highest line mass of ∼600 M⊙/pc.

Figure 26 .
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.

Figure 30 .
Figure30.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 −2 .The p-value from the Rayleigh's test indicates the similarity between an observed and a uniform distribution, and p< 0.05 favors a non-uniform distribution.PRS Zx indicates whether the relative orientations tend to be parallel (Zx > 0) or perpendicular (Zx < 0).

Figure 32 .
Figure 32.Identified clumps in NGC 2264.Contours are the boundaries of the clumps identified using Dendrogram.The green, cyan, and red labels indicate whether the clumps correspond to trunks, branches, or leaves.

Table 1 .
Magnetic field strengths and mass-to-flux ratios in NGC 2264.