can be identified through extinction maps at optical and infrared wavelengths (Schneider & Elmegreen 1979 ;Hatchell et al. 2005 ;Myers 2009 ;Jackson et al. 2010 ), as well as through molecular line observations (Loren 1989 ;Mizuno et al. 1995 ;Schneider et al. 2010 ). Far-infrared and submillimeter dust emission maps have also shown the presence of filaments. In particular, the Herschel survey revealed the ubiquity of filaments and their importance in the star formation process with dense star-forming cores being distributed along filaments akin to beads strung along a thread (Schneider & Elmegreen 1979 ;Andr é et al. 2010 ;Molinari et al. 2010 ;Andr é et al. 2014 ). This illuminates a possible picture where the dense gas within molecular clouds first assembles into dense filaments and then forms dense star-forming cores through fragmentation, as already proposed by Schneider & Elmegreen ( 1979 ).

It has long been suggested that molecular clouds have a fractal nature, with a hierarchy of structures from clouds to clumps to cores (Scalo 1990 ;Falgarone, Phillips & Walker 1991 ;Williams, Blitz & McKee 1999 ;Heyer & Dame 2015 ). Similarly, observations show that filamentary structures span several orders of magnitude in scale (see e.g. the re vie w by Hacar et al. 2023 ), from large-scale giant MNRAS 528, 6370-6387 (2024) molecular filaments associated with spiral arms on scales of dozens to hundreds of pc (Jackson et al. 2010 ;Goodman et al. 2014 ;Wang et al. 2015Wang et al. , 2016 ; ;Zucker, Battersby & Goodman 2015 ;Zhang et al. 2019 ;Ge & Wang 2022 ), to filaments on scales of 1-10 pc (Kainulainen & Tan 2013 ;Li et al. 2016 ;Xiong et al. 2019 ;Schisano et al. 2020 ), and even down to fiber structures on scales of 0.1 pc (Hacar et al. 2013 ;Kainulainen et al. 2016 ;Chung et al. 2021 ). This continuous distribution o v er such a vast range of scales suggests that filaments possess a hierarchical nature. For example, Orion A, when observed at large scales, appears as an elongated filamentary cloud extending up to ∼90 pc (Großschedl et al. 2018 ). At intermediate scales, Orion A can be resolved into numerous parsec-size filaments (Johnstone & Bally 1998 ;Nagahama et al. 1998 ), while at higher resolutions, small-scale sub-parsec filaments can be identified (Wiseman & Ho 1994 ;Takahashi et al. 2013 ;Hacar et al. 2018 ;Suri et al. 2019 ).

Notably, the choice of scale may result in divergent conclusions regarding filament stability. Larger-scale filaments tend to be supercritical, meaning they have line masses exceeding the critical threshold for gravitational collapse, exhibiting greater line masses. In contrast, smaller-scale filaments are more inclined towards subcritical states, where their line masses fall below this threshold, thus exhibiting smaller line masses (Hacar et al. 2013 ;Henshaw et al. 2014 ;Chen et al. 2019 ). Filaments at various scales ef fecti vely sample distinct gas densities within the ISM. The mean gas density in molecular filaments has been observed to decrease with increasing length (Hacar et al. 2023 ). Parsec-scale filaments typically exhibit average densities of ∼10 3 -10 4 cm -3 , while sub-parsec filaments, commonly referred to as fibres, display densities of > 10 5 cm -3 .

Many theoretical works have sought to investigate the formation and evolution of filaments at different scales. Turbulence plays a crucial role in the formation of large-scale filament networks. The origin of the intricate fibres within filaments is a subject of ongoing debate. One proposed scenario, known as the top-down approach, suggests that these structures may arise from their parental filaments, as indicated by the large-scale coherence and parallel organization of fibres with respect to the main filament, as seen in regions like B213-L1495 (Tafalla & Hacar 2015 ;Clarke et al. 2017 ). Ho we ver, some simulations support an alternative bottom-up mechanism, in which small subsonic filaments initially form within a turbulent medium and are subsequently gathered through collapse and shear flows (Smith et al. 2016 ). Hacar et al. ( 2023 ) integrated observations from the past decade to provide an updated description of filaments at different scales and environments, and identified a scaling relation for filament mass and length L ∝ M 0.5 , which is derived from the 1st and 2nd Larson relations (Larson 1981 ). Hacar et al. ( 2023 ) suggest that the distribution of filaments in the mass-length ( M -L ) diagram might be controlled by how they evolve in terms of a balance between accretion, fragmentation, collapse and destruction. Various studies show velocity gradients tend to be perpendicular to the main axis of the filament, suggesting ongoing gas accretion on to filaments (Schneider et al. 2010 ;Beuther et al. 2015 ;Dhabal et al. 2018 ;Williams et al. 2018 ;Shimajiri et al. 2019 ;Chen et al. 2020 ). The estimated accretion rates, denoted as ṁ , typically range from a few 10 to a few 100 M Myr -1 pc -1 (Kirk et al. 2013 ;Palmeirim et al. 2013 ;Schisano et al. 2014 ;Bonne et al. 2020 ;Gong et al. 2021 ). Environmental effects play a role in influencing the observed accretion rates, as higher mass filaments tend to be embedded in regions with higher background column density, which is expect to obtain a higher accretion rate o v er 500 M Myr -1 pc -1 (Heitsch 2013 ;G ómez & V ázquez-Semadeni 2014 ;Rivera-Ingraham et al. 2016 ). Moreo v er, studies suggest that accretion may not stabilize the filaments against gravitational collapse, it can induce fragmentation leading to the formation of dense cores and possibly star formation (Mac Low et al. 1998 ;Hacar & Tafalla 2011 ;Arzoumanian et al. 2013 ;Hacar et al. 2013 ;Seifried & Walch 2015 ;Clarke et al. 2017 ;Heigl, Gritschneder & Burkert 2020 ).

Accretion along filamentary structures heightens gas density, thereby influencing instantaneous fragmentation properties by reducing the Jeans length of these structures, especially at filament intersection points or along the filament's entire length. Largescale filaments are inherently turbulent, with supersonic compression creating local o v erdensities and thus instigating fragmentation before the decay of turbulent energy (Mac Low et al. 1998 ;Mac Low 1999 ). Numerical simulations further reveal that ongoing material accretion on to the filament significantly impacts gravity-induced fragmentation, deviating notably from the fragmentation of initially quiescent, pre-existing filaments (Inutsuka & Miyama 1992 ;Clarke, Whitworth & Hubber 2016 ). Due to the turbulent energy cascade, accreted material is expected to be moderately turbulent, with fragment locations determined by turbulent motions in the transto mildly supersonic regime where fragmentation is shifted from a gravity-dominated to turbulence-dominated process, highlighting the significance of turbulent motions and moderate density enhancements during filament formation (Seifried & Walch 2015 ;Clarke et al. 2017 ;Heigl, Gritschneder & Burkert 2020 ). Chira et al. ( 2018 ) suggest that filament fragmentation must occur before reaching the critical line mass. This is because once the filament line mass exceed the critical line mass, the filament will collapse radially in a free-fall time.

