Bimodal orientation distribution and head-tail asymmetry of a sample of filamentary molecular clouds

The morphology of molecular clouds is crucial for understanding their origin and evolution. In this work, we investigate the morphology of the filamentary molecular clouds (filaments for short) using a portion of the $^{12}\text{CO} (J=1-0)$ data from the Milky Way Imaging Scroll Painting (MWISP) project. The data cover an area spanning $104.75^\circ<l<150.25^\circ , \vert b\vert<5.25^\circ$ in Galactic coordinates, with $V_\text{LSR}$ ranging from $-95$ to 25 $\text{km s}^{-1}$. Our primary focus is on the orientation and morphological asymmetry of the filaments. To achieve this, we apply several criteria on the data to create a sample of filaments with well-defined straight shape, and we use elliptical fitting to obtain the orientation of each filament, with an estimated error of $\sim1.6^\circ$ for the orientation. We find that the filament orientation with respect to the Galactic plane exhibits a bimodal distribution, a double-Gaussian fitting of which has two centres located at $-38.1^\circ $ and $42.0^\circ $, with 1$\sigma$ of the two Gaussian functions being $35.4^\circ$ and $27.4^\circ$. We do not find significant correlation between the orientation and other parameters, including the Galactic coordinates, radial velocity, velocity width, and physical scale. A considerable fraction of filaments ($\gtrsim 40$ per cent) display head-tail asymmetry, which suggests that mass concentration tends to occur at one end of the filaments.


INTRODUCTION
Recent studies have demonstrated that molecular clouds display a variety of structural forms, one of the most prominent being filamentary molecular clouds (Arzoumanian et al. 2011;Palmeirim et al. 2013;Könyves et al. 2015).The results from the Gould Belt Survey with Herschel have revealed numerous filamentary structures in Galactic molecular clouds, and a connection between the filamentary structure and the formation of dense cloud cores has been made (Pilbratt et al. 2010;Jackson et al. 2010;Molinari et al. 2010;André et al. 2010).
Filaments have been observed at different scales from ≲1 pc (Hacar et al. 2013(Hacar et al. , 2018) ) to hundreds of pc (Li et al. 2013b;Ragan et al. 2014;Pineda et al. 2022).Various investigations using diverse tracers have significantly enhanced our understanding of filaments, such as extinction (Cambrésy 1999;Myers 2009), dust emission (Abergel et al. 1994), H I (McClure-Griffiths et al. 2006), and CO emission from diffuse molecular gas (Falgarone et al. 2001;Hily-Blant & Falgarone 2009).For instance, Wang et al. (2016) present a study with the dense gas tracers HCO + (3−2) and N 2 H + (3−2), which includes a comprehensive catalog and physical parameters of 54 large-scale filaments.They find that the filaments are widely distributed across the Galactic disk, and most of the filaments are associated with major spiral arms.Yuan et al. (2021) use 12 CO as the tracer.They classified the morphology of 18190 molecular clouds ⋆ E-mail: fjdu@pmo.ac.cn visually.Despite the lower fraction of filaments versus nonfilaments (∼11 per cent vs ∼89 per cent), filaments contribute ∼90 per cent to the total flux of the whole sample.In addition, filaments tend to have larger spatial scales than nonfilaments.
There are several possible explanations for the formation of filaments (Hacar et al. 2022), necessarily including some anisotropic factors such as large-scale gas motions sweeping up material (Burkert & Hartmann 2004), shock interactions induced by turbulence (Federrath 2016), magnetic fields (Soler et al. 2013), stellar feedback (Peretto et al. 2012), or Galactic shear (Duarte-Cabral & Dobbs 2017).The plurality of these physical processes may explain why filaments are ubiquitous in ISM.
Morphological analysis is crucial for understanding the formation of filamentary molecular clouds, because the shaping of them by large-scale processes of a galaxy may be reflected in their morphological characteristics.Some studies found that orientations of filaments tend to be either parallel or perpendicular to the directions of the local magnetic field (Li et al. 2013a;Heyer et al. 2020;Planck Collaboration et al. 2016a), while other studies also showed that the magnetic field is usually aligned in parallel with the filaments in the molecular clouds (McClure-Griffiths et al. 2006;Clark et al. 2014).The numerical simulation of Barreto-Mota et al. (2021) showed a tendency for magnetic fields to be aligned perpendicularly to dense filaments, and parallel to low density filaments, while Dobbs & Wurster (2021) found in their simulations that filaments are perpendicular to the magnetic field in sub-Alfvénic systems, whereas the opposite is true in super-Alfvénic models.Besides, the cold interstel-lar medium is observed to be organized in bubbles and filamentary structures (Pineda et al. 2022), and some investigations also indicate the expanding nature of these bubbles shapes the surrounding medium and possibly plays a role in the formation and evolution of interstellar filaments (Ntormousi et al. 2011;Inutsuka et al. 2015).Palmeirim et al. (2013) study the orientation of substructures within the dense, star-forming filament B211 in the Taurus molecular cloud.They find that the low-density striations are almost orthogonal to the B211 main filament, and very close to the local direction of the magnetic field.Soler et al. (2021) perform a statistical study of the filamentary structure orientation in the CO emission data from the MWISP project.They find that while most of the filamentary structures in CO do not show a preference to be parallel or perpendicular to the Galactic plane, a subset of them are indeed aligned with the Galactic plane.
In this work, we study the orientation and the head-tail asymmetry of filamentary molecular clouds.In Section 2, we introduce the data we use, and we describe the procedures we follow to construct a sample of filaments with well-defined straight morphology.In Section 3, we use the direction of the semi-major axis obtained from the elliptical fitting as the orientation of the filament, then we study the distribution of the orientation.In Section 4, we study the integrated intensity along the major axis of the filaments.In addition to skewness, we also define a new parameter p to measure the head-tail asymmetry of the filaments.We present a summary of our findings in Section 5.

