Decaying turbulence in molecular clouds: how does it affect filament networks and star formation?

The fragmentation of gas to form stars in molecular clouds is intrinsically linked to the turbulence within them. These internal motions are set at the birth of the cloud and may vary with galactic environment and as the cloud evolves. In this paper, we introduce a new suite of 15 high-resolution 3D molecular cloud simulations using the moving mesh code AREPO to investigate the role of different decaying turbulent modes (mixed, compressive and solenoidal) and virial ratios on the evolution of a $10^4\mathrm{M}_{\odot}$ molecular cloud. We find that diffuse regions maintain a strong relic of the initial turbulent mode, whereas the initial gravitational potential dominates dense regions. Solenoidal seeded models thus give rise to a diffuse cloud with filament-like morphology, and an excess of brown dwarf mass fragments. Compressive seeded models have an early onset of star-formation, centrally condensed morphologies and a higher accretion rate, along with overbound clouds. 3D filaments identified using DisPerSE and analyzed through a new Python toolkit we develop and make publicly available with this work called FIESTA, show no clear trend in lengths, masses and densities between initial turbulent modes. Overbound clouds, however, produce more filaments and thus have more mass in filaments. The hubs formed by converging filaments are found to favour star-formation, with surprisingly similar mass distributions independent of the number of filaments connecting the hub.