This paper is organized as follows. Section 2 introduces the methodology with our simulation, algorithms, and properties calculation. Section 3 presents the basic result of the simulated gas properties. We discuss the evolution scenarios and the hierarchical characteristics of filaments in Section 4. Our main conclusions are presented in Section 5.

The Cloud Factory simulations are a suite of high-resolution numerical simulations designed to study the behaviour of the interstellar medium (ISM) in a typical spiral galaxy across various scales, ranging from the entire galaxy down to individual filaments and clumps within specific molecular clouds (Smith et al. 2020 ). The Cloud Factory simulation suite utilizes a modified version of the AREPO code (Springel 2010 ;Pakmor et al. 2016 ), which includes various physical and chemical modules that account for different processes occurring in the cold molecular ISM. The model includes a galactic gravitational potential, time-dependent evolution of CO and hydrogen chemistry, ultraviolet extinction, dust absorption, star formation via sink particles, and feedback from supernova explosions.

For the galactic gravitational potential, to minimize computational effort, an analytical approach is adopted to model the large-scale galactic potential. The self-gravity of the gas is calculated using the standard AREPO gravitational tree (Springel 2010 ). The axisymmetric component of our analytical gravitational potential is based on the best-fitting potential described in McMillan ( 2016 ), which was developed to align with various observational and theoretical constraints for the Milky W ay. W e incorporated a spiral perturbation to the axisymmetric potential, following the method outlined in MNRAS 528, 6370-6387 (2024) Smith, Glo v er & Klessen ( 2014 ), incorporating a four-armed spiral component from Cox & G ómez ( 2002 ).

The chemical evolution of the gas in this study follows the model described in Smith, Glo v er & Klessen ( 2014 ) and incorporates the hydrogen chemistry framework proposed by Glo v er & Mac Low ( 2007a , b ). The formation and destruction of CO are treated in a simplified manner based on Nelson & Langer ( 1997 ). The hydrogen chemistry model encompasses processes such as H 2 formation on grains, H 2 destruction via photodissociation and collisional dissociation of atomic hydrogen, H + recombination in both gas and grain phases, and cosmic ray ionization Glo v er & Mac Low ( 2007a ). The evolution of the CO abundance is determined by assuming a limited CO formation rate through an initial radiative association step, and the CO destruction rate is primarily due to photodissociation. Full details of the combined network and a comparison of its accuracy and efficiency with other approaches can be found in Glo v er & Clark ( 2012 ).

We utilize the framework of sink particles, which are non-gaseous cells that serve as representations of star formation sites. Cells that surpass a critical density, ρ c , are considered potential candidates for conversion into sink particles. Ho we ver, these cells must successfully pass a series of energy checks first to confirm whether the gas is unambiguously gravitationally bound and has inwardly directed velocities and accelerations. Moreover, prior to the transformation into a sink particle, a cell must be situated at a local gravitational potential minimum and outside the accretion radius of any existing sink particle. For more details on the methodology used for the creation of sink particles, please refer to the work by Tress et al. ( 2020 ). The formation of sinks below protostellar densities implies that not all mass is converted into stars. The formation density and accretion radius in the simulation are determined based on specific criteria and physical constraints. Sink particles are created in cells with densities exceeding a critical density, ρ c . These cells also need to undergo energy checks to ensure they are gravitationally bound and have inward velocities. The accretion radius of the sink particles is initially set to correspond to the Jeans length at their creation density. Over time, the accretion radius increases in order to maintain a constant acceleration at the sink surface. This process ensures that the sinks can ef fecti vely accrete mass from neighbouring bound cells. The formation density and accretion radius are adjusted based on the target mass resolution, ensuring accurate representation of star-forming regions and preventing artificial fragmentation in the simulation (Greif et al. 2011 ). We use a star formation efficiency of 33 per cent in keeping with Matzner & McKee ( 2000 ) because we place our sinks at densities typical of protostellar cores. The stellar content of the sinks is obtained by multiplying the sink mass by the assumed star formation efficiency.

To estimate the supernova feedback rate, we use two methods: (1) random supernova explosions, and (2) supernovae tied to sinks. In the first approach, we randomly select points from the initial gas density profile chosen for the disc, assuming a rate of one supernova per 300 yr. This is typical for Type Ia supernovae in the Milky Way (Diehl et al. 2006 ), which should be decoupled from the gas distribution. Ho we ver, purely random feedback may result in unrealistic cloud properties when considering self-gravity, as it fails to destroy large molecular cloud comple x es (Gatto et al. 2015 ;Walch et al. 2015 ) and does not consider the causal link between massive stars and supernovae. Our secondary feedback method aims to account for this issue by incorporating a randomly distributed supernova component, with a frequency of one supernova every 300 yr, as proposed by Tsujimoto et al. ( 1995 ), to simulate Type Ia superno vae. Additionally, superno vae from sink particles are also included in the model, irrespective of the refinement level at which they were formed. We employ the stellar initial mass function (IMF) as outlined in Kroupa ( 2002 ) and determine the number of high mass stars within the stellar composition of sink particles based on the methodology described by Sormani et al. ( 2017 ). As the supernovae from each sink explode individually, this naturally results in clusters of supernova explosions. For further details on our supernovae model, see Tress et al. ( 2020 ) and Smith et al. ( 2020 ).

It also includes a four-armed spiral component and a consistent spiral perturbation to the potential. The gas chemistry description adopts a simplified treatment that assumes a direct conversion between the C + and CO abundances, following the approach of Nelson & Langer ( 1997 ). The non-equilibrium hydrogen chemistry involves reactions between molecular (H 2 ), atomic (H), and ionized (H + ) hydrogen, electrons, cosmic rays, dust grains, and the UV radiation field. The star formation model uses a hybrid approach based on sink particles, which can represent either individual stellar systems or clusters of stars.

In this study, three star-forming molecular complex regions are chosen from the feedback-dominated Cloud Factory model. We identify the filamentary structures and examine them from an e volving vie wpoint using time-e volved snapshots of these regions. These regions are selected based on their distinctive characteristics. Region A is a large molecular complex spanning over 100 pc in size. It is in the early stages of evolution, with barely any sinks. Region B is a relatively small molecular cloud with a low column density and low star formation. Region C contains multiple independent molecular clouds, and it is more evolved than the other two regions.

By examining these three regions, our objective is to span a range of scenarios that are likely to occur for molecular clouds of different sizes and evolutionary stages. The regions are tracked for varying time periods based on their ability to maintain consistency as objects in the time evolution.

To generate synthetic observations, we perform the radiative transfer simulation using the POLARIS code (Reissl, Wolf & Brauer 2016 ), which is publicly available. To ensure consistency in resolution and physical quantities, we utilized the same Voronoi grid as the AREPO simulations as in Izquierdo et al. ( 2021 ). We use the pre-calculated optical properties by Draine & Flatau ( 2013 ), which consist of a mixture of 62.5 per cent silicate and 37.5 per cent graphite. Ho we ver, this model does not account for the formation of ice mantles around cores, which lead to increased brightness within cores (Ossenkopf & Henning 1994 ). To address this, we increase the maximum radius to 8 × 10 -6 m to obtain a value of κ ∼ 1 cm 2 g -1 , consistent with the results from Ossenkopf & Henning ( 1994 ) (see also in prep).