DATA AND SAMPLE SELECTION
The MWISP project is a systematic survey of CO and its isotopologues on the Galactic plane with the 13.7 m single-dish telescope of the Purple Mountain Observatory.Three molecular spectral lines of 12 CO / 13 CO / C 18 O J = 1 − 0 are observed at the same time.
The typical system temperatures are ∼250 K for 12 CO and ∼140 K for 13 CO and C 18 O, respectively (Sun et al. 2018;Su et al. 2019).The data used in this work are in the form of a fits cube, where the first two dimensions represent the Galactic longitude and latitude, while the third dimension represents the local standard-of-rest velocity (V LSR ).The spatial resolution of the data is 50 ′′ with an approximately Nyquist-sampling pixel size of 30 ′′ , and the velocity resolution is 0.167 km s −1 .
In this work, the 18190 molecular clouds that we use are obtained from the 12 CO(1 − 0) spectral line data using the DBSCAN algorithm (Ester et al. 1996), with Galactic longitudes 104.75 • ≤ l ≤ 150.25 • , and Galactic latitudes |b| ⩽ 5.25 • (Yuan et al. 2021).DBSCAN has two parameters (ε and MinPts), where ε defines the neighborhood of voxels and MinPts defines core points, which are voxels whose numbers of neighboring voxels (including itself) are at least equal to MinPts.See Yan et al. (2021) for details about the parameters of the DBSCAN algorithm.The minimum cutoff on the PPV data cubes is set to 2σ (∼1 K).
To study the orientation of filaments, we construct a sample with a simple straight morphology from the 18190 clouds of Yuan et al. (2021), since only for such filaments a meaningful orientation can be defined.We first exclude those that have too large or too small angular scales; the former could be molecular cloud complexes, while the latter are limited by the angular resolution, thus an orientation is not well-defined for them.We also exclude the clouds with narrow velocity widths.If the velocity width is too narrow, the S/N of the image may not be high enough.Thus, we set 7 0.835 km s −1 ) and 80 × 80 × 40 (40 ′ × 40 ′ × 6.680 km s −1 ) as the minimum and maximum cloud sizes in voxel, respectively.
In the following two subsections, we describe the specific steps for selecting the samples through morphology.In Section 2.1, we use elliptical-fitting to exclude those with a small aspect ratio and define a 'preliminary' orientation of each cloud with the major axis.In Section 2.2, we further exclude those that are too curved with respect to this orientation.Fig. 1 shows integrated intensity maps of some clouds that are included in or excluded from our sample.

