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 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^h41^m1197, 9°29'445) and (6^h41^m015, 9°35'395). We performed 21 sets of 40 minute observations for each pointing under a weather condition with τ_{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 141 (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 ⁹⁸ in the smurf package ⁹⁹ (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⁻² using an aperture flux conversion factor (FCF) of 2.79 Jy pW⁻¹ arcsec⁻² (including a factor of 1.35 for POL-2, equivalent to 630 Jy pW⁻¹ beam⁻¹) 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⁻² at the map center and gradually increases to ∼0.04 mJy arcsec⁻² near the edge of the mapped region. The Stokes I image has a higher rms noise of up to ∼0.07 mJy arcsec⁻² for pixels with large intensities (I > 6 mJy arcsec⁻²), 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

$$\begin{eqnarray}&&P=\displaystyle \frac{1}{I}\sqrt{\left({Q}^{2}+{U}^{2})-\tfrac{1}{2}({\sigma }_{Q}^{2}+{\sigma }_{U}^{2}\right)}\end{eqnarray} \tag{ 1 }$$

with a polarization uncertainty σ_P calculated as

$$\begin{eqnarray}&&{\sigma }_{P}=\sqrt{\tfrac{{Q}^{2}{\sigma }_{Q}^{2}+{U}^{2}{\sigma }_{U}^{2}}{({Q}^{2}+{U}^{2}){I}^{2}}+\tfrac{{\sigma }_{I}^{2}({Q}^{2}+{U}^{2})}{{I}^{4}}},\end{eqnarray} \tag{ 2 }$$

where σ_I , σ_Q , and σ_U are the uncertainties in the I, Q, and U Stokes parameters. The polarization position angle (PA) is

$$\begin{eqnarray}&&\mathrm{PA}=\displaystyle \frac{1}{2}{\tan }^{-1}\left(\displaystyle \frac{U}{Q}\right),\end{eqnarray} \tag{ 3 }$$

and the corresponding uncertainty is estimated as

$$\begin{eqnarray}&&\delta \mathrm{PA}=\displaystyle \frac{1}{2}\sqrt{\tfrac{({Q}^{2}\delta {U}^{2}+{U}^{2}\delta {Q}^{2})}{{\left({Q}^{2}+{U}^{2}\right)}^{2}}}.\end{eqnarray} \tag{ 4 }$$

The magnetic field orientations in this paper are inferred to be PA + 90°.

¹³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 85 × 55 map centered on (R.A., decl.) = (6^h41^m020, 9°34'432). The observations with a half-power beamwidth of 141 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 85 × 55 area centered on (R.A., decl.) = (6^h41^m10^s, 9°29'402). 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.

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

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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 uncertainties in the polarization PA are 954 and 45. 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 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.

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⁻¹.

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σ and T_L = 3σ. Full details of the technique may be found in the code documentation. ¹⁰⁰ 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).

Zoom In Zoom Out Reset image size

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).

Zoom In Zoom Out Reset image size

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 _G,i ) of all gravitational forces from masses at all pixels j over a map ¹⁰¹ (Wang et al. 2020a), acting at a pixel location i as

$$\begin{eqnarray}&&{{\boldsymbol{F}}}_{G,i}={{kI}}_{i}\displaystyle \sum _{j=1}^{n}\displaystyle \frac{{I}_{j}}{{r}_{i,j}^{2}}\hat{r},\end{eqnarray} \tag{ 5 }$$

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 positions i and j, and n is the total number of pixels within the area of relevant gravitational influence. r_i,j is the plane-of-sky projected distance between 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 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.

Zoom In Zoom Out Reset image size

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.