We use a fixed dust temperature, T d , as an input parameter for our radiative transfer simulations. This is because the preceding Cloud Factory simulation does not include accretion luminosity or photoionizing feedback, resulting in a dust temperature that was relati vely lo w compared to real observ ations. To address this, we fix T d to a pumping value of 20 K, which is similar to the peak of the temperature distribution in the probability distribution function of actual observations (Lin et al. 2016 ). However, this is still a simplified assumption for the dust model. In future work, we plan to impro v e the sink model to further enhance the accuracy of the dust temperature.

For each snapshot of the target regions, we extract Voronoi grid data from the Cloud Factory, using a default box size of 200 pc as input. We then model dust emission at 70, 160, 250, 350, and 500 µm by following the Herschel PACS instrument (Poglitsch et al. MNRAS 528, 6370-6387 (2024) 2010 ) and the SPIRE instrument (Griffin et al. 2010 ) guidelines. As input parameters to POLARIS , the specified distance and grid sizes can be easily changed, which provides us with a method to model the same region at different distances, resulting in different spatial resolutions. Specifically, we adopt a fixed angular resolution of ∼20 arcsec, which is comparable to Herschel observations (Hennemann et al. 2012 ;Lin et al. 2016 ;Arzoumanian et al. 2019 ). We choose three distinct distances for the observer detector parameter: 2.5, 5, and 10 kpc. This allows us to examine the same region at three distinct physical resolutions, simulating the effects of distance.

The output data is a uniform grid with different space resolutions of 800 × 800, 400 × 400, and 200 × 200 pixels, corresponding to 0.25, 0.5, and 1 pc pix -1 , respectively.

In order to derive the properties of the clouds in a manner analogous to observations, we need to perform spectral energy distribution (SED) fitting. SED fitting involves comparing the observed flux densities of an object at different wavelengths with model SEDs that assume a single blackbody component in each line of sight. By fitting the observed SED with a model SED, the physical properties of the object, such as temperature and mass, can be determined. Generally, the SED-fitting method requires adjusting the observed flux density in different bands to a common angular resolution that matches the longest wavelength band. In our case, since we use the same resolution for all bands during the radiative transfer simulation, we do not need to perform any smoothing prior to SED fitting.

We use a dust opacity law that is mentioned by Hildebrand ( 1983 ). This law expresses the flux density S ν at a given observing frequency ν as

and,

B ν is the Planck function for a given temperature T d , derived by

where

is the dust opacity, β is the dust opacity index, is the considered solid angle, µ = 2.8 is the mean molecular weight, m H is the mass of a hydrogen atom, and κ 230 = 0 . 09 cm 2 g -1 is the dust opacity per unit mass at 230 GHz. We interpolate the opacity κ ν from Ossenkopf & Henning ( 1994 ). A gas-to-dust mass ratio of 100 is assumed. We use the method described in Lin et al. ( 2016 ) to fit the molecular hydrogen column density N H 2 with all the radiative transfer simulation bands (70, 160, 250, 350, and 500 µm). We fix the value of T d equal to 20 K as discussed in Section 2.2 , therefore, the T d has no physical meaning, and we only consider the N H 2 result for subsequent analysis.

After deriving the maps of N H 2 , we define the boundaries of the clouds by creating masks for each region where N H 2 is greater than 10 21 cm -2 . We then calculate the gas mass for each pixel using the H 2 column density maps, based on the formula

where µ is the mean molecular weight per hydrogen molecule, which is assumed to be 2.8 in this work. The mass of atomic hydrogen is represented by m H , and D represents the distance to the object. The solid angle element is denoted as d . Finally, we sum all the pixels within the mask to obtain the total 'observed' gas mass of the molecular clouds. Additionally, we calculate the mean N H 2 within the boundary mask. The calculated gas mass is denoted as M H 2 , dust for synthetic dust data and M H 2 , raw for original projection data, while M H 2 , raw is calculated with the projection H 2 column density map from the original simulation data. In the Cloud Factory simulations, dense star-forming cores are represented by sink particles. These particles are associated with both a stellar mass M stellar (33 per cent of the sink mass) and a gaseous mass M sinkgas (67 per cent of the sink mass). Consequently, the H 2 mass derived from radiative transfer calculations misses a significant portion of the total cloud mass. To address this issue, when calculating the total mass for the filament, we consider the 'missing' gaseous component of the corresponding sink mass, in addition to the molecular gas mass derived from radiative transfer column density maps. Additionally, we examine the discrepancies between the raw data sets and the data produced by POLARIS . The detailed results can be found in Appendix A .

Notably, there is a good alignment in column density between the two data sets within the dense gas regions.

In this work, we primarily use the FILFINDER algorithm (Koch & Rosolowsky 2015 ) to identify structures. This extracts filaments from a column density map by creating a mask beyond an adaptive column density threshold, derived from the intensity image, and then reducing the mask to a skeleton via the Medial Axis Transform. It returns all the possible filamentary structures across the input map and is able to extract hierarchical filament structures uniformly (Liu & Li 2018 ;Zucker, Battersby & Goodman 2018 ;Xiong et al. 2019 ;Zhang et al. 2021 ). Ho we ver, FILFINDER can only be applied to 2D diagrams, although Kim et al. ( 2023 ) recently presented a method for identifying filaments in 3D position-position-velocity data based on applying FilFinder to velocity channel maps. We apply FILFINDER on the H 2 column density maps originating from the RT simulated dust emission, as described in Section 2.3 . We first define a setting of 'flatten percent = 99.9' to emphasize the main structure. This means that during image contrast adjustment, intensities up to the 99.9th percentile of the data are retained. This approach helps to minimize the impact of extreme brightness values in the image, which can interfere with subsequent image processing steps like filtering and feature detection. As the clouds are in different evolutionary states, their density varies significantly, so rather than choosing a single threshold, we aim to create a mask with a threshold o v er the top 1 per cent of the data for each region such that only the dense structure is included. This is similar to real observations, where no universal threshold defines a filament. Instead, different thresholds are applied to various regions based on the column density contrast with their environment, aiming to highlight the morphological features.

The value of the threshold for each mask is shown in Table 1 . The column density threshold is kept constant for each region for different snapshots. It should be noted that this threshold is relatively low and is only used for a pre-process to maintain most of the structures. Using the mask, the algorithm generates skeletons and branches. The shortest lengths of these are set by tw o variables: 'sk el thresh' and 'branch thresh'. 'skel thresh' is set to 8 pc, which corresponds to 8, 16, and 32 pixels in three respective resolutions. On the other hand, MNRAS 528, 6370-6387 ( 2024) 5 3,5,7,9,11,13,15 a The number of the snapshots series indicate the evolving time with a unit of 10 5 yr.