Elliptical fitting and aspect ratio selection
We use elliptical fitting method to retrieve the size and shape of the 18190 molecular clouds of the Yuan et al. (2021) sample.The elliptical fitting results are used for further filtering of the sample.
In this method, the integrated intensity of each pixel within a region is treated as the density of each point on a rigid plate, and the moment of inertia of this plate with respect to its centre of mass can be represented by an elliptical shape.
We first calculate the position of the centre of mass: x c = xρ(x, y)dxdy ρ(x, y)dxdy , y c = yρ(x, y)dxdy ρ(x, y)dxdy . (1) Then we calculate the moment of inertia tensor: with (5) The direction of the two eigenvectors of J corresponds to the direction of the major and minor axes of the two-dimensional intensity distribution under study.The two eigenvalues describe the length scales in these two directions.
We only include clouds with aspect ratio (the ratio between the lengths of the major and minor axes) greater than 10 in our sample, to ensure that they are visually 'slender'.The orientation of a filament is defined as the major axis direction of the elliptical fitting.

Removal of curved clouds
After the previous two steps of selection based on size and aspect ratio, we get a sample with 817 sources.Some of the selected clouds may still appear curved, for which an orientation is not well-defined.Therefore we take a further step to exclude these curved ones.
We make use of the boundary morphology of a cloud to decide whether it is curved or not (the minimum cut-off is set to 2σ ).We obtain the pixel coordinates of the boundary points from the integrated intensity map of the sample, then we calculate the distance (d) between each boundary point and the major axis of the elliptical fitting.Intuitively, if the sample is more curved, the variance of this distance d will be larger.Therefore, we calculate the variance of the distance over the boundary: where n is the number of boundary points of a cloud.Fig. 2 shows examples of straight and curved sample.Calculations of many other filaments also show that Var(d) of curved clouds are greater than those of straight clouds.To eliminate the influence of the filament's overall size, we define a dimensionless parameter where S is the number of pixels within a filament.We expect that a filament with a small value of s has a straight shape.Since we are unable to define a threshold of parameter s from 'first principles' to obtain a suitable sample, we try with three thresholds for it to obtain the three different samples.Then we compare the samples to select one suitable threshold to be used.To see the effect of using different thresholds, we scale all the filaments in each sample to the same spatial scale and co-add their intensities.The three overlay maps obtained from the three different thresholds are shown in panels (a-c) of Fig. 3.
We also randomly select 310 clouds and find no evidence of a bimodal orientation distribution after stacking them (panel (d) of Fig. 3), thereby confirming the necessity of the step to exclude curved clouds.This test also supports our view that the orientation of curved molecular clouds cannot be accurately described.
With the decrease in the threshold for parameter s, the number of sources in the sample also decreases, and we observe that the bimodal orientation distribution of the overlay becomes more pronounced.This confirms our expectation that a lower threshold yields a higher level of straightness and morphological simplicity.
As a compromise between the sample size (for statistical significance) and sample 'sharpness' (for morphological straightness), for subsequent studies, we use the sample comprising 522 filaments that correspond to s < 0.105.
In Appendix A, we also conduct an experiment to confirm that this step does not introduce bias to the sample selection.

Estimation of physical parameters of filament
We employ the Parallax-Based Distance Calculator V2 (Reid et al. 2016) to estimate the distances of filaments.This tool uses a source's Galactic coordinates and V LSR to assign it to a portion of a spiral arm based on its proximity to spiral arms seen in CO or H I surveys.As the distances to a large portion of the spiral arms are now known, we can obtain reliable distance estimations through this model.
We further calculate the physical lengths and widths of the filaments using the distances, angular lengths and angular widths.As shown in Fig. 4, we create a series of lines along the major axis, perpendicular to the orientation of the filament, and employ each line for width calculation.For example, we take the red line in Fig. 4 and identify the two closest points to the line.Subsequently, we compute the separation between these two points for each line parallel to the red line.Any separation that is less than or equal to 2 pixels is excluded.Finally, we calculate the average value of all the remaining separations to represent the average width of the filament.We apply the same method to calculate the average length of the filaments.
The velocity width of each filament is obtained through Gaussian fitting to the integrated line profile.The full-width at half maximum (FWHM = 2 √ ln 2σ ) of the fitting results are later used to investigate the possible correlation between velocity width and orientation in Section 3.2.
In Table B1 of Appendix B, we provide a list of relevant parameters, including the source name, orientation, coordinates, velocity, estimated distance, physical scale, velocity width, mass, maximum column density, and average column density for the 522 filaments in our sample.The red line is used as an example to illustrate the calculation process, its direction is perpendicular to the orientation of the filament.The pixel size here is 30 ′′ .

