SILCC-Zoom: the dynamic balance in molecular cloud substructures

How molecular clouds fragment and create the dense structures which go on to form stars is an open question. We investigate the relative importance of different energy terms (kinetic, thermal, magnetic, and gravity - both self-gravity and tidal forces) for the formation and evolution of molecular clouds and their sub-structures based on the SILCC-Zoom simulations. These simulations follow the self-consistent formation of cold molecular clouds down to scales of 0.1 pc from the diffuse supernova-driven interstellar medium in a stratified galactic disc. We study the time evolution of seven molecular clouds (five with magnetic fields and two without) for 1.5-2 Myr. Using a dendrogram, we identify hierarchical 3D sub-structures inside the clouds with the aim to understand their dynamics and distinguish between the theories of gravo-turbulent fragmentation and global hierarchical collapse. The virial analysis shows that the dense gas is indeed dominated by the interplay of gravity and turbulence, while magnetic fields and thermal pressure are only important for fluffy, atomic structures. Over time, gravitationally bound sub-structures emerge from a marginally bound medium (viral ratio $1 \leq \alpha_{\rm vir}^{\rm vol}<2$) as a result of large-scale supernova-driven inflows rather than global collapse. A detailed tidal analysis shows that the tidal tensor is highly anisotropic. Yet the tidal forces are generally not strong enough to disrupt either large-scale or dense sub-structures but cause their deformation. By comparing tidal and crossing time scales, we find that tidal forces do not seem to be the main driver of turbulence within the molecular clouds.