'branch thresh' is set to 4 pc, corresponding to 4, 8, and 16 pixels in the three resolutions, respectively. The skeleton and the branches trace the bones of the filament. The skeleton is the longest part of the main structure, and the branches are shorter parts connected to the skeleton. Here, we only focus on the skeleton part of the filaments when defining the length. The algorithm also derives the widths of the filaments, which are fitted by a Gaussian with a mean fixed to 0 and a constant background.

With the length and width determined, the algorithm also creates a boundary mask of filaments containing both skeleton and branches with a threshold larger than 10 21 cm -2 , as shown in Fig. 1 . This mask enables the calculation of the total mass of the filaments from the input image. The length of the filament is computed as the product of the total pixel number of the identified skeleton and the physical length of a pixel. Using the modelled skeleton and width, we generate the filament mask and calculate the mass of each pixel from the column density map. Subsequently, we sum up the mass of the pixels within the mask to obtain the mass of the filament. Finally, we extract the sink data from the original Cloud Factory data set, where each sink comprises 67 per cent of the gas portion and 33 per cent of the stellar mass. If the sink position is within the filament mask, we then sum up the mass derived from the radiative transfer dust image and 67 per cent of the sink mass to obtain the total mass of each filament.

Before identifying filaments within our data set, it is important to determine how the gas around the filaments affects them. Fig. 2 displays the synthetic H 2 column density maps for the first and last snapshots of each of the three regions. Table 2 summarizes the column densities and masses for the synthetic observations with a resolution of 0.25 pc pix -1 , which were calculated using the methodology described in Section 2.3 . Notably, the molecular gas mass M H 2 , dust in Table 2 excludes the portion of gas within the sinks, as the radiative transfer simulation is restricted to modelling gas that has not transitioned into a sink.

We trace different durations and use varying time steps for the snapshots of these three regions, due to their distinct stellar evolutionary stages. For Region A, as indicated in Table 2 , the stellar mass, M stellar , in the initial snapshot 3 is significantly lower than in the other regions, with almost no star formation. As a result, the evolutionary duration of Region A is the longest, spanning 2 Myrs. In contrast, Region C has a substantial amount of star formation right from its initial snapshot, leading us to track it for the shortest duration of 0.4 Myrs. This is because of the omission in our simulation of early stellar feedback components, including jets, outflows, and photoionization. Consequently, filament evolution becomes distorted after significant star formation, making further evolutionary tracking less meaningful. Region B is intermediate between the two, with a tracked duration of 1.2 Myr. In this case, we opted for different time-steps to ensure a comparable number of snapshots across the three regions.

Region A comprises a complex molecular system with a mass of molecular gas M H 2 , dust ≈ 2 × 10 5 M in the masked region. We examine 5 snapshots, 3, 8, 13, 18, and 23, representing the time (in increments of 10 5 yr) since we zoom into the region in the Cloud Factory simulation. Ho we ver, due to its complex molecular structure, it is difficult to track the evolution of the individual identified filaments in Region A from beginning to end, as these filaments interact and even intertwine with each other. Table 2. The gas properties in three regions, as derived from all snapshots of the highest resolution (0.25 pc pix -1 ) synthetic observational data and the original simulation data.

Snapshot

A 3 5.91 680 96 64 A 8 6.44 339 1014 676 A 13 6.57 008 1500 1000 A 18 6.28 333 1796 1197 A 23 5.77 101 2028 1352 B 3 2.61 337 1706 1137 B 5 2.61 105 2414 1609 B 7 2.59 355 3620 2413 B 9 2.65 440 4572 3048 B 11 2.65 590 5249 3499 B 13 2.60 409 5868 3912 B 15 2.61 445 6302 4201 C 3 4.52 442 18 024 12 016 C 4 4.54 925 21 027 14 018 C 5 4.61 708 23 501 15 667 C 6 4.80 836 28 485 18 990 Note . (1) Re gion of clouds (2) Evolving time snapshot, in units of 10 5 yr (3) The mean H 2 column density of clouds derived by SED fitting via dust emission (4) Total mass of all clouds in this re gion, deriv ed from dust emission (5) Total mass of the gas portion in sinks (6) Total mass of the stellar portion in sinks.

Region B consists of a single cloud that we analyse from snapshot 3 (0.3 Myrs) to snapshot 15 (1.5 Myrs) with a time-step of 0.2 Myrs. The total gas mass in this region is M H 2 , dust ≈ 3 × 10 3 M . The column density of Region B is lower compared to the other regions, and we can track one particular filament o v er a longer time period because of the region's simplicity.

Region C contains multiple isolated filamentary objects, making it a suitable case for tracing the average properties of multiple filaments evolving o v er time. The total mass of Region C is M H 2 , dust ≈ 3 × 10 5 M . We analyse four snapshots, from snapshot 3 (0.3 Myrs) to snapshot 6 (0.6 Myrs).

In this work, we identify filaments across the three regions at three distinct resolutions. Ho we ver, the number of filaments identified varied depending on the resolution and the specific snapshot in time. Initially, we find a total of 189 filaments across all snapshots and resolutions combined. The filament mass M fil , length L fil , and the stellar mass M stellar of these filaments are derived using the method mentioned in Section 2.3 . The line mass (mass per unit length) of the filaments is defined as m = M fil / L fil .

One of the chief findings of Hacar et al. ( 2023 ) is that filament structures are hierarchical, with nested smaller sub-structures contained within larger filaments. Thus, a filament that appears to be a single object may reveal smaller substructures when observed at higher resolution. Fig. 3 shows the identified filament skeletons in Region A at three different resolutions for snapshot 3. As expected, for the higher resolution, more filaments are identified, while for the lowest resolution, only two filaments are identified. The filaments seen at higher resolution are not nested within the structures identified at lower resolution but rather are new objects that have not been previously identified. At lower resolutions, these new objects cannot be identified because their morphology remains unresolved, causing them to appear as 'blobs' with a skeleton length shorter than the threshold. Specifically, we find 82 filaments at a resolution of 0.25 pc pix -1 in all snapshots across the three regions, 61 filaments at 0.5 pc pix -1 , and 46 filaments at 1 pc pix -1 .

The mass probability distribution functions (PDF) of the filaments identified at different resolutions are shown in Fig. 4 . Panel (a) indicates that there is a different distribution of mass among the filaments at different resolutions. For example, the 0.25 pc pix -1 resolution filaments exhibit numerous low-mass filaments, with two peaks at ∼ 200 M and ∼ 5000 M , whereas the 1 pc pix -1 resolution filaments have a mass distribution with only one peak at ∼ 7000 M . Meanwhile, the length PDF shown in Fig. 4 (b) reveals that the FILFINDER filaments distribute from ∼8 pc up to ∼90 pc, with 0.25 pc pix -1 resolution filaments exhibiting two peaks at 10 1.7 pc ( ∼60 pc) and at 10 1.1 pc ( ∼12 pc). 0.5 and 1 pc pix -1 resolution filaments also exhibit a two-peak distribution with similar peak v alues. Ho we ver, as the resolution becomes lower, the proportion of longer filaments increases while the proportion of shorter filaments decreases. The mass and length of filaments display a distinct bimodal distribution in the PDF, without any indications of hierarchy. Hierarchical structures should manifest as self-similar patterns characterised by a power-law PDF. This might be due to the resolution variation being insufficient to show the nested structures and the sample numbers being inadequate. The bimodal distribution is due to the limitations of our sample, which only showcases typical samples that the algorithm can identify within these three regions at specific resolutions.