BIMODAL ORIENTATION DISTRIBUTION OF THE SAMPLE
We use the major axis direction of elliptical fitting to define the orientation of filaments in our final sample and study the distribution of the orientation.To estimate the error in the obtained orientation angle, we add Gaussian noise with a standard deviation of √ nσ ∆v (where the 'σ ' is the RMS noise of the MWISP data, '∆v' is velocity resolution of the data, and 'n' is the number of velocity channels (for calculating integrated intensity) of the filaments) to each pixel of the filament, then use elliptical fitting to calculate the orientation.This process is repeated 1000 times to obtain the standard deviation of these orientation data.This procedure is performed for each filament, from which the standard deviation of the orientation angle of each filament is calculated.The errors so calculated for each filament are different, and typically lies in the range 0.2 • → 7.1 • , and we report 1.6 • as the most probable value.Additionally, we investigate the correlation between the orientation and other physical parameters of the filaments.As an intuitive way to show the orientation distribution and an alternative to Fig. 5, we make an overlay map in Fig. 6, where the background, like in Fig. 3, is obtained using the method described in Section 2.2.The histogram shows the orientation distribution in polar coordinates.Both the overlay map and the polar histogram reveal a bimodal orientation distribution.

Bimodal orientation distribution
To see whether the relative orientation between two filaments depends on the their spatial separations, we define the following function where the symbol # means to take the number of elements of a set, ∧ means 'logical and', d i j is the (angular) distance between filament i and j, and ∠ i j is the relative angle between the two.Thus f (ℓ, θ )dℓdθ quantifies the fraction of filament pairs with angular separations between ℓ and ℓ + dℓ having relative orientation between θ and θ + dθ .We calculate f (ℓ, θ )dℓdθ as follows.We first calculate the relative orientations between each pair of filaments and obtain a total of 135981 (= 522×521/2, where 522 is the number of filaments in our sample) relative orientations.We then divide the pairs into intervals of 5 • based on their angular separations, and divide the relative orientations (0 • to 90 • ) into six intervals.Finally, we calculate the fraction of the relative orientations in each separation interval.Fig. 7 shows the fraction of different relative orientations for each angular separation interval.We observe no correlation between the relative orientation and the angular separation of the filaments.
The distribution of all the filaments in our final sample in plane of the sky is shown in Fig. 8.The background colour scale is a density map of filament count, generated using the scipy function gaussian_kde1 , where darker colours indicate higher density.The bimodal orientation distribution is already noticeable 'by eye' from this figure, and apparently the orientation does not appear to be correlated with spatial location.

Correlation of the orientation with the other parameters
As an attempt to understand the physical origin of the bimodal distribution, we look for any possible correlation between the orientation and other physical parameters of the filaments, such as Galactic longitude and latitude, radial velocity, physical scale, and velocity width (see Section 2.3).
The scatter plots of these parameters versus orientation angle are shown in Fig. 9, where in each panel the parameters R 1 and R 2 are the correlation coefficients for sources with negative or positive orientation angle, respectively.The maximum absolute value of these correlation coefficients never exceeds 0.15, thus we do not find any obvious correlation between orientation and other physical parameters, and the origin of the bimodal orientation distribution remains a puzzle at this moment.

HEAD-TAIL ASYMMETRY
For the filaments in our sample, we observe that some have a morphological asymmetry, that is, one end is 'thicker' than the other (See, e.g., panel (a) in Fig. 1).One possible explanation for the observed accumulation of matter at the end of filaments is the 'edge effect' (Bastien 1983), caused by the filament's self-gravity.This effect has been predicted by theory, and has been observationally found (Zernickel et al. 2013;Bhadari et al. 2020;Yuan et al. 2020).In this work, we use a large sample from the MWISP survey to statistically analyse the head-tail asymmetry of filamentary molecular clouds.