Zoom In Zoom Out Reset image size

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.

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'' ($\sim \sqrt{2}/2\,\times$ 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 Rayleigh's test (p < 0.05). ¹⁰² 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

$$\begin{eqnarray}&&{Z}_{x}=\displaystyle \frac{{\sum }_{i}^{n}\cos {\theta }_{i}}{\sqrt{n/2}},\end{eqnarray} \tag{ 6 }$$

and the related uncertainties determined by

$$\begin{eqnarray}&&{\sigma }_{{Z}_{x}}^{2}=\displaystyle \frac{2{\sum }_{i}^{n}{\left(\cos {\theta }_{i}\right)}^{2}-{\left({Z}_{x}\right)}^{2}}{n},\end{eqnarray} \tag{ 7 }$$

where θ_i = 2ΔPA_i , leading to θ_i ∈ [0, 180°] for our measured PA_i ∈ [0, 90°]. 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.

Zoom In Zoom Out Reset image size

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.

Zoom In Zoom Out Reset image size

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.

Zoom In Zoom Out Reset image size

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.

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

$$\begin{eqnarray}&&\xi =\displaystyle \frac{{A}_{c}-{A}_{e}}{{A}_{c}+{A}_{e}},\end{eqnarray} \tag{ 8 }$$

where A_c and A_e are the areas of a histogram around its center (0° < ∣ΔPA∣ < 225) and at its extreme (675 < ∣ΔPA∣ < 90°). A positive ξ indicates relative orientations close to 0, while a negative ξ indicates perpendicular orientations. The uncertainty σ_ξ is obtained as

$$\begin{eqnarray}&&{\sigma }_{\xi }=\displaystyle \frac{4({A}_{e}^{2}{\sigma }_{{A}_{c}}^{2}-{A}_{c}^{2}{\sigma }_{{A}_{e}}^{2})}{{\left({A}_{c}+{A}_{e}\right)}^{4}},\end{eqnarray} \tag{ 9 }$$

with the variances ${\sigma }_{{A}_{c}}^{2}$ and ${\sigma }_{{A}_{e}}^{2}$ of A_c and A_e as

$$\begin{eqnarray}&&{\sigma }_{{A}_{c,e}}^{2}={h}_{k}(1-{h}_{k}/{h}_{\mathrm{tot}}),\end{eqnarray} \tag{ 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 ξ as a function of local intensity. In both 2264C and 2264D, the ξ 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⁻² (${N}_{{{\rm{H}}}_{2}}$ = 3 × 10²²–10²³ cm⁻²) for different pairs, and they occur at higher intensities in 2264C than in 2264D. To examine whether the turnover of ξ in these pairs is occurring for the same filaments, we plot the 2D kernel density distributions (KDEs) 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⁻² to exclude diffuse outskirts regions. The ΔPAs 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 the type I configuration in Section 5.2) and one in the lower left corner (B∥F and G∥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.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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 nonthermal velocity dispersion σ_v,NT, and assuming equipartition of kinetic and magnetic energy, the DCF method yields a plane-of-sky magnetic field strength B_pos as

$$\begin{eqnarray}&&{B}_{\mathrm{pos}}=Q\,\sqrt{4\pi \rho }\displaystyle \frac{{\sigma }_{v,\mathrm{NT}}}{\delta \phi },\end{eqnarray} \tag{ 11 }$$

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 2264C and D (Figure 3) and calculate the angular dispersion using pairs from only adjacent segments as

$$\begin{eqnarray}&&\delta \phi =\displaystyle \frac{{\rm{\Sigma }}({\mathrm{PA}}_{i}-{\mathrm{PA}}_{j})}{{N}_{\mathrm{pairs}}},\end{eqnarray} \tag{ 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 κ of 0.0125 cm² 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 $\sqrt{A/\pi }$, where A is the clump area. From the observed C¹⁸O velocity dispersion (Section 3.3), we remove the thermal velocity component and derive the nonthermal velocity dispersion as

$$\begin{eqnarray}&&{\sigma }_{v,\mathrm{NT}}^{2}={\sigma }_{\mathrm{obs}}^{2}-\displaystyle \frac{{k}_{{\rm{B}}}{T}_{\mathrm{kin}}}{{m}_{{{\rm{C}}}^{18}{\rm{O}}}},\end{eqnarray} \tag{ 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 ± 0.06 and 0.22 ± 0.05 mG 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

$$\begin{eqnarray}&&{\lambda }_{\mathrm{obs}}=2\pi \sqrt{G}\displaystyle \frac{\mu {m}_{{\rm{H}}}{N}_{{{\rm{H}}}_{2}}}{{B}_{\mathrm{pos}}},\end{eqnarray} \tag{ 14 }$$

where μ = 2.8 is the mean molecular weight per H₂ molecule, G is the gravitational constant, and ${N}_{{{\rm{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 λ becomes

$$\begin{eqnarray}&&\lambda =\displaystyle \frac{{\lambda }_{\mathrm{obs}}}{3}.\end{eqnarray} \tag{ 15 }$$

Since the large uncertainties in σ_v,NT 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 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 λ are ${1.4}_{-0.4}^{+0.4}$ and ${0.6}_{-0.2}^{+0.1}$ for clouds C and D.

Table 1. Magnetic Field Strengths and Mass-to-flux Ratios in NGC 2264

Region	${n}_{{{\rm{H}}}_{2}}$	σv,NT	δ ϕ	Bpos	λ	T/W	M/W
	(cm−3)	(km s−1)	(deg)	(mG)
NGC 2264C	(1.5 ± 0.1) × 105	0.94 ± 0.24	30.2 ± 0.3	0.24 ± 0.06	${1.4}_{-0.4}^{+0.4}$	${0.3}_{-0.2}^{+0.2}$	${0.2}_{-0.1}^{+0.1}$
NGC 2264D	(5.1 ± 0.4) × 104	1.04 ± 0.23	21.5 ± 0.4	0.22 ± 0.05	${0.6}_{-0.2}^{+0.1}$	${0.9}_{-0.4}^{+0.4}$	${0.8}_{-0.4}^{+0.4}$

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.

Download table as:  ASCIITypeset image

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

$$\begin{eqnarray}&&\displaystyle \frac{1}{2}\ddot{I}=2({ \mathcal T }-{{ \mathcal T }}_{s})+{ \mathcal M }+{ \mathcal W }\end{eqnarray} \tag{ 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

$$\begin{eqnarray}&&{ \mathcal T }=\displaystyle \frac{3}{2}M{\sigma }_{\mathrm{tot}}^{2}\end{eqnarray} \tag{ 17 }$$

is the total kinetic energy, where M is the total mass and σ_tot is the total velocity dispersion in a clump, estimated as ${\sigma }_{\mathrm{tot}}^{2}$ = ${\sigma }_{v,\mathrm{NT}}^{2}$ + ${\sigma }_{\mathrm{thermal}}^{2}$. The thermal velocity dispersion ${\sigma }_{\mathrm{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 ${{ \mathcal 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

$$\begin{eqnarray}&&{ \mathcal M }=\displaystyle \frac{1}{2}{{Mv}}_{{\rm{A}}}^{2},\end{eqnarray} \tag{ 18 }$$

where ${v}_{{\rm{A}}}=B/\sqrt{4\pi \rho }$ 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,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/π)B_pos (Crutcher et al. 2004). The term

$$\begin{eqnarray}&&{ \mathcal W }=-\displaystyle \frac{3}{5}\displaystyle \frac{{{GM}}^{2}}{R}\end{eqnarray} \tag{ 19 }$$

is the gravitational potential of a sphere with a uniform density ρ and a radius R. The derived ${ \mathcal T }/{ \mathcal W }$ and ${ \mathcal M }/{ \mathcal 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 (${ \mathcal T }$: ${ \mathcal M }$: ${ \mathcal W }$ = 0.3^±0.2:0.2^±0.1:1, normalized to ${ \mathcal W }$). In contrast to this, the constituents are comparably strong within the uncertainties in 2264D (${ \mathcal T }$: ${ \mathcal M }$: ${ \mathcal 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, 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 ¹³CO. On the other hand, C¹⁸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 ¹³CO is systematically larger than the one from C¹⁸O by 30%–50%. Additionally, the velocity dispersion estimated from N₂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 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¹⁸O, ¹³CO, N₂H⁺) in Appendix F. Nevertheless, regardless of whether the ¹³CO or C¹⁸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.

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}_{{{\rm{H}}}_{2}}={10}^{22}$–10²⁴ cm⁻²) and more compact (length = 0.2–1.0 pc) than the typical filamentary networks discovered by Herschel, i.e., ${N}_{{{\rm{H}}}_{2}}={10}^{21}$–10²² cm⁻² 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.

Zoom In Zoom Out Reset image size

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 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 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 Δ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.

Zoom In Zoom Out Reset image size

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 ω 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 ω, 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 panels in Figure 14) will initially experience fast accretion from the outside (small ω; 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 ω). The type II filament (right panels) will favor 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 an 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 ($\tfrac{\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 ($\tfrac{\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 approaching the filament ridge, where the radial gravitational potential approaches its bottom and $\tfrac{\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 ω 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.

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. = 6^h40^m36^s–6^h41^m37^s and decl. = 9°22'23''–9°42'11''.

Zoom In Zoom Out Reset image size

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.

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.

Zoom In Zoom Out Reset image size

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.