INTRODUCTION
Star formation occurs inside cold, dense and mostly molecular gas in the interstellar medium (ISM) -in so-called molecular clouds (MCs).They are typically thought of as the largest structure that can decouple from galactic dynamics (Chevance et al. 2020).The nature of MCs, their lifetime, formation mechanism, and whether they are bound or not, all have direct consequences on our understanding of the possible ways in which dense structures and eventually stars can form inside these clouds.
Several key discoveries over the past decades have facilitated and developed our understanding of molecular clouds.We know that they are in rough energy equipartition between kinetic and potential energy (e.g.see Dame et al. 1986;Solomon et al. 1987;Blitz et al. 2007).Only a small fraction of gas inside these clouds eventually ends up in stars (Kennicutt 1998;Genzel et al. 2010;Krumholz et al. 2012), which leads to a low star formation efficiency of the order of a few percent, although with some variation (Murray 2011).Molecular clouds are known to exhibit supersonic linewidths (Larson 1981;Orkisz et al. 2017), but generally lack a clear collapse signature such as redshifted lines (Zuckerman & Evans 1974).While this has recently been challenged with observations of gravity driven flows (Kirk et al. 2013;Peretto et al. 2014;Chen et al. 2019;Shimajiri et al. 2019), it is unclear if such flows extend to entire clouds (Chevance ★ E-mail: ganguly@ph1.uni-koeln.deet al. 2022).The linewidth is related to the mass and size of the molecular clouds themselves (Larson 1981;Solomon et al. 1987;Brunt 2003).The original power-law exponents proposed by Larson (1981) have since been revised and challenged.For example, Heyer et al. (2009) show that the Larson coefficient has a dependency on the MC surface density, and there have been studies pointing out that the Larson relation may not hold inside a single cloud (Wu et al. 2010;Traficante et al. 2018).
There are several possible ways to interpret these key observations.The traditional view, originally proposed by Zuckerman & Evans (1974), is that the supersonic linewidths in observed clouds are caused by small scale turbulence that provides overall support and prevents the entire cloud from freely collapsing under gravity.In this so-called gravo-turbulent (GT) scenario, clumps and cores inside molecular clouds form as a result of compression in the highly turbulent gas (see e.g.Klessen 2001;Padoan & Nordlund 2011).Only a few of the compressed regions become gravitationally bound structures, which collapse and form stars.In this scenario, the Larson-like power law connecting the turbulent velocity dispersion,   , and the linear size of the cloud, , i.e.   ∝  0.5 , originates from a turbulent cascade of large-scale motions (see e.g.excellent reviews by Mac Low & Klessen 2004;Ballesteros-Paredes et al. 2007;McKee & Ostriker 2007;Hennebelle & Falgarone 2012).This view has since been challenged.Magnetic fields have been proposed as a possible mechanism to provide support against gravitational collapse (Mouschovias 1976a,b;Padoan & Nordlund 2011;Federrath & Klessen 2012;Ibáñez-Mejía et al. 2022).However, after the emergence of key observations suggesting that most cloud structures are magnetically supercritical , it is generally accepted that magnetic fields cannot be solely responsible for the observed, low star formation efficiency (Crutcher 1999;Bourke et al. 2001;Troland & Crutcher 2008;Falgarone et al. 2008;Crutcher et al. 2010).
It is unclear if MCs need any support at all.Several studies suggest that MCs are transient entities, and that their lifetimes are short enough that they may not need any support (Bash et al. 1977;Leisawitz et al. 1989;Elmegreen 2000;Hartmann et al. 2001;Engargiola et al. 2003;Kawamura et al. 2009;Meidt et al. 2015;Mac Low et al. 2017).
In recent years, chaotic, hierarchical gravitational collapse has emerged as another possible alternative to the gravo-turbulent scenario.The hierarchical collapse scenario, originally developed by Hoyle (1953) and suggested as a mechanism for the collapse of a cloud by Goldreich & Kwan (1974), has been considered briefly over the past decades (e.g.Elmegreen 2000;Hartmann et al. 2001;Vázquez-Semadeni et al. 2007).In its most recent form, the socalled global hierarchical collapse (GHC) scenario (see e.g.Vázquez-Semadeni et al. 2019), the interstellar turbulence produces nonlinear density fluctuations, which facilitate a multi-scale hierarchical collapse of the entire molecular cloud through a mass cascade across scales.The chaotic nature of the collapse masks any obvious signatures of global collapse, and the low efficiency of star formation is explained by the fact that stellar feedback from the first generation of massive stars disperses the cloud before it can convert more of its gas mass into stars (see e.g.Vázquez-Semadeni et al. 2009, 2017;Ballesteros-Paredes et al. 2011;Kuznetsova et al. 2015Kuznetsova et al. , 2018)).
Most analyses of the evolution of MCs focus on the self-gravity of the gas itself (Dobbs et al. 2011).The complex nature of the medium, however, implies that tidal forces can have a potentially significant sway over the gas dynamics.For example, stability produced by tidal forces has been invoked to explain the peak of the initial mass function on core-scales (Lee & Hennebelle 2018;Colman & Teyssier 2020).For star clusters, it has been shown that fully compressive tides can slow down the dissolution of clusters (Renaud et al. 2011), and that young clusters are preferentially located in fully compressive regions of galaxies (Renaud et al. 2009).We aim to consider tidal forces here in the context of the (internal) dynamics of MCs.
In the presented work, we perform a detailed energy analysis of self-consistently formed MCs and their sub-structures based on the SILCC-Zoom simulations (Seifried et al. 2017).The MCs are situated inside a larger scale multi-phase ISM environment (Walch et al. 2015;Girichidis et al. 2016).This allows us to do a rigorous analysis of the MC dynamics and to test whether the GT or the GHC scenario are in control.We apply a three-dimensional dendrogram algorithm to define hierarchical structures (Rosolowsky et al. 2008).
The paper is structured as follows.We first present the details of the numerical setup in section 2. In section 3, we highlight the different scaling relations (Larson's velocity-size and mass-size relations, and the Heyer relation).We then discuss the energetic balance between self-gravity, kinetic, magnetic, and thermal energy in section 4. By studying the dynamic balance of these structures, we find that larger scale and denser, molecular structures are primarily governed by the interaction of the turbulent kinetic energy and self-gravity, with magnetic and thermal energies playing a subservient role.In section 5, we discuss the relative importance of self-gravity and the gravity from the surrounding medium, by both a direct comparison of acceleration terms as well as a more detailed tidal analysis.From the direct comparison of the acceleration vectors, we find that bound and marginally bound structures are dominated by self-gravity over external gravity.We further assess the nature and extent of deformation caused by the gravity of the surrounding clumpy medium by performing a tidal analysis, and find that the overall gravitational field is highly anisotropic but still mostly fully compressive.Finally, we present the summary of our findings in section 6.

NUMERICAL METHODS AND SIMULATION
We present results based on the SILCC-Zoom simulations (Seifried et al. 2017).The SILCC-Zoom simulations are zoom-in simulations of MCs forming self-consistently within the larger-scale, multi-phase ISM SILCC simulations (Walch et al. 2015;Girichidis et al. 2016).In this section, we highlight some key features of the simulations.Further details regarding the simulations can be found in Seifried et al. (2017).For description of the hydrodynamic (HD) and magnetohydrodynamic (MHD) clouds used for the present analysis, see Seifried et al. (2017) and Seifried et al. (2019), respectively.
All simulations were performed using the adaptive mesh refinement (AMR) code FLASH, version 4 (Fryxell et al. 2000;Dubey et al. 2008).We present results from runs with and without magnetic fields.The hydrodynamic cloud simulations have been performed using the MHD 'Bouchut 5-wave solver' (Bouchut et al. 2007;Waagan 2009), with magnetic field strength set to zero.The MHD simulations have been run using an entropy-stable solver which guarantees positive pressure and minimizes dissipation (Derigs et al. 2016(Derigs et al. , 2017(Derigs et al. , 2018)).
The entire simulation domain consists of a box of size 500 pc × 500 pc × ± 5 kpc, with the long axis representing the vertical −direction of the stratified galactic disc.The box is set with periodic boundary conditions in the − and −directions, and outflow boundary conditions in −direction.The initial gas surface density is set to Σ gas = 10 M pc −2 which corresponds to solar neighbourhood conditions.The vertical distribution of the gas is modelled as a Gaussian, i.e.  =  0 exp(− 2 /2ℎ 2  ), where ℎ  = 30 pc is the scale height and  0 = 9 × 10 −24 g cm −3 is the midplane gas density.The initial gas temperature is set to 4500 K.For runs with magnetic fields, the magnetic field is initialized along the −direction, i.e.B = (  , 0, 0) with   =  ,0 √︁ ()/ 0 and the magnetic field strength in the midplane  ,0 = 3.The field strength is chosen to be in accordance with recent observations (e.g.Beck & Wielebinski 2013).
All simulations include self-gravity as well as an external galactic potential due to old stars.The external potential is calculated using the assumption of a stellar surface density of Σ star = 30 M pc −2 with a scale height of 100 pc.The self-gravity of the gas is computed using the Octtree-based algorithm by Wünsch et al. (2018).
Apart from the dynamics of the gas, we also model its chemical evolution using a simple non-equilibrium chemical network based on hydrogen and carbon chemistry (Nelson & Langer 1997;Glover & Mac Low 2007;Glover et al. 2010).For this purpose, we follow the abundance of ionized, atomic, and molecular hydrogen (H + , H, H 2 ), carbon-monoxide (CO), ionized carbon (C + ), free electrons (e − ), and atomic oxygen (O).At the beginning of the simulation, hydrogen is fully atomic within the disc midplane and carbon is fully ionized.The simulations also include an interstellar radiation field (ISRF), with strength of  0 = 1.7 in Habing units (Habing 1968;Draine 1978).The ISRF is attenuated depending on the column density of the gas, thereby accounting for dust shielding as well as (self-) shielding of H 2 and CO.This is accounted for using the T R O D module developed by Wünsch et al. (2018) and takes into account the 3D geometry of the embedded region.
The turbulence in the simulations is generated by supernova ex-plosions (SNe).The explosion rate is set to 15 SNe per Myr, which is consistent with the star formation rate surface density expected for the given gas surface density from the Kennicutt-Schmidt relation (Schmidt 1959;Kennicutt 1998).50% of the SNe are placed following a Gaussian random distribution with a scale height of 50 pc along the −direction , while the other 50% are placed at density peaks of the gas.This prescription of mixed supernova driving creates a multiphase turbulent ISM which can be used as initial condition for the zoom-in simulations (Walch et al. 2015;Girichidis et al. 2016).
For the SILCC-Zoom simulations, we drive the medium with SNe up to a certain time  0 .Up to  0 , the uniform grid resolution is 3.9 pc.At  0 , the SNe are stopped and multiple regions are identified for the zoom-in process, primarily due to their molecular gas content at a later simulation time, when the simulation has been continued on the coarse resolution.This corresponds to the start of the evolution of the zoom-in regions, and we refer to this point as  evol = 0.The total simulation time  is related to the cloud evolution time  evol as (1) The different simulations are compared at similar  evol .For tracing the evolution of the clouds, in the selected regions, the AMR grid is allowed to refine to a higher resolution in order to capture the internal structure of the forming MCs.These regions are referred to as "zoom-in regions".In each SILCC box we run two zoom-in regions simultaneously, such that we follow eight zoom-in regions in total.The different MHD runs are seeded with a different set of of random supernovae.All other initial conditions are the same.All the runs presented here have an effective spatial resolution of 0.125 pc.
For more details on the zoom-in process see Seifried et al. (2017).

Cloud selection
For the analysis presented in this work, we look at eight different (62.5 pc) 3 zoom-in regions.Each zoom-in region contains one molecular cloud (called 'MC1', 'MC2', and so on).Six simulations include magnetic fields (run names with '-MHD'), while the other two are of hydrodynamic nature ( = 0; run names with '-HD').We present some basic information of the different MCs in Table 1, and a projected view of all eight clouds at time  evol = 3.5 Myr in Fig. 1.The clouds we analyze have total gas masses of 5.4 − 10.1 × 10 4 M , and H 2 masses between 0.9 − 2.1 × 10 4 M at t evol = 2 Myr.Note that the time at which we begin the zoom-in,  0 , varies for the different simulations (see Table 1).The evolution time  evol , which we refer to throughout the paper when comparing the results for the different MCs, represents the "age" of the molecular cloud and corresponds to the time that has passed since the start of the zoom-in.We perform a detailed analysis of the various clouds at different times between  evol = 1.5 − 3.5 Myr for the two hydrodynamic clouds, as well as for  evol = 2 − 3.5 Myr for the MHD clouds.We include the earlier time step for the hydrodynamic clouds as they evolve more rapidly compared to their MHD counterparts and show discernible structures early on.We do not include sink particles in these simulations, and hence the end time is chosen such that the cloud sub-structures are still well resolved with respect to the local Jeans length (Truelove et al. 1997).Of the eight simulations, we find that even at  evol = 3.5 Myr, MC3-MHD has a very low molecular gas fraction (∼1900 M molecular H 2 of a total mass of ∼20,000 M , see Table 1).The resulting dendrogram analysis is rather incomplete because this cloud lacks discernible dense structures.Therefore, we omit MC3-MHD from our further analysis.1. Basic information on the eight analyzed simulations.From left to right we list the run name, whether it is a magnetized run or not, the time when the AMR zoom-in started, as well as the total mass and the molecular hydrogen mass at  evol = 2 Myr. a We discard MC3-MHD from our further analysis because of its low molecular gas content and lack of dense structures (see also Fig. 1).

Sub-structure identification and modelling
To identify structures in the MCs, we use a dendrogram algorithm (Rosolowsky et al. 2008).Dendrograms are a model independent method to determine hierarchical structures.We specifically use the python package astrodendro (Robitaille et al. 2019) and apply the method to find density (sub-)structures in 3D.For this purpose, we convert the AMR data to a uniform grid at our effective spatial resolution of 0.125 pc.Given the initial density cubes, our dendrogram analysis essentially depends on three free parameters: the initial starting threshold  0 , the density jump Δ, and the minimum number of cells included in a structure,  cells .For the analysis presented here, unless otherwise stated, we use  0 = 10 −22 g cm −3 , Δ(log 10 ) = 0.1 and  cells = 100.Defining Δ(log 10 ) implies that the dendrogram is built on the logarithm of the density.In order to reduce the number of random fluctuations, we only retain those structures in the dendrogram tree, which eventually reach a density of  peak > 10 −21 g cm −3 .A view of the resulting dendrogram tree can be seen in the left panel of Fig. 2 for MC1-MHD at  evol = 2 Myr.We highlight two branches of the dendrogram tree to show the hierarchical nature of the structures, a larger scale structure in blue, and a contained denser filamentary structure in red.The projections of these two structures can be seen in the right panel of Fig. 2, as contours over column density maps projected along the −, −, and −direction, respectively.We note that we varied the three free parameters within a reasonable range in order to test the robustness of our results.Increasing  cells washes out the smallest sub-structures, while decreasing  0 results in larger structures on the largest scales, but leaves the smaller scales unaffected.Increasing the minimum value of Δ causes fewer structures to be identified but does not affect the general trends for different sub-structures.We have also performed the dendrogram analysis by using linear density jumps (i.e.dendrogram built on the density field instead of the logarithmic density field) of Δ = 10 −22 g cm −3 and Δ = 10 −23 g cm −3 , in addition to using logarithmic density bins.This also leaves the statistical properties of the resulting structures basically unaffected.
Examples of the resulting dendrogram structures can be seen in Fig. 3, for MC5-MHD at  evol = 3.5 Myr.Projections of different dendrogram leaf structures (i.e.strutures that contain no further sub-structures inside them) are plotted over the column density as contours.We differentiate between sub-structures based on their molecular content.Molecular structures (H 2 mass fraction > 0.5) are plotted in solid lines, while atomic (H 2 mass fraction < 0.5) struc- The side length of each plot is 62.5 pc.The MHD clouds (top right, and bottom row) have elongated filamentary structures with extended diffuse envelopes, while the hydrodynamic clouds (top left) appear to be more clumpy and fragmented.Note that MC3-MHD has been excluded in further analysis due to its lack of discerning structures and low molecular gas content (see Table 1).
Figure 2. Left: Dendrogram tree with two sub-structures highlighted in blue and red for MC1-MHD at  evol = 2 Myr, plotted using astrodendro.Right: We depict these two sub-structures as contours on the column density maps of MC1-MHD projected along the -, -, and -axis, respectively.The blue and red contours correspond to the branches of the same colour from the left panel.It is apparent that the structures look filamentary on the -, -projections, but rather clumpy on the -projection.
tures are plotted in dashed lines.Visually, we see that the molecular structures trace the denser spine of the cloud, while the atomic structures represent more the surrounding envelope.The number of total, as well as molecular structures for the various analyzed MCs can be seen in Table 2 for  evol = 3.5 Myr.From Table 2, we can see that apart from the more diffuse MC2-MHD, almost all the clouds are primarily dominated by molecular structures at the latter stage.The exact number of structures depends on the dendrogram parameters chosen, but the fraction of molecular and atomic structures is nonetheless indicative of the chemical state of the clouds.
Once we obtain the dendrogram, we aim to approximate the density sub-structures as ellipsoids.In a companion paper (Ganguly et al., in prep.), we use this method to classify the morphology of the given structures.For the present paper, however, it allows us to compute a number of physical quantities dependent on the shortest axis of a given structure (such as the crossing time, see Fig. 16).
For each structure, we compute an equivalent ellipsoid that has the same mass and the same moments of inertia as the original structure.Let us consider a uniform density ellipsoid of mass  and semi-axis lengths , ,  where  ≥  ≥ .The moments of inertia along the three principal axes will be given as follows: where   ≥   ≥   .If we now compute the principal moments of inertia of a given dendrogram structure to be , , and , respectively, then the ellipsoid has an equivalent moment of inertia if  From left to right are listed the run name, the number of total analyzed structures in each region, the number of structures which have >50% of their hydrogen mass in molecular form, and the fraction of such molecular structures at  evol = 3.5 Myr.The exact number of structures depends on the dendrogram parameters.The high fraction of molecular structures is partially a result of specifying that all structures should have a peak density  peak > 10 −21 g cm −3 .
This leads to the following equation for computing , , and  of the equivalent ellipsoids: We emphasize that the (sub-)structures identified with our dendrogram analysis represent volumina that are enclosed by the defined density contours.This does not imply that the structures are longlived entities.They may either condense or diffuse and disperse completely.They may also accrete gas, grow in volume, merge with each other, or fragment.In the following, we analyse the dynamical state of the density (sub-)structures in detail and relate them with observed relations.

SCALING RELATIONS
In his seminal work, Larson (1981) showed that MCs follow a relation between their linewidth, representing the 1D gas velocity dispersion, and their linear size, (5) While Larson (1981) originally found  = 0.38, this has later been revised to  = 0.5 (Solomon et al. (1987), or more recently, Brunt This points towards a change in the cloud dynamics once self-gravity becomes dominant (Wu et al. 2010;Traficante et al. 2018).Reexamining Larson's laws, Heyer et al. (2009) show that this can be explained when taking into account the surface density of the structures.Here, we investigate the Larson and Heyer scaling relations of our dendrogram structures.We compute the mass-weighted 1D velocity dispersion,  1 , from where  is the total mass of the structure and v 0 is the velocity of its centre of mass, computed as The integration is performed over the entire volume  of the structure.We choose a mass-weighted average because it roughly represents an intensity-weighted average, which is often used in observations.We estimate the linear size  from the volume of the structure as In Fig. 4, we plot  1D vs.  for different times (from left to right).The top row shows a compilation of all structures extracted from MC1-HD and MC2-HD, while the bottom row shows all structures found in the five MHD clouds.Note that we show slightly different times for the hydrodynamic ( evol = 1.5, 2.0, 3.5 Myr) and the MHD simulations ( evol = 2.0, 2.5, 3.5 Myr), owing to the fact that the former show discernible structures comparatively quicker (see below, also see Seifried et al. (2020)).We distinguish between mostly atomic structures (empty symbols) and mostly molecular structures (filled symbols) as defined in section 2.2 (see also Fig. 3).The colour bar shows the density threshold,  thr , of each structure, defined as the minimum density value contained in the structure.
For early times, all structures seem to follow a power law relation for both HD and MHD clouds.The exponents were derived from fitting the logarithm of the data points with a non-linear least squares method.For the hydrodynamic clouds we find  = 0.61 ± 0.03, which is slightly steeper than Larson's relation, while the best fit for the MHD clouds gives  = 0.54 ± 0.03 in quite good agreement with the observed scaling relation (Eq.5), although the scatter is significant.
At later times, high density structures tend to trail off to the left of the power law relation, as these sub-structures have decreasing size but retain a high (mass-weighted) velocity dispersion.We will later see that this is a sign of gravity becoming dynamically important.The obvious branch-like collections of points leading towards small scales are actually identified as those dendrogram branches which lead to the densest sub-structures.Due to the presence of these branches, the fitted power law slope is much steeper than the expected  = 0.5 for later times.From these plots one can also clearly see that the MHD clouds evolve more slowly than the non-magnetized clouds, as the MHD clouds only show dense, flat branches at  = 3.5 Myr, while these are already present at  evol = 2.0 Myr for the HD runs.
A similar behaviour is seen in the mass-size relation shown in Fig. 5. Apart from a best-fit power law, based on a non-linear least squares fit on the logarithm of the values (similar to discussed above), we also plot here lines corresponding to  ∝  and  ∝  2 in dash-dotted and dotted lines, respectively.The overall fitted slope is generally steeper compared to the  ∝  2 expected from the Larson relation.This steeper slope rather describes the inability of a single exponent to explain the behaviour of all the structures.Visually, the less dense and mostly atomic structures seem to fit better with a  ∝  2 type power law.The denser structures on the other hand clearly branch off to follow a much shallower power law more similar to  ∝ .This is reminiscent of more centrally condensed profiles such as Bonnor-Ebert spheres (Ebert 1955;Bonnor 1956), although in our case the structures are far from spherically symmetric.
Larson's velocity dispersion-size relation was later modified by Heyer et al. (2009) to show that the Larson coefficient  1D / 1/2 can depend on the surface density, Σ, of the structure.We plot the Heyer relation, i.e.  1D / 1/2 vs. Σ in figure 6.The surface density of a given structure is calculated as Here  avg is the volume-average density of the structure: The dash-dotted diagonal line represents the  1D = (/5) 1/2 Σ 1/2  1/2 line from Heyer et al. (2009) representing virial equilibrium between kinetic energy and self-gravity ( = 1, where  is the kinetic virial parameter).Other values of  correspond to lines parallel to the Heyer line in log-log space and are plotted with dotted lines.The virial parameter  obtained in this manner is only a rough representation of the true virial equilibrium, and includes only the kinetic and potential energy terms.We find that at earlier times there is very weak, if any, trend along the Heyer line of  = 1 for both hydro and MHD structures.This picture evolves over time.At  evol = 3.5 Myr (figure 6, right column), we find two groups of structures: mostly molecular and dense structures that follow the Heyer line, for both HD and MHD clouds and less dense, mostly (but not exclusively) atomic structures that show no correlation at all.Camacho et al. (2016) have found similar differing behaviour when they look at structures at different density thresholds in turbulent box simulations.Similar results have been also been found by Weis et al., in prep., when analyzing structures forming in colliding flow simulations.We additionally find that the molecular content of the gas is also tracing these two families of structures relatively well.
There are various possible interpretations for structures following the Heyer relation.The original interpretation of Heyer et al. (2009) was that they represent gravitational equilibrium.Ballesteros-Paredes et al. (2011) point out that, given the typical error estimates, they could also be consistent with gravitational infall, as the infall velocity and virial equilibrium velocity only differ by a factor of √ 2. A similar discussion was also presented in Larson (1981).In addition, the estimate of the surface density here is rather crude, and does not incorporate the fractal nature of the complicated medium.To investigate different possible scenarios, we perform an analysis on the different energy terms in the next section.

ENERGETICS
The interplay of various forces determines the dynamics of the structures involved.For a structure embedded in a more diffuse medium, we consider the self-gravity, kinetic energy, as well as thermal and magnetic energy (for the MHD runs) of the structure.We limit ourselves to the discussion of the various volume terms and do not consider the time dependent or the surface terms of the full virial theorem.Below, we briefly discuss how we determine the various terms.

Energy terms
The gravitational energy of a structure can be separated in two parts -the self-gravity due to its own mass, and the gravitational energy due to the fact that the structure is embedded in an asymmetric gravitational potential caused by the surrounding gas.For the analysis of this section, we have limited ourselves to considering only selfgravity for a more traditional virial analysis.In section 5, we analyse the relative importance of the external medium.
The self gravitational energy of the structure is defined as We compute this term by using a KD tree with a geometric opening angle of 0.5 radians.By comparing with the exact potential computed from a direct sum, we find that the errors introduced by the tree method are typically ∼ 1%.We also compute the volume term of the magnetic energy, which can contribute to the important magnetic pressure term.This pressure term can act against compression and prevent a structure from collapsing.The magnetic energy of a given structure is computed as where |B| is the magnitude of the magnetic field.We compute the thermal energy of a structure as where  is the thermal pressure of the gas.The computation of the kinetic energy,  KE , of a structure is done as where v 0 is the bulk velocity computed from Eq. 7.
From the above quantities, we estimate a virial ratio of the volume energy terms: The magnetic energy is absent for purely hydrodynamic simulations.Eq. 15 is a simplification of the full virial theorem (McKee & Zweibel 1992;Dib et al. 2007), which includes both time-dependent as well as surface terms, which is beyond the scope of this paper.However, we discuss the surface terms briefly at the end of Section 4.2.We define bound structures as those with a virial ratio  vol vir < 1 and marginally bound structures as 1 ≤  vol vir < 2. The definition of marginally bound structures here is motivated by the fact that, in presence of only self-gravity and kinetic energy, a structure with  vol vir = 2 would have a net zero total energy,  tot , i.e.
and as a consequence, is reduced to  = 2.

Virial analysis
In Fig. 7, we show the time evolution of the ratio of different energy terms and the self gravitational energy, as well as  vol vir (bottom row) as a function of their size  (see Eq. 8) for the two hydrodynamic clouds and in Fig. 8 for the different MHD clouds.We plot the substructures at 2 Myr (left column) and at 3.5 Myr (right column).For the HD clouds (Fig. 7), the top two rows present the evolution of thermal and kinetic energy over  PE as a function of size .The colors represent the density threshold of each analyzed structure.For the MHD clouds (Fig. 8), we also plot the magnetic energy over  PE in the top row.The red shaded region shows the marginally bound region.According to Eq. 15, for the ratio of magnetic to potential energy, marginally bound is where 1 ≤  B /| PE | < 2; for the kinetic and thermal energies, this represents a region where 0.5 ≤  KE /| PE | < 1 and 0.5 ≤  TE /| PE | < 1, respectively.This results 1 ≤  vol vir < 2. For the sub-structures in the HD clouds (Fig. 7), we find that the thermal energy is mostly irrelevant for the large scale, as well as for the denser and molecular structures.However, it seems to be important for atomic structures on all scales and for some smallscale molecular structures, particularly at later times.The  KE /| PE | ratio is far closer to unity for large scale (>10 pc) and denser ( thr 10 −20 g cm −3 ) structures.Due to the significant contribution of  KE , we find that only a few structures seem to be bound, although the number of such structures increases from 2 to 3.5 Myr.When we consider the full virial parameter, the behaviour is quite similar to that of the kinetic energy, with the addition that even fewer structures are bound at all times.The emergence of gravitationally bound structures over time suggests a growing importance of gravity as the cloud becomes more evolved.This complements the picture we see in Larson's  1D - relation (Fig. 4), where a flat line-width size relation emerges for the densest branches over time.
For the MHD clouds (Fig. 8), the magnetic energy is relevant only for less dense and mostly atomic structures, and does not appear to be significant for the larger scale structures or the denser and potentially star forming structures.This is true at earlier as well as later times.This behaviour is echoed in the evolution of  TE /| PE |, which is quite similar to the behaviour in the hydrodynamic clouds.As in the hydrodynamic case, the dynamics of larger scale, as well as dense structures is mostly dominated by the interplay of  KE and  PE .The behaviour of  vol vir for HD and MHD clouds is also similar and follows the overall trend of  KE /| PE |.We note one exception to the purely hydrodynamic structures here though.While at earlier times, we barely have any bound structures at 2 Myr, this number is much higher at 3.5 Myr.This suggests that, with magnetic fields, bound structures form more slowly from a marginally bound environment.
The emergence of gravitationally bound structures over time from a marginally bound environment, suggests that gravity only gradually becomes important for the densest branches, but is not the primary driving force for the formation of the structures.Our findings are more in line with the GT scenario of structure formation, where gravity takes over only once the turbulent structures are sufficiently compressed.
We stress that the aspect of gradually increasing boundness is essential here.If we look at only the relatively late stages at  evol = 3.5 Myr for the large-scale cloud structures, which are all molecular and lie close to  vol vir = 1, it is difficult to distinguish between GT and GHC.This can be resolved only by looking at earlier times.The behaviour of the kinetic virial ratio, as well as  vol vir , for the largestscale structure in each cloud can be seen in Table 3.We find that at 2 Myr, almost all largest-scale structures, with the exception of MC6-MHD, have  vol vir > 1.8, with three of the clouds having  vol vir > 2. This is dramatically reduced at later times, with only MC2-MHD being unbound ( vol vir > 2) at  evol = 3.5 Myr.The decrease in  vol vir for large scale structures which are unbound or marginally bound suggests that the clouds are being compressed, likely by larger scale flows, towards becoming gravitationally bound (see also Section 4.3).While we explicitly show this here only for  the largest structures, the dependence of  vol vir on  is relatively flat at large scales ( > 10 pc, Figs.7 and 8, bottom panels) and we can therefore conclude that this trend is not a consequence of e.g. the choice of our initial dendrogram threshold.
While over time we do see the emergence of individual structures below the bound line, as well as a decrease in the value of  vol vir at large scales, the overall statistical behaviour of most of the cloud substructures remains more or less unchanged.This can be seen in Fig. 9, where we show the time evolution of the average of the logarithm of  vol vir , log 10  vol vir , for atomic, molecular, and dense molecular structures (structures which are molecular and have a threshold density  thr > 10 −20 g cm −3 ), respectively.Here, we sum over all the different cloud structures, differentiating only between hydrodynamic (Fig. 9, left panel) and MHD clouds (Fig. 9, right panel).The vertical error bar shows the standard error on the mean of log 10  vol vir .Fig. 9 illustrates that the atomic and molecular structures have clearly different energetic behaviour.We can attribute this to the fact that they trace different parts of the parent cloud, as can be visually seen in Fig. 3.The chemical evolution followed in the simulations seems to trace the dynamical difference of the cloud sub-structures in different environments.We further observe that, as we trace denser gas, the average value of  vol vir decreases, suggesting that we trace gas which is close to being gravitationally bound.
As a conclusion to the virial analysis, we briefly remark on the expected behaviour of the surface terms which also contribute to the full virial ratio.Since we find that for most of the structures whose dynamics is important for eventual star formation (both large scale, as well as denser), the kinetic and the potential energy dominate their dynamics, we expect that their corresponding surface terms could be important as well.We perform a detailed anaylsis on the tidal effect of the external medium in Section 5.

Inflow or infall
The GHC scenario suggests that the high line-widths seen in MCs are a result of infalling gas.On the other hand, if gas is swept up due to larger scale flows, it could also result in possible signs of inflowing gas.To investigate this, we compute the radial velocity as where r 0 is the position of the centre of mass of a given structure.We show the ratio of the radial velocity of the structures over their 1-dimensional velocity dispersion against  in Fig. 10.The color bar here represents  vol vir , with bound structures ( vol vir < 1) marked with an additional black outline.The top row corresponds to all sub-structures contained in the two hydrodynamic clouds at  evol = 2 (left) and 3.5 Myr (right), while the bottom row represents the MHD sub-structures.A negative value of   suggests that the gas inside the structure has motion towards the centre of mass while a positive value suggests the reverse.If the high kinetic energy in our dendrogram structures originates from inward motions due to gravitational collapse, we would expect more bound structures to have lower (more negative)   / 1D values.From Fig. 10, we do indeed see that the marginally bound points (reddish) tend to have a negative velocity, especially at larger scales, although the trend is far from clear.In order to investigate which effect is dominating (gravitational infall or large-scale, compressive inflows), we perform a correlation analysis between  vol vir and   / 1D .This can be seen in Table 4, where we compute the Kendall  parameter at 2 and 3.5 Myr for all (both HD and MHD) large-scale ( > 10 pc) as well as bound and marginally bound ( vol vir < 2) sub-structures.In each case, we give the Kendall  coefficient, k  , and the p-value.k  quantifies the rank correlation between two quantities: +1 denotes perfect positive correlation and -1 perfect negative correlation (Kendall 1938).
Firstly, we note that for all structures, there is a statistically significant, mild positive correlation (0.19 at 2 Myr and 0.21 at 3.5 Myr) between   / 1D and  vol vir .Since most of the sub-structures are unbound, this reflects rather the fact that structures with higher  vol vir tend to have higher dispersive motions.
For large-scale structures, at 2 Myr, we do not find a statistically significant correlation (p-value of 0.46).However, at 3.5 Myr, there is a relatively strong anti-correlation (k  = −0.43).This implies that a more negative   / 1D is associated with a higher  vol vir , i.e. at large scales the more unbound structures tend to feel a stronger inwards motion.This is contrary to our expectations of a gravitational infall, and rather suggest an inflow, where the stronger compression is associated with a higher kinetic energy, and consequently, a higher virial ratio.
Over all structures with  vol vir < 2, we find an even milder or almost no correlation at both earlier and later times (k  =-0.14 at 2 Myr and -0.07 at 3.5 Myr).The mild anti-correlation here is likely a result of the anti-correlation seen for the larger scale structures, which are mostly marginally bound.

Gravo-turbulence vs global hierarchical collapse
The gradual formation of dense, gravity dominated branches in Larson's  1D - relation (Fig. 4, see also Appendix A), the emergence of gravitationally bound structures over time (Figs.7 and 8, bottom panel), the decrease in  vol vir for large-scale unbound and marginally bound structures (Table 3), and the significant negative correlation between   / 1D and  vol vir at large scales (Table 4) suggest that the dense molecular gas forms from the rarefied medium being swept up by large scale flows, likely originating from the supernovae being driven in the original SILCC simulations.We therefore conclude that our findings support the GT over GHC scenario of structure formation in our simulated clouds.Our results are in line with the new review by Hacar et al. (2022) on filamentary structures and their fragmentation, which provides independent evidence against GHC.In our companion paper (Ganguly et al., in prep.), we suggest that the MCs in our simulations rather form as sheet-like structures tracing the shells of expanding or interacting supernova bubbles, in agreement with the bubble-driven filament formation scenario of MCs The points represent the mean of log 10  vol vir , while the error bars represent the error on the mean.There is a clear distinction between the average behaviour of these ratios between atomic structures on the one side, and molecular as well as dense structures on the other side.This suggests that the molecular content helps us distinguish between two kinds of structures -unbound and mostly atomic, and marginally bound and mostly molecular.(Koyama & Inutsuka 2002;Inoue & Inutsuka 2009;Inutsuka et al. 2015).
In similar simulations of forming MCs in a supernova-driven, stratified, turbulent medium, Ibáñez-Mejía et al. ( 2017) have found that their clouds remain overall gravitationally bound, and compres-sions due to supernovae are insufficient to drive the turbulence in the dense molecular gas.This is in contrast to our results.While it is difficult to compare the simulations directly, we hypothesise that the difference between the results could lie in the technique of forming dense structures.Ibáñez-Mejía et al. ( 2017 for all, large scale, and bound structures at different times.In each case, the  statistic,   , and the corresponding p-value is shown.Overall, the structures show a weak correlation between   / 1D and  vol vir , i.e. a larger  vol vir leads to a larger outwards velocity.However, the large scale structures experiencing inflow show an anti-correlation at  evol = 3.5 Myr, i.e. less bound structures on the large scales have a higher inflow velocity, suggesting compression. the beginning and model gravitational interactions after the inception of the clouds.In contrast, our clouds include self-gravity from the very beginning and therefore form more self-consistently.This perhaps suggests that modelling gravitational interactions, including tidal effects, in the rare gas is important for consistently modelling the velocity structures and dynamics of the denser sub-structures that form later on.Once the dense, molecular gas has been formed, additional external nearby supernovae (at distances less than 50 pc) only drive turbulence for a short time (Seifried et al. 2018).
The role of supernovae in creating a turbulent cascade has been extensively discussed as a possible origin of the Larson-like scaling relations in simulated clouds (e.g.Kritsuk et al. 2013;Padoan et al. 2016).The role of the combined effect of supernova-driven turbulence and self-gravity has also been suggested with relation to the observed relations (see e.g.Izquierdo et al. 2021).Our findings are in agreement with these results.

Comparison of acceleration terms
In the analysis done so far, we have concentrated only on the selfgravity of individual structures.Due to the chaotic nature of the medium, the clumpy surrounding medium contributes to an asymmetric gravitational pull on any given structure, and thereby influences its dynamics.The gravity of the surrounding medium causes tidal effects within each structure, which can result in both, shear motions to tear a structure apart, as well as tidal compression to enhance its collapse along certain directions.Note that tidal compression due to the external medium can never be compressive from all directions, and rather deforms the mass distribution.
We highlight the relation between the acceleration due to selfgravity and the surrounding medium in Fig. 11.We consider a structure embedded in a more diffuse medium.For each point inside this structure, the total gravitational acceleration, g tot , can be computed as where Φ tot is the global gravitational potential due to all of the gas in the simulation and is obtained directly from the simulation data.
We estimate the acceleration due to self-gravity, g int , by a KD tree algorithm using only the mass contained inside the structure.We can then obtain the acceleration solely due to the surrounding medium by the vector subtraction of these two quantities.
Note that g tot (and as a consequence g ext ) here is caused by the gas distribution, and does not include the galactic potential.The scale height of the galactic potential due to old stars in the simulations is A given structure mass M and centre of mass velocity v 0 is embedded in a more diffuse surrounding medium.For each point inside the structure, we compute the net gravitational acceleration g tot as the negative gradient of the global potential Φ tot , and the self gravitational acceleration g int based on the mass distribution of the structure itself using a KD tree.The effect of solely the surrounding medium is then the vectorial subtraction of these two terms according to Eq. 20.
100 pc, much larger than our largest structures.We therefore neglect its effect.If self-gravity is the only important term, then the g int and g tot terms should be comparable, and the angle between them should peak at a value close to 0. The large-scale structures should also peak near zero, by virtue of the fact that there simply is little gas outside the structures to exert significant tidal forces.
We show the distribution of angles between these two terms, denoted by cos −1 ( ĝtot • ĝint ), as PDF computed over all cells for each dendrogram structure of one example cloud (MC1-HD at 3.5 Myr) in Fig. 12.Here, ĝtot and ĝint are the unit vectors of g tot and g int , respectively, i.e.
The histogram of each individual dendrogram structure is shown by a grey line.If the angle between these two vectors is completely random, as could be the case if self-gravity is not at all important, then the histogram of angles should follow a distribution P() ∝ sin().This is denoted by the dashed cyan line.The average behaviour for atomic, molecular, and dense molecular ( thr > 10 −20 g cm −3 ) substructures for this particular cloud is shown in red, blue and yellow, respectively.
From the individual grey outlines, we see that for many structures, The average for the atomic, molecular, and dense molecular structures is shown in red, blue, and yellow, respectively.We see that a preferential alignment between the two vectors is seen prominently for molecular structures.
the angle between the two vectors is not random, and there is some preferential alignment.A few of the structures do indeed peak near 0, while there are also a significant number of structures where the angle distribution is relatively random.A clear distinction is seen when calculating the average behaviour for the molecular structures and the atomic structures, separately.The atomic structures show a behaviour which is in good agreement with the random line, suggesting that they are only marginally affected by their self-gravity.In contrast, the average behaviour of the molecular structures clearly shows that self-gravity plays a more pivotal role in their evolution.We find consistent results when looking at the median angle distribution against  vol vir , for all the different cloud structures in Fig. 13 (left panel).The median angle for a given sub-structure here is computed by obtaining the median for cos −1 ( ĝtot • ĝint ) over all cells in the structure.The horizontal red band represents the marginally bound region (1 ≤  vol vir < 2), while the vertical dotted line represents the expected median for a completely random alignment, i.e. 90 degrees.We see a clear correlation between the median angle and the virial ratio, with smaller median angles corresponding to smaller virial ratios.Interestingly, there seem to be some structures, which are marginally bound but still have a relatively high median angle.These could possibly represent structures that are at least partially compressed due to tidal forces.
We find a similar picture if we compare the density weighted mean of the ratio of magnitudes of g int and g ext for each structure.This is shown in the right panel of Fig. 13, where we plot the log of the average of We compute this by: The colorbar here represents  vol vir , with bound structures ( vol vir < 1) marked with an additional black outline.We find that large scale (>10 pc) structures experience only a small effect from external gravity, owing to the fact that they represent almost the entire cloud and there is little mass outside to impart significant tidal influence.At smaller scales, we do continue to have structures where self-gravity is far more important compared to the gravity of the surrounding medium, but we also find a large number of (predominantly atomic, but also molecular) structures where the external gravity can play an important role.There is also a clear trend that more bound structures tend to have stronger self-gravity, and therefore experience less influence due to the external medium.It is important to note that this analysis compares only gravity terms, and therefore does not tell us about the overall importance of gravity compared to other forces (see Section 4).
Finally, we attempt to directly compute the energy associated with the external gravitational field g ext .We estimate this from The comparison of  ext PE with the self-gravitating potential energy and kinetic energy at  evol = 3.5 Myr is seen in Fig. 14 for all sub-structures, both HD and MHD.The x-axis represents the size  of structure, while the colorbar plots the density threshold.The horizontal dotted line represents a value of 1 in the ratio for both panels.We see that  ext PE is smaller than the kinetic (right panel) or the self-gravitating potential energy (left panel) for almost all structures.This, however, hides the potentially dynamic effect tidal forces can have.If a structure is being partially compressed on one side and torn apart on the opposite side, in terms of energy these two effects would tend to cancel each other but the structure is actually being deformed.This can be captured by performing a more nuanced tidal analysis, which we do in the following.

Tidal analysis
The gravity of the external medium is important for a significant number of structures in our clouds.It is, however, unclear what exact effect the external gravity has.The structures in the simulations live in a chaotic environment, and therefore the net effect of the tides can cause shear motions, disrupting their possible collapse or hindering compression.Conversely, tidal effects may cause at least partial compression.To quantify the effect of the gravity of the medium surrounding the structure, we perform an analysis based on the tidal tensor (see e.g.Renaud et al. 2009Renaud et al. , 2011)).
For a gravitational potential field Φ, the tidal tensor T is defined such that The eigenvalues of the tidal tensor encode the information related to deformation (compression/extension) occurring due to gravity at a certain point due to the local gravitational field.Because the tidal tensor is symmetric and real-valued, it can be written in orthogonal form.The diagonalized tidal tensor is given as, where   are its eigenvalues.The trace of the tidal tensor Tr(T) = 3 =1   contains the local density information in the Poisson equation: This implies that the trace is 0 if the point at which the tidal tensor is evaluated is outside the mass distribution.The virial parameter  vol vir of the structures, plotted against the median angle distribution between the self gravitational acceleration g int and the total gravitational acceleration g tot , for all clouds, both HD and MHD at  evol =3.5 Myr.The red shaded region shows the region between marginally bound and virialized, while the vertical dotted line shows the expected value for a completely random alignment between g int and g tot .More bound structures have a preferential alignment between g int and g tot .Right: Average ratio of g int to the external gravitational acceleration g ext , plotted for each structure against its size, for all clouds at  evol =3.5 Myr.The colorbar shows the virial ratio  vol vir , and gravitationally bound structures are marked with a black circle.For the unbound, and mostly but not exclusively atomic structures, the external gravity seems to be more important compared to the self-gravity of the structures.PE to the self gravitation potential energy (left) and to the kinetic energy (right) for all clouds, both HD and MHD, at  evol =3.5 Myr.The potential energy due to the external medium is generally less for almost all structures compared to both potential and kinetic energy.This, however, masks compression and extension happening to the same structure due to the external medium, which would cancel out in the energy term but still deform the structure.
For example, for a point of mass , the tidal tensor at distance  is given as: The sign of the eigenvalues indicates whether the given mode is compressive or extensive, while the magnitude represents the strength of the respective compressive/extensive mode (see e.g.Renaud et al. 2009).A positive eigenvalue   > 0 implies that a clump of gas will expand along that particular eigen-direction due to the local gravitational field, while a negative eigenvalue of   < 0 implies the opposite.This does not of course automatically imply that the structure will successfully collapse or disperse along the given eigendirection.The present analysis is done only on the gravitational field, and reflects the nature of gravitational deformation, in absence of all other resisting forces such as thermal or magnetic pressure.
For the purpose of our analysis, we are interested if an entire structure will compress/extend due to the gravitational field.The average tidal tensor for an entire structure can be computed as the volume average of the tidal tensor at each point: Similar to splitting the gravitational acceleration vector in the previous section, we can split the tidal tensor also into three parts: the tidal tensor due to only the matter inside the structure, T int , due to only matter distribution external to the structure, T ext , and their sum due to the entire matter distribution: Here, T tot represents the net deformation due to both, the structure itself and its surroundings, while T ext represents the deformation introduced solely due to the external medium.The trace of the three different tidal tensors are as follows: Tr( T tot ) = −4  avg ( 30) where  avg is the average density computed in Eq. 10.We show in the Appendix B that the traces are retrieved.An important consequence of Eqs.30-32 is that T ext must contain both compressive and extensive (disruptive) modes, corresponding to negative and positive   respectively.In contrast, T int and T tot must contain at least one compressive mode, but might also be fully compressive.
In Fig. 15, we plot the ratio of the maximum to the minimum absolute eigenvalue of the total tidal tensor | max tot |/| min tot | against the virial ratio of the structures.These are computed as Since the eigenvalues can become negative, we consider their absolute values in order to disentangle the relative importance of tidal compression ( < 0) vs extension ( > 0).Fully compressive deformations (all three   < 0) are shown in red, while those with at least one extensive mode are shown in cyan.The two vertical dashed lines represent  vol vir = 1 and  vol vir = 2.The top and right side-panels show the fraction f fraction of structures against  vol vir and | max tot |/| min tot |, respectively, for both fully compressive and partially extensive structures.
We first note that most structures have three compressive modes, suggesting that the inclusion of the surrounding clumpy medium, while modifying the nature of gravitational interaction, still results in compressive deformation due to the self-gravity.Note that this analysis is done on the gravity alone, and does not therefore imply that the structure is actually collapsing, only that this is the net effect the overall gravitational field will achieve.
We further find that the compression due to gravity is highly anisotropic, and on average appears to be more anisotropic (higher | max tot / min tot |) for structures with at least one extensive mode (Fig. 15, right panel histogram).A higher anisotropy indicates that gravity is attempting to flatten or elongate these structures.Partially extensive structures also seem to experience such an elongating effect more strongly compared to fully compressive structures.
From Fig. 15, top panel histogram, we also see that structures which have at least one extensive mode, also tend to be more unbound.This raises the intriguing possibility that perhaps their high  vol vir values are even generated by the highly anisotropic deformation introduced due to tides, by converting tidal energy into kinetic energy.
To investigate this possible scenario, and to quantify the relevance of the tidal force compared to other terms (such as self-gravity and turbulence), we perform a further timescale analysis based on the tidal tensor.For this purpose, we are interested in quantifying the relevant timescales and energies introduced by the external medium, and therefore use T ext for our analysis.
The eigenvalues of the tidal tensor have units of [time] −2 , and can therefore directly be converted into a timescale.Assuming that the largest eigenvalue dominates the deformation of a given structure, we can define a tidal timescale which represents the typical deformation timescale of the given structure solely due to the external gravitational field as follows: where | max ext | is the maximum of the absolute eigenvalues | ,ext |, similar to Eq. 33.
To compare this deformation timescale with the typical timescale of the gravity and turbulence, we use the free-fall time,  ff , and the crossing time,  crossing , respectively.The free fall time represents the timescale over which a uniform density spherical structure would collapse in the absence of any pressure forces, solely due to its own self-gravity: If the ratio  ext tidal / ff >> 1, this implies that the self-gravity of the structure dominates and acts on a much shorter timescale compared to any deformation due to external tidal forces.In contrast, if  ext tidal / ff << 1, then tidal deformation is important in terms of modifying any possible collapse of a given structure.
Similarly, the crossing timescale is computed as where  is the shortest axis of the fitted ellipsoid (Eq.4).The crossing timescale represents the timescale over which supersonic turbulence is established throughout the medium.In the absence of driving forces, supersonic turbulence is also expected to decay over  crossing (Mac Low et al. 1998;Stone et al. 1998).If  ext tidal / crossing >> 1, the turbulence dissipates on a much shorter timescale compared to the tidal deformation timescale.On the other hand, if  ext tidal / crossing << 1, the crossing timescale is much longer compared to the tidal timescale, and gravitational deformation occurs on a timescale short enough to possibly generate the high amount of kinetic energy seen in our energy analysis.
The behaviour of the two ratios against the virial ratio of the different structures is shown in Fig. 16 (top:  ext tidal / ff , middle:  ext tidal / crossing ).In both plots, the colorbar represents the size of the structures.The vertical dotted lines show the boundary between virialized and marginally bound.The circles with a black outline indicate structures which have at least one extensive  tot  (corresponding to the cyan points in Fig. 15).
We see that the ratio  ext tidal / ff becomes slightly higher with a lower  vol vir , and therefore self-gravity becomes more and more important as we go to more bound structures.This agrees with the previous picture in Fig. 13, where the external gravity was magnitude-wise less important for larger scale, as well as more bound structures.Structures with at least one extensive mode have tidal timescales comparable to or slightly shorter than the free fall time.
The ratio  ext tidal / crossing has quite a different behaviour (Fig. 16, middle panel).For almost all structures, but especially for unbound structures, the crossing timescale is much shorter compared to the tidal timescale.This suggests that turbulence dissipates on a much shorter timescale compared to the gravitational deformation timescale introduced by tides and as such, tides cannot be the principal contributor to the high kinetic energy in the structures.
The difference in this average behaviour can be seen in a histogram of these two ratios in Fig. 16, bottom panel.Here, the distribution of  ext tidal / ff and  ext tidal / crossing are plotted as a PDF in red and cyan, respectively.The vertical dotted lines represent the mean of each distribution.We see that the distribution of  ext tidal / crossing is shifted to the right compared to the distribution of  ext tidal / ff , and in both cases the mean and most of the structures lie above a ratio of 1. Overall we find a scenario where  crossing <  ext tidal for most of the structures, while  ff is more comparable to  ext tidal .This suggests that turbulence acts on a short enough timescale to modify the structures as they evolve.

CONCLUSIONS
We perform a detailed energetic analysis of seven different simulated molecular clouds (5 with magnetic fields and 2 without) from the SILCC-Zoom simulations (Seifried et al. 2017).In these simulations, we study the evolution of the multi-phase interstellar medium in a supernova-driven, stratified galactic disc environment.We identify structures forming inside each cloud using a dendrogram algorithm, and trace the evolution of their statistical properties over 1.5 Myr during the early stages of cloud evolution before stellar feedback complicates the picture.We include a simple chemical network which allows us to follow the formation of H 2 as the cloud assembles and as such, distinguish between mostly atomic (H 2 mass fraction < 50%) and mostly molecular (H 2 mass fraction > 50%) structures.
• We find that our clouds show a Larson-like power-law behaviour with significant systematic deviations for the densest dendrogram branches in both  1D -size and mass-size relations (Fig. 4).We further see that these deviations evolve over time from a roughly Larson-like law, suggesting that they emerge as gravity becomes more and more important.We observe that all sub-structures can be divided into roughly two kinds depending on their behaviour in the Heyer plot (Fig. 6); denser and molecular structures that roughly follow the Heyer line, and less dense and mostly atomic structures that show no trend with surface density.
• In terms of energetics, we find that the dynamics of large scale (≥ 10 pc) structures, as well as smaller scale denser structures is primarily governed by the interplay between turbulence and gravity (Figs. 7 and 8).Thermal and magnetic energies are dynamically important only for more diffuse, mostly atomic structures.Our results are based on an analysis of the volume energy terms.
• From a virial analysis, we find that on the large scales, the dendrogram structures are unbound or marginally bound at earlier times, but become much closer to virialized over time.Denser, potentially star-forming structures emerge over time, suggesting that they are created by compression in a turbulent, marginally bound environment.
By performing a correlation analysis between the virial parameter and radial velocity tracing the nfall or inflow for each structure, we find that on the larger scales a more unbound structure is likely to have a stronger inflowing motion.This suggests that we see signs of supernova driven compression, rather than gravitational infall and is in agreement with the bubble-driven filament formation scenario proposed by Inutsuka et al. (2015).
• Finally, we attempt to assess the importance of gravitational tidal forces and their role in the dynamics of cloud sub-structures.We compare a given structure's self-gravity and the gravity due to its inhomogeneous surrounding medium by comparing the magnitude of, as well as the angle between, the two respective acceleration terms (Fig. 13).We find that the surrounding medium has some influence on a number of smaller scale unbound structures.Energetically, we find that the external gravitational energy is smaller in magnitude compared to both self-gravity and kinetic energy (Fig. 14).
To quantitatively evaluate the nature of the tidal field, we perform an analysis based on the tidal tensor.We find that the tidal effects are mostly compressive (Fig. 15).However, the gravitational field is highly anisotropic and even more so for structures which have extensive modes.Based on a timescale analysis, we find that timescales related to tidal deformation are generally longer compared to the turbulence decay timescale, suggesting that tidal deformations cannot be the primary source of generating kinetic energy in the structures (Fig. 16).Howeverm the tidal timescale is rather comparable to the freefall time, implying that the structures will deform as they collapse.
Overall, we find a structure formation scenario consistent with the gravo-turbulent scenario of structure formation, with bound structures emerging over time from a largely unbound or marginally bound medium.Magnetic and thermal energy seem to play a subservient role compared to turbulence and self-gravity, with gravitational tides modifying the nature of gravitational compression and leading to formation of more anisotropic, elongated structures.

Figure 1 .
Figure 1.Column densities of the analyzed clouds, each projected along the −axis at  evol =3.5 Myr.The side length of each plot is 62.5 pc.The MHD clouds (top right, and bottom row) have elongated filamentary structures with extended diffuse envelopes, while the hydrodynamic clouds (top left) appear to be more clumpy and fragmented.Note that MC3-MHD has been excluded in further analysis due to its lack of discerning structures and low molecular gas content (see Table1).

Figure 3 .
Figure 3. Left to right: Projections of MC5-MHD at  evol =3.5 Myr along the -, -, and -axis, respectively.The contours show the projections along the same axis of the leaf dendrogram structures.Molecular structures (> 50% H 2 ) are plotted with solid lines, and atomic structures (< 50% H 2 ) are plotted with dashed lines.The molecular structures nicely trace the dense spine of the cloud, while the atomic structures mostly represent the envelope.
) -roughly consistent with virialization(Solomon et al. 1987), turbulent cascade(Padoan et al. 2016), as well as gravitational collapse (Ballesteros-Paredes et al. 2011).More recently, for a single molecular cloud, a much flatter velocity dispersion-size relation has been observed when studying structures with high column density.

Figure 4 .
Figure 4. Time evolution of the velocity dispersion-size relation for HD clouds at 1.5, 2 and 3.5 Myr (top row, from left) and for MHD clouds at 2, 2.5 and 3.5 Myr (bottom row, from left).The colour bar represents the threshold (minimum) density inside each structure.Plotted as red, dashed and black, dotted lines are the best fit and the  ∝  1/2 obtained bySolomon et al. (1987), respectively.We show the hydrodynamic clouds at an earlier time to highlight that they behave similar to the MHD clouds.Some denser branches trail off to the left, showing complete departure from a Larson-like power law.These dense branches are absent or less pronounced at earlier times.

Figure 5 .
Figure 5. Same as in Fig. 4 but now for the mass-size relation.The red, dashed line represents the best-fit, while the black, dash-dotted and the black, dotted lines represent  ∝  and  ∝  2 , respectively.The denser branches follow roughly  ∝ , while the atomic structures seem to be more consistent with  ∝  2 .

Figure 6 .
Figure 6.Same as in Fig. 4 but now for the Heyer relation.The dash-dotted line represents the ( Σ/5) 1/2 line from Heyer et al. (2009) (corresponding to the virial parameter  = 1), while the parallel dotted lines represent other possible values of  in a purely kinetic virial analysis.We find two populations of structures: denser and mostly molecular structures following the Heyer line towards higher Σ, and mostly atomic structures that show no correlation with surface density.

Figure 7 .
Figure7.Energetics of the two HD clouds at  evol = 2 Myr (left column) and  evol = 3.5 Myr (right column).Each plot is a ratio of a different energy to the potential energy.The red shaded region represents the region between marginally bound and virialized (1 ≤  < 2) for that particular energy term.Large scale structures are marginally bound in the virial analysis.At smaller scales, most of the structures are kinetically dominated with only a few marginally bound structures, as well as a small number of bound structures.

Figure 8 .
Figure 8. Same as in Fig 7 but now for the MHD clouds.Additionally  B / | PE | is plotted in the top row.The red shaded region represents the region between marginally bound and virialized (1 ≤  < 2) for that particular energy term.Note that due to Eq. 15, this implies a factor of 2 difference between the energy ratios and  vol vir , except for  B / | PE |.Large scale structures are mostly marginally bound in the virial analysis.Gravitationally bound structures are few and emerge over time, suggesting their origin is from the dynamics of the marginally bound gas.

Figure 9 .
Figure 9.Time evolution of the average behaviour of  vol vir for HD (left) and MHD clouds (right).The points represent the mean of log 10  vol vir , while the error bars represent the error on the mean.There is a clear distinction between the average behaviour of these ratios between atomic structures on the one side, and molecular as well as dense structures on the other side.This suggests that the molecular content helps us distinguish between two kinds of structures -unbound and mostly atomic, and marginally bound and mostly molecular.

Figure 10 .
Figure 10.Radial velocity over 1D velocity dispersion plotted against the size for HD clouds at  evol =2 Myr (top left) and  evol =3.5 Myr (top right), and for MHD clouds at  evol =2 Myr (bottom left) and  evol =3.5 Myr (bottom right).A negative value means that the structure has overall a radial velocity towards its centre of mass (i.e.inflowing).The color-bar represents the boundness of the structures, with reddish structures representing  vol vir < 2. The bound structures ( vol vir < 1), corresponding to points below the horizontal black bar in the color bar) are marked with black outlines.Large scale marginally bound structures show signs of inflowing velocity,

Figure 11 .
Figure11.Sketch representing the internal and external gravitational acceleration.A given structure mass M and centre of mass velocity v 0 is embedded in a more diffuse surrounding medium.For each point inside the structure, we compute the net gravitational acceleration g tot as the negative gradient of the global potential Φ tot , and the self gravitational acceleration g int based on the mass distribution of the structure itself using a KD tree.The effect of solely the surrounding medium is then the vectorial subtraction of these two terms according to Eq. 20.

Figure 12 .
Figure12.Histogram of angles between the self gravitational acceleration g int and the total gravitational acceleration g tot due to gas both inside and outside the structure, for all different dendrogram structures of the cloud HD 1 at  evol =3.5 Myr.The histogram for each individual structure is plotted with a grey line.The cyan dashed line shows how the distribution should be if these two vectors were randomly aligned (i.e. if self-gravity plays no role).The average for the atomic, molecular, and dense molecular structures is shown in red, blue, and yellow, respectively.We see that a preferential alignment between the two vectors is seen prominently for molecular structures.

Figure 13 .
Figure13.Left: The virial parameter  vol vir of the structures, plotted against the median angle distribution between the self gravitational acceleration g int and the total gravitational acceleration g tot , for all clouds, both HD and MHD at  evol =3.5 Myr.The red shaded region shows the region between marginally bound and virialized, while the vertical dotted line shows the expected value for a completely random alignment between g int and g tot .More bound structures have a preferential alignment between g int and g tot .Right: Average ratio of g int to the external gravitational acceleration g ext , plotted for each structure against its size, for all clouds at  evol =3.5 Myr.The colorbar shows the virial ratio  vol vir , and gravitationally bound structures are marked with a black circle.For the unbound, and mostly but not exclusively atomic structures, the external gravity seems to be more important compared to the self-gravity of the structures.

Figure 14 .
Figure 14.Ratio of the magnitudes of external gravitational energy  extPE to the self gravitation potential energy (left) and to the kinetic energy (right) for all clouds, both HD and MHD, at  evol =3.5 Myr.The potential energy due to the external medium is generally less for almost all structures compared to both potential and kinetic energy.This, however, masks compression and extension happening to the same structure due to the external medium, which would cancel out in the energy term but still deform the structure.

Figure 15 .
Figure 15.The ratio of the largest to the smallest eigenvalue of the tidal tensor T tot plotted against  vol vir of the structures, for all clouds at  evol =3.5 Myr.The two vertical lines denote the marginally bound region (1 ≤  vol vir < 2).The red symbols represent fully compressive structures (all   < 0), while the cyan ones are structures which have at least extensive mode.The side-panes show the 1D distribution of the fraction of structures based on the distribution along that axis, and are normalized by the total number of structures.The  max tot / min tot ratio represents the degree of anisotropy in the gravitational field.For most structures, the overall deformation due to gravity is still fully compressive, suggesting that disruptive tidal effects are not dominant.Structures with at least one extensive mode seem to experience a more anisotropic (higher | max tot |/ | min tot |) gravitational field.They are also more likely to be more unbound, while marginally or gravitationally bound structures are almost always fully compressive, with a few exceptions.

Figure 16 .
Figure 16.Ratio of the tidal deformation timescale due to T ext , to the gravitational free fall timescale (top) and the crossing time (middle), representing the turbulence decay timescale, plotted against the virial parameter  vol vir of the substructures, for all clouds at  evol =3.5 Myr.The colorbar represents the size of the structures.Structures which have at least one extensive mode (one  tot  > 0) are marked with a black circle.The vertical dotted lines represent the marginally bound region 1 ≤  volvir < 2 for the structures.The top plot shows that self-gravity is becoming more important compared to the tidal deformation timescale on average, for more bound structures.The bottom plot shows that the crossing time is shorter than the tidal deformation timescale for almost all structures, suggesting that tidal deformations due to the surrounding medium cannot be the primary source of turbulent energy.The histogram of the distribution of the two ratios,  ext tidal / crossing and  ext tidal / ff , with the means plotted in dotted vertical lines (bottom).

Figure A1 .
Figure A1. 1D - relation for all MHD sub-structures, at  evol = 2 (left) and 3.5 Myr (right).The colour bar represents  vol vir , with structures with  vol vir < 1 being marked with an additional black outline.

Figure B1 .
Figure B1.Trace of the averaged total tidal tensor T tot (left) and trace of the averaged external tidal tensor T ext (right), plotted against the average density, for all cloud sub-structures at  evol = 3.5 Myr.The dashed line represents a one to one line, and the green shaded region represents a factor of 2 in each direction.The trace of T tot well reflects the average density.The trace of T ext should ideally be 0, but is typically only 1% of the 1:1 line.

Table 2 .
Information on the number of molecular structures obtained in the dendrogram analysis.

Table 4 .
The Kendall  correlation values between   / 1D and  vol vir ) have no self-gravity in