One-dimensional intensity distribution along the major axis and the number of peaks
We quantify this head-tail asymmetry using the integrated intensity along the major axis of the filaments.Fig. 10 shows one example for the one-dimensional integrated intensity profile.
The number of peaks in the curve represents the concentration of integrated intensity, which can influence the asymmetry of the curve.Therefore, we perform statistical analysis to study the headtail asymmetry in filaments with an equal number of peaks.
To determine the number of peaks in the curve, we perform multi-Gaussian fitting on the integrated profile and extract a series of coordinates for candidate peaks from the centres of the multi-Gaussian functions.The peak number of the curve is defined based on the parameters of the Gaussian function.We exclude candidate peaks that have peak values less than 3 √ nσ , where the 'σ ' is the RMS noise of the MWISP data, and 'n' is the number of integration steps, which includes two components: the number of pixels in the velocity dimension and the number of pixels in the direction perpendicular to the filament's major axis before the interpolation.Next, we calculate the separation between the two closest candidate peaks along the major axis.We use 0.85FWHM of the Gaussian fitting to a peak to represent its width and compare it with the separation of the two peaks.Velocity width (km s 1 ) If both widths of the two peaks are less than the separation, we consider that these two candidate peaks are two independent peaks.In our sample, 11.7 per cent of the filaments' one-dimensional integrated intensity profiles are single-peaked, 31.8 per cent are double-peaked, 30.7 per cent have three peaks, and the fraction of filaments with four or more peaks is ≈26 per cent (see Fig. 11).
Subsequently, we analyse the filaments based on different peak numbers.Since the limited number of samples with five or six peaks, we only construct subsamples for filaments containing four peaks or fewer.

Head-tail asymmetry of the integrated intensity
We use two approaches to quantify the asymmetry of the integrated intensity curves of the filaments.We first use skewness to quantify the degree of asymmetry in the one-dimensional integrated intensity.Skewness is a numerical measure of the degree of asymmetry of a probability density distribution curve, with positive and negative skewness representing left and right deviation of the data, defined as the third-order standardized moment of the random variable: If we consider the value of each point on the integrated intensity profile as a probability density, the skewness can be used to describe the head-tail asymmetry of the curve.Alternatively, we define a parameter p to quantify the degree of asymmetry in a curve.We calculate two moments of inertia of the cloud relative to two axes parallel to the minor axis and located at the two endpoints of the major axis: and the parameter p is defined as: Thus a higher value for p means the one dimensional distribution along the major axis is more asymmetric.To demonstrate the morphological difference of filaments with different p, we present 12 CO (J = 1 − 0) velocity-integrated intensity maps of representative filaments in increasing order of p in Appendix C. A visual inspection there indicates that filaments with p ≳ 0.35 may be reasonably considered to be 'obviously head-tail asymmetric'.After calculating the parameter p for all the filaments in the sample, we sort them in descending order of p.We then normalize the length and one-dimensional integrated intensity of all filaments and adjust their orientation so that the 'heavier' end is aligned.Using the serial number of filaments as the x-axis and the normalized length as the y-axis, we draw the intensity distribution and co-add the normalized intensities of all filaments to obtain the one-dimensional superimposed integrated intensity curve (Fig. 12).It can be seen that as the parameter p decreases, the intensity of the filaments gradually transitions from asymmetrical to symmetrical, and the overall intensity distribution along the major axis indeed tends to be one-sided.
The complementary cumulative distribution function (CCDF) of skewness and p are shown in panels (a) and (b) of Fig. 13.Panel (c) of that figure shows the correlation between skewness and p, which demonstrates that both can be used to characterize the asymmetric distribution.The CCDF plot shows that about 40 per cent of all the filaments in our sample have p > 0.35, so we may say that about 40 per cent of our sample are 'obviously' head-tail asymmetric.Filaments with a single peak tend to be more asymmetric (having a higher value of skewness or p) than those with more peaks, which is not unexpected.2

Correlation of the parameter p with other physical parameters
To understand whether the degree of head-tail asymmetry is related to the evolution of filaments, we analyse the relationship between parameter p and other parameters of filaments, including physical scale, velocity width, mass, maximum column density, and orientation.
The column densities and masses are simply calculated using the X-factor approach, with X CO = 2×10 20 cm −2 (K km s −1 ) −1 and a ±30 per cent uncertainty, as recommended by Bolatto et al. (2013): where the sum is over all the pixels of a filament, W(CO) is the integrated line intensity of 12 C 16 O (J = 1 − 0), ∆Ω is the angular area of a pixel, and d is the distance to the filament.
We do not see any significant correlation between parameter p and other physical parameters in the resulting scatter plot Fig. 14, while we do notice that the maximum column density, N(H 2 ) max , seems to be weakly correlated with p.This could be understood if the very head-tail asymmetric filaments have been undergoing mass concentration towards one end, leading to the elevated column density.