The line mass of FILFINDER filaments varies o v er a large range, from ∼5 to ∼ 500 M pc -1 . Filaments with 0.25 and 0.5 pc pix -1 resolutions have a peak at ∼ 100 M pc -1 , while those with 1 pc pix -1 resolution have a higher peak at ∼ 200 M pc -1 . Unlike the other quantities, the line mass distribution does not show a clear bimodal profile, indicating a narrower range of line mass. This may suggest that underlying mechanisms are constraining the distribution.

Fig. 5 presents the mass and length properties of all our samples (including different resolutions) in log-space.

We follow the approach of Hacar et al. ( 2023 ) to compare the distribution of our filament sample to equilibrium filaments. We first compare our results with the critical line mass, which is analogous to the isothermal Jeans mass of molecular clouds (e.g. Klessen & Burkert 2000, 2001 ;Larson 2005 ). The grey solid line in Fig. 5 represents the critical line mass of a hydrostatic, isothermal cylinder (Stod ólkiewicz 1963 ;Ostriker 1964 ):

where c s is the sound speed, which is given by

Here k B is the Boltzmann constant, m p is the proton mass, the molecular weight µ = 2.8, and T is the gas kinetic temperature. Filaments with a line mass m > m crit (supercritical) are unstable and undergo collapse due to their own gravity, while those with m < m crit (subcritical) can remain in hydrostatic equilibrium.

The red solid line represents the virial line mass m vir , which is given by

It is important to note that this virial mass includes an extra non-thermal pressure contribution based on the emission line-width scaling seen in Milky Way filaments (see fig. 2 in Hacar et al. 2023 ). The thermal-only case is shown by the grey line, and thus, between the grey and red lines, the filament is only stable if the observed supersonic velocity dispersion is interpreted as a supportive pressure force. Below the red line, gravity dominates, resulting in collapse and fragmentation. Abo v e this line, a combination of thermal and nonthermal pressure dominates, potentially attributable to turbulence.

In our case, as shown in Fig. 5 , most of the filaments are distributed below the thermally critical line mass.This might be because the model for thermally critical line mass is based on several simple assumptions, including that the filament is in hydrostatic equilibrium, isolated, and isothermally. Ho we ver, filaments are generally more complex than the ideal model, which can result in differences from the theoretical predictions.

The yellow solid line indicates an upper limit for the filament distribution to remain molecular, driven by self-shielding against where N shield = 10 21 cm -2 (Van Dishoeck & Black 1988 ). It is based on the assumption that the central density is 10 times the average density and the filament is in hydrostatic equilibrium. Ho we ver, the upper limit relation is only a rough estimate because materials in the surrounding environment also contribute to the shielding effect. Most of the filaments in our sample lie in the regime between the yellow and red lines (with only one exception, which is slightly below the virial mass line). This is consistent with the observed distribution MNRAS 528, 6370-6387 (2024) of filaments in the Milky Way (Hacar et al. 2023 ), and suggests that the distribution of our simulated filaments aligns well with the expectations from these simple models. Ho we ver, this does not imply that the filaments are static objects because they can still grow (by moving rightward), shrink (by moving leftward), and fragment to form stars as the evolutionary paths shown in Fig. 5 . We will discuss this evolutionary path in more detail in Section 3.4 . Hacar et al. ( 2023 ) show that molecular filaments follow a continuous distribution across a large range of mass and lengths, corresponding to an approximate scaling relation L ∝ M 0.5 . Over this range, filaments possess different stability and dynamical properties. For example, filaments that exceed the expected m vir for a hydrostatic filament would be more susceptible to collapse and fragmentation. On the other hand, many of the shorter filaments in Hacar et al. ( 2023 ) display a sub-critical line mass, suggesting these filaments might be transient structures or evolving (Hacar et al. 2023 ). In this case, the distribution of filaments in the M-L diagram might be able to trace the evolving status of filaments.

Here, we use the highest resolution (0.25 pc pix -1 ) when tracking the evolutionary path, as this identifies the greatest number of filaments. We examine all of these filaments and select only those that were clearly identifiable in at least three snapshots for our subsequent analysis of their evolution, which we refer to as tracked-filaments. This is because the small structure undergoes significant morphology changes during evolution, particularly in complex surrounding gas environments, which make it difficult to track with an automatic algorithm. The properties of these tracked-filaments are listed in Table 3 . The ID number of the tracked-filaments corresponds to the re gion the y belong to. F or e xample, 'RBF1' stands for Region B Filament 1. We finally track three filaments in Region A, one filament in Region B, and five filaments in Region C. Not all trackedfilaments can be identified in every snapshot because of the changing morphology. RAF2 and RCF4 are lost from our tracking before the final snapshot. RAF1 is the only filament that is lost and re-found due to the algorithm bias. After we trace the time evolution of filaments, we can now test how filaments evolve in the M-L diagram.

Fig. 5 shows that 7 of 9 filaments undergo a pronounced accretion process as they move from left to right on the M -L plot. This is the dominant evolutionary path for all the filaments in our samples. This process can result in either an increase in both mass and length or a significant increase in mass with little change in length. In both scenarios, the filamentary structures accrete material from their surroundings. Filaments RAF1, RBF1, RCF2, and RCF3 all exhibit a significant accretion path with continuous mass increases.

We consider the filament growth rate by comparing the line mass of two adjacent snapshots to assess the rate of line mass increase. The equation employed for calculating the growth rate is presented below:

where m is the line mass of the filament, m = M / L , while M refers to the mass of the filament and L represents its length. The time between the two snapshots is denoted by t with a unit of Myr. A positi ve filament gro wth rate might be due to accretion or length shortening. We remo v e the ne gativ e values caused by segmentation. Note.

(1) The ID number of the filament (2) Time snapshot (3) The total gas mass, including the gas portion within the sink mass (4) The length of the skeleton identified by FILFINDER (5)

The line mass (6) The mass of the stellar portion of the sinks, which is 33 per cent of the sink mass.

This is because the line mass of filaments is not uniform, so it might lead to a decrease in line mass after segmentation. Since we are only interested in studying the relationship between accretion and line mass, such drastic changes are not considered. RBF1 is the longest-tracked filament in our samples, and it is relatively independent, simple, and at an early stage of evolution with little star formation. Therefore, we chose RBF1 for further study to analyse the long-term evolution process of the filament (10 6 yr). Fig. 6 shows a simplified version of the M -L plot containing only RBF1. As shown in Table 3 , RBF1 increases in mass and length in the earlier part of the evolution and then drops in the later part. Meanwhile, the line mass increases the whole time.