INTRODUCTION
The study of star formation spans several centuries, but the basic physical processes governing the earliest stages still remain poorly understood.The current paradigm of present-day star-formation suggests that the cold, dense phase of the interstellar medium (ISM), in the form of molecular clouds, undergo runaway gravitational collapse and birth stars in the process.As early as 1907, it was hypothesised that these clouds were not homogeneous, but contained substructure in the form of filaments connecting dense cores (Barnard 1907).Over the decades, the literature has vastly expanded to include modelling and characterisation of these sub-structures -from protostars, to dense cores and clumps, to parsec-scale filaments, to several hundred parsec-scale molecular clouds.In particular, the discovery of the omni-presence of these filaments in the ISM by Herschel (e.g.Men'shchikov et al. 2010) has brought a flurry of excitement and research on connecting the ISM geometry to the complex multi-scale physics of star formation, especially that of massive stars (see Hacar et al. 2022, and references therein).
The typical characteristics of molecular clouds are generally understood: they are dense (≳ 100 cm −3 ), and largely composed of ★ E-mail: jvd29@cam.ac.uk † E-mail: zoe.faes@esa.int‡ E-mail: rjs22@st-andrews.ac.ukH 2 , HI, H + , CO, and other molecules of C, H, O and N (Bergin et al. 1997).Efficient line cooling by CO and a high population of absorptive dust grains impart a characteristically low temperature of 10 K -20 K (see, e.g., reviews Bergin & Tafalla 2007;Dobbs et al. 2014) -conditions that are perfect for gravitational collapse of the cloud.The collapse is, however, not linear; the clouds are permeated by supersonic, compressible turbulence, as demonstrated in the seminal work by Larson (1981), which play a crucial role in the evolution of the cloud (see, e.g., reviews Mac Low & Klessen 2004;McKee & Ostriker 2007).Fully characterising the collapse of molecular clouds is a hard task due to the complexity of the scales at which they occur -from kiloparsec-scale in the galactic interstellar medium, down to less than 0.01 parsec for single proto-stellar objects.Idealised simulations reveal that the gas develops a lognormal density probability function (PDF) which relates to the properties of the turbulence and possibly the subsequent stellar initial mass function (IMF) (e.g.Klessen & Burkert 2000;Krumholz & McKee 2005;Burkhart 2018;Federrath & Klessen 2013).Generally, turbulence is seen as playing a dual role: it supports the cloud against global collapse, and also accelerates it through local compressions which give rise to a complex filamentary networks of dense gas.The mass, size and density of filaments can span several orders of magnitude, with lengths ranging from of 0.01 pc to several hundreds of pc, and masses from 0.01 M ⊙ -10 5 M ⊙ (Hacar et al. 2022).They generally have central number densities higher than 1000 cm −3 .These filaments fragment into denser clumps and cores, where stars and stellar clusters are born (Hartmann 2002;André et al. 2010;Myers 2011;Hacar et al. 2013;Könyves et al. 2015) like beads on a string.This makes filaments a key ingredient in any star formation theory, especially in the post-Herschel era.Numerical simulations have shown that contracting sheet-like clouds accrete material onto filaments with velocity flows perpendicular to them; and once channelled onto the filament, material flows longitudinally towards cores (e.g.Smith et al. 2011;Gómez & Vázquez-Semadeni 2014;Smith et al. 2016).This is not just inferred using simulations, but also seen in observations (Kirk et al. 2013;Peretto et al. 2013;Hacar et al. 2017;Liu et al. 2021).Hence, filaments form the bridge for collapsing gas to be transferred from large scales into dense cores where stars may form.Modern theoretical works have explored the many different mechanisms that give rise to these filaments (Hennebelle 2013;Federrath 2016;Smith et al. 2016;Abe et al. 2021).
Despite their importance, there is no universal, physically motivated definition for a filament in both observations and numerical simulations.Early pioneering works used an idealised filament geometry of an infinite isothermal cylinder (Ostriker 1964a,b;Larson 1985;Inutsuka & Miyama 1992), and found there existed a criticalline mass above which it becomes gravitationally unstable and fragments into cores.However, modern studies find filaments to be far from ideal, with complex substructure.They are dynamic objects affected by magnetic fields (Nagasawa 1987;Nakamura et al. 1993;Adam et al. 2016;Ade et al. 2016;Li & Klein 2019), external pressures (Fiege & Pudritz 2000), turbulence (Hartmann 2002), filament shape (Gritschneder et al. 2017), and evolve through interactions with their environment (Smith et al. 2014b).There has been significant discourse on the possibility of a universal filament width of 0.1 pc, first suggested by Arzoumanian et al. (2011).This remains contested with some suggesting that it may be an effect of the methodology used to derive the width (Smith et al. 2014b) or an effect of distance and resolution (Panopoulou et al. 2022) .For a full discussion on filaments, see the review by Hacar et al. (2022).
The evolution of filaments is inextricably linked to that of massive stars.Massive stars play a key role in galaxy evolution as they regulate the gas-star cycle through wind and radiation, and eventual supernovae explosions (Mac Low & Klessen 2004).Studies of massive star formation are also important to constrain the IMF and characterise star-formation rates in galaxies.Despite this, a full understanding of their formation remains elusive.Two modern theories exist: core accretion and competitive accretion.In the core accretion model, massive self-gravitating pre-stellar cores that condense from a fragmenting clump (protostellar region) give rise to massive stars (Shu et al. 1987;McKee & Tan 2003).In the competitive accretion model, the formation of dense clumps occurs simultaneously with stellar clusters.Initially, this theory proposed that stars grew by Bondi-Hoyle type accretion (Zinnecker 1982;Bonnell et al. 2001a).However modern versions of the theory focus on accretion during the large scale collapse of a parent clump of gas whereby protostars at the centre of a potential accrete mass from large radii and become massive over time (Bonnell & Bate 2006;Smith et al. 2009).For a comprehensive review on this topic we refer to Tan et al. (2014) and references therein.Regardless of the mechanism, high-resolution observations from ALMA have hinted at a relation between complex networks of filaments and massive-star forming clumps (e.g.Peretto et al. 2013;Hacar et al. 2018;Motte et al. 2018;Chen et al. 2019;Dewangan et al. 2020).Naturally, intersections of filaments, or 'hubs', are the primary suspect for these clumps, having very high densities that can support massive star formation (Myers 2009;Peretto et al. 2012;Schneider et al. 2012;Kirk et al. 2013;Rayner et al. 2017).
The complex dynamics of turbulence and global collapse in molecular clouds, along with the formation and evolution of filamentary networks and hubs, form the perfect basis for numerical simulations -to constrain the recipe for star-formation.The use of magnetohydrodynamical (MHD) codes is widespread for this purpose (Fryxell et al. 2000;Springel 2010).The theory, and subsequently simulations, are divided into two competing scenarios on this topic: the turbulent scenario and the global hierarchical and chaotic collapse scenario.
In the turbulent scenario, molecular clouds are in virial equilibrium between self-gravity and strongly supersonic turbulence and evolve over many free-fall times, maintaining a low rate of star formation throughout (see reviews by Mac Low & Klessen 2004;Ballesteros-Paredes et al. 2007;McKee & Ostriker 2007;Bergin & Tafalla 2007;Gnedin et al. 2006, and references therein).In order for this scenario to play out, there must be a continuous source of turbulence that supports these clouds.Proposed mechanisms are (i) internal feedback from stellar sources (e.g.Norman & Silk 1980;Franco & Cox 1983;McKee 1989;Krumholz et al. 2006;Goldbaum et al. 2011), or (ii) external driving by converging flows in the warm neutral medium (e.g.Vázquez-Semadeni et al. 2007, 2010;Heitsch & Hartmann 2008;Klessen & Hennebelle 2010;Micic et al. 2013).However, the former fails to explain the non-thermal velocities observed in clouds without star-formation, while the latter drivers are too weak to support the cloud.Simulations of this scenario drive turbulence at large scales (e.g.Federrath et al. 2008;Heyer & Brunt 2004;Federrath 2016), akin to 'shaking' the simulation box periodically (called turbulent forcing).
The opposing view of global hierarchical collapse proposes that molecular clouds are transient, irregular features born from largescale HI flows (see review Vázquez-Semadeni et al. 2019, and references therein).Turbulence is generated by the converging flows and amplified by the collapse itself (so-called 'chaotic infall') giving rise to a moderately supersonic cloud.Most of the turbulent energy decays through a large number of weak shocks (Smith et al. 2000a,b), and the evolution of the cloud takes place in the form of a mass cascade -collapse occurs on all scales with smaller scales accreting from their parental structure (e.g.Hopkins 2013).Small, dense regions contract later, but on shorter timescales than large, diffuse ones.Overall, the global contraction onto local centres of collapse occurs on ∼ one free-fall timescale.The complex inner structure of molecular clouds is explained by the amplification of density fluctuations through self-gravity, magneto-hydrodynamic instabilities (Heitsch et al. 2005(Heitsch et al. , 2006(Heitsch et al. , 2007) ) and thermal instabilities (Hennebelle & Pérault 1999;Audit & Hennebelle 2005).As transient objects, the gravitational potential energy of the cloud starts low and increases monotonically until it catches up with the kinetic and thermal energies.Thus, rather than evolving in true virial equilibrium, the cloud evolves far from it until gravity dominates the energy balance, causing the near equipartition condition  grav ∼ 2 kin (Vázquez-Semadeni et al. 2007;Bonilla-Barroso et al. 2022;Vázquez-Semadeni et al. 2019).Due to the successes of this model, it forms the basis of our numerical simulations.This scenario can be simulated by feeding a turbulent velocity field to the molecular cloud, and letting selfgravity take over the evolution of the cloud leading to a 'decay' in the initial conditions.
Depending on the source, the mode of the turbulence in clouds can be compressive (curl-free) or solenoidal (divergence-free).Dynamical mechanisms (galactic spiral shocks and converging flows) or stellar feedback (supernovae and radiation-pressure-driven shells) excite compressive modes in clouds, while high shear environments (at galactic centres, for example) or magneto-rotational instabilities induce solenoidal modes (Federrath et al. 2016).
In this work, through a suite of 15 gravoturbulent simulations, we explore how freely decaying turbulence seeded with different initial turbulent modes (mixed, compressive and solenoidal) and dynamical states (characterised by their virial ratios) affect the evolution of a molecular cloud.The paper is broadly divided into two main parts.In Section 2, we outline the numerical model and methodology used to identify and characterise filamentary structures in the evolved clouds.In Section 3, we present and discuss an array of results ranging from cloud morphologies and density distributions, to star-formation statistics, IMFs and filament statistics.Finally, we outline some caveats and future work in Section 4 and summarise our conclusions in Section 5.

Numerical model
We perform our simulations on the moving mesh code Arepo (Springel 2010;Weinberger et al. 2020), a multi-physics cosmological (M)HD code that uses a finite volume approach to solve the MHD equations, and performs calculations on an unstructured mesh formed by the Voronoi tessellation of a set of points that move freely with the fluid flow.This results in a Galilean-invariant solution accurate to second order in all four dimensions.The freedom of geometry allows for complex dynamics to be simulated.The smooth fluid flow automatically adjusts resolution at localised clusters (i.e., higher resolution at high density regions), while the ability to perform adaptive mesh refinement on a grid captures the effects of shocks and fluid instabilities well.Furthermore, all calculations are made on a hierarchical time-step that is optimised for parallelisation.In our setup, the refinement of the Voronoi mesh is set such that it applies Jeans refinement criteria of 8 cells per local Jeans length to regions with number density  > 100 cm −3 , to avoid artificial fragmentation (Truelove et al. 1997).Furthermore, to avoid discontinuous jumps in the cell size, adjacent cells may not differ in their radii by a factor larger than two throughout the simulation volume.The self-gravity of the gas is calculated using the standard Arepo gravitational tree (Springel 2010).The effect of magnetic fields and stellar feedback processes such as supernovae explosions, stellar winds, jets and outflows are absent in our models.
The chemical network in the simulations consist of five species: H 2 , HI, H + , CO and C + , with the inclusion of dust grains.The hydrogen network is based on Glover & Low (2007a,b) and the CO network is that of Nelson & Langer (1997).The treatment includes formation and destruction of molecules through collisions and chemical reactions, and ionisation by interstellar radiation field (ISRF) and cosmic rays.The strength of the ultraviolet ISRF field is set to  0 = 1.7 Habing units, as calculated by Draine (1978) assuming same strength and spectrum of the field as the solar neighbourhood.The cosmic ray ionisation rate is set to  H = 3 × 10 −17 s −1 for atomic hydrogen, as calculated using observations of 7 protostars by van der Tak & van Dishoeck (2000), and twice the amount for molecular hydrogen.This model was first implemented in Arepo by Smith et al. (2014a).The effect of ISRF attenuation due to H 2 selfshielding, CO self-shielding, CO shielding by H 2 and dust absorption is modelled using the TreeCol algorithm developed by Clark et al. (2012).Heating and cooling of the gas from radiative processes is computed alongside the chemistry using the atomic and molecular cooling function described in Clark et al. (2019).

Star formation characterisation
We characterise star formation in the simulations through the framework of sinks, first introduced in Arepo by Greif et al. (2011).We use a modified version detailed in Treß et al. (2020), briefly discussed here.A sink is formed when a gas cell passes a density threshold of  sink (equivalently number density of  sink ) within an accretion radius  acc .Once a cell or group of cells passes this threshold, the code checks if the region is indeed gravitationally bound -i.e.potential energy  > 2 ( kinetic +  thermal ) -and has a negative divergence of velocity and acceleration to ensure gas is unambiguously collapsing.On passing the checks, the cells are replaced with collisionless 'sink' particles which inherit the combined mass and velocity of the parent gas cells and are located at the centre of mass of the region.The sink particles henceforth interact with their environment solely through gravity and may accrete mass (above the density threshold) from adjacent gas cells.We set  sink to ≈ 4 × 10 −16 g cm −3 (or  sink ≈ 1.7 × 10 8 cm −3 ), which is orders of magnitude higher than the average molecular cloud density, and close to those seen in protostellar regions (e.g.Whitworth & Ward-Thompson 2001).To ensure the entire instability is contained within the sink particle, the sink accretion radius  acc is set to 110% of the Jeans length at  sink , yielding a value of ≈ 200 AU.At this length scale, sinks may represent individual stellar systems (in case of low mass sinks) or small stellar sub-clusters (in case of high mass sinks).
Due to the absence of complex feedback processes like stellar winds and jets, sinks may be better thought of as 'accretion reservoirs' -they do not represent the mass in stars, but the mass available for the formation of stars.In actuality, a variable percentage of the sink mass would go towards the star formation, due to varying efficiencies depending on stellar environments and length scales.

Initial conditions
The initial molecular cloud used across all simulations is a spherical, uniform cloud of mass 10 4 M ⊙ and radius 8.83 pc, embedded centrally in a periodic bounding cubic box of side 64.82 pc, as shown in

Mixed
Compressive Solenoidal  the cloud and 10,000 to the surrounding ISM, which changes as the mesh evolves.The cloud has a number density of ≈ 120 cm −3 and temperature of 27 K.The surrounding ISM environment has temperature of more than 7500 K and a much lower number density of 10 −5 -10 −4 cm −3 .The global free-fall time of the cloud is  ff ≈ 6.5 Myrs.
The cloud is initially composed entirely of atomic hydrogen with a helium abundance (by number) of 0.1, H + abundance of 0.01, carbon abundance of 1.4 × 10 −4 , oxygen abundance of 3.2 × 10 −4 and a metallic abundance of 1.0 × 10 −7 , all relative to hydrogen abundance.The dust-to-gas ratio is 1/100.The gas cells in the cloud are yet to be assigned velocities, which is discussed in the next section.
We generate homogeneous, normalised turbulent velocity fields using the code described in Girichidis et al. (2011).The algorithm creates a Fourier (k-space) grid and calculates the energies from an energy spectrum  (), which can be set to either Kolmogorov spectrum  () ∼  − 5 3 (Kolmogorov 1941) or Burger's turbulence  () ∼  −2 (Burgers 1939).We choose the latter since it is consistent with observations and simulations of supersonic, compressible turbulence in molecular clouds (Kritsuk et al. 2007;Schmidt et al. 2009;Federrath 2013, albeit centred around turbulent forcing).The code then generates amplitudes and random phases for each k-space grid point.The amplitude is just the energy  (), while the phase is produced through random number generation (RNG) on a uniform distribution from a specified seed parameter.This gives û(k), the Fourier transform of the turbulent velocity field u(x).The field is optionally Helmholtz-decomposed into longitudinal (compressive) and transversal (solenoidal) fields, or left unchanged (mixed turbulence, which is 2 parts compressive to 1 part solenoidal) depending on the mode of turbulence chosen.A final inverse-Fourier transform gives the turbulent velocity field u(x), which is then normalised before outputting.The velocity field generated is defined on a Euclidean grid of size  3 where  is the number of grid points along an axis and can be varied as a parameter.
To simulate a decaying turbulent velocity field, we feed the output generated by this code to the molecular cloud at the start of the simulation.Since these are generated on a cubic grid of dimensions  3 and the geometry of Arepo is a Voronoi mesh, the velocities have to be interpolated and scaled to fit the spherical molecular cloud.The  3 velocity grid is first centred on the cloud, and all gas cells belonging to the molecular cloud in Arepo's ROI are assigned a velocity by linear interpolation between the cubic grid points.The computational intensity of interpolation sets the limit on how finely the velocities are resolved, i.e. how big  is.In our setup, we use  = 750 which corresponds to a length of ≈ 0.025 pc in the simulation.
The velocities   , once assigned to the cloud, are then scaled in accordance to the virial ratio: where  iterates over the cloud cells in the simulation,  is the total mass of the cloud,  is its radius, and  is the gravitational constant.For virial-type equipartition, the cloud must be  vir = 1.Since the gravitational potential energy of the molecular cloud in our simulation remains unchanged, i.e. the cloud mass is not varied, an under-bound cloud is one which has large velocities leading to  vir > 1, while an over-bound cloud is one which has a weak velocity field, i.e.  vir < 1. Figure 2 shows the   component of virial balanced clouds for different turbulent modes.
Having discussed the simulation technique, we now present the suite of fifteen different simulations with varying turbulent velocity modes and virial ratios.Nine of these simulations explore the statistics of different turbulent velocity modes at virial-type equipartiation by using different random seeds, while the other six explore different virial ratios for the three modes.Table 1 summarises the different initial velocity conditions.These simulations are all available publicly on Zenodo 1 .

Filament identification
As stated in the introduction, filaments are not strictly defined features.They can be diffuse, clumpy, dense, smooth, stub-like or long and tortuous, and often blend into their background.This makes it hard to assign a length, width or mass to them, whether in observational images or in simulations.However, their undebatable ubiquity at all scales in our Universe warrants a method of characterising them.In order to identify filamentary structures in our molecular clouds, we use the open-source tool DisPerSE (DIScrete PERsistent Structures Extractor; Sousbie 2011). 2 Although the code was originally written for the identification of web-like structures in cosmology (e.g.Kleiner et al. 2017;Laigle et al. 2018;Kraljic et al. 2018;Libeskind et al. 2018), the so-called 'cosmic web', it has found wide use in molecular cloud studies for identifying filaments (e.g.Arzoumanian et al. 2011;Federrath 2016;Smith et al. 2016Smith et al. , 2020)).The algorithm employs Morse theory to identify structures like voids, sheets, filaments and peaks in an -dimensional topological space.It first constructs a Morse-Smale complex of an input density distribution (or field), and identifies two key features in it: critical points -points where the gradient of the distribution is zero, i.e. ∇(, , ) = 0 -1 See Section 5 for data availability.and integral lines -curves that are tangent to the gradient at every point, and start and end at critical points.The critical points may be maxima, minima or saddle-points.The set of integral lines that join maxima to saddle points, called ascending 1-manifolds, delineate filaments in the input field, and form the 'skeleton' network of filaments.Each filament is defined by a number of sampling points it passes through, and ends in critical points.Since DisPerSE requires uniform cubic grids for identifying filaments, the density distribution of Arepo simulations within the ROI are cast from Voronoi geometry to a uniform cubic grid of size  3 using Arepo's in-built ray-tracing algorithm.Due to the computational intensity of running DisPerSE on 3D density cubes, we use  = 200 which sets an artificial resolution limit of ≈ 0.1 pc for our filament networks.Thus, in the hierarchical substructures of filaments, our access is limited to a specific lowest scale.This is not unlike observations, where the resolution limit is set by the telescope.
DisPerSE has several free parameters for identifying filaments.The persistence threshold parameter is the most crucial and largely determines what the identified network looks like.A pair of critical points, joined by a filament, is defined to have a numerical value called 'persistence'.This quantity has significance in defining the robustness of a topological feature, and is the main focus of the persistence theory.An intuitive explanation for this is offered in Sousbie (2011, Section 4), which we briefly summarise here.Consider a function used to model the density of a region.The relative 'height' of the adjacent local maxima and minima in the density field contains information about how tenuous or significant that topological feature is.Formally, this is known as the lifetime of the component in the excursion set.Small relative heights are akin to 'noise' (or small scale perturbations) in the local density field.The higher the threshold, the denser the regions being traced, while low values normally generate many spurious, short filaments.
Most often in observation-based literature, this value is provided in the context of filamentary structures in 2D column density maps, where it is set to ≈ 10 21 cm −2 (Arzoumanian et al. 2011;Federrath 2016).In 3D simulations, Chira et al. (2018) have used two values, 100 cm −3 and 5000 cm −3 , for their analysis.In this work, we pick  thresh = 8500 cm −3 by analysing the regions traced for various persistence thresholds, as done in Figure 3.A more detailed description of our methodology is provided in Appendix A, where we also quantify the errors introduced by this parameter.We also fix a threshold angle, set to 90 • , above which longer filaments are split into shorter pieces.

A Fiesta of filaments
In order to study the nature of these filaments in our simulations, we have developed a set of tools to work with outputs from Arepo and DisPerSE, and analyse them.This includes reading and writing data for both softwares, visualisations of the simulations and filaments, and also functionality for characterising their properties such as lengths, masses, density profiles and more, which we describe in the next section.These tools are released as an open-source Python toolkit called FIlamEntary STructure Analysis (Fiesta), available on GitHub3 .Although primarily used for the aforementioned software, the algorithms developed can be applied more generally to study any filamentary networks defined on structured grids or unstructured meshes.

Filament characterisation
The filament network identified by DisPerSE only returns the coordinates of the filament spines.The length, width, density profile, and mass have to be determined through other means.The conversion of a 1D spine to a 3D tortuous cylinder is not a trivial one.
The sampling points for each filament are first fit by a 1D curve passing through them, by means of cubic B-spline interpolation4 , hence giving each filament a parameterised characteristic function that fully defines it:  () = (, , ).Here,  ∈ [0, 1] is the normalised distance of a point along the filament and (, , ) are the co-ordinates of the same point.The utility of this function is that it allows the filament to be sampled as finely or coarsely as needed.The length is then trivially calculated by integrating the Euclidean length along the characteristic function.
Determining the mass is a much harder task since filaments are not rigidly defined structures, but diffuse and merge with their background.To assign mass, we use an algorithm which relies on splitting the filament into longitudinal cylinders and querying the mass of all Arepo cells within these cylinders.The step-wise algorithm is as follows: consecutive sampling points along a filament are joined by vectors.Each vector is rotated and translated to the origin, along with all Arepo cells in a cubic region surrounding the base of the vector.This region selection avoids re-centering and rotation of all cells in the simulation box, which would be too computationally heavy.A cylindrical region is then defined with a fixed radius of  cyl around the vector.In the origin-centred vector frame of reference, this constrains the - values, while the  axis is constrained by the ends of the vector.Arepo cells within this cylinder are queried and the process is repeated for all sampling points in the filament.The total mass contained in these Arepo cells is defined as the mass of the filament, correcting for double-counting between cylinders.
The top panel of Figure 4 shows an example of the spine of the filament and the cylinders queried around it.

Caveats to the characterisation
Filament identification by DisPerSE is not perfect -the main issue lies with the 0.1 pc artificial resolution introduced by ray casting the finely resolved densities on the Voronoi mesh to a coarsely resolved cubic grid.The grid casting washes out local overdensitieshence, irregular clumps on length scales shorter than the resolution are not resolved properly.Filaments passing through or near these dense clumps may not be identified reliably by DisPerSE, and are noted to have inflated masses.To partially mitigate this, we apply two corrections: first, all filaments shorter than 0.1 pc, the minimum resolution of the cubic grid, are considered unreliable detections and discarded.Second, we apply a spine correction to the filament where the spine is re-adjusted using the resolved Voronoi mesh.This is done by querying all Arepo cells in a spherical region of radius  sph = 0.1 pc around each sampling point of the filament, and then shifting the point to either the densest point in the sphere ('max density' method) or its mean ('centre of density' method).The former method introduces more tortuosity in the filament due to strong local density fluctuations in the resolved mesh, as seen in Fig- ure 4. Hence, we adopt the 'center of density' method in this paper.This gives smoother filaments, and doesnt bias the resultant filament completely towards dense cores.The mode of correction noticeably affects the filament profiles as discussed in the next section.

On filament widths and profiles
As prefaced in the introduction, filament widths have been intensely studied in literature.A 'universal' width of 0.1 pc has been suggested by many (Arzoumanian et al. 2011;Malinen et al. 2012;Juvela et al. 2012;Kirk et al. 2015;Federrath 2016), but remains contested (Smith et al. 2014b;Panopoulou et al. 2022).In this work, we adopt a fixed  cyl = 0.1 pc radius for each filament since it is the resolution of the grid used for DisPerSE, hence twice the suggested universal width.For the purpose of this work, a fixed radius can shed light on how much mass is contained in and around these filaments as a comparative metric, rather than a definitive property of filaments in these turbulent clouds.
In observational data or 2D simulation projections, a Plummer fit is often used for filament density profiles (see, e.g., Federrath 2016) to calculate its width.We find that such a fit is not well suited for 3D filaments.Figure 4 shows an example of the density profile obtained for a randomly chosen filament, for both types of spine correction mentioned in section 2.6.2.Indeed, there is a filament profile in the density, but it is embedded within a background of both diffuse and clumpy regions scattered across the profile.The sharp peaks are likely due to dense cores (as suggested by the magnitude of their density) offset from the filament spine.The hierarchical nature of filaments also means that there are mini-filaments (or 'ribbons') that thread the main filament.Thus the profiles vary from filament to filament, with some less messy than others, making it hard to define a single fit that would work for all.A similar conclusion is reached in Smith et al. (2014b) where sub-filaments contribute to a non-smooth 3D density profile.In our case, there is an additional layer of subjectivity introduced by the form of spine correction, which affects the density profile as seen in Figure 4.

RESULTS
Throughout this work, we shall focus our analysis on four epochs of interest based on the total mass in sinks: 250 M ⊙ , 500 M ⊙ , 750 M ⊙ and 1000 M ⊙ .Since the total mass of the molecular cloud is 10 4 M ⊙ , these epochs correspond to 2.5%, 5%, 7.5% and 10% global starformation efficiency of the cloud (or more accurately, sink formation efficiency).As the cloud evolves differently in each simulation, these thresholds are achieved at different times.Hence, instead of comparing clouds at fixed time-steps, we compare them at similar stages in their evolutionary track.In the former approach, any differences in filamentary structures or density distributions between the simulations would be an artefact of the stage of cloud collapse; intuitively, one would expect compressive and overbound clouds to collapse faster than solenoidal and underbound clouds.The latter approach gives a truer picture of how the modes and virial ratios differ fundamentally once seeded into the cloud, and is closer to how an observer may interpret the cloud without knowledge of its absolute time evolution and only current star formation efficiency. Figure 5 shows the total sink mass evolution as a function of time, which we shall discuss in detail in Section 3.3.

Cloud morphology
We first investigate the morphology of the clouds.Figure 6 shows a 3D view of the dense regions in each cloud (> 1000 cm −3 ) at the final epoch, along with their 2D projections.We see that despite being at the same stage of collapse (i.e., same mass in sinks), solenoidal seeded turbulence tends to form elongated 'filament-like' clouds at the global scale, while compressive seeded turbulence forms 'centrally condensed' clouds.Looking at the evolution 5 , there is a definite rapid convergent flow in the compressive seeded clouds which forms an intermediate filament-like structure, before ending with a centrally condensed morphology.In comparison, solenoidal seeded clouds have a diffuse and slow collapse, decaying into compressive sub-modes through self-gravitation at much later times.There is a similar trend seen across underbound and overbound clouds, where the underbound models generate a more 'filament-like' morphology, while overbound models generate a more 'centrally condensed' morphology for the same initial turbulence seed.To quantify this observation better, we calculate the point cloud spread of all cells above 1000 cm −3 as the mass-weighted standard deviation: 5 see Zenodo repository, available in Section 5, for animations.where (  ,   ,   ) are the co-ordinates of cell ,   is the mass of the cell, and ( x * , ȳ * , z * ) is the mass-weighted centroid of all cells.
Table 2 shows this quantity for all simulations.There is a clear large dispersion of dense gas in the solenoidal seeded clouds, while mixed turbulence and compressive seeded turbulence forms more compact clouds.This is also seen in the underbound versus overbound cases, where S57U seems to have the largest spread of any simulation, almost twice as much as some compressive clouds.This naturally raises the question, are the solenoidal clouds really diffuse and less dense than others?-a question readily answered by looking at the probability density functions (PDFs) of these clouds.

Probability distribution functions
PDFs are an invaluable tool in probing the complex processes that affect cloud evolution -gravitational contraction, magnetic field effects (Molina et  give rise to a log-normal distribution while high density regions, characterised by gravitationally unstable, free-falling material and star-forming cores, give rise to a power-law distribution (Burkhart 2018).
Figure 7a shows the mass-weighted PDF for all simulations at the final epoch.The solenoidal seeded models show a significantly larger fraction of diffuse regions, followed by mixed and compressive seeded turbulence.To support this and further de-clutter the figure, we plot the cumulative distribution function (CDF) and group different turbulent seeds and different virial ratios for comparison in Figure 7b,c,d,e.As the cloud evolves between the two epochs being compared, the increase in diffuse regions is manifested in the CDF as a shift upwards for  < 100 cm −3 across all simulations.
In Figure 7b,c, diffuse regions dominate solenoidal seeded clouds.The different modes bring a very clear separation in the distribution during the early stages of the cloud, which is further deepened as the cloud evolves.Compressive seeded clouds are characterised by fewer low density regions while mixed seeded clouds lies between the two extremes.At the high density end, the difference between turbulent modes is unnotable due to local dynamics taking over the decayed turbulence.
In Figure 7d,e, comparing initial virial ratios, solenoidal seeded models already show a lag in their evolution at early stages.Across all modes of turbulence, we see a common trend at the low-density end: underbound > virial balanced > overbound, with the latter two following each other closely.S57U again remains an extreme case, showing a very large fraction of low density regions.We interpret this as the combination of the solenoidal turbulence mode and the initial underbound state of the cloud both impeding the local and global collapse, resulting in an excess of diffuse regions as compared to the other clouds.At the high density end, which is dominated by the strength of local (and to some extent global) gravitational collapse, we see dense regions preferentially forming in overbound clouds (see Figure 7a).
To test these claims more rigorously, we perform Kolmogorov-  apparent in the distributions for the underbound clouds, which are most likely to differ from the initial virial equilibrium and overbound cloud distributions.
It is clear, however, that the initial turbulence modes have a much stronger effect on the lower end of the density distributions.The most pronounced difference in the underlying distribution as measured by the KS test can be seen for the solenoidal underbound (S57U) distribution compared to all other CDFs.We also find a stark difference between C57O and all the other underbound cases, with the -values indicating that the underlying CDF for C57O is different from the other underbound distributions up to almost 3 for C57U and beyond for M57U and S57U, again suggesting a combined effect between the turbulence mode and initial virial ratio.
Thus, the low density end of distributions is governed by a decreased efficiency of collapse which can be brought about by two mechanisms: having dominant solenoidal modes of turbulence, or by a less gravitationally bound cloud.The effect of either is the same, and in combination can give an extremely diffuse cloud.We conclude that this low density regime, dominated by turbulence, maintains a strong relic of the initial turbulent mode and is weakly sensitive to the initial gravitational potential of the cloud.

Mass growth
We now look at the evolution of sink mass in the cloud.Figure 5 shows that compressive seeded turbulence is characterised by a very early onset of star formation, followed by mixed and lastly solenoidal seeded turbulence.In fact, by the time solenoidal seeded clouds reach 250 M ⊙ in sinks, the compressive clouds are ≈ 10 times higher in their stellar mass.This is consistent with previous studies that note similar star formation suppression in solenoidal clouds (Federrath & Klessen 2012;Padoan et al. 2014).The mixed turbulence seeded clouds lag by about 0.5 Myr (≈ 7.5% of free-fall time  ff ) compared to their compressive counterparts; in turn, the solenoidal seeded clouds lag by almost 1 Myr (≈ 15% of  ff ) to the mixed modes.The compressive seeded clouds also have a higher accretion rate onto sinks when comparing each RNG seed.In general, the effect of RNG seed on the total sink mass at a given timestep is very large, with variations of > 50% at 0.5  ff , akin to the findings of Jaffa et al. (2022).
Despite the mixed and compressive seeded underbound clouds commencing star formation earlier than their overbound counterparts, the rate of accretion is higher for overbound clouds and they catch up around the second epoch, when  sink = 500 M ⊙ .The same is true for the solenoidal case, although S57U shows a remarkably delayed collapse.Howard et al. (2016) studied variation in boundedness for a larger  vir range and 10 times more massive molecular cloud, and found a similar delay in cluster formation but higher star formation rates for overbound clouds.Their work also suggests, although at larger scales than ours, that initial gravitational boundedness has greater effect than radiative feedback during the first 5 Myrs, with the latter only affecting star-formation efficiency by a few percent.

Mass distribution and brown dwarf mass excess
We now look at the distribution of the sinks masses, i.e. the sink mass function, for the first and final epoch in Figure 9.Note that this is truncated to only the star-forming regime > 0.1 M ⊙ .An understanding of the origin of stellar IMF is crucial for the microphysics of galactic evolution and formation of planetary systems, and has been the subject of a great deal of numerical simulations.Some of the physical processes studied in this context include fragmentation (e.g.Zinnecker 1984;Larson 1985;Klessen et al. 1998), turbulence (Padoan & Nordlund 2002), accretion (e.g.Bonnell et al. 1997Bonnell et al. , 2001b;;Bate & Bonnell 2005), mergers (e.g.Bonnell & Bate 2002), or combination of many processes (Adams & Fatuzzo 1996).See Offner et al. (2014) for a review on stellar IMF studies.
At the first epoch, the number and distribution of sinks is seemingly independent of the turbulent mode.The median sink mass lies somewhere in the 0.3 -0.5 M ⊙ range, and and is consistently higher for overbound clouds than underbound clouds which have a more bottom-heavy distribution.This effect is less pronounced in the compressive case (C57U) where the competing effects of convergent flow and diffuse collapse balance each other.The solenoidal underbound cloud (S57U) shows an interesting feature -the cloud has an abundance of sink particles, roughly ×3 more than the other simulations, but the distribution is dominated by low-mass sinks.The variation between the turbulence modes in virial balance is not significant enough just yet.
At the final epoch, the initial distributions are maintained at the lower end by continual production of low mass sinks and at the higher end by the accretion of mass onto heavy sinks.As the clouds evolve, solenoidal seeded turbulence shows a characteristic abundance of sinks, dominating in number, but with a bottom-heavy distribution so that each individual sink mass is very low.The trend seen for S57U at the first epoch is now also seen across the other solenoidal modes at the final epoch.
To compare distributions, now also including sinks below the 0.1 M ⊙ threshold, we plot the mass function of the sinks in Figure 10a.As suggested by Figure 9, the characteristic mass for underbound clouds is lower than overbound clouds.Furthermore, across all simulations, we see that the mass functions roughly follow the Chabrier (2005) IMF model but is steeper at the low mass end, and shallower at the high mass end (the exact IMF shape is sensitive to binning at high masses.The shape and characteristic masses also vary largely with the turbulent seed.Most interestingly, an abundance of sinks in the solenoidal seeded clouds are of sub-stellar mass, with distributions peaking below 0.1 M ⊙ -equating to an excess at brown dwarf masses!This leads us to question -where are these low mass sinks forming?To investigate this, we plot the spatial distribution of sinks, split into low mass (< 0.1 M ⊙ ) and high mass (> 0.1 M ⊙ ) mass, in S26 in Figure 10b.The two mass regimes show similar spatial clustering, which is also seen in the other solenoidal simulations.This suggests that the brown dwarf mass sinks are not formed in unique regions that prevent accretion.
Mixed and compressive clouds attain massive sinks by efficient accretion and have a larger fraction of sinks at the high mass end, but dont show any noticeable difference between the two turbulent modes.At the low mass end, underbound clouds produce more low mass sinks compared to their overbound counterparts, while high mass sinks are favoured by the overbound clouds.This is in agreement with Howard et al. (2016) who see the same effect in cluster formation for a 10 6 M ⊙ cloud.In the case of S57O, the competing effects of the less efficient solenoidal turbulence and compressive gravitational energy counteract each other, producing a mass distribution similar to other mixed and compressive seeded simulations.
Figure 11 shows the number of sinks which are candidates for massive star formation, with masses greater than 16 M ⊙ (yielding >8 M ⊙ stars assuming 50% efficiency).The number of massive sinks is comparable for all simulations, with the exception of S57U as it is an extreme case, meaning that the most massive clumps form regardless of initial conditions.However, the relative fraction of massive sinks is greater for overbound clouds (except in case of C57O) compared to underbound clouds, and mixed/compressive seeded clouds compared to solenoidal seeded clouds, seen in both the mass function (Figure 10a) and count histogram (Figure 11).The abundance of sinks in the solenoidal seeded clouds likely means that they are able to coalesce some together, making it possible to have a comparable number of massive sinks.We do not include any radiative feedback in our simulations, instead choosing to focus only on the first few Myrs of collapse.The presence of protostellar heating could potentially change the abundance of brown dwarfs (see, e.g., Bate 2009a) or affect sink-formation in the immediate vicinity of massive protostars.Barring this caveat, we arrive at a similar result to the PDF investigation, where the low mass end of the distribution is controlled by the presence of turbulent solenoidal modes in the cloud or low gravitational potential energy.At the high mass end, the distribution is dominated by strong gravitational collapse and mixed/compressive turbulent seeds.

Observational aspects
On the side of observations, studies of the Central Molecular Zone (CMZ, Federrath et al. 2016;Rani et al. 2022) and the Orion B cloud complex (Orkisz et al. 2017) have found that the low star formation efficiency in these clouds is likely due to stronger solenoidal modes in turbulence.The CMZ, being at the Galactic centre, is certainly subject to a high shear environment exciting these solenoidal modes which decrease in fraction with Galactocentric distance.Our results here are consistent with this idea of star formation inhibition.
The presence of moderately large numbers of low-mass brown dwarfs in clusters like Trapezium, IC 348 and  Oph (Luhman et al. 2000) could perhaps be explained by the solenoidal mechanism.Although these clusters have not been studied in this context, the velocity PDFs of the Huygens region in Orion Nebula (home to the Trapezium cluster) show characteristics of solenoidal turbulence (Anorve-Zeferino 2019).Furthermore, studies of spatial distribution of sub-stellar objects in Taurus (Luhman 2006;Parker et al. 2011),  Oph (Parker et al. 2012) and IC 348 (Parker & Alves de Oliveira 2017) reveal that they are not clustered differently than higher mass stars, consistent with our simulations and hint at a common formation mechanism across the mass range.

Filament networks
We now look at the nature of filamentary structures produced in the simulations, identified and characterised using the methodology detailed in Section 2.5 and 2.6.The total number of filaments at each epoch of interest is shown in Figure 12.Over time, there is a marginal increase in the number of filaments in most simulations, however some show greater variability such as M90, S57 and S57O.Furthermore, although overbound and virial balanced clouds don't show significant differences, they produce more filaments than their underbound counterparts.

Filament length
Figure 13a shows the distribution of filament lengths in the simulations.The lower limit of 0.1 pc is artificially set by the resolution for DisPerSE, and tenuous filaments of length > 5 pc are shown as outliers.Remarkably, there is very little differentiation between the initial turbulent mode or virial ratio in the distributions.There is a marginal evolution over time towards shorter lengths which perhaps suggests fragmentation.One global distinction comes from the total length of the network in the last epoch which follows the trend: overbound > virial balanced > underbound, likely because the overbound and virial balanced clouds produce more filaments.Figure 14 highlights this trend clearly.Overall, the last epoch shows an interesting consistency in median length of ≈ 0.5 pc, and mean length of ≈ 0.7 -0.8 pc across all simulations.

Filament mass
Figure 13b shows the distribution of filament masses in the simulations.As with the length distributions, there is no clear distinction between the initial turbulent modes, and no clear evolution in their statistics, other than accretion as all simulations gather more mass onto filaments.The mean filament mass of ≈ 10 M ⊙ and median of ≈ 5 -6 M ⊙ remains consistent across simulations and time.The trend of overbound > virial balanced > underbound seen in the lengths at the final epoch is also reflected in the masses, so they are roughly proportional, as seen in Figure 14.In fact, filaments in overbound clouds carry twice as much mass in total as underbound clouds.Simulations with more high density regions in their PDFs (Figure 7a), such as C26, M57, S57, M57O, S57O and C57O, also favour more mass in filaments.

Filament line density
We now look at how the mass per unit length evolves in the clouds.Note that the mass-to-length ratio here is defined as / where  is the total mass of the filament and  is its total length; hence, the quantity referred to here is the average mass-to-length ratio rather than an integrated quantity.Early works (Ostriker 1964a) used idealised filament geometry to study fragmentation into cores/clumps, suggesting a critical line density of  line,cr = 2  / ≈ 16.7 As prefaced in the introduction, and examined in detail by the authors in Chira et al. (2018), this simplistic model does not account for the complexity of fragmentation, with sub-critical filaments able to fragment and super-critical filaments seen in both observations and simulations.Figure 13c shows the distribution of filament line densities seen in our simulations.We see a few extreme super-critical filaments, but a majority of them (∼ 50% -75%) sit just below the critical line.This is also seen in Figure 14 where almost all filament networks (on average, in their entirety) are subcritical with respect to thermal pressure.Since we are looking at fixed epochs rather than  evolving time-steps, it is likely that the extreme cases have very low survivability, and may be too fragmented to identify as a filaments.
As seen with the filament masses, overbound clouds form more supercritical filaments than underbound clouds.There are marginally more super-critical filaments in compressive clouds in the first epoch, but this is washed out in the final epoch, and the only correlation that remains is those simulations producing more high density regions having more super-critical filaments.As the cloud evolves, we see the expected increase in number of filaments in the super-critical regime.

Comparison with previous work
Statistical analysis of filaments in 3D is sparse in literature.The closest work to ours is that of Smith et al. (2014bSmith et al. ( , 2016)).In Smith et al. (2016), the authors identify sub-filaments in a cloud with very similar initial conditions and mixed turbulence, and characterize velocities within them.After 4 Myrs, they find a mean sub-filament mass of ∼ 10 M ⊙ consistent with our work.The longer filaments (≳ 1.5 pc) identified in Smith et al. (2014b) with masses of 40 -100 M ⊙ are roughly in line with the higher end of simulations M26, M57 and M90 seen in Figures 13a and 13b.Clarke et al. (2017) study accretiondriven turbulence in a 3 pc cylinder and find greater heirarchical and filamentary substructure in mixed turbulence compared to compressive; although we do not study hierarchy specifically, we find no evidence for it in the number of filaments.Chira et al. (2018)  A direct comparison to observations proves difficult since they identify filaments in 2D column density maps.Zamora-Avilés et al. ( 2017) investigated 3D filaments in clouds with similar initial cloud mass, and found that these fibres (the filament family for ≲ 1 pc length scales) in one 2D projection are not necessarily fibres in other projections.Suri et al. (2019) use position-position-velocity (PPV) cubes to study such fibres/short filaments in 3D in Orion A, but focus on filament widths which we fix in our work.

Hubs
As discussed in the introduction, there exists a possible link between massive-star formation and 'hubs', regions where distinct filaments meet.We present a simple model here to characterise them in our simulations: hubs are identified in a filamentary network by constructing the set of all points in filaments, and checking the number of repetitions (if any) of each point in different filaments, which we define as its rank.Rank- points, for  > 1, are classified as hubs where  filaments converge.
To study these hubs, we pick the  sink = 1000 M ⊙ epoch since the sink particles have had more time to accrete mass and form potential sites of massive stars.Figure 15 shows the mean number density in a radius  = 0.05 pc around a hub, for different hub ranks across all simulations.All hubs lie above  > 1000 cm −3 , as expected from dense regions.Since rank-2 hubs are just where two filaments meet, they are far more abundant than rank-3 and rank-4 hubs.Rank-4 hubs are primarily high density spots, followed by rank-3 and rank-2 hubs which suggests that the geometry of these junctions is indeed correlated to their density.A possible scenario is that as dense regions within a filament fragment, they create shorter filaments and naturally a junction between them.There also exists an interesting bi-modality in the distribution, which is less prominent prior to filament spine correction (see Figure B1 in Appendix B), although the correlation between hub rank and density still remains.It is possible that longer filaments, just above the DisPerSE identification threshold, produce the first mode whereas shorter, fragmented filaments around dense clumps produce the second mode.
We also investigate, for the first time on these scales, the link be-tween massive star formation and the hubs.Figure 16 shows the distance to the nearest hub for all sinks, as a function of sink mass across the simulations.Low mass sinks, in the range 0.01 M ⊙ -1 M ⊙ , show some preference for hub activity but are highly scattered.They are more likely to be strongly affected by dynamics of their surroundings, being potentially slung away from hubs.High mass sinks (> 1 M ⊙ ) show less scatter and a preference for hub activity.Hub ranks, despite being correlated to local density, have very similar distribution in the sink masses they host and their distance to these sinks.
Recent observations of hub-filament systems are limited to highmass, supercritical filaments/large cores (several hundred M ⊙ , e.g.Williams et al. 2018;Kumar et al. 2020Kumar et al. , 2022)).As more of these systems are found on lower mass scales, it will shed light on how protostars interact with these dense hubs, and compare to simulations.

CAVEATS AND FUTURE WORK
Although care has been taken to select a reasonable persistence threshold for filament identification in DisPerSE, the properties of filaments will always depend on the method of identification, and the particular set of parameters for that method -do filaments that bend at angles > 90 • count as two merging filaments or one?How does one account for discontinuities in filaments?An irregular clump or just a short filament?Applying the same methodology across simulations (such as fixing the filament identification parameters and filament width in our case) provides a means for comparative work, but cannot produce definitive conclusions.
Furthermore, our simulations also omit important feedback processes such as supernovae explosions, stellar winds, jets and outflows.These would become crucial in the late-stage evolution of the clouds, hence why our focus lies in the early stages to reduce this effect.We also lack magnetic fields, which may play a role particularly in the formation of filaments.This is an issue that would be worth exploring in future work.
Our simulations and Python library have opened up the possibility for a lot more to be studied.The statistics of the turbulent velocity field and its evolution on local and global scales are particularly interesting, and the relation between velocity flows and densities along the filament could reveal a lot about accretion mechanisms.

CONCLUSIONS
In this paper, we study the effects of large-scale turbulence in molecular clouds and its links to dense filamentary networks and starforming regions.We perform 15 high-resolution Arepo simulations of a spherical, molecular cloud of mass 10 4 M ⊙ and radius ≈ 9 pc, seeded with different types of decaying turbulence.Nine of these simulations focus on comparing mixed, purely compressive and purely solenoidal modes of turbulence in virial balance, while the other six focus on comparing different virial ratios, producing gravitationally overbound and underbound clouds.The suite of simulations is publicly available on Zenodo.Dense cores and clumps are formed in the simulations, and sink particles are used to characterise star-formation.Filaments are identified using DisPerSE, and characterised using novel algorithms released as a Python toolkit called Fiesta.Our main conclusions from the simulations are as follows: (i) Clouds seeded with different initial turbulent modes collapse differently, and a different timescale on which star-formation occurs.Compressive seeded turbulence results in an early onset of collapse into stars, approximately 1.5 Myrs -2 Myrs after the initial collapse (30% of 1 free-fall timescale  ff ), followed by mixed turbulence.Solenoidal seeded turbulence is the slowest, with star-formation beginning only 3 Myrs -4 Myrs (60% of  ff ) after collapse.The mass accretion rate is higher for compressive modes, and also for overbound clouds.
(ii) Compressive seeded and overbound clouds give rise to a centrally condensed morphology, with more compact global structure, while solenoidal seeded and underbound clouds give rise to more filament-like morphologies.The latter is characterised by a more diffuse cloud, spread over a larger volume.
(iii) The low mass end of the gas density distribution (PDFs) maintains an imprint of the initial turbulent mode.Solenoidal seeded turbulence produces more low-density regions, as do underbound clouds.However, the high-mass end of the PDF is dominated by gravitational collapse (both local and global) and hence is favoured in overbound clouds.
(iv) The sink mass functions (representing the mass going into stellar systems) are dominated by brown dwarf mass sinks in the case of solenoidal seeded clouds.These clouds generate an abundance of sinks, far more than the other turbulent modes, but are indeed bottomheavy in their mass distribution.The high-mass regime generally constitute a larger number fraction in mixed/compressive seeded and overbound clouds, but the most massive sinks (> 16 M ⊙ ) form equally across all simulations.
(v) The distribution of filament lengths and masses is seemingly unaffected by the mode of turbulence seeded, since they trace dense regions where the initial turbulence is decayed and potentially washed over by local collapse dynamics.However, clouds with higher gravitational potential energy produce more filaments, and hence longer networks and more total mass in filaments.Across all simulations, the filaments accrete mass and fragment over time, increasing the abundance of super-critical filaments.In the later stages of collapse (at  sink = 1000 M ⊙ in our simulations), the mean filament length is ≈ 0.75 pc while the mean mass is ≈ 10 M ⊙ .
(vi) Hubs/junctions where filaments converge show proportionally higher densities with increasing number of converging filaments.Sinks have a preference for hub activity, with low mass sinks showing larger scatter away from hubs.The distribution of sinks is agnostic to the number of filaments converging at the hub.
To succinctly summarise the role of initial conditions in molecular clouds: the initial turbulence imprints diffuse regions in the clouds, seen through the different diffuse morphologies, brown dwarf mass excess and low density end of the PDFs; while initial gravity dominates dense regions, seen through the high density end of PDFs, abundance of mass in filaments, and increased hub activity.conclusion of a strong correlation between hub rank and density remains.Even in the uncorrected filaments, rank-4 hubs form the densest regions followed by rank-3 and rank-2.The clumpy nature and complex hierarchy of structure within filaments adds difficulty in rigorously defining hubs and making simple assertions regarding their properties.
This paper has been typeset from a T E X/L A T E X file prepared by the author.Mean number density at hubs (in a sphere of radius  = 0.05 pc), for different hub ranks before correcting filament spines (c.f. Figure 15).

Figure 1 .
Figure 1.Column density projection in the -direction of our initial condition for the molecular clouds.The cloud is spherically symmetric and centred in the simulation box.The region enclosed by dashed grey lines, spanning ≈ 19.4 pc 3 delineates the region of interest (ROI).Our analysis is restricted to this region throughout the paper.

Figure 2 .
Figure 2. Example of turbulent velocity field with virial-type equipartition, generated using RNG seed = 57.Figures show the heat-map of   component of the velocity field for mixed, compressive and solenoidal modes respectively.The velocities shown are on the surface on the cloud, for a portion of the cloud facing the viewer; the clouds are indeed spherical, but appear ellipsoidal due to skewed axes.Note the large-scale nature of fluctuations and in particular, the striated pattern typical of compressive turbulence.
2 For more on filament identification codes, Chira et al. (2018) provides an excellent comparison of the several available tools.

Figure 3 .
Figure 3. Example of filamentary structures identified using DisPerSE in simulation C57.Top panel shows a projection of the entire central region of the simulation box with filaments (white).Middle panel is similar but shows the projection of dense regions > 5000 cm −3 only.Bottom panel shows a 3D plot of these dense regions (colour), their axial projections (grey shadows), and filaments (red).

Figure 4 .
Figure 4. Example of the two filament spine correction algorithms implemented -'centre of density' and 'max density' -for a randomly chosen filament, and their corresponding radial density profile.Top row: 3D plots showing the filament (and projection) pre-correction in dark grey (light grey) and post-correction in red (light red), with the Arepo cells (colour) belonging to each longitudnal cylinders of radius  cyl = 0.1 pc along the filament spine.Bottom row: Plots shows the density of Arepo cells versus their distance from filament spine (black scatter), along with contours (colour).

Figure 5 .
Figure 5. Top: Evolution of total sink mass as a function of time for all 15 simulations.The horizontal dashed grey lines demarcate the epochs of interest used in our analysis, corresponding to a fixed total sink mass -250 M ⊙ , 500 M ⊙ , 750 M ⊙ and 1000 M ⊙ -achieved at different times for each simulation.The vertical grey lines show quarters of 1 free-fall timescale ( ff ≈ 6.5 Myrs).Bottom: Accretion rate onto sinks as a function of time.

Figure 6 .
Figure 6.Morphology of all clouds at the last epoch,  sink = 1000⊙.The 3D plots show dense regions > 1000 cm −3 (colour), and their axial projections (grey).Compressive seeded clouds show more compact centrally condensed structures, while solenoidal seeded clouds are more filament-like.
Figure 7. (a) Mass-weighted PDFs of number density for all simulations at the  sink = 1000 M ⊙ epoch, and the mean  ∼ 100 cm −3 (dotted black line).(b) CDFs of number density for clouds in initial virial balance, but different seeds, at  sink = 250 M ⊙ epoch.(c) Same as (b) but at  sink = 1000 M ⊙ epoch.(d) CDFs of number density for clouds with different virial ratios at  sink = 250 M ⊙ epoch.(e) Same as (d) but at  sink = 1000 M ⊙ epoch.

Figure 8 .
Figure 8. Two-sampled Kolmogorov-Smirnov (KS) test -values for each combination of log-scaled mass-weighted CDFs in Figure 7.We focus on a single RNG seed = 57 to compare the effects of turbulent mode and boundedness.

Figure 9 .
Figure 9. Violin plots of sink masses for all simulations at  sink = 250 M ⊙ and  sink = 1000 M ⊙ epochs, truncated to those > 0.1 M ⊙ since they are likely to form stars.The plots show 25%-75% quartiles (coloured bars), along with the median (black dot) and mean (black star), and total number of sinks along the -axis.

Figure 11 .
Figure 10.(a) Probability density function of sink masses for all simulations at the final epoch, taking the form of an IMF (note, however, that these mass functions are for sinks, and not stars, so it is strictly not an IMF).The Chabrier (2005) model is generated using github.com/keflavich/imfand normalized to the displayed mass range.Most notably, the solenoidal clouds show an excess of low mass sinks (brown dwarfs) < 0.1 M ⊙ .(b) 2D projection of simulation S57 (grey), overlayed with the distribution of low mass sinks < 0.1 M ⊙ (red scatter) and high mass sinks > 0.1 M ⊙ (blue scatter).They appear to be similarly clustered, indicating that these brown dwarfs are not forming in unique spatial environments.

Figure 12 .
Figure 12.Number of filaments at the four epochs of interest for all simulations, showing a consistent or increasing trend.

Figure 16 .
Figure 16.Distance of individual sinks from the closest hub as a function of their sink mass for all simulations at  sink = 1000 M ⊙ , with different colors indicating the rank of the closest hub.The distribution along  and  axes for each hub rank are plotted as histograms with median (solid line).The grey region demarcates sink masses below 0.1 M ⊙ , which are unlikely to form stars.

Figure A1 .
Figure A1.Total length of filaments as a function of DisPerSE persistence threshold for the two epochs  sink = 250 M ⊙ and  sink = 1000 M ⊙ .The dashed line indicates the ideal persistence threshold, used for analysis in this paper, while the dotted lines indicate 3 confidence region around this value.
Figure B1.Mean number density at hubs (in a sphere of radius  = 0.05 pc), for different hub ranks before correcting filament spines (c.f.Figure15).

Table 1 .
(Vázquez-Semadeni et al. 2019)the suite of 15 simulations performed.For each turbulent mode, three seeds (26, 57 and 90) were used to generate turbulent velocity fields.Once generated, they are scaled to either virial balance (simulations 1 -9) or non-virial ratios (simulations 10 -15).The virial ratios  vir = 2.0 and  vir = 0.6 correspond to underbound and overbound clouds respectively.A Mach number M ≈ 3 -4, defined as the ratio of velocity dispersion and sound speed   /  , is typical of the GHCC model(Vázquez-Semadeni et al. 2019).As expected, the overbound and underbound Mach numbers are

Table 2 .
Point cloud spread of cells above 1000 cm −3 , as defined in Equation 2, at the epoch  sink = 1000 M ⊙ .
Violin plots showing the filament length statistics at  sink = 250 M ⊙ and  sink = 1000 M ⊙ epochs.The plots show 25%-75% quartiles (coloured bars), along with the median (black dot) and mean (black star).Empty circles are outlier filaments > 5 pc.
(b) Same as (a) but for filament mass.
simulate a much larger 40 kpc long vertical box with a minimum 1.7 pc filament length, beyond the length scale relevant to this work.Figure 15.Mean number density at hubs (in a sphere of radius  = 0.05 pc), for different hub ranks.It can be seen that higher rank hubs exist in denser regions.

Table A1 .
Example of 1 errors introduced in filament properties at the final epoch due to persistence threshold parameter in DisPerSE.