CONCLUSIONS AND DISCUSSIONS
In this work, we use part of the 12 CO (J = 1 − 0) data from the MWISP project to study the orientation and head-tail asymmetry of filaments with a simple straight morphology.We construct a series of samples during the process of excluding curved clouds.From the results, a distinct bimodal orientation distribution is seen (Figs. 5  and 6), and the bimodal orientation becomes more pronounced as the screening criteria become more stringent (while still ensuring the inclusion of a certain amount of clouds).We also find that certain filamentary molecular clouds exhibit distinct head-tail asymmetry, with one end being wider or brighter than the other.The main conclusions are as follows: (i) We find a bimodal orientation distribution of the sample.By employing a double Gaussian fitting, we obtain that the centres of the two Gaussian functions are −38.1 • and 42.0 • , and the 1σ of the two Gaussian functions are 35.4• and 27.4 • , respectively.
(ii) Most filaments exhibit a small number of peaks.Specifically, 43.5 per cent of the filaments' one-dimensional integrated intensity profiles show either a single peak or double peaks, while ∼ 26 per cent of filaments have four or more peaks.We note that the number of peaks could be affected by the spatial resolution of the data.
(iii) Based on the complementary cumulative distribution function (CCDF) of the two asymmetric parameters, absolute skewness and p, we find that filaments with a single peak tend to be more asymmetrical than those with more peaks, and about 40 per cent of the whole sample have obvious head-tail asymmetry.
The orientation of filaments can be influenced by various factors, including the local environment such as the magnetic field, physical processes within the cloud, and their formation and evolution.Since turbulence can compress gas in any direction, we may use the bimodal orientation distribution of filaments to exclude the possibility turbulence solely dominates the formation of filaments.
The data used in this work are in the form of a data cube in position-position-velocity (ppv) space, rather than position-positionposition (ppp) space.The impact of projection effects cannot be eliminated.If the three-dimensional orientations of filaments in the sample could be obtained, their distribution might indeed be entirely different.The filaments might not even exhibit a filamentary structure.However, this does not negate the fact that the two-dimensional projection of orientations exhibits bimodality.
Several studies on filaments have concluded that filaments are oriented along the Galactic plane.However, it is noteworthy that these studies focus on giant molecular filaments (GMFs).Goodman et al. (2014) find that the very long and thin infrared dark cloud 'Nessie' (Jackson et al. 2010) lies not just parallel to the Galactic plane, but  2021) found that most of the filamentary structures in the 12 CO and 13 CO emission do not show a global preferential orientation either parallel or perpendicular to the Galactic plane, and they propose that a selection bias may be introduced by the GMF identification in the general trends found for the GMFs, which appear to be aligned with the Galactic plane.
Our work differs from these previous studies in a few aspects.Firstly, in the sample selection, we focus on small-scale filaments with simple straight morphology, as opposed to large-scale filaments (GMFs).Secondly, in the orientation calculation, we derive orientations from the major axis direction through elliptical fitting, while previous studies typically employed methods based on ideas from computer graphics.For example, Soler et al. (2013) and Soler et al. (2021) employ methods based on the direction perpendicular to the local intensity gradient and the direction characterized by the eigenvector of the Hessian matrix, respectively, to statistically analyze the orientation of the density structures.In Appendix D, we also applied the method of Soler et al. (2013) to our data to cross-check our results.Finally, the distribution of filament orientation is expected to vary across different regions of the Galaxy, and our work only focuses on the molecular clouds in the second Galactic quadrant.All these factors contribute to the disparities between our and previous results.Li et al. (2013a) found that filaments in the Gould Belt tend to align either parallel or perpendicular to the local mean magnetic field directions.They interpreted this as evidence that the magnetic fields within the intercloud medium are strong enough to influence gravitational contraction, resulting in the formation of flattened condensations oriented perpendicular to the magnetic fields, and channeling the turbulence and leading to the alignment of filaments along the magnetic fields.Planck Collaboration et al. (2016b) performed a statistical analysis to assess the relationship between the magnetic field projected onto the plane of the sky (B B B ⊥ ) and the gas column density (N H ) structures.They observed a systematic change in the relative orientation between B B B ⊥ and N H structure across different regions, with the orientation being parallel in areas with the lowest column density and perpendicular in areas with the highest column density.
While in the current paper we did not notice any correlation between orientation and other physical parameters, we do notice hints of presence of a global pattern of the filament orientations (Fig. 8).
In Appendix E we show that a global sinusoidal arrangement of the filaments could lead to a bimodal orientation distribution.In our next study, we will investigate the potential correlation between the local magnetic field and the orientation of the filaments.
The head-tail asymmetry of filaments may be explained by the 'edge effect' that has been previously studied.Yuan et al. (2020) conducted an investigation on filament S242, which contains concentrated massive clumps and YSO clusters at its end regions.They proposed a two-stage fragmentation scenario for filament S242.Firstly, gravitational focusing leads to the preferential concentration of material at the end regions, resulting in a significant acquisition of mass by the end-clumps.Secondly, longitudinal accretion becomes dominant in the vicinity of the high-mass clumps.This 'edge effect' explains the accumulation of material at both ends of a filament and has been observed in several studies.However, our current study emphasizes the concentration of matter towards one end rather than both ends as described by the 'edge effect'.The head-tail asymmetry of filaments may indicate different physical conditions at the two ends of a filament.By further investigating this asymmetry and its correlation with the local environment, we can deepen our understanding of the structure and evolution of molecular clouds.