In Fig. 7 (a), we present the evolution of the filament growth rate across the snapshots. For RBF1, the mean filament growth rate in dust is 94 M pc -1 Myr -1 . The filament growth rate ho v ers around a few dozen before the last snapshot and then increases to 300 M pc -1 Myr -1 in the final snapshot. Ho we ver, this final increase is driven more by a decrease in the length of the filament, as we will discuss in the next section.

The mass and length of the filaments do not uniformly increase at all times. Fragmentation often occurs when star formation begins during the latter stages of filament development, causing the filament to segment into multiple smaller pieces. RCF5 is one of a filament sample that undergoes such a segmentation process and then resumes growth, as is shown in Fig. 8 .

At snapshot 4, the filament RCF5 breaks into two parts and combines with RCF4 at snapshot 6. These processes lead to a 'U'turn trend in the M -L diagram (see the purple path 5 in Fig. 5 ). The segmentation of RCF5 splits the filament into an upper and bottom part, with the upper part having a mass of approximately 900 M and a length of 16 pc, while the bottom part is the main part with a mass of approximately 4600 M and a length of 18 pc. When gas is consumed by star formation in the filament, gaps appear along its length, and it is broken into smaller pieces by the filament-finding algorithm. In the last snapshot, the upper segment splits from the bottom and is no longer identified as a single structure by the algorithm. During this process, the gas gathers into dense cores, eventually transforming into stellar mass.

The skeleton determined by the algorithm is sensitive to the column density threshold in these regions. This sensitivity can significantly alter the skeleton path since the algorithm draws the path based on the surrounding pix els. F or instance, the choice of threshold is likely responsible for the segmentation and merging of RCF5.

In the M -L diagram, filaments tend to mo v e to the right due to accretion until they are disrupted by segmentation (shift towards the bottom-left), and then they resume moving rightward. For example, although RAF3 and RCF5 show a significant decrease in mass in the last snapshot, it is probably due to the segmentation mechanism changing the structure of the filament identified rather than a rapid change in the gas density.

As suggested by Hacar et al. ( 2023 ), the balance between fragmentation and accretion also determines whether a filament can survive. If fragmentation is too rapid, all the gas in the filament will be in the fragments rather than a continuous density structure, and so a filament will not be identified. Fig. 7 panel (b) show the accretion and fragmentation time of filament RBF1 calculated using the below relations (Hacar et al. 2023 ): 11)

The mean accretion time-scale is 1.8 Myrs, similar to that found in the Galactic Cold Core Herschel sample (1-2 Myrs; Rivera-Ingraham et al. 2017 ), while the mean fragmentation time-scale is 2.5 Myrs. The fragmentation time-scale τ frag is consistently larger than τ acc during the evolution period for Filament RBF1. The longer time-scale of fragmentation shows that filament RBF1 more rapidly accretes gas from its surroundings than it can fragment into smaller structures. This is also shown by the increasing line mass. This accounts for why filament RBF1 can be followed for o v er 1 Myr in our simulations, as the gas is replenished as the filament evolves despite the ongoing star formation.

Fig. 7 panel (c) presents the conversion time-scale for all the gas in the filament gas to be transformed into stars, derived by

where Ṁ Stellar = M Stellar / t is the stellar mass conversion rate.

The conversion time-scale increases from ∼ 1.7 to ∼ 6 Myrs, with an average value of 3.3 Myrs. Compared to the mean accretion timescale (1.8 Myrs), the conversion time-scale is similar at the early stage of the evolution and much larger at the end. The mass conversion rate Ṁ Stellar also supports this conclusion, as it decreases from 2800 to 1000 M Myr -1 o v er the evolution, representing a decline of nearly 2/3. In this specific case, the reduction in star formation is due to the densest part of the filament being converted into stars and the remainder of the filament being less dense and, consequently, less actively star forming.

The final mechanism we see is 'dispersal', where rather than the filament being consumed by internal star formation, it is disrupted due to external forces. Sources of external pressure may come from mechanical and radiation feedback or compression by the galactic spiral potential, acting at different scales. When an external force interacts with a filament, it may undergo a morphological change resulting in a shorter length due to a combination of compression of the dense gas and removal of the more diffuse gas. This is e x emplified well by Filament RBF1, which, as seen in Fig. 6 , undergoes a sharp decrease in length at the end of its evolution while retaining most of its mass. Alternatively, filaments may also be stretched by shear Smith et al. ( 2016 ) or by differential rotation (Smith, Glo v er & Klessen 2014 ; Duarte-Cabral & Dobbs 2017 ). The surrounding subplots in MNRAS 528, 6370-6387 (2024) Fig. 6 show the morphology change of REF1, with the contour of H 2 column density. It can be clearly seen that the upper gas structure is moving downward, causing a dispersal with length shortening.

To determine the reason for this, we first inspect the location of RBF1 in the galaxy. RBF1 is located in between two spiral arms, as shown in Fig. 9 . We then examine the supernovae locations from the Cloud Factory model. Ho we ver, we find no supernova explosion within a distance of 50 pc from RBF1 during the time frame under consideration. Considering RBF1 is not in the spiral arm, we suggest the dispersal of RBF1 is probably due to the differential rotation of gas moving between the arms. Such dispersal can be visually observed in Fig. 2 . In the last snapshot 15, the bottom part of RBF1, which shows active star formation, appears denser than in the first snapshot 3. This indicates that diffuse gas from the upper part has mo v ed to the bottom part during this time period.

Ho we ver, we only find one filament associated with this mechanism. This is likely due to the fact that only Region B is located in the interarm region, and we stop tracing the evolution before many possible supernova explosions could occur after active star formation.

We therefore propose that the evolutionary progression of filaments in our models in the M -L phase diagram can be reduced to three physical processes, which are illustrated in Fig. 10 . Fig. 10 (a) shows the directionality of the evolutionary path, while the corresponding four physical mechanisms are displayed in the sub-figures below in (a). Fig. 10 (b) displays accretion, which is characterized by an increase in filament mass o v er time, while the length may either increase or remain steady. The entire gas filamentary structure is undergoing accretion from its surrounding material, leading to significant increases in mass.

The remaining two mechanisms result in a decrease in mass and length. Segmentation of filaments due to the conversion of gas into stars, as shown in Fig. 10 (c), causes both a decrease in filament mass and length. Meanwhile, Fig. 10 (d) illustrates dispersal, characterized by a significant decrease in length with little change in mass. This mechanism is attributable to external forces acting on the filament, such as shocks generated by supernovae or shear resulting from larger-scale galactic motions.