APPENDIX A: TEST OF THE METHOD OF REMOVAL OF CURVED CLOUDS
To ensure our sample selection is unbiased, here we test whether the method used in Section 2.2 for removal of curved clouds would artificially introduce any pattern to the orientation distribution.We generate 10,000 ellipses with varying aspect ratios and uniformly distributed orientations (from −90 • to 90 • ).Then we calculate their s parameter using equation 7. It is expected that this s parameter is related with the aspect ratio, but not the orientation.We then apply different thresholds of s to select different sets of subsample to see whether a preference for certain orientation angle would emerge.The results, depicted in Fig. (A1), indicate that the method introduces no orientation preference, thus we conclude that no bias is introduced during sample construction.

APPENDIX B: PARAMETERS OF FILAMENTS WE USED
For the 522 filaments used for the analysis in this work, we list their relevant parameters in Table B1.The complete table is available online.These parameters include the source name, orientation, coordinates, velocity, estimated distance, physical scale, velocity width, mass, maximum column density, and average column density.

APPENDIX C: EXAMPLE IMAGES OF FILAMENTS WITH VARYING DEGREES OF ASYMMETRY
In Fig. (C1), we present the 12 CO (J = 1 − 0) velocity-integrated intensity maps of some filaments, sorted in ascending order of parameter p.

APPENDIX D: CROSS-CHECK WITH DIFFERENT ANALYSIS METHOD
Since our analysis involves a sample selection step and thus cannot be directly applied to other dataset, we opt to cross-check our results using alternative methods.Specifically, we employ the orientation description method used by Soler et al. (2013), to conduct a crosscheck on our data.
This method primarily relies on the gradient of the density to characterize the directionality of the density structures.We employ Gaussian first derivative kernels in both horizontal and vertical directions to convolve the strength data separately, obtaining gradients in both directions.We use the Sobel operator (an approximate Gaussian derivative kernel) to perform this calculation: Then convolve the 2D integrated intensity with the Sobel operator: where I is the 2D integrated intensity.
According to equation (10) in Soler et al. (2013), the 'orientation description' of each pixel is given by: The weakening of the bimodal feature could be related with that the Sobel operator is applied to all the pixels of a filament.To test this, we generate a series of ellipses with well-defined orientations, while randomly setting the aspect ratio between 2 and 5.The intensity of these ellipses follows a two-dimensional Gaussian distribution overlaid with some Gaussian noise.Fig. (D2) displays the results for these ellipses, where we observe that the orientation of the pixels exhibits a certain degree of dispersion.This observation helps explain the dilution of the bimodal distribution seen in Fig. (D1).

Figure 1 .
Figure 1.Integrated intensity map of some clouds of the full sample that is included in or excluded from our sample.Panels (a) and (b) show examples of the filamentary molecular clouds that are included in our final sample.Panel (c) shows the curved cloud that we exclude from the final sample, and panel (d) shows the example that is excluded because of large angular scale and network-like substructures.

Figure 2 .
Figure 2.These two cloud silhouettes are examples selected from the filaments obtained by elliptical fitting.The blue line represents the major axes obtained from elliptical fitting, which indicate the orientation of the filament.The top panel shows a straight cloud (included in our final sample), which Var(d)=1.069pixel 2 , and the bottom panel shows a curved cloud (not included in our final sample), with Var(d)=6.021pixel 2 .The pixel size here is 30 ′′ .

Figure 3 .
Figure 3.The overlay map of filaments obtained by different thresholds.Panels (a), (b), and (c) show the subsamples obtained using thresholds of 0.095, 0.1, and 0.105, respectively, with the number of samples marked in each panel.Panel (d) presents the overlay map of 310 randomly selected samples.

Figure 4 .
Figure 4. Schematic diagram illustrating the calculation of the average width.The red line is used as an example to illustrate the calculation process, its direction is perpendicular to the orientation of the filament.The pixel size here is 30 ′′ .

Fig. 5 Figure 5 .
Fig. 5 shows the distribution of the angle between the orientation of filaments in our final sample and the Galactic plane.We can see that the orientation of the filaments mainly concentrates on two ranges.To obtain the parameters of the bimodal orientation distribution, we use a double Gaussian function to fit the histogram.The centres of the two Gaussian functions are −38.1 • and 42.0 • , and the 1σ of the two Gaussian functions are 35.4• and 27.4 • .

Figure 7 .
Figure 7.The fraction of different relative orientations for each angular separation interval.The six lines represent the relative orientations corresponding to the six intervals indicated by the labels, and the horizontal axis represents the corresponding angular separation.

Figure 8 .
Figure 8. Distribution of all filaments in our final sample in Galactic coordinates.Each filament is coloured in black (for orientation angle > 0) or blue (for orientation angle < 0).The background colour indicates the number density of filaments in our sample.The bimodal orientation distribution is already visually perceivable in this figure.

Figure 9 .Figure 10 .
Figure9.The (lack of) correlation between orientation and other physical parameters of the filaments in our final sample.R 1 is the correlation coefficient obtained from sources with nagative orientation angles, and R 2 from sources with positive orientation angles.

Figure 11 .
Figure 11.Distribution of the number of peaks of the filaments in our final sample.

Figure 12 .Figure 13 .
Figure 12.The colour shows the intensity distribution along the major axis (y axis, with length normalized to 1) of the filaments in our sample.The filaments are sorted by the parameter p (top panel), and are oriented such that the 'heavier' half are on the upper side.The right panel shows a sum of the intensity profile of all the filaments.

Figure 14 .
Figure 14.The (lack of) correlation between parameter p and other physical parameters of the filaments in our final sample.The background colour indicates the number density of the scatter points.

Figure A1 .
Figure A1.Orientation distribution for samples of mock filamentary clouds with different threshold for the s parameter.Panel (a): orientation distribution of the full 10000 mock filaments; here the histogram is based on the given value of the orientation angles.Panel (b): orientation of these filaments calculated using the elliptical fitting method.Panels (c, d, e, f): orientation distribution of subsamples obtained with different thresholds of the s parameter for the removal of curved clouds (thresholds are 0.09, 0.07, 0.05, 0.03, respectively).The numbers in each panel are the size of each subsample.
Fig. (D1) shows the distribution of pixel-based orientations in the final samples calculated using the method proposed by Soler et al.We can see that the pixel-based orientation exhibits a weak bimodal distribution.The weakening of the bimodal feature could be related with that the Sobel operator is applied to all the pixels of a filament.To test this, we generate a series of ellipses with well-defined orientations, while randomly setting the aspect ratio between 2 and 5.The intensity of these ellipses follows a two-dimensional Gaussian distribution overlaid with some Gaussian noise.Fig.(D2) displays the results for these ellipses, where we observe that the orientation of the pixels exhibits a certain degree of dispersion.This observation helps explain the dilution of the bimodal distribution seen in Fig.(D1). p=0

Figure C1 .
Figure C1.Example images of filaments with different value of parameter p.

Figure E1 .
Figure E1.Orientation distribution of the tangents of a sinusoidal curve.Panel (a): an example sinusoidal curve.The yellow arrows illustrate the filaments that are aligned with the curve.Panel (b): distribution of the angle of tangential directions (relative to the x axis).The blue histogram is obtained from numerically differentiating the curve in panel (a) and then calculating the orientation angle of the tangents, while the orange histogram is obtained the same way excepts that noise is artificially added to the two components of the tangent vector.The magenta curve is calculated using equation (E2).

Table B1 .
The parameters of the sources we select in our sample.The complete table is available online.