In our samples, 7 out of 9 tracked-filaments exhibit significant accretion paths, except RCF1 and RCF4. Both RCF1 and RCF4 are activ e star-forming re gions and hav e only evolv ed for a short period of time. Three filaments (RAF3, RCF4, and RCF5) experience a segmentation process, while only one filament, REF1, undergoes a dispersal process. These results indicate that accretion is the dominant mechanism in filament e volution. Ho we ver, it is important to note that these mechanisms are merely observed in our samples and do not represent all possible scenarios. All these mechanisms control the distribution of filaments in the M -L diagram, maintaining the hierarchical nature of filaments in a large spread range of scale.

We want to emphasize that we do not incorporate internal dispersal mechanisms in our study because our simulation lacks early stellar feedback components such as jets, outflows, and photoionization. The absence of these factors in our paper should not be interpreted as an indication that filament dispersal through these mechanisms is unimportant; it simply means that we have not included them in our current models. We plan to investigate these internal processes in future research.

In Section 3.4 , a rightward shift in the M -L diagram for most filaments pointed to accretion as a probable cause. To validate this, we first focus on the line mass m , a parameter sensitive to accretion activity. Chira et al. ( 2018 ) find that the average line masses of filaments in their simulation al w ays increase with time. We also notice the same trend as shown in Fig. 11 . Most of our tracked-filaments samples (7 of 9) end up with a higher line mass compared to their first snapshot. The increasing line mass is reflected in the M -L diagram as filaments tend to shift towards the right, as shown in Fig. 5 . In the maximum case, the line mass can increase by o v er 100 M pc -1 Myr -1 compared to the initial line mass (RBF1).

Ho we v er, the observ ed increase in line mass may also be influenced by algorithmic identification, as structural variations in the same filament across different snapshots could affect the line mass calculations. To eliminate this potential confounding factor, an analysis of the raw simulation data was conducted for the filaments.

In real observation, the kinematics of molecular gas can be inferred MNRAS 528, 6370-6387 (2024) Figure 12. The figure provides a detailed visualization of a filament spine and its associated velocity components. In the top panel, the column density map of Region A at snapshot 3 is displayed, with a designated rectangle highlighting the area expanded upon in the bottom panels. The bottom panels present the velocity structures from the Cloud Factory simulation raw data within the selected region. The left-hand, centre, and right-hand panels represent the XY , XZ , and YZ plane projections, respectively. The background colourmap shows the standardized velocity components of V z , V y , and V x . The colourbar representing the standardization velocity component is displayed in the centre of the picture. The purple contours indicate the dense gas structures. The black arrows represent the velocity vector in its respective projection. All three bottom panels show an ongoing accretion process.

In Fig. 12 , we present the velocity structure of a segment of the spine of a filament in Region A snapshot 3 to further inspect the ongoing accretion process within the filament. Across the three different projected view angles, the velocity components exhibit symmetrical distribution, indicating gas convergence towards the center of the filament from all directions. It should be noted that the velocities used here are derived directly from simulation raw data, and thus cannot be simply compared with V LSR from real observations.

In addition, although there is a broadly symmetrical distribution of velocities, there are also complex patterns present. The XZ plane projection reveals that the filament is not an ideal cylinder object, but instead exhibits a curvature. The projected velocity components (shown by black vectors) do not al w ays follow a perpendicular trend to the filament spine, but might also be parallel to it. In fact, the v elocity v ectors appear to be oriented towards the densest re gion. This indicates that the gas accretion is shifting towards the dense core, either from the surroundings or along the filament itself.

One might expect accretion to have a more substantial influence in regions where the line mass is high. Hacar et al. ( 2023 ), for example, proposes that the normalization of filaments in the M -L plot is set by accretion with high line masses seen in high accretion environments. Fig. 13 illustrates the relationship between the filament growth rate and the line mass and shows there is a potential relationship between the two in our models. We performed a linear regression on the logarithmic transformation of m and ṁ to obtain the fit parameters. The standard error of the fit was calculated from the mean squared error, and the 99 per cent confidence intervals were derived by exponentiating the fitted values ± 2.58 times the standard error, providing an error band around our linear model. As discussed in Section 3.4.1 , a positive filament growth rate can relate to both mass increasing and/or length shortening. To further investigate the role of accretion, we plotted a bar for each point, corresponding to ṁ multiplied by a length correction factor L m + 1 L , to demonstrate the impact of length variation on ṁ . We observe that the corrected values tend to be closer to the fitted linear relationship. This suggests that

In the preceding sections, we have delved into the evolution pathways of filaments in the M -L diagram and discussed the probable role of accretion. Ho we ver, it is important to note that filament identification largely depends on 2D projections of inherently 3D structures. Consequently, this projection effect may introduce systematic errors in our understanding of filament properties and their evolution. In this section, we aim to discuss the implications of these projection effects on our analysis.

In our previous radiative transfer realisations, the detector's lineof-sight is aligned along the Z -axis. As a result, the dust continuum images represent a projection on to the X -Y plane. We now simulate projections on to the X -Z plane by selecting the Y -axis as the new lineof-sight direction with POLARIS and applying the same analytical procedures. Fig. 14 displays the evolutionary paths of filaments in different projections in Region C and shows that the length and mass of the filaments can vary due to projection. However, despite these variations, a similar trend in the evolutionary paths was observed for most filaments, except for RCF5. This consistency is largely due to the dispersed distribution of filaments in Region C, allowing the algorithm to identify individual filaments with ease. For RCF5, the example we use for the segmentation process in Section 3.4 , the filament also became segmented in the X-Z plane.

Ho we ver, RCF5 appears as multiple fragmented segments in the X-Z plane, making it difficult to correlate these parts with the more unified structure observed in the X -Y plane. Filaments in Region A are also difficult to correlate to individual structures across different projection planes. This is primarily because the molecular gas environment in Region A forms a continuous molecular complex. Areas of relati vely lo w-density dif fuse gas interconnect dense regions, significantly affecting how the algorithm identifies distinct filament structures in different projections. In other words, filaments of Region A identified in different projections are barely the same objects. These findings underscore the critical need to consider projection effects and the limitations of our current identification algorithm when analysing filamentary structures in diverse environments.

The Mass-Radius relationship of molecular clouds has been extensively studied and is found to be linear, as evidenced by both observations and simulations (Larson 1981 ;Stutzki et al. 1998 ;Simon et al. 2001 ;Kauffmann et al. 2010 ;Lombardi, Alves & Lada 2010 ;Roman-Duval et al. 2010 ;Qian, Li & Goldsmith 2012 ). The origin of this linear relationship may be attributed to the hierarchical structure induced by turbulence and fragmentation (Kauffmann et al. 2010 ). Similarly, based on e xtensiv e observational data, Hacar et al. ( 2023) hav e pro vided the mass-length relationship of filaments with approximately L ∝ M 0.5 and explored their hierarchical structure.

There are generally two mechanisms for the hierarchical structure of filaments. One mode of filament formation is the top-down process, where larger filaments fragment to form smaller filaments (Hacar et al. 2013 ). The other mode is the bottom-up process, where small filaments are assembled onto pre-existing larger filaments under the influence of gravity to make them larger (Smith, Glo v er & Klessen 2014 ). It is worth noting that both mechanisms may coexist. As we mentioned in Section 3.2 , the bimodal distribution of filament sizes is clearly seen in the PDF distribution in Fig. 4 and can also be observed in the M-L diagram in Fig. 5 . The grey triangles in Fig. 5 can be clearly divided into two parts, one with a mass of up to 10 4 M and the other with a size of around 10 pc and a mass ranging from 10 2 to 10 3 M . Ho we ver, due to the resolution limitation, we are not able to identify those smallest fibre-like structures with a size of around 0.1 pc. In most cases, especially in Regions B and C, the filaments are distinct and separate. This means we cannot link these filaments to the two mechanisms because we do not find any cases where filaments are within filaments.

Fig. 15 shows the tracked-filaments paths in different resolutions in Regions A and C. The filaments in Region C are clearly separated and distinct, pro viding a relativ ely simple environment for comparing the impact of resolution on their identification. Generally, the filaments identified at the lowest resolution have the highest masses, while those identified at the highest resolution have lower masses. The evolutionary paths show similarities across different resolutions, except for filament RCF1, where there is a large change in the mask. Fig. 16 illustrates the impact of resolution on the observed properties of filaments ( M gas and m ). As the resolution changes from 0.25 to 1 pc pix -1 , there is a noticeable trend of increasing mean values for both M gas and m , with the predominant distribution of the violin plot contours shifting upward. This trend is consistent MNRAS 528, 6370-6387 (2024) with the findings presented in Fig. 15 , corroborating the influence of spatial resolution on the characterization of filaments. In conclusion, the small structures are difficult for the algorithm to discern at low resolutions. This explains how it is challenging to identify filaments with a mass distribution spanning several orders of magnitude in a single observation. Still, when the statistics of all the filaments from dif ferent observ ations are combined and placed in the log-axes M -L diagram, these filaments exhibit a linear relationship. We conducted a linear regression on the log-transformed mass and length. In our case, we identified a fitted relation of L ∝ M 0.45 as shown in Fig. 5 , which is indistinguishable from observations. The four mechanisms in Section 3.4.4 maintain this linear relation with a dispersion that falls within the upper and lower limits predicted by the theoretical models.

We examined three molecular cloud regions in the Cloud Factory galactic-scale ISM suite, using synthetic dust observations using the radiative transfer code POLARIS . We observed the evolution of filaments, which are identified in the observational plane, with FILFINDER .

We identified 189 filaments across the three regions at different resolutions in all snapshots and analysed their properties. We plotted their mass and length on an M -L plot and found a scaling relation L ∝ M 0.45 similar to that seen in observations. The distribution of our simulated filaments aligns well with that presented by Hacar et al. ( 2023 ). Comparisons with equilibrium models (virial line mass m vir and critical line mass m crit ) suggest that most filaments are supercritical when considering only thermal support, making them susceptible to collapse and fragmentation. Ho we ver, when turbulence is included, these filaments are found to be abo v e the virial line, indicating a more complex stability balance.

We identified 9 filaments that could be tracked o v er an e xtended time period and studied their evolution process in a galactic-scale environment. Filaments are dynamic objects that can grow and fragment to form stars, as in the evolutionary paths shown in Fig. 5 . By studying their evolution, we found three main filament evolution mechanisms that dictate the evolutionary progression of filaments in the M -L phase diagram. Accretion leads to a mass increase, while length might increase or remain steady. Segmentation causes significant mass and length decreases, while dispersal leads to reduced length with little change in mass.

Most tracked filaments in our samples exhibit paths from smaller masses to larger masses in the M -L plot, indicating the dominant role of accretion in filament evolution. We also found a potential linear trend between filament growth rate and line mass, suggesting gravity-dominated accretion.

These three mechanisms determine the distribution of filaments in the M -L diagram. Ho we ver, we ackno wledge that internal dispersal mechanisms, such as stellar feedback from jets, outflows, and photoionization, are not included in our current simulation. Magnetic fields, which might help keep filaments more coherent and linear, are also not included. We plan to explore their impact in future research.

We investigate the effects of projection on filament evolution in the M -L diagram and find that it introduces systematic errors in the identification and analysis of filamentary structures. The general trend of the evolutionary paths of individual filaments in a simple environment is preserved in different projections. However, measurements like length and mass may vary. In more complex regions, projection affects the accuracy of our current identification algorithms, meaning the structures identified at different viewing angles do not correspond to each other. These findings emphasize the need for caution when interpreting results based solely on 2D projections.

Finally, we inspect the hierarchical nature of the filaments. We compare the filaments identified at different resolutions. The higher resolution images reveal more filaments with smaller masses, shorter lengths, and lower line masses. However, due to the resolution limit and algorithmic limitations, we could not find any instances of filaments nested within their parental filaments.

We gratefully acknowledge helpful discussions with Andres Izquierdo, Jouni Kainulainen, Gina Panopoulou, Fabian Heitsch, and Xuepeng Chen. JF acknowledges support of the National Natural Science Foundation of China (grant No. 12041305) and the CAS International Cooperation Program (grant No. 114332KYSB20190009), and grants from the STFC and CSC 201904910935, without which, this work would not have been possible. RJS gratefully acknowledges an STFC Ernest Rutherford fellowship (grant ST/N00485X/1). AH acknowledges support and funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (grant agreement No. 851435). SEC acknowledges support from the National Science Foundation under grant No. AST-2106607. DS acknowledges support of the Bonn-Cologne Graduate School, which is funded through the German Excellence Initiative as well as funding by the Deutsche Forschungsgemeinschaft (DFG) via the Collaborative Research Center SFB 956 'Conditions and Impact of Star Formation' (subproject C6) and the SFB 1601 'Habitats of massive stars across cosmic time' (subprojects B1 and B4). Furthermore, DS received funding from the programme 'Profilbildung 2020', an initiative of the Ministry of Culture and Science of the State of Northrhine Westphalia.

This work used the DiRAC@Durham facility managed by the Institute for Computational Cosmology on behalf of the STFC DiRAC HPC Facility ( www.dirac.ac.uk). The equipment was funded by BEIS capital funding via STFC capital grants ST/P002293/1, ST/R002371/1, and ST/S002502/1, Durham University and STFC operations grant ST/R000832/1. DiRAC is part of the National e-Infrastructure.Author: In order to validate your funding information prior to publication, please check and confirm whether the name of the funding body given in your manuscript is complete and correct. If any edits are required please mark them on the text. Please also e xpand an y acron yms used in this section. If multiple grants are cited, please ensure the text of your funding statement clearly indicates which grant applies to which funding body.Author: In order to validate your funding information prior to publication, please check and confirm whether the name of the funding body given in your manuscript is complete and correct. If any edits are required please mark them on the text. Please also expand any acronyms used in this section. If multiple grants are cited, please ensure the text of your funding statement clearly indicates which grant applies to which funding body.

Data supporting this study is available in the Zenodo repository at the following link: https://zenodo.org/r ecor ds/10117334 . The data set encompasses the column density maps derived from the synthetic observations for all the snapshots across the three resolutions in three regions. Further data products are available upon request.