# Cloud properties across spatial scales in simulations of the interstellar medium

Tine Colman, Noé Brucy, Philipp Girichidis, Simon C. O. Glover, Milena Benedettini, Juan D. Soler, Robin G. Tress, Alessio Traficante, Patrick Hennebelle, Ralf S. Klessen, et al.

► To cite this version:

Tine Colman, Noé Brucy, Philipp Girichidis, Simon C. O. Glover, Milena Benedettini, et al. Cloud properties across spatial scales in simulations of the interstellar medium. Astronomy and Astrophysics - A&A, 2024, 686, pp.A155. 10.1051/0004-6361/202348983. cea-04605689

HAL Id: cea-04605689

https://cea.hal.science/cea-04605689

Submitted on 7 Jun 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Cloud properties across spatial scales in simulations of the interstellar medium

Tine Colman, Noé Brucy, Philipp Girichidis, Simon C. O. Glover, Milena Benedettini, Juan D. Soler, Robin G. Tress, Alessio Traficante, Patrick Hennebelle, Ralf S. Klessen, et al.

► To cite this version:

Tine Colman, Noé Brucy, Philipp Girichidis, Simon C. O. Glover, Milena Benedettini, et al. Cloud properties across spatial scales in simulations of the interstellar medium. Astronomy & Astrophysics - A&A, 2024, 686, pp.A155. 10.1051/0004-6361/202348983. cea-04605689

HAL Id: cea-04605689

https://cea.hal.science/cea-04605689v1

Submitted on 7 Jun 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Cloud properties across spatial scales in simulations of the interstellar medium

Tine Colman<fnref n="1" />, Noé Brucy<fnref n="2" />, Philipp Girichidis<fnref n="2" />, Simon C. O. Glover<fnref n="2" />, Milena Benedettini<fnref n="3" />
Juan D. Soler<fnref n="3" />, Robin G. Tress<fnref n="4" />
Alessio Traficante<fnref n="3" />, Patrick Hennebelle<fnref n="1" />, Ralf S. Klessen<fnref n="2" />,<fnref n="5" />
Sergio Molinari<fnref n="3" />, and Marc-Antoine Miville-Deschênes<fnref n="6" />

Received 17 December 2023 / Accepted 01 March 2024

## ABSTRACT

**Context.** Molecular clouds (MCs) are structures of dense gas in the interstellar medium (ISM) that extend from ten to a few hundred parsecs and form the main gas reservoir available for star formation. Hydrodynamical simulations of a varying complexity are a promising way to investigate MCs evolution and their properties. However, each simulation typically has a limited range in resolution and different cloud extraction algorithms are used, which complicates the comparison between simulations.

**Aims.** In this work, we aim to extract clouds from different simulations covering a wide range of spatial scales. We compare their properties, such as size, shape, mass, internal velocity dispersion, and virial state.

**Methods.** We applied the HOP cloud detection algorithm on (M)HD numerical simulations of stratified ISM boxes and isolated galactic disk simulations that were produced using FLASH, RAMSES, and AREPO.

**Results.** We find that the extracted clouds are complex in shape, ranging from round objects to complex filamentary networks in all setups. Despite the wide range of scales, resolution, and sub-grid physics, we observe surprisingly robust trends in the investigated metrics. The mass spectrum matches in the overlap between simulations without rescaling and with a high-mass power-law index of -1 for logarithmic bins of mass, in accordance with theoretical predictions. The internal velocity dispersion scales with the size of the cloud as $\sigma \propto R^{0.75}$ for large clouds ($R \ge 3$ pc). For small clouds we find larger $\sigma$ compared to the power-law scaling, as seen in observations, which is due to supernova-driven turbulence. Almost all clouds are gravitationally unbound with the virial parameter scaling as $\alpha_{vir} \propto M^{-0.4}$, which is slightly flatter compared to observed scaling but in agreement given the large scatter. We note that the cloud distribution towards the low-mass end is only complete if the more dilute gas is also refined, rather than only the collapsing regions.

**Key words.** methods: numerical – ISM: clouds

## 1. Introduction

Historically, (molecular) clouds have been considered important building blocks in the interstellar medium (ISM). The warm and diffuse gas condenses to form cloud-like structures with masses of ~$10^4–10^6$ M$_{\odot}$, typical sizes of ~10 – 100 pc, mean densities of $n_{H_2}$ ~ 100 cm$^{-3}$, and temperatures of T ~ 10 K (e.g. Solomon et al. 1987; Scoville et al. 1987; Dame et al. 1987, 2001, or Klessen & Glover 2016 for a review). Most of the mass in these conditions is in the form of molecular hydrogen, and hence the term molecular clouds (MCs). However, most of the dynamical properties and the connection to hydrodynamical evolution does not require the gas being molecular. One particularly important

aspect connected to clouds is that they harbour the coldest (~10 K) regions of the ISM including the collapsing cores that eventually form stars.

For a long time, clouds were considered to be stable and long-lived structures that are gravitationally bound and supported against collapse by magnetic fields, for example (e.g. Mouschovias & Spitzer 1976; Shu et al. 1987; Mouschovias & Ciolek 1999). This paradigm has changed. In particular, the idea that clouds are magnetically supported has given way to a picture in which they are formed by a combination of gravity and turbulence (e.g. Ballesteros-Paredes et al. 1999; Ballesteros-Paredes & Mac Low 2002; Mac Low & Klessen 2004; Krumholz & McKee 2005; Ballesteros-Paredes 2006; Hennebelle & Chabrier

<footnote n="1">¹ AIM, CEA, CNRS, Université Paris-Saclay, Université Paris Diderot, Sorbonne Paris Cité, 91191 Gif-sur-Yvette, France
e-mail: tine.colman@cea.fr</footnote>

<footnote n="2">² Universität Heidelberg, Zentrum für Astronomie, Institut für Theoretische Astrophysik, Albert-UEberle-Str. 2, 69120 Heidelberg,
Germany
e-mail: noe.brucy@uni-heidelberg.de; philipp@girichidis.com</footnote>

<footnote n="3">³ INAF - Istituto di Astrofisica e Planetologia Spaziali, via Fosso del Cavaliere 100, 00133 Roma, Italy</footnote>

<footnote n="4">⁴ Institute of Physics, Laboratory for Galaxy Evolution and Spectral Modelling, EPFL, Observatoire de Sauverny, Chemin Pegasi 51, 1290 Versoix, Switzerland</footnote>

<footnote n="5">⁵ Universität Heidelberg, Interdisziplinäres Zentrum für Wissenschaftliches Rechnen, Im Neuenheimer Feld 205, 69120 Heidelberg, Germany</footnote>

<footnote n="6">⁶ Laboratoire de Physique de l'École Normale Supérieure, ENS, Université PSL, CNRS, Sorbonne Université, Université de Paris, 75005 Paris, France</footnote>
# Cloud properties across spatial scales in simulations of the interstellar medium

Tine Colman<fnref n="1" />, Noé Brucy<fnref n="2" />, Philipp Girichidis<fnref n="2" />, Simon C. O. Glover<fnref n="2" />, Milena Benedettini<fnref n="3" />
Juan D. Soler<fnref n="3" />, Robin G. Tress<fnref n="4" />
Alessio Traficante<fnref n="3" />, Patrick Hennebelle<fnref n="1" />, Ralf S. Klessen<fnref n="2" />,<fnref n="5" />
Sergio Molinari<fnref n="3" />, and Marc-Antoine Miville-Deschênes<fnref n="6" />

Received 17 December 2023 / Accepted 01 March 2024

## ABSTRACT

**Context.** Molecular clouds (MCs) are structures of dense gas in the interstellar medium (ISM) that extend from ten to a few hundred parsecs and form the main gas reservoir available for star formation. Hydrodynamical simulations of a varying complexity are a promising way to investigate MCs evolution and their properties. However, each simulation typically has a limited range in resolution and different cloud extraction algorithms are used, which complicates the comparison between simulations.

**Aims.** In this work, we aim to extract clouds from different simulations covering a wide range of spatial scales. We compare their properties, such as size, shape, mass, internal velocity dispersion, and virial state.

**Methods.** We applied the HOP cloud detection algorithm on (M)HD numerical simulations of stratified ISM boxes and isolated galactic disk simulations that were produced using FLASH, RAMSES, and AREPO.

**Results.** We find that the extracted clouds are complex in shape, ranging from round objects to complex filamentary networks in all setups. Despite the wide range of scales, resolution, and sub-grid physics, we observe surprisingly robust trends in the investigated metrics. The mass spectrum matches in the overlap between simulations without rescaling and with a high-mass power-law index of -1 for logarithmic bins of mass, in accordance with theoretical predictions. The internal velocity dispersion scales with the size of the cloud as $\sigma \propto R^{0.75}$ for large clouds ($R \ge 3$ pc). For small clouds we find larger $\sigma$ compared to the power-law scaling, as seen in observations, which is due to supernova-driven turbulence. Almost all clouds are gravitationally unbound with the virial parameter scaling as $\alpha_{vir} \propto M^{-0.4}$, which is slightly flatter compared to observed scaling but in agreement given the large scatter. We note that the cloud distribution towards the low-mass end is only complete if the more dilute gas is also refined, rather than only the collapsing regions.

**Key words.** methods: numerical – ISM: clouds

## 1. Introduction

Historically, (molecular) clouds have been considered important building blocks in the interstellar medium (ISM). The warm and diffuse gas condenses to form cloud-like structures with masses of ~$10^4–10^6$ M$_{\odot}$, typical sizes of ~10 – 100 pc, mean densities of $n_{H_2}$ ~ 100 cm$^{-3}$, and temperatures of T ~ 10 K (e.g. Solomon et al. 1987; Scoville et al. 1987; Dame et al. 1987, 2001, or Klessen & Glover 2016 for a review). Most of the mass in these conditions is in the form of molecular hydrogen, and hence the term molecular clouds (MCs). However, most of the dynamical properties and the connection to hydrodynamical evolution does not require the gas being molecular. One particularly important

aspect connected to clouds is that they harbour the coldest (~10 K) regions of the ISM including the collapsing cores that eventually form stars.

For a long time, clouds were considered to be stable and long-lived structures that are gravitationally bound and supported against collapse by magnetic fields, for example (e.g. Mouschovias & Spitzer 1976; Shu et al. 1987; Mouschovias & Ciolek 1999). This paradigm has changed. In particular, the idea that clouds are magnetically supported has given way to a picture in which they are formed by a combination of gravity and turbulence (e.g. Ballesteros-Paredes et al. 1999; Ballesteros-Paredes & Mac Low 2002; Mac Low & Klessen 2004; Krumholz & McKee 2005; Ballesteros-Paredes 2006; Hennebelle & Chabrier

<footnote n="1">¹ AIM, CEA, CNRS, Université Paris-Saclay, Université Paris Diderot, Sorbonne Paris Cité, 91191 Gif-sur-Yvette, France
e-mail: tine.colman@cea.fr</footnote>

<footnote n="2">² Universität Heidelberg, Zentrum für Astronomie, Institut für Theoretische Astrophysik, Albert-UEberle-Str. 2, 69120 Heidelberg,
Germany
e-mail: noe.brucy@uni-heidelberg.de; philipp@girichidis.com</footnote>

<footnote n="3">³ INAF – Istituto di Astrofisica e Planetologia Spaziali, via Fosso del Cavaliere 100, 00133 Roma, Italy</footnote>

<footnote n="4">⁴ Institute of Physics, Laboratory for Galaxy Evolution and Spectral Modelling, EPFL, Observatoire de Sauverny, Chemin Pegasi 51, 1290 Versailles, Switzerland</footnote>

<footnote n="5">⁵ Universität Heidelberg, Interdisziplinäres Zentrum für Wissenschaftliches Rechnen, Im Neuenheimer Feld 205, 69120 Heidelberg,
Germany</footnote>

<footnote n="6">⁶ Laboratoire de Physique de l'École Normale Supérieure, ENS, Université PSL, CNRS, Sorbonne Université, Université de Paris,
75005 Paris, France</footnote>
2008; Ballesteros-Paredes et al. 2011; Vázquez-Semadeni et al. 2017, 2019), with the magnetic field playing an important but not dominant role (e.g. Crutcher 2012; Girichidis 2021; Whitworth et al. 2023). Turbulence creates over-densities and large-scale coherent structures, whose shapes are similar to the complex shapes of molecular clouds (e.g. Schneider et al. 2011; Ebagezio et al. 2023). As a result, the gas in the ISM is constantly reshaped by turbulence without permanent structures on a crossing time $t_{cross} = L/v$, where L is the cloud size and v the local velocity dispersion on a scale L. In the solar neighbourhood, representative values for L and v are 30 pc and 3 km s⁻¹, respectively, yielding $t_{cross} \sim 10$ Myr, a small fraction of the orbital period of the Galaxy. Clouds can form in regions of converging flows and are dispersed by shear and/or feedback from stars (e.g. Chevance et al. 2023, for a recent review). Therefore, molecular clouds should not be regarded as well-defined discrete entities with clearly identifiable boundaries.

Larson (1981) was the first to establish power-law scaling relations between molecular cloud properties. He observed that cloud mass and size followed the relation $M \propto R^{1.90}$ while the velocity dispersion scaled as $\sigma_v \propto R^{0.38}$. Both of these relations showed significant scatter. In the Milky Way and neighbouring galaxies, MCs seem to adhere to the Larson relations, where the non-thermal line-width increases with cloud size following an approximate power law of index one half (Larson 1981; Solomon et al. 1987; Bolatto et al. 2008; Fukui et al. 2008; Heyer et al. 2009; Roman-Duval et al. 2010; Wong et al. 2011). This relation extends within clouds (Heyer & Brunt 2004) and may be attributed to the power-law scaling expected for compressible turbulence (McKee & Ostriker 2007). Larson (1981) also argued that MCs exhibited similar levels of kinetic and gravitational energy and hence were approximately in virial equilibrium (see also Blitz 1993). This can be quantified through the use of the virial parameter (Bertoldi & McKee 1992)

$$ \alpha_{vir} \equiv \frac{5\sigma_v^2 R}{GM}, \qquad (1) $$

in which $\sigma_v$ is the one-dimensional velocity dispersion, G is the gravitational constant, and R and M are the size and mass of the cloud, respectively. Virial equilibrium corresponds to $\alpha_{vir} = 1$. Although early studies of massive clouds found that $\alpha_{vir} \sim 1$, more recent surveys that are sensitive to a much broader range of cloud masses find a more complex picture (Heyer et al. 2001; Oka et al. 2001; Gratier et al. 2012; Donovan Meyer et al. 2013; Rice et al. 2016; Miville-Deschênes et al. 2017; Colombo et al. 2019; Rosolowsky et al. 2021). These surveys show that low-mass clouds are unbound with virial ratios $\gg 1$ and only massive clouds with masses of order $10^6 M_{\odot}$ or more are marginally bound. However, there are large uncertainties involved in determining the mass of MCs from observations (Szűcs et al. 2016), and so the characteristic mass at which MCs become gravitationally bound has been difficult to determine accurately. Clouds with virial parameters exceeding unity mostly follow a relation with $\alpha_{vir} \propto M^{-1/2}$ (Chevance et al. 2023).

The formation of MCs and the dynamical evolution of the gas within them has been addressed in hydrodynamical simulations with a range of different complexities. Early studies either focussed on simulating isolated clouds (Vazquez-Semadeni et al. 1995; Stone et al. 1998; Klessen & Burkert 2000; Ostriker et al. 2001; Bate et al. 2003; Padoan et al. 2007) or covered large fractions of the interstellar medium without resolving the interior of clouds (de Avillez 2000; de Avillez & Breitschwerdt 2005;

Joung & Mac Low 2006; Kim & Ostriker 2007; Hill et al. 2012; Gent et al. 2013; Walch et al. 2012). However, increasing computational power now allows us to cover both the large-scale environment, including the hot and dilute gas from which clouds condense as well as the coldest collapsing phases. Seifried et al. (2017) showed the importance of embedding the clouds into the turbulent large-scale (~100 pc) environment, in order to accurately follow the turbulent cascade and the formation of cold gas. They adopted a zoom-in approach in which they followed the evolution of individual clouds selected from a large-scale simulation of the ISM with successively improving resolution. This approach allowed them to reach resolutions as small as 0.1 pc, but limited them to modelling only a few clouds, thereby preventing them from drawing robust statistical conclusions. A companion study at a lower resolution of 0.25 pc by Girichidis (2021) complements these efforts. Both models only consider environmental conditions similar to the solar neighbourhood. Similar simulations of stratified boxes but with different turbulent driving recipes and total gas masses have been performed by other groups (e.g. Joung et al. 2009; Simpson et al. 2016; Martizzi et al. 2016; Kim & Ostriker 2017; Brucy et al. 2020, 2023; Colman et al. 2022). These simulations cover slightly larger scales and employ different cooling recipes. In addition, it has now become possible to carry out simulations of entire disk galaxies with a resolution sufficient to resolve individual clouds, allowing one to study the roles of global galactic rotation and differential shear and to cover a large range of local surface densities (e.g. Tress et al. 2020, 2021; Jeffreson et al. 2020). All of these models have their own strengths and weaknesses in terms of physical processes included, environmental properties of molecular clouds and maximum resolution. Therefore, an important question to ask prior to comparing the results of these simulations to observations is whether the predictions of different simulations for the properties of the MCs distribution agree with one other, that is whether these predictions are robust to changes in the numerical approach. Because the definition of a cloud is not trivial, it is important to use the same cloud extraction method with identical parameters for cloud identification, so that we can be sure that any differences we find are due to differences in the simulations and not in the cloud identification algorithm. A thorough comparison between numerical models in overlapping ranges of the cloud masses is crucial for verifying that the full range of cloud masses and sizes can be compared to observations without numerically induced breaks in the scaling relations. Furthermore, a comparison between simulations can be considerably more detailed than a comparison with observations, as we have the full 3D distribution of all relevant physical variables available for our simulated clouds, which is not the case for real, observed MCs.

In this paper, we investigate the molecular cloud properties in a variety of numerical simulations with different resolutions and different numerical recipes for cooling, star formation and stellar feedback. Our goal is to verify the existence of smooth transitions between the simulations and quantify the clouds in terms of global properties that can be transferred to the observational domain. Our study is structured as follows: in Sect. 2 we discuss the algorithm used to extract the clouds and the numerical simulations on which the algorithm was applied. Then, in Sect. 3, we compare the cloud population's masses, sizes and shapes. In Sect. 4 we focus on the internal properties of the clouds, looking at the thermal state, the turbulence and the gravitational stability. We discuss the importance of resolution and of the specifics of the cloud extraction process in Sect. 5 and conclude in Sect. 6.
2008; Ballesteros-Paredes et al. 2011; Vázquez-Semadeni et al. 2017, 2019), with the magnetic field playing an important but not dominant role (e.g. Crutcher 2012; Girichidis 2021; Whitworth et al. 2023). Turbulence creates over-densities and large-scale coherent structures, whose shapes are similar to the complex shapes of molecular clouds (e.g. Schneider et al. 2011; Ebagezio et al. 2023). As a result, the gas in the ISM is constantly reshaped by turbulence without permanent structures on a crossing time $t_{cross} = L/v$, where L is the cloud size and v the local velocity dispersion on a scale L. In the solar neighbourhood, representative values for L and v are 30 pc and 3 km s⁻¹, respectively, yielding $t_{cross} \sim 10$ Myr, a small fraction of the orbital period of the Galaxy. Clouds can form in regions of converging flows and are dispersed by shear and/or feedback from stars (e.g. Chevance et al. 2023, for a recent review). Therefore, molecular clouds should not be regarded as well-defined discrete entities with clearly identifiable boundaries.

Larson (1981) was the first to establish power-law scaling relations between molecular cloud properties. He observed that cloud mass and size followed the relation $M \propto R^{1.90}$ while the velocity dispersion scaled as $\sigma_v \propto R^{0.38}$. Both of these relations showed significant scatter. In the Milky Way and neighbouring galaxies, MCs seem to adhere to the Larson relations, where the non-thermal line-width increases with cloud size following an approximate power law of index one half (Larson 1981; Solomon et al. 1987; Bolatto et al. 2008; Fukui et al. 2008; Heyer et al. 2009; Roman-Duval et al. 2010; Wong et al. 2011). This relation extends within clouds (Heyer & Brunt 2004) and may be attributed to the power-law scaling expected for compressible turbulence (McKee & Ostriker 2007). Larson (1981) also argued that MCs exhibited similar levels of kinetic and gravitational energy and hence were approximately in virial equilibrium (see also Blitz 1993). This can be quantified through the use of the virial parameter (Bertoldi & McKee 1992)

$$ \alpha_{vir} \equiv \frac{5\sigma_v^2 R}{GM}, \qquad (1) $$

in which $\sigma_v$ is the one-dimensional velocity dispersion, G is the gravitational constant, and R and M are the size and mass of the cloud, respectively. Virial equilibrium corresponds to $\alpha_{vir} = 1$. Although early studies of massive clouds found that $\alpha_{vir} \sim 1$, more recent surveys that are sensitive to a much broader range of cloud masses find a more complex picture (Heyer et al. 2001; Oka et al. 2001; Gratier et al. 2012; Donovan Meyer et al. 2013; Rice et al. 2016; Miville-Deschênes et al. 2017; Colombo et al. 2019; Rosolowsky et al. 2021). These surveys show that low-mass clouds are unbound with virial ratios $\gg 1$ and only massive clouds with masses of order $10^6 M_{\odot}$ or more are marginally bound. However, there are large uncertainties involved in determining the mass of MCs from observations (Szűcs et al. 2016), and so the characteristic mass at which MCs become gravitationally bound has been difficult to determine accurately. Clouds with virial parameters exceeding unity mostly follow a relation with $\alpha_{vir} \propto M^{-1/2}$ (Chevance et al. 2023).

The formation of MCs and the dynamical evolution of the gas within them has been addressed in hydrodynamical simulations with a range of different complexities. Early studies either focussed on simulating isolated clouds (Vazquez-Semadeni et al. 1995; Stone et al. 1998; Klessen & Burkert 2000; Ostriker et al. 2001; Bate et al. 2003; Padoan et al. 2007) or covered large fractions of the interstellar medium without resolving the interior of clouds (de Avillez 2000; de Avillez & Breitschwerdt 2005;

Joung & Mac Low 2006; Kim & Ostriker 2007; Hill et al. 2012; Gent et al. 2013; Walch et al. 2012). However, increasing computational power now allows us to cover both the large-scale environment, including the hot and dilute gas from which clouds condense as well as the coldest collapsing phases. Seifried et al. (2017) showed the importance of embedding the clouds into the turbulent large-scale (~100 pc) environment, in order to accurately follow the turbulent cascade and the formation of cold gas. They adopted a zoom-in approach in which they followed the evolution of individual clouds selected from a large-scale simulation of the ISM with successively improving resolution. This approach allowed them to reach resolutions as small as 0.1 pc, but limited them to modelling only a few clouds, thereby preventing them from drawing robust statistical conclusions. A companion study at a lower resolution of 0.25 pc by Girichidis (2021) complements these efforts. Both models only consider environmental conditions similar to the solar neighbourhood. Similar simulations of stratified boxes but with different turbulent driving recipes and total gas masses have been performed by other groups (e.g. Joung et al. 2009; Simpson et al. 2016; Martizzi et al. 2016; Kim & Ostriker 2017; Brucy et al. 2020, 2023; Colman et al. 2022). These simulations cover slightly larger scales and employ different cooling recipes. In addition, it has now become possible to carry out simulations of entire disk galaxies with a resolution sufficient to resolve individual clouds, allowing one to study the roles of global galactic rotation and differential shear and to cover a large range of local surface densities (e.g. Tress et al. 2020, 2021; Jeffreson et al. 2020). All of these models have their own strengths and weaknesses in terms of physical processes included, environmental properties of molecular clouds and maximum resolution. Therefore, an important question to ask prior to comparing the results of these simulations to observations is whether the predictions of different simulations for the properties of the MCs distribution agree with one other, that is whether these predictions are robust to changes in the numerical approach. Because the definition of a cloud is not trivial, it is important to use the same cloud extraction method with identical parameters for cloud identification, so that we can be sure that any differences we find are due to differences in the simulations and not in the cloud identification algorithm. A thorough comparison between numerical models in overlapping ranges of the cloud masses is crucial for verifying that the full range of cloud masses and sizes can be compared to observations without numerically induced breaks in the scaling relations. Furthermore, a comparison between simulations can be considerably more detailed than a comparison with observations, as we have the full 3D distribution of all relevant physical variables available for our simulated clouds, which is not the case for real, observed MCs.

In this paper, we investigate the molecular cloud properties in a variety of numerical simulations with different resolutions and different numerical recipes for cooling, star formation and stellar feedback. Our goal is to verify the existence of smooth transitions between the simulations and quantify the clouds in terms of global properties that can be transferred to the observational domain. Our study is structured as follows: in Sect. 2 we discuss the algorithm used to extract the clouds and the numerical simulations on which the algorithm was applied. Then, in Sect. 3, we compare the cloud population's masses, sizes and shapes. In Sect. 4 we focus on the internal properties of the clouds, looking at the thermal state, the turbulence and the gravitational stability. We discuss the importance of resolution and of the specifics of the cloud extraction process in Sect. 5 and conclude in Sect. 6.
<table><caption>**Table 1.** Set of simulations.</caption><thead><tr><th>Name</th><th>Setup</th><th>Size<br>(kpc)</th><th>Δx<sub>min</sub><br>(pc)</th><th>B<sub>0</sub><br>(μG)</th><th>Code</th><th>References</th><th>No. of<br>snapshots</th></tr></thead><tbody><tr><td>SILCC-1pc-3μG</td><td>Stratified box</td><td>0.5</td><td>0.98</td><td>3</td><td>FLASH</td><td>Girichidis et al. (2018)</td><td>4</td></tr><tr><td>SILCC-1pc-6μG</td><td>Stratified box</td><td>0.5</td><td>0.98</td><td>6</td><td>FLASH</td><td>Girichidis et al. (2018)</td><td>4</td></tr><tr><td>SILCC-0.5pc-3μG</td><td>Stratified box</td><td>0.5</td><td>0.49</td><td>3</td><td>FLASH</td><td>Girichidis (2021)</td><td>6</td></tr><tr><td>SILCC-0.25pc-3μG</td><td>Stratified box</td><td>0.5</td><td>0.24</td><td>3</td><td>FLASH</td><td>Girichidis (2021)</td><td>2</td></tr><tr><td>LS-no-driving</td><td>Stratified box</td><td>1</td><td>0.24</td><td>7.6</td><td>RAMSES</td><td>Colman et al. (2022)</td><td>1</td></tr><tr><td>LS-weak-driving</td><td>Stratified box</td><td>1</td><td>0.24</td><td>7.6</td><td>RAMSES</td><td>Colman et al. (2022)</td><td>1</td></tr><tr><td>LS-medium-driving</td><td>Stratified box</td><td>1</td><td>0.24</td><td>7.6</td><td>RAMSES</td><td>Colman et al. (2022)</td><td>1</td></tr><tr><td>LS-strong-driving</td><td>Stratified box</td><td>1</td><td>0.24</td><td>7.6</td><td>RAMSES</td><td>Colman et al. (2022)</td><td>1</td></tr><tr><td>M51</td><td>Full galaxy</td><td>—</td><td>≈1.4</td><td>0</td><td>AREPO</td><td>Tress et al. (2020)</td><td>1</td></tr><tr><td>Ramses-F20</td><td>Full galaxy</td><td>—</td><td>0.92</td><td>0</td><td>RAMSES</td><td>Brucy (2022, Chapter 11)</td><td>1</td></tr></tbody></table>

**Notes.** We list the type of setup, the size of the computational domain, the spatial resolution, the initial mid-plane magnetic field strength, the simulation code used to run the simulation, the reference in which the simulation was first described and the number of snapshots analysed.

## 2. Numerical methods, simulation data

### 2.1. Structure finding algorithm

Our structure finding algorithm is built around the HOP clump finder (Eisenstein & Hut 1998). While originally it was developed for finding groups of particles in N-body simulations, the algorithm and its original implementation are quite general and can be used in many more applications. One advantage HOP has is that it simply requires a list of particles as input. While this may seem sub-optimal for analysing regular grid-based hydrodynamics simulations, it allows for a direct comparison between the results of particle-based codes, unstructured meshes as well as regular grid-based (AMR) codes. To achieve this, we transform the AMR grid into a list of cells.

Overall, there are three steps towards obtaining a structure catalogue, which are described in depth in Appendix A. First, we use HOP to identify peak patches around local density maxima. Then, we use the regrouping tool from HOP to merge peaks based on a number of specified criteria. For clouds, we adopt a density threshold of 30 cm⁻³. Everything below this threshold is not considered to be part of a cloud. Then we require a cloud to have a peak density of at least 60 cm⁻³, which corresponds to a minimum peak-to-background ratio of 2. When the saddle density between two structures is larger than 300 cm⁻³, that is to say ten times the density threshold, we consider them as being substructures within the same cloud and therefore merge them. Lastly, we require a structure to be made up of at least 100 cells. After this merging step, we know for each cell or particle the structure to which it belongs, or whether it does not belong to any structure. Based on that, we then finally calculate the properties of all identified structures as detailed in Appendix A.

### 2.2. Numerical simulations of the ISM

We apply the HOP algorithm to several (M)HD simulations of the ISM. We combine simulations with a wide range of simulation techniques, resolution, box size and physics included. An overview is given in Table 1. A face-on projection of one snapshot of each simulation type is featured in Fig. 1, where we over-plot the positions of the extracted clouds colour-coded by their derived velocity dispersion for illustration.

#### 2.2.1. SILCC: Stratified boxes with FLASH

We include high-resolution magnetised simulations from the SILCC collaboration (Walch et al. 2015; Girichidis et al. 2016), which are described in detail in Girichidis et al. (2018) and Girichidis (2021). The setup consists of a stratified box with a volume of 0.5 × 0.5 × 0.5 kpc³. The gas is initially at rest and in hydrostatic equilibrium with a total gas mass surface density of 10 M☉ pc⁻². The two different magnetic field models considered here have central magnetic field strengths at z = 0 of Bₓ,₀ = 3 µG and Bₓ,₀ = 6 µG and are oriented along the x direction. The field strength scales vertically with the gas density as Bₓ(z) = Bₓ,₀ [ρ(z)/ρ(z = 0)]¹/². The initial setup is magnetically supercritical, that is the field cannot support the gas against collapse.

The MHD equations are solved using the HLLR5 solver (Bouchut et al. 2007, 2010; Waagan 2009; Waagan et al. 2011), which is implemented in FLASH in Version 4 (Fryxell et al. 2000; Dubey et al. 2008). Heating of the gas includes spatially clustered supernovae (SNe), as well as cosmic ray (Goldsmith & Langer 1978) and X-ray heating (Wolfire et al. 1995). The CR ionisation and heating rates are ζCR = 3 × 10⁻¹⁷ s⁻¹ and ΓCR = 3.2 × 10⁻¹¹ ζCR n erg s⁻¹ cm⁻³, respectively. Photoelectric heating follows Bakes & Tielens (1994), Bergin et al. (2004) and Wolfire et al. (2003). We assume a spatially constant interstellar radiation field (ISRF) with a strength of 1.7 in the units of the Habing field G₀ (Habing 1968; Draine 1978) that is then locally attenuated in dense, shielded gas using the TreeCol algorithm (Clark et al. 2012; Wünsch et al. 2018). We fix the dust-to-gas mass ratio to 0.01, and adopt dust opacities from Mathis et al. (1983) for wavelengths shorter than 1 µm and Ossenkopf & Henning (1994) for longer wavelengths.

Heating and radiative cooling processes are computed using a chemical network that tracks the non-equilibrium concentrations of ionised hydrogen (H⁺), atomic hydrogen (H), molecular hydrogen (H₂), as well as singly ionised carbon (C⁺) and carbon monoxide (CO). The hydrogen chemistry is based on the models by Hollenbach & McKee (1989); Glover & Mac Low (2007a,b) and Micic et al. (2012). The reactions connected to CO follow the model developed by Nelson & Langer (1997). Radiative cooling incorporates contributions from the fine structure lines of C⁺, O, and Si⁺. Furthermore, we consider rotational and vibrational
<table><caption>**Table 1.** Set of simulations.</caption><thead><tr><th>Name</th><th>Setup</th><th>Size<br>(kpc)</th><th>Δx<sub>min</sub><br>(pc)</th><th>B<sub>0</sub><br>(μG)</th><th>Code</th><th>References</th><th>No. of<br>snapshots</th></tr></thead><tbody><tr><td>SILCC-1pc-3μG</td><td>Stratified box</td><td>0.5</td><td>0.98</td><td>3</td><td>FLASH</td><td>Girichidis et al. (2018)</td><td>4</td></tr><tr><td>SILCC-1pc-6μG</td><td>Stratified box</td><td>0.5</td><td>0.98</td><td>6</td><td>FLASH</td><td>Girichidis et al. (2018)</td><td>4</td></tr><tr><td>SILCC-0.5pc-3μG</td><td>Stratified box</td><td>0.5</td><td>0.49</td><td>3</td><td>FLASH</td><td>Girichidis (2021)</td><td>6</td></tr><tr><td>SILCC-0.25pc-3μG</td><td>Stratified box</td><td>0.5</td><td>0.24</td><td>3</td><td>FLASH</td><td>Girichidis (2021)</td><td>2</td></tr><tr><td>LS-no-driving</td><td>Stratified box</td><td>1</td><td>0.24</td><td>7.6</td><td>RAMSES</td><td>Colman et al. (2022)</td><td>1</td></tr><tr><td>LS-weak-driving</td><td>Stratified box</td><td>1</td><td>0.24</td><td>7.6</td><td>RAMSES</td><td>Colman et al. (2022)</td><td>1</td></tr><tr><td>LS-medium-driving</td><td>Stratified box</td><td>1</td><td>0.24</td><td>7.6</td><td>RAMSES</td><td>Colman et al. (2022)</td><td>1</td></tr><tr><td>LS-strong-driving</td><td>Stratified box</td><td>1</td><td>0.24</td><td>7.6</td><td>RAMSES</td><td>Colman et al. (2022)</td><td>1</td></tr><tr><td>M51</td><td>Full galaxy</td><td>—</td><td>≈1.4</td><td>0</td><td>AREPO</td><td>Tress et al. (2020)</td><td>1</td></tr><tr><td>Ramses-F20</td><td>Full galaxy</td><td>—</td><td>0.92</td><td>0</td><td>RAMSES</td><td>Brucy (2022, Chapter 11)</td><td>1</td></tr></tbody></table>

**Notes.** We list the type of setup, the size of the computational domain, the spatial resolution, the initial mid-plane magnetic field strength, the simulation code used to run the simulation, the reference in which the simulation was first described and the number of snapshots analysed.

## 2. Numerical methods, simulation data

### 2.1. Structure finding algorithm

Our structure finding algorithm is built around the HOP clump finder (Eisenstein & Hut 1998). While originally it was developed for finding groups of particles in N-body simulations, the algorithm and its original implementation are quite general and can be used in many more applications. One advantage HOP has is that it simply requires a list of particles as input. While this may seem sub-optimal for analysing regular grid-based hydrodynamics simulations, it allows for a direct comparison between the results of particle-based codes, unstructured meshes as well as regular grid-based (AMR) codes. To achieve this, we transform the AMR grid into a list of cells.

Overall, there are three steps towards obtaining a structure catalogue, which are described in depth in Appendix A. First, we use HOP to identify peak patches around local density maxima. Then, we use the regrouping tool from HOP to merge peaks based on a number of specified criteria. For clouds, we adopt a density threshold of 30 cm⁻³. Everything below this threshold is not considered to be part of a cloud. Then we require a cloud to have a peak density of at least 60 cm⁻³, which corresponds to a minimum peak-to-background ratio of 2. When the saddle density between two structures is larger than 300 cm⁻³, that is to say ten times the density threshold, we consider them as being substructures within the same cloud and therefore merge them. Lastly, we require a structure to be made up of at least 100 cells. After this merging step, we know for each cell or particle the structure to which it belongs, or whether it does not belong to any structure. Based on that, we then finally calculate the properties of all identified structures as detailed in Appendix A.

### 2.2. Numerical simulations of the ISM

We apply the HOP algorithm to several (M)HD simulations of the ISM. We combine simulations with a wide range of simulation techniques, resolution, box size and physics included. An overview is given in Table 1. A face-on projection of one snapshot of each simulation type is featured in Fig. 1, where we over-plot the positions of the extracted clouds colour-coded by their derived velocity dispersion for illustration.

#### 2.2.1. SILCC: Stratified boxes with FLASH

We include high-resolution magnetised simulations from the SILCC collaboration (Walch et al. 2015; Girichidis et al. 2016), which are described in detail in Girichidis et al. (2018) and Girichidis (2021). The setup consists of a stratified box with a volume of 0.5 × 0.5 × 0.5 kpc³. The gas is initially at rest and in hydrostatic equilibrium with a total gas mass surface density of 10 M☉ pc⁻². The two different magnetic field models considered here have central magnetic field strengths at z = 0 of Bₓ,₀ = 3 µG and Bₓ,₀ = 6 µG and are oriented along the x direction. The field strength scales vertically with the gas density as Bₓ(z) = Bₓ,₀ [ρ(z)/ρ(z = 0)]¹/². The initial setup is magnetically supercritical, that is the field cannot support the gas against collapse.

The MHD equations are solved using the HLLR5 solver (Bouchut et al. 2007, 2010; Waagan 2009; Waagan et al. 2011), which is implemented in FLASH in Version 4 (Fryxell et al. 2000; Dubey et al. 2008). Heating of the gas includes spatially clustered supernovae (SNe), as well as cosmic ray (Goldsmith & Langer 1978) and X-ray heating (Wolfire et al. 1995). The CR ionisation and heating rates are ζCR = 3 × 10⁻¹⁷ s⁻¹ and ΓCR = 3.2 × 10⁻¹¹ ζCR n erg s⁻¹ cm⁻³, respectively. Photoelectric heating follows Bakes & Tielens (1994), Bergin et al. (2004) and Wolfire et al. (2003). We assume a spatially constant interstellar radiation field (ISRF) with a strength of 1.7 in the units of the Habing field G₀ (Habing 1968; Draine 1978) that is then locally attenuated in dense, shielded gas using the TreeCol algorithm (Clark et al. 2012; Wünsch et al. 2018). We fix the dust-to-gas mass ratio to 0.01, and adopt dust opacities from Mathis et al. (1983) for wavelengths shorter than 1 µm and Ossenkopf & Henning (1994) for longer wavelengths.

Heating and radiative cooling processes are computed using a chemical network that tracks the non-equilibrium concentrations of ionised hydrogen (H⁺), atomic hydrogen (H), molecular hydrogen (H₂), as well as singly ionised carbon (C⁺) and carbon monoxide (CO). The hydrogen chemistry is based on the models by Hollenbach & McKee (1989); Glover & Mac Low (2007a,b) and Micic et al. (2012). The reactions connected to CO follow the model developed by Nelson & Langer (1997). Radiative cooling incorporates contributions from the fine structure lines of C⁺, O, and Si⁺. Furthermore, we consider rotational and vibrational
<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Column density map for a single snapshot from each of the four listed groups of simulations. The panels show SILCC-0.5pc-3µG (a), LS-weak-driving (b), M51 (c), Ramses-F20 (d). Details of these simulations can be found in Table 1. The coloured symbols indicate the locations of the clouds identified by the HOP algorithm, with the colour corresponding to the internal velocity dispersion.</figcaption></figure>

lines of H₂ and CO, as well as Lyman-α cooling. The energy transfer from the gas to the dust follows Glover et al. (2010) and Glover & Clark (2012a). Above $10^4$ K, we assume collisional ionisation equilibrium (CIE) for helium and all heavy elements, and adopt the corresponding CIE cooling rates from Gnat & Ferland (2012). Hydrogen is not assumed to be in collisional ionisation equilibrium, and its contribution to the cooling process is calculated self-consistently at all temperatures, accounting for any deviations from chemical equilibrium.

Gravitational forces include self-gravity as well as an external potential based on an isothermal sheet (Spitzer 1942) with a surface density of $30 M_{\odot} pc^{-2}$, representing the stellar mass distribution. The vertical scale height of this sheet is set to 100 pc. The Poisson equation is solved using the tree-based method as described in Wünsch et al. (2018).

Star formation and the related supernova feedback are included at a fixed rate. We derive a star formation rate based on the Kennicutt-Schmidt relation (Kennicutt 1998) and convert it
<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Column density map for a single snapshot from each of the four listed groups of simulations. The panels show SILCC-0.5pc-3µG (a), LS-weak-driving (b), M51 (c), Ramses-F20 (d). Details of these simulations can be found in Table 1. The coloured symbols indicate the locations of the clouds identified by the HOP algorithm, with the colour corresponding to the internal velocity dispersion.</figcaption></figure>

lines of H₂ and CO, as well as Lyman-α cooling. The energy transfer from the gas to the dust follows Glover et al. (2010) and Glover & Clark (2012a). Above $10^4$ K, we assume collisional ionisation equilibrium (CIE) for helium and all heavy elements, and adopt the corresponding CIE cooling rates from Gnat & Ferland (2012). Hydrogen is not assumed to be in collisional ionisation equilibrium, and its contribution to the cooling process is calculated self-consistently at all temperatures, accounting for any deviations from chemical equilibrium.

Gravitational forces include self-gravity as well as an external potential based on an isothermal sheet (Spitzer 1942) with a surface density of $30 M_{\odot} pc^{-2}$, representing the stellar mass distribution. The vertical scale height of this sheet is set to 100 pc. The Poisson equation is solved using the tree-based method as described in Wünsch et al. (2018).

Star formation and the related supernova feedback are included at a fixed rate. We derive a star formation rate based on the Kennicutt-Schmidt relation (Kennicutt 1998) and convert it
to a SN rate using the Chabrier (2003) IMF. This yields 15 SNe per Myr for the volume simulated here. Each explosion injects $10^{51}$ erg of thermal energy. For the positioning of the SNe we distinguish between a type Ia (20%) and a type II component (80%). The type Ia SNe are individual explosions placed at random (uniformly chosen) positions for *x* and *y*. The *z* position is randomly sampled from a Gaussian distribution with a vertical scale height of 300 pc (Bahcall & Soneira 1980; Heiles 1987), where we ensure that SNe that are placed outside the box are re-drawn to be placed inside the box. The type II SNe are further split into isolated SNe (40% of the type II) and clustered counterparts (60%) (Heiles 1987; Kennicutt et al. 1989; McKee & Williams 1997; Clarke & Oey 2002). For the positioning of the individual type II SNe and the clusters we also chose random positions for *x* and *y*. The vertical distribution is drawn from a Gaussian with a smaller scale height of 90 pc. All SNe within a cluster explode at the same cluster position. The positions of the SNe are computed beforehand and stored in a table to ensure an identical feedback configuration for the different magnetic field strengths and resolutions.

We consider three different resolutions for the simulations. In all cases, the base resolution is $128^3$, corresponding to a maximum cubic cell with a side length of 3.9 pc. On top we add two additional levels of refinement for simulations SILCC-1pc-3 µG and SILCC-1pc-6 µG, which reach a minimum cell size of 0.98 pc. For simulations SILCC-0.5pc-3µG and SILCC-0.25pc-3µG we restart simulation SILCC-1pc-3µG at a time of 20 Myr and add one and two additional levels of refinement, which yields minimum cell sizes of 0.49 and 0.24 pc, respectively. The simulations with a resolution of 1 pc run for a total time of 60 Myr; the higher resolution ones were evolved for a shorter time of a few Myr (Girichidis 2021).

### 2.2.2. LS: Stratified boxes with RAMSES

For our comparison, we also use RAMSES (Teyssier 2002) simulations of a stratified piece of a galactic disk. The simulations are fully described in Colman et al. (2022) and references therein. The computational domain is cubic with a box size of 1 kpc. The base grid with a resolution of 3.9 pc is further refined using a mass-based criterion up to a minimum cell size of 0.24 pc in the densest regions. Compared to the SILCC simulations, the computational domain is larger, but the different refinement strategy can results in fewer AMR grids (see discussion in Sect. 5.1).

The initial density profile is Gaussian $n(z) = n_0 \exp[-0.5(z/z_0)^2]$ with a mid-plane particle density $n_0 = 1.5 \text{ cm}^{-3}$ and a thickness $z_0 = 150$ pc, corresponding to a column density of $\Sigma = 19.1 M_{\odot} \text{pc}^{-2}$. An initial level of turbulence is introduced by adding a turbulent velocity field with a root mean square dispersion of 5 km s$^{-1}$ and a Kolmogorov power spectrum $E(k) \propto k^{-5/3}$ with random phase. The initial temperature is 5333 K, which is a typical value for the warm neutral medium phase in the ISM. We also include a magnetic field oriented along the *x*-axis with $B_x(z) = B_0 \exp[-0.5(z/z_0)^2]$ initially, where $B_0 = 7.62 \mu\text{G}$ is the mid-plane field strength. The magnetic field is solved using the ideal MHD approximation, as implemented by Fromang et al. (2006). Aside from self-gravity, we also apply an external gravitational disk potential as prescribed by Kuijken & Gilmore (1989) and Joung & Mac Low (2006), which accounts for the profile of (old) stars and dark matter.

We use an ISM cooling model based on Audit & Hennebelle (2005) that includes the most important processes responsible for regulating the thermal balance of the atomic ISM. At low

temperatures, cooling is provided by the fine structure transitions of $C^+$ and $\text{O}^{\text{II}}$, while at high temperatures this is supplemented by contributions from grain surface recombination of atomic hydrogen and from Lyman-$\alpha$ cooling. Heating is provided primarily by the photoelectric effect, calculated assuming a constant uniform UV background with a strength equal to the solar neighbourhood value. In contrast to the SILCC and AREPO simulations, we do not account for local attenuation of this radiation field within dense clouds. We also account for cosmic ray heating, using the rate given in Goldsmith (2001). The ionisation of hydrogen is treated by the RT module from Rosdahl et al. (2013). We use sink particles to track star formation and stellar feedback self-consistently. Sinks represent newly formed star clusters (Bleuler & Teyssier 2014), which accrete gas which is within a 4 cell accretion radius and above the sink formation threshold of $10^4 \text{ cm}^{-3}$. Each time a sink has accreted enough mass, it forms a massive star particle with a mass between 8 and $120 M_{\odot}$ randomly drawn from the Salpeter IMF (Salpeter 1955). This massive star emits ionising radiation and explodes as a supernova at the end of its lifetime (Rosdahl et al. 2013; Geen et al. 2016; Colling et al. 2018). The radiation is treated using a moment-based method (Rosdahl et al. 2013). In addition to the turbulence generated by stellar feedback, we also include external driving on scales between 1 box length and 1/3 of the box length with solenoidal fraction 0.75. Several forcing strengths were explored: weak, medium and strong driving, as well as no driving, which result in mass weighted velocity dispersions of $\sigma_{3D} \approx 9, 12, 20$ and $8.5 \text{ km s}^{-1}$, respectively. The simulations were evolved for 60 Myr.

### 2.2.3. M51: Full galaxy with AREPO

Stratified box set-ups allow one to simulate the ISM at high resolution but miss important aspects linked to large scales: self-consistent large-scale driving, galactic dynamical effects such as rotation or the influence of spiral density waves, and large-scale instabilities such as the Toomre instability (Toomre 1964) or the wiggle instability (Sormani et al. 2017). As this may influence the properties of the clouds that form in the ISM, we included two full galaxy simulations in our study. The first one was carried out using AREPO, a moving-mesh hydrodynamic code (Springel 2010), and was fully presented in Tress et al. (2020). The clouds from this simulation were already extracted and analysed in Tress et al. (2021), but with a different cloud extraction algorithm. The modelled galaxy is interacting with a smaller companion and its properties are chosen to be similar to the M51 ‘Whirlpool’ galaxy. The characteristics of the different components (dark matter halo, stellar bulge, stellar disk and gaseous disk) are summarised in Table B.1. Gravitational interactions between all those components were self-consistently accounted for.

Resolution in AREPO is defined by the target mass, which is the typical mass of a cell. This is analogous to the particle mass in a smoothed particle hydrodynamics (SPH) computation. For the M51 simulation, the target mass is 300 $M_{\odot}$, with additional refinement in the denser part of the ISM such that the Jeans length is resolved by at least four cells up to a density of $10^{-21} \text{ g cm}^{-3}$. As a result of this refinement scheme, most of the cells above the HOP density threshold have a radius of 1 pc or less, with the smallest cells having a radius of only 0.4 pc. For comparison with the minimum $\Delta x$ in the AMR simulations, we quote in Table 1 the typical diameter of a cell at the sink creation density threshold, which is ~1.4 pc. The effective resolution of this simulation in dense regions is thus slightly worse than the low resolution SILCC runs.
to a SN rate using the Chabrier (2003) IMF. This yields 15 SNe per Myr for the volume simulated here. Each explosion injects $10^{51}$ erg of thermal energy. For the positioning of the SNe we distinguish between a type Ia (20%) and a type II component (80%). The type Ia SNe are individual explosions placed at random (uniformly chosen) positions for *x* and *y*. The *z* position is randomly sampled from a Gaussian distribution with a vertical scale height of 300 pc (Bahcall & Soneira 1980; Heiles 1987), where we ensure that SNe that are placed outside the box are re-drawn to be placed inside the box. The type II SNe are further split into isolated SNe (40% of the type II) and clustered counterparts (60%) (Heiles 1987; Kennicutt et al. 1989; McKee & Williams 1997; Clarke & Oey 2002). For the positioning of the individual type II SNe and the clusters we also chose random positions for *x* and *y*. The vertical distribution is drawn from a Gaussian with a smaller scale height of 90 pc. All SNe within a cluster explode at the same cluster position. The positions of the SNe are computed beforehand and stored in a table to ensure an identical feedback configuration for the different magnetic field strengths and resolutions.

We consider three different resolutions for the simulations. In all cases, the base resolution is $128^3$, corresponding to a maximum cubic cell with a side length of 3.9 pc. On top we add two additional levels of refinement for simulations SILCC-1pc-3 µG and SILCC-1pc-6 µG, which reach a minimum cell size of 0.98 pc. For simulations SILCC-0.5pc-3µG and SILCC-0.25pc-3µG we restart simulation SILCC-1pc-3µG at a time of 20 Myr and add one and two additional levels of refinement, which yields minimum cell sizes of 0.49 and 0.24 pc, respectively. The simulations with a resolution of 1 pc run for a total time of 60 Myr; the higher resolution ones were evolved for a shorter time of a few Myr (Girichidis 2021).

### 2.2.2. LS: Stratified boxes with RAMSES

For our comparison, we also use RAMSES (Teyssier 2002) simulations of a stratified piece of a galactic disk. The simulations are fully described in Colman et al. (2022) and references therein. The computational domain is cubic with a box size of 1 kpc. The base grid with a resolution of 3.9 pc is further refined using a mass-based criterion up to a minimum cell size of 0.24 pc in the densest regions. Compared to the SILCC simulations, the computational domain is larger, but the different refinement strategy can results in fewer AMR grids (see discussion in Sect. 5.1).

The initial density profile is Gaussian $n(z) = n_0 \exp[-0.5(z/z_0)^2]$ with a mid-plane particle density $n_0 = 1.5 \text{ cm}^{-3}$ and a thickness $z_0 = 150$ pc, corresponding to a column density of $\Sigma = 19.1 M_{\odot} \text{pc}^{-2}$. An initial level of turbulence is introduced by adding a turbulent velocity field with a root mean square dispersion of 5 km s$^{-1}$ and a Kolmogorov power spectrum $E(k) \propto k^{-5/3}$ with random phase. The initial temperature is 5333 K, which is a typical value for the warm neutral medium phase in the ISM. We also include a magnetic field oriented along the *x*-axis with $B_x(z) = B_0 \exp[-0.5(z/z_0)^2]$ initially, where $B_0 = 7.62 \mu\text{G}$ is the mid-plane field strength. The magnetic field is solved using the ideal MHD approximation, as implemented by Fromang et al. (2006). Aside from self-gravity, we also apply an external gravitational disk potential as prescribed by Kuijken & Gilmore (1989) and Joung & Mac Low (2006), which accounts for the profile of (old) stars and dark matter.

We use an ISM cooling model based on Audit & Hennebelle (2005) that includes the most important processes responsible for regulating the thermal balance of the atomic ISM. At low

temperatures, cooling is provided by the fine structure transitions of $C^+$ and $\text{O}^{\text{II}}$, while at high temperatures this is supplemented by contributions from grain surface recombination of atomic hydrogen and from Lyman-$\alpha$ cooling. Heating is provided primarily by the photoelectric effect, calculated assuming a constant uniform UV background with a strength equal to the solar neighbourhood value. In contrast to the SILCC and AREPO simulations, we do not account for local attenuation of this radiation field within dense clouds. We also account for cosmic ray heating, using the rate given in Goldsmith (2001). The ionisation of hydrogen is treated by the RT module from Rosdahl et al. (2013). We use sink particles to track star formation and stellar feedback self-consistently. Sinks represent newly formed star clusters (Bleuler & Teyssier 2014), which accrete gas which is within a 4 cell accretion radius and above the sink formation threshold of $10^4 \text{ cm}^{-3}$. Each time a sink has accreted enough mass, it forms a massive star particle with a mass between 8 and $120 M_{\odot}$ randomly drawn from the Salpeter IMF (Salpeter 1955). This massive star emits ionising radiation and explodes as a supernova at the end of its lifetime (Rosdahl et al. 2013; Geen et al. 2016; Colling et al. 2018). The radiation is treated using a moment-based method (Rosdahl et al. 2013). In addition to the turbulence generated by stellar feedback, we also include external driving on scales between 1 box length and 1/3 of the box length with solenoidal fraction 0.75. Several forcing strengths were explored: weak, medium and strong driving, as well as no driving, which result in mass weighted velocity dispersions of $\sigma_{3D} \approx 9, 12, 20$ and 8.5 km s$^{-1}$, respectively. The simulations were evolved for 60 Myr.

### 2.2.3. M51: Full galaxy with AREPO

Stratified box set-ups allow one to simulate the ISM at high resolution but miss important aspects linked to large scales: self-consistent large-scale driving, galactic dynamical effects such as rotation or the influence of spiral density waves, and large-scale instabilities such as the Toomre instability (Toomre 1964) or the wiggle instability (Sormani et al. 2017). As this may influence the properties of the clouds that form in the ISM, we included two full galaxy simulations in our study. The first one was carried out using AREPO, a moving-mesh hydrodynamic code (Springel 2010), and was fully presented in Tress et al. (2020). The clouds from this simulation were already extracted and analysed in Tress et al. (2021), but with a different cloud extraction algorithm. The modelled galaxy is interacting with a smaller companion and its properties are chosen to be similar to the M51 ‘Whirlpool’ galaxy. The characteristics of the different components (dark matter halo, stellar bulge, stellar disk and gaseous disk) are summarised in Table B.1. Gravitational interactions between all those components were self-consistently accounted for.

Resolution in AREPO is defined by the target mass, which is the typical mass of a cell. This is analogous to the particle mass in a smoothed particle hydrodynamics (SPH) computation. For the M51 simulation, the target mass is 300 $M_{\odot}$, with additional refinement in the denser part of the ISM such that the Jeans length is resolved by at least four cells up to a density of $10^{-21} \text{ g cm}^{-3}$. As a result of this refinement scheme, most of the cells above the HOP density threshold have a radius of 1 pc or less, with the smallest cells having a radius of only 0.4 pc. For comparison with the minimum $\Delta x$ in the AMR simulations, we quote in Table 1 the typical diameter of a cell at the sink creation density threshold, which is ~1.4 pc. The effective resolution of this simulation in dense regions is thus slightly worse than the low resolution SILCC runs.
Unlike the previous sets of simulations, the magnetic field is not modelled. The thermal and chemical evolution of the gas are solved simultaneously, using the NL97 chemical network described in Glover & Clark (2012b) and the atomic and molecular cooling function described in Clark et al. (2019). Star formation is modelled using sink particles that are created within gravitationally bound regions of radius 2.5 pc when the gas density exceeds a threshold of $10^{-21}$ g cm$^{-3}$, if the region satisfies several additional criteria: it must be a local minimum of the gravitational potential, and both the velocity and the acceleration must be converging. Once formed, sink particles can accrete gas from within an accretion radius of 2.5 pc, with the accretion limited to gas with a density exceeding the sink creation density. Two feedback processes are modelled: type II and type Ia SNe. The type II SNe are attached to the sinks, via a process similar to the one described in Sect. 2.2.2. Type Ia SNe are exploded at the position of a randomly selected old star with a rate of one explosion every 250 yr. These are the stellar particles of the disk and bulge components of the N-body model of the galaxy. Supernovae were modelled by injecting energy in a region containing 40 resolution elements surrounding the injection region. Momentum or thermal energy is injected based on whether the Sedov-Taylor phase is resolved.

### 2.2.4. Ramses-F20: Full galaxy with RAMSES

To complete the study, we also ran the algorithm on an isolated disk simulation performed with RAMSES that was initially presented in Brucy (2022, Chapter 11). The model is an isolated disk with a radius of 15 kpc, within a box with a side length of 120 kpc. Initial conditions are generated with the code DICE, which allow one to specify the different components of a galaxy. DICE generates particles representing dark matter, stars and gas. The two first are directly used in RAMSES while the gas distribution was mapped from the particles onto the AMR grid. The parameters for the generation of the initial conditions are summed up in Table B.1.

Using the adaptive mesh refinement capabilities of RAMSES, the maximal cell size is 1 kpc outside the disk and 7 pc within the disk. The minimal cell size is 0.92 pc within the disk. Dark matter and initial stars are modelled via particles that only undergo gravity. The gas follows the law of hydrodynamics (without magnetic field) and gravity. The simulation is run over a period of 300 Myr.

The simulation also includes a sub-grid model for star formation and stellar feedback, here limited to SN explosions and photo-ionisation-triggered heating. These models are those developed and described in Kretschmer & Teyssier (2020). The star formation rate is computed for each cell from its properties. First the Mach number $M$ within the cell is computed thanks to a model of the evolution of the sub-grid turbulent velocity. This is then used to estimate the sub-grid density PDF, assuming that it is a log-normal (Vazquez-Semadeni 1994; Federrath & Klessen 2012), as well as the virial parameter of the cell. The expected star formation rate within the cell is then computed using the multi-free-fall analytical model for the star formation rate adapted from Krumholz et al. (2005) by Hennebelle & Chabrier (2011) and popularised by Federrath & Klessen (2012). At each time-step the number of new star particles created is randomly drawn from a Poisson distribution, calibrated so that the mean star formation rate over time is equal to the one given by the analytical model. The SN injection scheme depends whether the cooling radius the SN is resolved by at least one grid cell. If yes, the energy from the explosion is directly injected in the form

of thermal energy. Otherwise, the terminal momentum of the SN is also injected. Photo-ionisation is modelled by maintaining the cell where the star sits at a temperature of $10^4$ K until the last SN explodes. The cooling method is the same as the one used in the RAMSES stratified box simulation presented in Sect. 2.2.2.

### 2.3. Examples of extracted clouds

Figure 2 displays some examples of extracted objects in the various simulations. The images show the projected column density in a region around the clouds, so we can also see their environment. The white contours show the projected cloud boundaries. The mass and velocity dispersion are listed at the top of the image. A coloured symbol is assigned to each example cloud, which will be used to identify them in further figures.

The column density projections illustrate the wide variety in sizes and shapes we find. The clouds range from simple roundish shapes to very complex structures, which cannot be captured with simple measures such as the sphericity or the triaxiality. This is particularly clear when comparing the projected images for the LS clouds in Fig. 2b with a 3D rendering of the same clouds in Fig. 3, where we also show the ellipsoid approximation to the cloud from which its size is estimated (see Sect. 3 and Appendix A.3). In most cases, the cloud's shape is far from ellipsoid (see also Ebagezio et al. 2023). However the corresponding sizes give a reasonable approximation to the extent of the objects. We also note that projection effects are significant. Especially large clouds can appear quite different when viewing them from different angles. Furthermore, there are often multiple clouds along the same line of sight. Consequently, the projections along the line of sight are limited to the cloud under consideration and its direct environment rather than being the integration throughout the entire simulation box.

In the stratified box simulations, large clouds are filamentary in nature and their internal structure is resolved. Small clouds typically do not have a resolved substructure and appear more regular in shape, which is true for all simulations. In the galaxy simulations, some of the largest clouds have disk-like shapes as illustrated in the left panel of Figs. 2c and 2d. These are unrealistic and result from a lack of resolution.

### 2.4. Notes on comparing with cloud catalogues from observations

Many catalogues of ISM structures have been compiled from observations using various techniques. A complete and robust comparison of simulation results to these catalogues requires the creation of synthetic observations, which could then be processed by the same pipeline as the observations to extract the clouds in the exact same way. This is clearly beyond the scope of this work. Nevertheless, it is useful to compare the trends we find in the simulations, which will be discussed in the following sections, to the ones observed in the Milky Way and nearby galaxies.

The classic molecular cloud catalogues are obtained from CO observations. The main CO isotopologue, $^{12}$CO, is detectable at relatively low densities (Snow & McCall 2006), potentially as low as our cloud extraction density threshold of $30$ cm$^{-3}$, provided there is molecular gas at these densities. Emission lines from the $^{13}$CO isotopologue trace somewhat higher densities. The spectral information can be used to try to disentangle overlapping projected clouds. The observed linewidth of the transition is assumed to be set mainly by turbulence. While the dimension of the observed clouds is typically derived
Unlike the previous sets of simulations, the magnetic field is not modelled. The thermal and chemical evolution of the gas are solved simultaneously, using the NL97 chemical network described in Glover & Clark (2012b) and the atomic and molecular cooling function described in Clark et al. (2019). Star formation is modelled using sink particles that are created within gravitationally bound regions of radius 2.5 pc when the gas density exceeds a threshold of $10^{-21}$ g cm$^{-3}$, if the region satisfies several additional criteria: it must be a local minimum of the gravitational potential, and both the velocity and the acceleration must be converging. Once formed, sink particles can accrete gas from within an accretion radius of 2.5 pc, with the accretion limited to gas with a density exceeding the sink creation density. Two feedback processes are modelled: type II and type Ia SNe. The type II SNe are attached to the sinks, via a process similar to the one described in Sect. 2.2.2. Type Ia SNe are exploded at the position of a randomly selected old star with a rate of one explosion every 250 yr. These are the stellar particles of the disk and bulge components of the N-body model of the galaxy. Supernovae were modelled by injecting energy in a region containing 40 resolution elements surrounding the injection region. Momentum or thermal energy is injected based on whether the Sedov-Taylor phase is resolved.

### 2.2.4. Ramses-F20: Full galaxy with RAMSES

To complete the study, we also ran the algorithm on an isolated disk simulation performed with RAMSES that was initially presented in Brucy (2022, Chapter 11). The model is an isolated disk with a radius of 15 kpc, within a box with a side length of 120 kpc. Initial conditions are generated with the code DICE, which allow one to specify the different components of a galaxy. DICE generates particles representing dark matter, stars and gas. The two first are directly used in RAMSES while the gas distribution was mapped from the particles onto the AMR grid. The parameters for the generation of the initial conditions are summed up in Table B.1.

Using the adaptive mesh refinement capabilities of RAMSES, the maximal cell size is 1 kpc outside the disk and 7 pc within the disk. The minimal cell size is 0.92 pc within the disk. Dark matter and initial stars are modelled via particles that only undergo gravity. The gas follows the law of hydrodynamics (without magnetic field) and gravity. The simulation is run over a period of 300 Myr.

The simulation also includes a sub-grid model for star formation and stellar feedback, here limited to SN explosions and photo-ionisation-triggered heating. These models are those developed and described in Kretschmer & Teyssier (2020). The star formation rate is computed for each cell from its properties. First the Mach number $M$ within the cell is computed thanks to a model of the evolution of the sub-grid turbulent velocity. This is then used to estimate the sub-grid density PDF, assuming that it is a log-normal (Vazquez-Semadeni 1994; Federrath & Klessen 2012), as well as the virial parameter of the cell. The expected star formation rate within the cell is then computed using the multi-free-fall analytical model for the star formation rate adapted from Krumholz et al. (2005) by Hennebelle & Chabrier (2011) and popularised by Federrath & Klessen (2012). At each time-step the number of new star particles created is randomly drawn from a Poisson distribution, calibrated so that the mean star formation rate over time is equal to the one given by the analytical model. The SN injection scheme depends whether the cooling radius the SN is resolved by a least one grid cell. If yes, the energy from the explosion is directly injected in the form

of thermal energy. Otherwise, the terminal momentum of the SN is also injected. Photo-ionisation is modelled by maintaining the cell where the star sits at a temperature of $10^4$ K until the last SN explodes. The cooling method is the same as the one used in the RAMSES stratified box simulation presented in Sect. 2.2.2.

### 2.3. Examples of extracted clouds

Figure 2 displays some examples of extracted objects in the various simulations. The images show the projected column density in a region around the clouds, so we can also see their environment. The white contours show the projected cloud boundaries. The mass and velocity dispersion are listed at the top of the image. A coloured symbol is assigned to each example cloud, which will be used to identify them in further figures.

The column density projections illustrate the wide variety in sizes and shapes we find. The clouds range from simple roundish shapes to very complex structures, which cannot be captured with simple measures such as the sphericity or the triaxiality. This is particularly clear when comparing the projected images for the LS clouds in Fig. 2b with a 3D rendering of the same clouds in Fig. 3, where we also show the ellipsoid approximation to the cloud from which its size is estimated (see Sect. 3 and Appendix A.3). In most cases, the cloud's shape is far from ellipsoid (see also Ebagezio et al. 2023). However the corresponding sizes give a reasonable approximation to the extent of the objects. We also note that projection effects are significant. Especially large clouds can appear quite different when viewing them from different angles. Furthermore, there are often multiple clouds along the same line of sight. Consequently, the projections along the line of sight are limited to the cloud under consideration and its direct environment rather than being the integration throughout the entire simulation box.

In the stratified box simulations, large clouds are filamentary in nature and their internal structure is resolved. Small clouds typically do not have a resolved substructure and appear more regular in shape, which is true for all simulations. In the galaxy simulations, some of the largest clouds have disk-like shapes as illustrated in the left panel of Figs. 2c and 2d. These are unrealistic and result from a lack of resolution.

### 2.4. Notes on comparing with cloud catalogues from observations

Many catalogues of ISM structures have been compiled from observations using various techniques. A complete and robust comparison of simulation results to these catalogues requires the creation of synthetic observations, which could then be processed by the same pipeline as the observations to extract the clouds in the exact same way. This is clearly beyond the scope of this work. Nevertheless, it is useful to compare the trends we find in the simulations, which will be discussed in the following sections, to the ones observed in the Milky Way and nearby galaxies.

The classic molecular cloud catalogues are obtained from CO observations. The main CO isotopologue, $^{12}$CO, is detectable at relatively low densities (Snow & McCall 2006), potentially as low as our cloud extraction density threshold of $30$ cm$^{-3}$, provided there is molecular gas at these densities. Emission lines from the $^{13}$CO isotopologue trace somewhat higher densities. The spectral information can be used to try to disentangle overlapping projected clouds. The observed linewidth of the transition is assumed to be set mainly by turbulence. While the dimension of the observed clouds is typically derived
<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Column density maps of selected clouds extracted from the simulations. The different rows show examples from SILCC-0.5-pc (a), LS-no-driving (b), M51 (c), Ramses-F20 (d). Column density is calculated in a box of three times the maximum extent of the object. Labels are added showing the cloud mass and velocity dispersion. The coloured symbol is used to refer to the cloud in the following figures in the paper.</figcaption></figure>
<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Column density maps of selected clouds extracted from the simulations. The different rows show examples from SILCC-0.5-pc (a), LS-no-driving (b), M51 (c), Ramses-F20 (d). Column density is calculated in a box of three times the maximum extent of the object. Labels are added showing the cloud mass and velocity dispersion. The coloured symbol is used to refer to the cloud in the following figures in the paper.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>**Fig. 3.** 3D visualisations of the clouds from Fig. 2b, extracted from the LS stratified box simulations. The density scale is linear. The ellipsoid approximation (in cyan) and its half-axes (in black) are also shown.</figcaption></figure>

from the circularised radius of the cloud area in all the catalogues, the mass is derived with different methods that include different assumptions and approximations. A key parameter in the mass and radius estimate is the distance of the cloud. In the Milky Way, this is typically derived using a model for the rotation curve of the Galaxy, a method which is prone to large uncertainties (see e.g. Reid 2022). Thanks to the recent effort to create 3D dust maps (e.g. Chen et al. 2019; Leike et al. 2020), it has become possible to accurately estimate the distances to clouds detected in dust extinction for structures up to a distance of ~400 pc, offering an alternative to CO cloud catalogues. Chen et al. (2020) detected clouds directly in PPP space, but cloud properties were still derived from 2D projections. Very recently, Dharmawardena et al. (2023) and Cahlon et al. (2024) calculated the cloud mass and size from the reconstructed 3D dust density distributions. 3D dust maps do not provide kinematic information and so one cannot recover the velocity dispersion with this data alone. So far, no one has attempted to combine spectral line and 3D dust information to derive cloud velocity dispersion. Extra-galactic surveys of MCs in nearby galaxies do not suffer from the same distance uncertainties, but are typically sensitive to only the most massive MCs in all but the closest Local Group galaxies (Rosolowsky et al. 2021).

All current Milky Way molecular cloud catalogues have specific lower limits in the recoverable masses and radii. In particular, the smallest recoverable radius depends on the spatial resolution of the observations. The smallest observable mass depends on the sensitivity. Both resolution and sensitivity decrease with increasing cloud distance. As a consequence, catalogues such as those of Rice et al. (2016) and Miville-Deschênes et al. (2017) that were derived from the first generations of CO surveys (Dame et al. 2001) with spatial resolution of 8'5 and spectral resolution of 1.3 km s⁻¹ contain on average larger and more massive clouds than more recent compilations of molecular clouds (e.g. Benedettini et al. 2020, 2021; Duarte-Cabral et al. 2021) that are extracted from the new generation of CO surveys with sub-arcminute spatial resolution and spectral resolution well below 1 km s⁻¹. Moreover, in general all the catalogues of the observed molecular clouds are incomplete at the lower masses and radii. The threshold for incompleteness depends on the distance and thus can be variable for catalogues of the entire Milky Way.

It must also be noted that CO, especially ¹²CO, becomes optically thick in the denser parts of clouds (see e.g. Tielens 2010; Draine 2011). This means that we cannot see beyond a certain maximum surface density, leading us to underestimate of the mass if not properly accounted for, which is hard to do. Dust, on the other hand, is an optically thin tracer and thus does not have this problem. This, and the other remarks made here, must

be kept in mind when putting the results described in the next sections into an observational context.

### 3. Comparison of cloud mass, size, and shape

We start our comparison of the cloud catalogues by looking at the masses, sizes and shapes of the clouds. The cloud mass is simply the sum of the masses of all of the cells/particles it contains. The size is determined through the diagonalisation of the inertia matrix. Using the resulting eigenvalues, we construct an ellipsoid of uniform density which has the same moments of inertia as the cloud. The size is then taken to be the geometric average of the axis lengths of this ellipsoid. The corresponding equations can be found in Appendix A. Figure 3 illustrates that the ellipsoid approximation traces well the extent of a cloud.

An alternative definition of the size is the cubic root of the total volume of the cloud, which neglects the shape information. In the following, we use both definitions as they provide different insights.

#### 3.1. Mass and size distribution

The mass and size distributions are shown in Fig. 4. To facilitate the comparison, we normalise the histograms according to the surface area covered by the simulation. In other words, the figure shows the number of clouds per unit area (here kpc²). For the stratified box simulations (SILCC and LS), this is straightforward as their computational domain defines the area of the region. For the stratified boxes we combine the catalogues of several individual simulations and/or snapshots into one histogram to keep the plot readable. This is justified because the distributions do not differ significantly. The distribution for SILCC-1pc contains the two SILCC runs with 1 pc resolution and different magnetic field strengths (SILCC-1pc-3 µG and SILCC-1pc-6 µG). The label LS includes all four LS runs with different driving strengths. The number of snapshots used for this compilation is listed in Table 1. The total surface area is then the box area multiplied by the number of snapshots. For the galaxy simulations, the relevant area is less straightforward to define. To get an estimate, we divided the face-on column density image into 1 kpc² tiles and counted how many of them have an average column density larger than 3 M⊙ pc⁻².

##### 3.1.1. Power-law tail

The high-mass end of the mass distribution of clouds is typically described by a power law $dN/d\log M \propto M^\alpha$, or alternatively $dN/dM \propto M^\gamma$ with $\gamma = \alpha - 1$, optionally with a truncation. Unless stated differently, listed exponents correspond to $\alpha$. With
<figure><img src="image_4.png" /><figcaption>**Fig. 3.** 3D visualisations of the clouds from Fig. 2b, extracted from the LS stratified box simulations. The density scale is linear. The ellipsoid approximation (in cyan) and its half-axes (in black) are also shown.</figcaption></figure>

from the circularised radius of the cloud area in all the catalogues, the mass is derived with different methods that include different assumptions and approximations. A key parameter in the mass and radius estimate is the distance of the cloud. In the Milky Way, this is typically derived using a model for the rotation curve of the Galaxy, a method which is prone to large uncertainties (see e.g. Reid 2022). Thanks to the recent effort to create 3D dust maps (e.g. Chen et al. 2019; Leike et al. 2020), it has become possible to accurately estimate the distances to clouds detected in dust extinction for structures up to a distance of ~400 pc, offering an alternative to CO cloud catalogues. Chen et al. (2020) detected clouds directly in PPP space, but cloud properties were still derived from 2D projections. Very recently, Dharmawardena et al. (2023) and Cahlon et al. (2024) calculated the cloud mass and size from the reconstructed 3D dust density distributions. 3D dust maps do not provide kinematic information and so one cannot recover the velocity dispersion with this data alone. So far, no one has attempted to combine spectral line and 3D dust information to derive cloud velocity dispersion. Extra-galactic surveys of MCs in nearby galaxies do not suffer from the same distance uncertainties, but are typically sensitive to only the most massive MCs in all but the closest Local Group galaxies (Rosolowsky et al. 2021).

All current Milky Way molecular cloud catalogues have specific lower limits in the recoverable masses and radii. In particular, the smallest recoverable radius depends on the spatial resolution of the observations. The smallest observable mass depends on the sensitivity. Both resolution and sensitivity decrease with increasing cloud distance. As a consequence, catalogues such as those of Rice et al. (2016) and Miville-Deschênes et al. (2017) that were derived from the first generations of CO surveys (Dame et al. 2001) with spatial resolution of 8'5 and spectral resolution of 1.3 km s⁻¹ contain on average larger and more massive clouds than more recent compilations of molecular clouds (e.g. Benedettini et al. 2020, 2021; Duarte-Cabral et al. 2021) that are extracted from the new generation of CO surveys with sub-arcminute spatial resolution and spectral resolution well below 1 km s⁻¹. Moreover, in general all the catalogues of the observed molecular clouds are incomplete at the lower masses and radii. The threshold for incompleteness depends on the distance and thus can be variable for catalogues of the entire Milky Way.

It must also be noted that CO, especially ¹²CO, becomes optically thick in the denser parts of clouds (see e.g. Tielens 2010; Draine 2011). This means that we cannot see beyond a certain maximum surface density, leading us to underestimate of the mass if not properly accounted for, which is hard to do. Dust, on the other hand, is an optically thin tracer and thus does not have this problem. This, and the other remarks made here, must

be kept in mind when putting the results described in the next sections into an observational context.

### 3. Comparison of cloud mass, size, and shape

We start our comparison of the cloud catalogues by looking at the masses, sizes and shapes of the clouds. The cloud mass is simply the sum of the masses of all of the cells/particles it contains. The size is determined through the diagonalisation of the inertia matrix. Using the resulting eigenvalues, we construct an ellipsoid of uniform density which has the same moments of inertia as the cloud. The size is then taken to be the geometric average of the axis lengths of this ellipsoid. The corresponding equations can be found in Appendix A. Figure 3 illustrates that the ellipsoid approximation traces well the extent of a cloud.

An alternative definition of the size is the cubic root of the total volume of the cloud, which neglects the shape information. In the following, we use both definitions as they provide different insights.

#### 3.1. Mass and size distribution

The mass and size distributions are shown in Fig. 4. To facilitate the comparison, we normalise the histograms according to the surface area covered by the simulation. In other words, the figure shows the number of clouds per unit area (here kpc²). For the stratified box simulations (SILCC and LS), this is straightforward as their computational domain defines the area of the region. For the stratified boxes we combine the catalogues of several individual simulations and/or snapshots into one histogram to keep the plot readable. This is justified because the distributions do not differ significantly. The distribution for SILCC-1pc contains the two SILCC runs with 1 pc resolution and different magnetic field strengths (SILCC-1pc-3 µG and SILCC-1pc-6 µG). The label LS includes all four LS runs with different driving strengths. The number of snapshots used for this compilation is listed in Table 1. The total surface area is then the box area multiplied by the number of snapshots. For the galaxy simulations, the relevant area is less straightforward to define. To get an estimate, we divided the face-on column density image into 1 kpc² tiles and counted how many of them have an average column density larger than 3 M⊙ pc⁻².

##### 3.1.1. Power-law tail

The high-mass end of the mass distribution of clouds is typically described by a power law $dN/d\log M \propto M^\alpha$, or alternatively $dN/dM \propto M^\gamma$ with $\gamma = \alpha - 1$, optionally with a truncation. Unless stated differently, listed exponents correspond to $\alpha$. With
<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Comparison of the mass (top) and size (bottom) distribution of clouds in the different simulations, normalised according to the relevant surface area of the simulation. To guide the eye, we add a grey line indicating an appropriate power-law slopes of −1 for the mass spectrum in the top panel and −3 for the size distribution in the bottom panel. The dashed lines correspond to multiples of the spatial resolution in the SILCC and LS simulations. The dotted lines show multiples of the resolution for the Ramses galaxy.</figcaption></figure>

the appropriate normalisation, the agreement of the slope of both the mass and size distribution is remarkable. Especially for the mass spectrum which spans six orders of magnitude ranging from giant cloud complexes of $10^6–10^7 M_\odot$ down to tiny cloudlets of barely $10 M_\odot$. The value of the power-law exponent is close to −1, as indicated by the grey line. The exponent of the size distribution is around −3. For large masses and sizes, we typically see a cut-off in the distributions for the stratified box simulations due to the limited box size and the resulting limited total mass. Since we analyse only a handful of snapshots, values of the histogram below 1 cloud per $kpc^2$ correspond to rare events and are not statistically significant.

The only massive clouds which do not lie on the general mass spectrum are the most massive clouds in the Ramses-F20 simulation, which show an excess compared to the overall trend. We found 286 clouds with a mass above $10^6 M_\odot$, indicating that this is not a statistical fluctuation. Furthermore, the cloud mass distribution in this simulation seems to consist of multiple components. Upon closer inspection (cf. Figs. 5 and C.1), we see a secondary peak around $5 \times 10^5 M_\odot$, indicating an excess of clouds with high masses, in both the Ramses-F20 galaxy and M51. These object are spurious disk-shaped clouds like the examples shown in the left column of Fig. 2c, and thus not physical because of limited resolution. Once formed, they are very difficult to destroy by regular stellar feedback, which is limited as well by the overall worse resolution.

We fitted the mass spectrum slope in the appropriate range for individual simulations (Appendix C). The resulting slopes for the stratified boxed and M51 agree well with one another, ranging from −1.16 to −1.35. Ramses-F20 has a shallower slope of −0.76. This difference is likely related to differences in the numerical feedback recipes used in each simulation. The SN recipe used in Ramses-F20 injects energy or momentum only in one of the cell neighbouring the exploding star, which makes it hard to destroy large structures. Meanwhile, because the injection radius is linked to a minimal number of cells, the recipe used in the AREPO simulation tends to overestimate the feedback effect when the resolution is poor. This difference alone could explain why clouds in Ramses-F20 are typically larger than those in M51.

### 3.1.2. Turn-over and low-mass end

The low-mass end of the spectra are shaped by resolution effects. In all cases, the mass and size distribution is incomplete for values below the peak. The SILCC setups show a relatively sharp cut-off slightly above four times the resolution limit, an effect of our requirement that a cloud has to contain at least 100 cells<fnref n="1" />. We note that RAMSES simulations (LS and Ramses-F20) have different shapes at the low-mass end, with a much shallower transition between a peak and a low-mass cut-off. We discuss in Sect. 5.1 that this stems from a difference in the grid refinement strategy used in RAMSES and FLASH. For LS, the peak is observed at a value slightly larger than the spatial resolution of the coarse grid (refinement level 8, corresponding to a cell size of 4 pc). For the Ramses-F20 galaxy, it is somewhere between twice and four times the coarse grid resolution. For M51 the cut-off is less straightforward to determine, since both the cell size and the cell mass vary in AREPO simulations. This results in a relatively wide mass and size distribution for the cells. We also typically have overall fewer cells in AREPO runs compared to AMR simulations (FLASH and RAMSES). The combination of these two aspects result in a low cloud mass cut-off that is largely determined by the minimum number of cells for a structure used by HOP, which is 100. At the highest densities reached in the M51 simulation, the typical particle mass is around 10–20 $M_\odot$, and so we see a cut-off in the cloud mass distribution at 100 times this value, namely 1000–2000 $M_\odot$. We note that this cut-off is not as sharp as in the SILCC runs because there are some cells in the M51 simulation with masses below 10 $M_\odot$ and hence HOP identifies a few clouds containing these cells with total masses below 1000 $M_\odot$.

<footnote n="1">¹ The actual minimal size depends on the shape, but we note that $\sqrt[3]{100} \approx 4.64$.</footnote>
<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Comparison of the mass (top) and size (bottom) distribution of clouds in the different simulations, normalised according to the relevant surface area of the simulation. To guide the eye, we add a grey line indicating an appropriate power-law slopes of −1 for the mass spectrum in the top panel and −3 for the size distribution in the bottom panel. The dashed lines correspond to multiples of the spatial resolution in the SILCC and LS simulations. The dotted lines show multiples of the resolution for the Ramses galaxy.</figcaption></figure>

the appropriate normalisation, the agreement of the slope of both the mass and size distribution is remarkable. Especially for the mass spectrum which spans six orders of magnitude ranging from giant cloud complexes of $10^6–10^7 M_\odot$ down to tiny cloudlets of barely $10 M_\odot$. The value of the power-law exponent is close to −1, as indicated by the grey line. The exponent of the size distribution is around −3. For large masses and sizes, we typically see a cut-off in the distributions for the stratified box simulations due to the limited box size and the resulting limited total mass. Since we analyse only a handful of snapshots, values of the histogram below 1 cloud per $kpc^2$ correspond to rare events and are not statistically significant.

The only massive clouds which do not lie on the general mass spectrum are the most massive clouds in the Ramses-F20 simulation, which show an excess compared to the overall trend. We found 286 clouds with a mass above $10^6 M_\odot$, indicating that this is not a statistical fluctuation. Furthermore, the cloud mass distribution in this simulation seems to consist of multiple components. Upon closer inspection (cf. Figs. 5 and C.1), we see a secondary peak around $5 \times 10^5 M_\odot$, indicating an excess of clouds with high masses, in both the Ramses-F20 galaxy and M51. These object are spurious disk-shaped clouds like the examples shown in the left column of Fig. 2c, and thus not physical because of limited resolution. Once formed, they are very difficult to destroy by regular stellar feedback, which is limited as well by the overall worse resolution.

We fitted the mass spectrum slope in the appropriate range for individual simulations (Appendix C). The resulting slopes for the stratified boxed and M51 agree well with one another, ranging from −1.16 to −1.35. Ramses-F20 has a shallower slope of −0.76. This difference is likely related to differences in the numerical feedback recipes used in each simulation. The SN recipe used in Ramses-F20 injects energy or momentum only in one of the cell neighbouring the exploding star, which makes it hard to destroy large structures. Meanwhile, because the injection radius is linked to a minimal number of cells, the recipe used in the AREPO simulation tends to overestimate the feedback effect when the resolution is poor. This difference alone could explain why clouds in Ramses-F20 are typically larger than those in M51.

### 3.1.2. Turn-over and low-mass end

The low-mass end of the spectra are shaped by resolution effects. In all cases, the mass and size distribution is incomplete for values below the peak. The SILCC setups show a relatively sharp cut-off slightly above four times the resolution limit, an effect of our requirement that a cloud has to contain at least 100 cells<fnref n="1" />. We note that RAMSES simulations (LS and Ramses-F20) have different shapes at the low-mass end, with a much shallower transition between a peak and a low-mass cut-off. We discuss in Sect. 5.1 that this stems from a difference in the grid refinement strategy used in RAMSES and FLASH. For LS, the peak is observed at a value slightly larger than the spatial resolution of the coarse grid (refinement level 8, corresponding to a cell size of 4 pc). For the Ramses-F20 galaxy, it is somewhere between twice and four times the coarse grid resolution. For M51 the cut-off is less straightforward to determine, since both the cell size and the cell mass vary in AREPO simulations. This results in a relatively wide mass and size distribution for the cells. We also typically have overall fewer cells in AREPO runs compared to AMR simulations (FLASH and RAMSES). The combination of these two aspects result in a low cloud mass cut-off that is largely determined by the minimum number of cells for a structure used by HOP, which is 100. At the highest densities reached in the M51 simulation, the typical particle mass is around 10–20 $M_\odot$, and so we see a cut-off in the cloud mass distribution at 100 times this value, namely 1000–2000 $M_\odot$. We note that this cut-off is not as sharp as in the SILCC runs because there are some cells in the M51 simulation with masses below 10 $M_\odot$ and hence HOP identifies a few clouds containing these cells with total masses below 1000 $M_\odot$.

<footnote n="1">¹ The actual minimal size depends on the shape, but we note that $\sqrt[3]{100} \approx 4.64$.</footnote>
<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Comparison of the mass distribution of clouds in the inner 7 kpc of a galaxy versus those in the outskirts. While Ramses-F20 overall has a shallower mass spectrum than M51 (and the kpc boxes), there is no difference in the measured slope between the inner and the outer galaxy for either simulation.</figcaption></figure>

It is worth noting that the ranges of masses and radii mea-
sured in observational data are in agreement with the ranges
derived from simulations with a minimum cell size comparab-
le to the spatial resolution of the observations. When quoting
median values for the property distributions, it thus only makes
sense to compare with studies which have comparable resolution.

3.1.3. Predictions from theory

Turbulent fragmentation theories predict a mass spectrum with a power law for massive objects on all scales from GMCs to cores. The self-similarity of gravity, which is argued to be the dominant force, results in $dN/d\log M \propto M^\alpha$ with $\alpha = -1$. Detailed derivations of turbulent fragmentation theories, which take into account several corrections, predict slopes which are slightly shallower (Hennebelle & Chabrier 2008; Hopkins 2012). Also accretion theories predict $\alpha = -1$ when the accretion rate scales as $\dot{M} \propto M^2$ (Kuznetsova et al. 2018). The global slope we find is indeed close to $-1$.

3.1.4. Slope variations in observations

Observationally, the obtained slope varies significantly between
studies. Roman-Duval et al. (2010) find $\alpha = -1.64 \pm 0.25$, con-
sistent with the value of -1.5 obtained in older work (Sanders
et al. 1985; Solomon et al. 1987; Williams & McKee 1997).
Miville-Deschênes et al. (2017) find $\alpha = -2.0 \pm 0.1$ for their
full Milky Way catalogue. Rice et al. (2016), who used the same
underlying CO data as Miville-Deschênes et al. (2017), report
a shallower slope and differences between the outer Galaxy
for which $\alpha = -1.2 \pm 0.1$, and the inner Galaxy where $\alpha =$
$-0.6 \pm 0.1$. A similar diversity was obtained by Rosolowsky
(2005) who found $\alpha = -0.5 \pm 0.1$ and $\alpha = -1.1 \pm 0.2$ for the
inner and outer Milky Way respectively, $\alpha = -1.9 \pm 0.4$ for the
M33 galaxy and $\alpha = -0.7 \pm 0.2$ for the Large Magellanic Cloud.
A study of the actual M51 galaxy also reports variations with
environments: $\alpha \approx -0.3$ for the centre and bar, $\alpha \approx -0.8$ to

-0.6 for the molecular ring and spiral arms and $\alpha \approx -1.5$ for the inter-arm and outer galaxy regions (Colombo et al. 2014). Moreover, the mass spectrum is not well-described by a power law in all regions, suggesting different dominant mechanisms for the formation and destruction of clouds in the different parts of the galaxy. In summary, the observed slopes range from as shallow as -0.3 to as steep as -2.0. This variety is in contradiction to the universality of our results which are all fairly close to -1. Figure 5 investigates more closely the mass spectrum in our two galaxy simulations, dividing the cloud population into two groups: clouds in the inner and outer galaxy with the boundary being at 7 kpc from the centre of mass. We find no significant difference in the slope. A small difference in typical cloud mass is obtained for clouds in the very centre of the galaxies, as discussed in Appendix D. A similar trend was already found by Tress et al. (2021), who extracted clouds from their M51 galaxy using a different algorithm than we apply here. Possibly the measured slope depends on the exact method with which clouds have been extracted from observations. The same holds for simulations as will be discussed in Sect. 5.2. We also note also that it can be difficult to estimate the completeness limit, which may also affect the fit of the slope. The situation is even more complex for catalogues for which the cloud distance is not uniform. Since the effective spatial resolution is lower at greater distances, small clouds appear blended into larger objects, which results in a shallower mass spectrum slope. A more rigorous comparison with observations through synthetic observations could provide a more direct conclusion about whether or not the simulation results are in agreement with the observations.

3.2. Mass–size relation

Observations show that cloud properties are not independent
from each other (Hennebelle & Falgarone 2012; Heyer & Dame
2015). A tight relation between cloud mass and size is reported
by many studies (e.g. Solomon et al. 1987; Roman-Duval et al.
2010; Kauffmann et al. 2013; Miville-Deschênes et al. 2017).
In Fig. 6a, we plot the relation between the mass and size of
each cloud for the various simulations of our comparison set.
We remind the reader that the size of the cloud is defined as
the average between the three principal axes, as identified by the
structure algorithm. The tight correlation we see with a power-
law slope of 3 is expected and is due to the way clouds are
extracted. By construction, the mean density of the cloud cannot
be below the chosen threshold of 30 cm⁻³ above which cells are
considered by the cloud detection algorithm. Also, it is rather
improbable to have an average density above the saddle den-
sity threshold. Unless the density gradient is very sharp, such
a high density region will be associated with a less dense enve-
lope which together form a bigger, more massive and on average
less dense cloud. In Fig. 6b, we replace the size by the cubic root
of the volume, and by doing so, remove the information about
the shape. We can then see clearly that almost all the clouds
have indeed an average density between the algorithm density
and saddle thresholds.

It is thus interesting to understand the spread and the outliers
in both of these figures. The spread in Fig. 6a is mainly due to
the shape of the clouds. Indeed, for a given mass and volume,
a spherical cloud will have a smaller effective size than a more
elongated one. We study in detail the distribution of shapes in
Sect. 3.3. The simulation LS 0.25 pc stands out from the others
with a significant number of small clouds with high densities.
These are dense clumps where the low density envelope has been
<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Comparison of the mass distribution of clouds in the inner 7 kpc of a galaxy versus those in the outskirts. While Ramses-F20 overall has a shallower mass spectrum than M51 (and the kpc boxes), there is no difference in the measured slope between the inner and the outer galaxy for either simulation.</figcaption></figure>

It is worth noting that the ranges of masses and radii mea-
sured in observational data are in agreement with the ranges
derived from simulations with a minimum cell size comparab-
le to the spatial resolution of the observations. When quoting
median values for the property distributions, it thus only makes
sense to compare with studies which have comparable resolution.

3.1.3. Predictions from theory

Turbulent fragmentation theories predict a mass spectrum with a power law for massive objects on all scales from GMCs to cores. The self-similarity of gravity, which is argued to be the dominant force, results in $dN/d\log M \propto M^\alpha$ with $\alpha = -1$. Detailed derivations of turbulent fragmentation theories, which take into account several corrections, predict slopes which are slightly shallower (Hennebelle & Chabrier 2008; Hopkins 2012). Also accretion theories predict $\alpha = -1$ when the accretion rate scales as $\dot{M} \propto M^2$ (Kuznetsova et al. 2018). The global slope we find is indeed close to $-1$.

3.1.4. Slope variations in observations

Observationally, the obtained slope varies significantly between
studies. Roman-Duval et al. (2010) find $\alpha = -1.64 \pm 0.25$, con-
sistent with the value of -1.5 obtained in older work (Sanders
et al. 1985; Solomon et al. 1987; Williams & McKee 1997).
Miville-Deschênes et al. (2017) find $\alpha = -2.0 \pm 0.1$ for their
full Milky Way catalogue. Rice et al. (2016), who used the same
underlying CO data as Miville-Deschênes et al. (2017), report
a shallower slope and differences between the outer Galaxy
for which $\alpha = -1.2 \pm 0.1$, and the inner Galaxy where $\alpha =$
$-0.6 \pm 0.1$. A similar diversity was obtained by Rosolowsky
(2005) who found $\alpha = -0.5 \pm 0.1$ and $\alpha = -1.1 \pm 0.2$ for the
inner and outer Milky Way respectively, $\alpha = -1.9 \pm 0.4$ for the
M33 galaxy and $\alpha = -0.7 \pm 0.2$ for the Large Magellanic Cloud.
A study of the actual M51 galaxy also reports variations with
environments: $\alpha \approx -0.3$ for the centre and bar, $\alpha \approx -0.8$ to

-0.6 for the molecular ring and spiral arms and $\alpha \approx -1.5$ for the inter-arm and outer galaxy regions (Colombo et al. 2014). Moreover, the mass spectrum is not well-described by a power law in all regions, suggesting different dominant mechanisms for the formation and destruction of clouds in the different parts of the galaxy. In summary, the observed slopes range from as shallow as -0.3 to as steep as -2.0. This variety is in contradiction to the universality of our results which are all fairly close to -1. Figure 5 investigates more closely the mass spectrum in our two galaxy simulations, dividing the cloud population into two groups: clouds in the inner and outer galaxy with the boundary being at 7 kpc from the centre of mass. We find no significant difference in the slope. A small difference in typical cloud mass is obtained for clouds in the very centre of the galaxies, as discussed in Appendix D. A similar trend was already found by Tress et al. (2021), who extracted clouds from their M51 galaxy using a different algorithm than we apply here. Possibly the measured slope depends on the exact method with which clouds have been extracted from observations. The same holds for simulations as will be discussed in Sect. 5.2. We also note also that it can be difficult to estimate the completeness limit, which may also affect the fit of the slope. The situation is even more complex for catalogues for which the cloud distance is not uniform. Since the effective spatial resolution is lower at greater distances, small clouds appear blended into larger objects, which results in a shallower mass spectrum slope. A more rigorous comparison with observations through synthetic observations could provide a more direct conclusion about whether or not the simulation results are in agreement with the observations.

3.2. Mass–size relation

Observations show that cloud properties are not independent
from each other (Hennebelle & Falgarone 2012; Heyer & Dame
2015). A tight relation between cloud mass and size is reported
by many studies (e.g. Solomon et al. 1987; Roman-Duval et al.
2010; Kauffmann et al. 2013; Miville-Deschênes et al. 2017).
In Fig. 6a, we plot the relation between the mass and size of
each cloud for the various simulations of our comparison set.
We remind the reader that the size of the cloud is defined as
the average between the three principal axes, as identified by the
structure algorithm. The tight correlation we see with a power-
law slope of 3 is expected and is due to the way clouds are
extracted. By construction, the mean density of the cloud cannot
be below the chosen threshold of 30 cm⁻³ above which cells are
considered by the cloud detection algorithm. Also, it is rather
improbable to have an average density above the saddle den-
sity threshold. Unless the density gradient is very sharp, such
a high density region will be associated with a less dense enve-
lope which together form a bigger, more massive and on average
less dense cloud. In Fig. 6b, we replace the size by the cubic root
of the volume, and by doing so, remove the information about
the shape. We can then see clearly that almost all the clouds
have indeed an average density between the algorithm density
and saddle thresholds.

It is thus interesting to understand the spread and the outliers
in both of these figures. The spread in Fig. 6a is mainly due to
the shape of the clouds. Indeed, for a given mass and volume,
a spherical cloud will have a smaller effective size than a more
elongated one. We study in detail the distribution of shapes in
Sect. 3.3. The simulation LS 0.25 pc stands out from the others
with a significant number of small clouds with high densities.
These are dense clumps where the low density envelope has been
<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Correlation between cloud size and mass, with the size computed in two different ways. In panel a, the size is defined as the average between the three half-axis of the approximating ellipse, as determined by the structure algorithm. In panel b, the size is computed as the cubic root of the volume. Each dot represents one cloud. The contour lines show a kernel density estimate of the distribution. Starting from the innermost line, approximately 10, 30, 50, and 70% of the distribution lie within the respective contours. The symbols highlight the example clouds from Fig. 2. The grey line in panel a indicates a slope of 3. The solid and dotted line in panel b are iso-density lines for the detection and saddle threshold respectively, as defined in Sect. 2.1 and Appendix A. Their slopes are 3.</figcaption></figure>

stripped by stellar feedback, specifically ionising radiation. This
is less destructive than SN feedback, and in a cloud with dense
substructure, it tends to preferentially blow away the diffuse
gas, leaving the isolated dense clumps, which are then picked
up by the cloud detection algorithm. This form of feedback is
not included in the other simulations, which either include only
SN feedback (SILCC, M51) or SN feedback plus localised pre-
SN heating (Ramses-F20), which explains why we only see this
population of objects in the LS simulations.

3.2.1. Relation between the number of dimensions and the M-R slope

In Table 2 we list some of the recent mass-size relation slopes
obtained from observed cloud catalogues. Both CO and dust
extinction studies report mass-size relations with an exponent
of about 2 when cloud properties are obtained from projec-
tions on the sky. This suggests a constant surface density
for clouds (Larson 1981; Lada & Dame 2020). Some studies
report a slightly steeper slope, but Ballesteros-Paredes et al.
(2019) demonstrated that this is likely an effect of overlapping
clouds. Interestingly, when cloud properties are calculated from
3D data, that is the radius from the volume and the mass from
the volume density, a scaling relation of $M \propto R^3$ is obtained,
in line with the global relation established in our simulations.

Cahlon et al. (2024) investigated exactly the effect of projection
on the mass-size relation in their study of clouds extracted
from 3D dust maps. They obtained $M \propto R^{2.9}$ for 3D clouds
but recovered a scaling of $M \propto R^{2.1}$ for projected clouds, in
agreement with the classical observed relation. This effect was
predicted by theory (Ballesteros-Paredes & Mac Low 2002;
Ballesteros-Paredes et al. 2012, 2019) and previous analysis of
hydrodynamics simulations (Shetty et al. 2010). These works
indicate that the slope is determined by the way the size is
measured and the underlying number of dimensions. Using the
area results in a slope of 2 while using the volume results in
a slope of 3. We note that clouds are thought to be fractal in
nature (Falgarone et al. 1991; Elmegreen & Falgarone 1996;
Elmegreen & Elmegreen 2001) and thus in principle neither of
these methods uses a good approximation for the shape. This
behaviour is due to the cloud mass being dominated by low den-
sity, volume-filling material. This can be seen clearly in the 3D
visualisations of our extracted clouds (Fig. 3). This makes the
mass sensitive to how the cloud boundaries are defined which
relates to the shape. We thus can expect the mass–size relation to
change when studying denser structures, such as dense clumps,
for which this condition may no longer hold. The LS simulations
do contain a the population of small objects which have been
stripped of their envelope. Interestingly, they indeed follows a
relation which is shallower than 3, as can be seen in Fig. 6a.
<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Correlation between cloud size and mass, with the size computed in two different ways. In panel a, the size is defined as the average between the three half-axis of the approximating ellipse, as determined by the structure algorithm. In panel b, the size is computed as the cubic root of the volume. Each dot represents one cloud. The contour lines show a kernel density estimate of the distribution. Starting from the innermost line, approximately 10, 30, 50, and 70% of the distribution lie within the respective contours. The symbols highlight the example clouds from Fig. 2. The grey line in panel a indicates a slope of 3. The solid and dotted line in panel b are iso-density lines for the detection and saddle threshold respectively, as defined in Sect. 2.1 and Appendix A. Their slopes are 3.</figcaption></figure>

stripped by stellar feedback, specifically ionising radiation. This is less destructive than SN feedback, and in a cloud with dense substructure, it tends to preferentially blow away the diffuse gas, leaving the isolated dense clumps, which are then picked up by the cloud detection algorithm. This form of feedback is not included in the other simulations, which either include only SN feedback (SILCC, M51) or SN feedback plus localised pre-SN heating (Ramses-F20), which explains why we only see this population of objects in the LS simulations.

### 3.2.1. Relation between the number of dimensions and the M-R slope

In Table 2 we list some of the recent mass-size relation slopes obtained from observed cloud catalogues. Both CO and dust extinction studies report mass-size relations with an exponent of about 2 when cloud properties are obtained from projections on the sky. This suggests a constant surface density for clouds (Larson 1981; Lada & Dame 2020). Some studies report a slightly steeper slope, but Ballesteros-Paredes et al. (2019) demonstrated that this is likely an effect of overlapping clouds. Interestingly, when cloud properties are calculated from 3D data, that is the radius from the volume and the mass from the volume density, a scaling relation of $M \propto R^3$ is obtained, in line with the global relation established in our simulations.

Cahlon et al. (2024) investigated exactly the effect of projection on the mass-size relation in their study of clouds extracted from 3D dust maps. They obtained $M \propto R^{2.9}$ for 3D clouds but recovered a scaling of $M \propto R^{2.1}$ for projected clouds, in agreement with the classical observed relation. This effect was predicted by theory (Ballesteros-Paredes & Mac Low 2002; Ballesteros-Paredes et al. 2012, 2019) and previous analysis of hydrodynamics simulations (Shetty et al. 2010). These works indicate that the slope is determined by the way the size is measured and the underlying number of dimensions. Using the area results in a slope of 2 while using the volume results in a slope of 3. We note that clouds are thought to be fractal in nature (Falgarone et al. 1991; Elmegreen & Falgarone 1996; Elmegreen & Elmegreen 2001) and thus in principle neither of these methods uses a good approximation for the shape. This behaviour is due to the cloud mass being dominated by low density, volume-filling material. This can be seen clearly in the 3D visualisations of our extracted clouds (Fig. 3). This makes the mass sensitive to how the cloud boundaries are defined which relates to the shape. We thus can expect the mass–size relation to change when studying denser structures, such as dense clumps, for which this condition may no longer hold. The LS simulations do contain a the population of small objects which have been stripped of their envelope. Interestingly, they indeed follows a relation which is shallower than 3, as can be seen in Fig. 6a.
<table><caption>**Table 2.** Compilation of observational studies, with summary of their underlying data used to produce the molecular cloud catalogue, and their reported mass-size relation.</caption><thead><tr><th>References</th><th>Data type</th><th>Milky Way region</th><th>M-R slope</th></tr></thead><tbody><tr><td>Rice et al. (2016)</td><td><sup>12</sup>CO(1–0)</td><td>Full</td><td>1.98</td></tr><tr><td>Miville-Deschênes et al. (2017)</td><td><sup>12</sup>CO(1–0)</td><td>Full</td><td>2.2</td></tr><tr><td>Roman-Duval et al. (2010)</td><td><sup>12/13</sup>CO(1–0)</td><td>Galactic ring</td><td>2.36</td></tr><tr><td>Benedettini et al. (2020)</td><td><sup>12</sup>CO(1–0)</td><td>Outer Galaxy</td><td>≈2</td></tr><tr><td>Benedettini et al. (2021)</td><td><sup>13</sup>CO(1–0)</td><td>Outer Galaxy</td><td>2.13</td></tr><tr><td>Duarte-Cabral et al. (2021)</td><td><sup>13</sup>CO(2–1)</td><td>Inner Galaxy</td><td>≈2.3</td></tr><tr><td>Ma et al. (2022)</td><td><sup>12/13</sup>CO(1–0)</td><td>Third quadrant</td><td>2.24</td></tr><tr><td>Chen et al. (2020)</td><td>3D dust (projected)</td><td>Solar neighbourhood</td><td>1.96</td></tr><tr><td>Dharmawardena et al. (2023)</td><td>3D dust</td><td>Solar neighbourhood</td><td>≈3</td></tr><tr><td>Cahlon et al. (2024)</td><td>3D dust</td><td>Solar neighbourhood</td><td>2.9</td></tr><tr><td>Cahlon et al. (2024)</td><td>3D dust (projected)</td><td>Solar neighbourhood</td><td>2.1</td></tr></tbody></table>

### 3.2.2. Normalisation of M-R relation

When $M = AR^2$ the normalisation A provides a roughly constant surface density for clouds, or alternatively a constant average volume density when $M = AR^3$. Lada & Dame (2020) reported a systematic shift of the M-R relation to higher surface densities when clouds have been extracted from dust extinction versus <sup>12</sup>CO, in other words the normalisation of the relation is different for different tracers. In fact, for any area in the sky the mass derived from CO is lower than the mass derived from cold dust. Small differences can also be observed between results derived from <sup>12</sup>CO and <sup>13</sup>CO (Benedettini et al. 2021). This discrepancy could be (partially) due to optical depth effects, which can be severe for <sup>12</sup>CO as mentioned already in Sect. 2.4. If the mass is mainly set by the low density material in the cloud envelop rather than the dense cores and filaments, the effect of optical depth on the mass determination should be limited. In this case, the normalisation is set by the density threshold of the cloud extraction algorithm or the equivalent sensitivity limit of the observations, as we concluded from Fig. 6b. The threshold sets a characteristic average density or column density which is a few times larger than the threshold value. Ballesteros-Paredes et al. (2012) demonstrate that this is due to the power-law density PDF of clouds. The amount of mass within a certain density bin decreases rapidly with increasing density. The cumulative mass above a certain density threshold is thus quite sensitive to where you make the cut in the PDF. This has indeed been demonstrated in observations (Lombardi et al. 2010; Beaumont et al. 2012).

### 3.3. Cloud shapes

We characterise the shape of a cloud by considering the ratio of the half-axes $a \le b \le c$ of the ellipsoid approximation of the cloud. We compute them from the matrix of inertia of the cloud, assuming the cloud is an ellipsoid with an uniform density (see Eq. (A.5)). The ratio of the shortest half-axis $a$ to the longest $c$ defines how spherical the cloud is, and is called the sphericity. The ratio $b/c$ tells us whether the cloud is oblate, that is pancake-shaped or disk-shaped ($a \ll b \sim c$) or prolate, that is filamentary shaped ($a \sim b \ll c$). This is illustrated in Fig. 7a where we apply the shape categorisation from van der Wel et al. (2014). We remark that, as we already discussed in Sect. 2.3, the 3D cloud structures are very complex, and the ellipsoidal approximation fails to account for the very different shapes that clouds

can actually have. However they are a useful simplification to obtain the extents of the clouds and help classifying them.

### 3.3.1. Comparison of the typical shapes between simulations

In Fig. 7 we compare the distributions of the shape quantities, which are very similar for all the simulations, except for Ramses-F20. The typical sphericity of the clouds is around 0.35 with values ranging from 0.05 to 0.8, and the middle half-axis is typically equal to 0.55 times the large one. This corresponds to rather prolate structures, such as filament hubs, bent filaments, merging disks and compound structures, or unresolved blobs. We verified that this similarity between all the simulations is not an artefact caused by the cloud finding algorithm: HOP is indeed able to extract objects with largely varying shapes (see Appendix A.7). If we look in more detail, we note that filaments, defined as structures with $a \approx b \ll c$, are common as can be seen in bottom left corner of the diagrams in Fig. 7a. In contrast, spherical clouds are rare as indicated by the emptiness in the top right corner of the diagrams. This indicates that most of the structures are connected with their environment, the filamentary ISM. Oblate disk or sheet shaped structures, that is structures with $a \ll b \approx c$, are also very rare. Ramses-F20 shows a clear excess of them compared to the other simulations. The presence of such disk-like structures is most likely linked with a lack of resolution that prevents further collapse in these regions as already discussed in the previous sections.

Overall, while there is some variety in sphericity within cloud populations, the different simulations feature very similar distributions of the shape parameter. Since here the shape of the clouds is approximated to match an ellipsoid, we can conclude that the global proportions of the clouds are not affected by the details of the physics and the geometry of the simulations. More elaborate measures, such as the fractal index, may provide more insights on the detailed structure of the cloud. This is however beyond the scope of this work.

### 3.3.2. Discussion of example clouds

To aid with the interpretation of the shape diagnostics and to improve our intuition, we discuss here the shapes of the example clouds displayed in Fig. 2. We start with the clouds extracted from the galaxies Ramses-F20 (fourth row) and M51 (third row) which are generally most massive and illustrate some extreme shapes, before moving to those extracted from the stratified boxes
<table><caption>**Table 2.** Compilation of observational studies, with summary of their underlying data used to produce the molecular cloud catalogue, and their reported mass-size relation.</caption><thead><tr><th>References</th><th>Data type</th><th>Milky Way region</th><th>M-R slope</th></tr></thead><tbody><tr><td>Rice et al. (2016)</td><td><sup>12</sup>CO(1–0)</td><td>Full</td><td>1.98</td></tr><tr><td>Miville-Deschênes et al. (2017)</td><td><sup>12</sup>CO(1–0)</td><td>Full</td><td>2.2</td></tr><tr><td>Roman-Duval et al. (2010)</td><td><sup>12/13</sup>CO(1–0)</td><td>Galactic ring</td><td>2.36</td></tr><tr><td>Benedettini et al. (2020)</td><td><sup>12</sup>CO(1–0)</td><td>Outer Galaxy</td><td>≈2</td></tr><tr><td>Benedettini et al. (2021)</td><td><sup>13</sup>CO(1–0)</td><td>Outer Galaxy</td><td>2.13</td></tr><tr><td>Duarte-Cabral et al. (2021)</td><td><sup>13</sup>CO(2–1)</td><td>Inner Galaxy</td><td>≈2.3</td></tr><tr><td>Ma et al. (2022)</td><td><sup>12/13</sup>CO(1–0)</td><td>Third quadrant</td><td>2.24</td></tr><tr><td>Chen et al. (2020)</td><td>3D dust (projected)</td><td>Solar neighbourhood</td><td>1.96</td></tr><tr><td>Dharmawardena et al. (2023)</td><td>3D dust</td><td>Solar neighbourhood</td><td>≈3</td></tr><tr><td>Cahlon et al. (2024)</td><td>3D dust</td><td>Solar neighbourhood</td><td>2.9</td></tr><tr><td>Cahlon et al. (2024)</td><td>3D dust (projected)</td><td>Solar neighbourhood</td><td>2.1</td></tr></tbody></table>

### 3.2.2. Normalisation of M-R relation

When $M = AR^2$ the normalisation A provides a roughly constant surface density for clouds, or alternatively a constant average volume density when $M = AR^3$. Lada & Dame (2020) reported a systematic shift of the M-R relation to higher surface densities when clouds have been extracted from dust extinction versus $^{12}$CO, in other words the normalisation of the relation is different for different tracers. In fact, for any area in the sky the mass derived from CO is lower than the mass derived from cold dust. Small differences can also be observed between results derived from $^{12}$CO and $^{13}$CO (Benedettini et al. 2021). This discrepancy could be (partially) due to optical depth effects, which can be severe for $^{12}$CO as mentioned already in Sect. 2.4. If the mass is mainly set by the low density material in the cloud envelop rather than the dense cores and filaments, the effect of optical depth on the mass determination should be limited. In this case, the normalisation is set by the density threshold of the cloud extraction algorithm or the equivalent sensitivity limit of the observations, as we concluded from Fig. 6b. The threshold sets a characteristic average density or column density which is a few times larger than the threshold value. Ballesteros-Paredes et al. (2012) demonstrate that this is due to the power-law density PDF of clouds. The amount of mass within a certain density bin decreases rapidly with increasing density. The cumulative mass above a certain density threshold is thus quite sensitive to where you make the cut in the PDF. This has indeed been demonstrated in observations (Lombardi et al. 2010; Beaumont et al. 2012).

### 3.3. Cloud shapes

We characterise the shape of a cloud by considering the ratio of the half-axes $a \le b \le c$ of the ellipsoid approximation of the cloud. We compute them from the matrix of inertia of the cloud, assuming the cloud is an ellipsoid with an uniform density (see Eq. (A.5)). The ratio of the shortest half-axis $a$ to the longest $c$ defines how spherical the cloud is, and is called the sphericity. The ratio $b/c$ tells us whether the cloud is oblate, that is pancake-shaped or disk-shaped ($a \ll b \sim c$) or prolate, that is filamentary shaped ($a \sim b \ll c$). This is illustrated in Fig. 7a where we apply the shape categorisation from van der Wel et al. (2014). We remark that, as we already discussed in Sect. 2.3, the 3D cloud structures are very complex, and the ellipsoidal approximation fails to account for the very different shapes that clouds

can actually have. However they are a useful simplification to obtain the extents of the clouds and help classifying them.

### 3.3.1. Comparison of the typical shapes between simulations

In Fig. 7 we compare the distributions of the shape quantities, which are very similar for all the simulations, except for Ramses-F20. The typical sphericity of the clouds is around 0.35 with values ranging from 0.05 to 0.8, and the middle half-axis is typically equal to 0.55 times the large one. This corresponds to rather prolate structures, such as filament hubs, bent filaments, merging disks and compound structures, or unresolved blobs. We verified that this similarity between all the simulations is not an artefact caused by the cloud finding algorithm: HOP is indeed able to extract objects with largely varying shapes (see Appendix A.7). If we look in more detail, we note that filaments, defined as structures with $a \approx b \ll c$, are common as can be seen in bottom left corner of the diagrams in Fig. 7a. In contrast, spherical clouds are rare as indicated by the emptiness in the top right corner of the diagrams. This indicates that most of the structures are connected with their environment, the filamentary ISM. Oblate disk or sheet shaped structures, that is structures with $a \ll b \approx c$, are also very rare. Ramses-F20 shows a clear excess of them compared to the other simulations. The presence of such disk-like structures is most likely linked with a lack of resolution that prevents further collapse in these regions as already discussed in the previous sections.

Overall, while there is some variety in sphericity within cloud populations, the different simulations feature very similar distributions of the shape parameter. Since here the shape of the clouds is approximated to match an ellipsoid, we can conclude that the global proportions of the clouds are not affected by the details of the physics and the geometry of the simulations. More elaborate measures, such as the fractal index, may provide more insights on the detailed structure of the cloud. This is however beyond the scope of this work.

### 3.3.2. Discussion of example clouds

To aid with the interpretation of the shape diagnostics and to improve our intuition, we discuss here the shapes of the example clouds displayed in Fig. 2. We start with the clouds extracted from the galaxies Ramses-F20 (fourth row) and M51 (third row) which are generally most massive and illustrate some extreme shapes, before moving to those extracted from the stratified boxes
<figure><img src="image_8.png" /><figcaption>**Fig. 7.** General shapes of the clouds within the simulations. Panel a: relation between the ratio of the middle half-axis over the largest *b/c* and the ratio of the smallest half-axis over the largest *a/c* assuming clouds have a ellipsoidal shape. The characterisation of the shape the one used in van der Wel et al. (2014). The contour lines and symbols are the same as in Fig. 6a. Panel b: marginal distributions of the ratio of the ellipsoidal half-axes of the clouds between the simulations.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>**Fig. 7.** General shapes of the clouds within the simulations. Panel a: relation between the ratio of the middle half-axis over the largest *b/c* and the ratio of the smallest half-axis over the largest *a/c* assuming clouds have a ellipsoidal shape. The characterisation of the shape the one used in van der Wel et al. (2014). The contour lines and symbols are the same as in Fig. 6a. Panel b: marginal distributions of the ratio of the ellipsoidal half-axes of the clouds between the simulations.</figcaption></figure>
LS (second row) and SILCC (first row) which have properties
that are more representative for the entire sample.

The most massive example for each galaxy (left panels, and
indicated as a flat vertical diamond in the scatter plots) is a spu-
rious disk with spiral arm pattern. These are amongst the most
massive clouds found in their respective simulation. Their disk-
like shape is reflected by the shape measures: the Ramses-F20
cloud has *a/c* as low as 0.1 with *b/c* ≈ 0.75, while the M51 cloud
is a little thicker with *a/c* ≈ 0.2 and *b/c* ≈ 0.95. This indeed clas-
sifies their shape as disky in Fig. 7a. Both these examples have a
velocity dispersion high above the global relation due to ordered
rotation.

Next, we take the second and third examples of each galaxy
simulations, resulting in a collection of four clouds with com-
parable mass but a variety of shapes. The Ramses-F20 cloud
indicated by a triangle appears very isolated and has a smooth
appearance with a density gradient towards the centre. It is an
categorised as an oblate spheroid with *a/c* ≈ 0.3 and *b/c* ≈ 0.85.
The second example from M51 looks completely different. It
consists of a filament hub system where one major filament is
perpendicular to the other. The ellipsoidal shape approximation
clearly breaks down in this case. While *a/c* is indeed very low,
*b/c* ≈ 0.55. This inflates the average size which is almost four
times larger than the smooth compact Ramses-F20 cloud of the
same mass. This illustrates how the shape alters the mass-size
relation in Fig. 6a. If, on the other hand, we estimate the size as
the cubic root of the volume, both clouds appear very similar (see
Fig. 6b). The third example from each galaxy is a filament which
is connected to a larger filament network. The example from
M51 is one of the purest filaments in our entire collection with
very low *a/c* and *b/c*. The Ramses-F20 filament is somewhat
thicker with *a/c* ≈ 0.15 and *b/c* ≈ 0.25. Both of these seem to
be very normal clouds when placed on the mass-size (cf. Fig 6a)
and size-velocity dispersion relation (cf. Fig. 9, which will be
discussed in Sect. 4.1).

The final example for each galaxy features a small round
cloud. For Ramses-F20 this is one of the smallest clouds recov-
ered in this simulation. With *a/c* ≈ 0.75 and *b/c* ≈ 0.95 it is
also one of the most spherical objects in the entire collection.
We note that these are also the least resolved clouds, in which
further substructure is subsequently harder to resolve. Extreme
non-spherical shapes are therefore less likely to be driven from a
natural cascade of turbulent motions.

The example clouds from the stratified box simulations have less extreme shapes. We see that clouds are shaped by a network of connecting filaments, which generally do not align, resulting in the range of *a/c* and *b/c* values shown in Fig. 7b. Of note is the fourth example of the LS runs: a comma shapes clump which has been stripped of its low density envelope, making it lie above both versions of the mass-size relation. Generally, a similar argument applies here for the resolution limit and the shape of the cloud. However, this effect is less severe compared to the two disk simulations for two reasons. First, in these simulations we find overall better resolved feedback physics and consequently better resolved turbulent motions. Second, the higher resolution and the overall lower masses for the clouds allows for a stronger impact of sudden external perturbations such as strong shock fronts, which might reshape the cloud to non-spherical objects.

3.3.3. Dependence of shape on size

From the examples discussed above, we might have gotten the
impression that small clouds tend to be spherical. This is how-
ever not generally true. Figure 8a shows the relation between the
A155, page 14 of 29

<figure><img src="image_9.png" /><figcaption>**Fig. 8.** Sphericity as a function of structure size, computed from the volume (a) or the number of cells (b). The solid line depicts the mean in a logarithmic bin of size, and is computed only if the bin contains at least 100 clouds. Contour lines and symbols are as in Fig. 6a.</figcaption></figure>

sphericity of a cloud and the typical size defined as the cubic
root of the volume. We do not use the average size defined in
the beginning of this section to avoid introducing any spurious
correlation. For all the stratified box simulations, the distribu-
tion stays unchanged for all size bins between 1 and 10 parsec.
Interestingly, the mean value of the sphericity decreases when
the typical size is around 10 pc. A possible interpretation is that
large clouds are more likely to undergo external influence, like
turbulent motion. Also, they may be younger objects which kept
their initial filamentary shape. It may also be just a resolution
effect: clouds with a smaller number of cells tend to be more
spherical while a higher number of cells allows for more com-
plex shapes. However, Fig. 8b demonstrates that this is not the
case: expect for M51, the distribution of the sphericity does not
vary much when looking at objects with a different number of
cells. Although clouds with very high sphericity are found only
when the number of cells is small, they represent only a tiny frac-
tion of the sphericity distribution, which is well sampled for all
cell numbers. In the stratified box simulations, there is no signif-
icant decrease of the sphericity with an increase of the number
LS (second row) and SILCC (first row) which have properties
that are more representative for the entire sample.

The most massive example for each galaxy (left panels, and
indicated as a flat vertical diamond in the scatter plots) is a spu-
rious disk with spiral arm pattern. These are amongst the most
massive clouds found in their respective simulation. Their disk-
like shape is reflected by the shape measures: the Ramses-F20
cloud has *a/c* as low as 0.1 with *b/c* ≈ 0.75, while the M51 cloud
is a little thicker with *a/c* ≈ 0.2 and *b/c* ≈ 0.95. This indeed clas-
sifies their shape as disky in Fig. 7a. Both these examples have a
velocity dispersion high above the global relation due to ordered
rotation.

Next, we take the second and third examples of each galaxy
simulations, resulting in a collection of four clouds with com-
parable mass but a variety of shapes. The Ramses-F20 cloud
indicated by a triangle appears very isolated and has a smooth
appearance with a density gradient towards the centre. It is an
categorised as an oblate spheroid with a/c ≈ 0.3 and b/c ≈ 0.85.
The second example from M51 looks completely different. It
consists of a filament hub system where one major filament is
perpendicular to the other. The ellipsoidal shape approximation
clearly breaks down in this case. While a/c is indeed very low,
b/c ≈ 0.55. This inflates the average size which is almost four
times larger than the smooth compact Ramses-F20 cloud of the
same mass. This illustrates how the shape alters the mass-size
relation in Fig. 6a. If, on the other hand, we estimate the size as
the cubic root of the volume, both clouds appear very similar (see
Fig. 6b). The third example from each galaxy is a filament which
is connected to a larger filament network. The example from
M51 is one of the purest filaments in our entire collection with
very low a/c and b/c. The Ramses-F20 filament is somewhat
thicker with a/c ≈ 0.15 and b/c ≈ 0.25. Both of these seem to
be very normal clouds when placed on the mass-size (cf. Fig 6a)
and size-velocity dispersion relation (cf. Fig. 9, which will be
discussed in Sect. 4.1).

The final example for each galaxy features a small round
cloud. For Ramses-F20 this is one of the smallest clouds recov-
ered in this simulation. With *a/c* ≈ 0.75 and *b/c* ≈ 0.95 it is
also one of the most spherical objects in the entire collection.
We note that these are also the least resolved clouds, in which
further substructure is subsequently harder to resolve. Extreme
non-spherical shapes are therefore less likely to be driven from a
natural cascade of turbulent motions.

The example clouds from the stratified box simulations have less extreme shapes. We see that clouds are shaped by a network of connecting filaments, which generally do not align, resulting in the range of *a/c* and *b/c* values shown in Fig. 7b. Of note is the fourth example of the LS runs: a comma shapes clump which has been stripped of its low density envelope, making it lie above both versions of the mass-size relation. Generally, a similar argument applies here for the resolution limit and the shape of the cloud. However, this effect is less severe compared to the two disk simulations for two reasons. First, in these simulations we find overall better resolved feedback physics and consequently better resolved turbulent motions. Second, the higher resolution and the overall lower masses for the clouds allows for a stronger impact of sudden external perturbations such as strong shock fronts, which might reshape the cloud to non-spherical objects.

3.3.3. Dependence of shape on size

From the examples discussed above, we might have gotten the
impression that small clouds tend to be spherical. This is how-
ever not generally true. Figure 8a shows the relation between the
A155, page 14 of 29

<figure><img src="image_9.png" /><figcaption>**Fig. 8.** Sphericity as a function of structure size, computed from the volume (a) or the number of cells (b). The solid line depicts the mean in a logarithmic bin of size, and is computed only if the bin contains at least 100 clouds. Contour lines and symbols are as in Fig. 6a.</figcaption></figure>

sphericity of a cloud and the typical size defined as the cubic
root of the volume. We do not use the average size defined in
the beginning of this section to avoid introducing any spurious
correlation. For all the stratified box simulations, the distribu-
tion stays unchanged for all size bins between 1 and 10 parsec.
Interestingly, the mean value of the sphericity decreases when
the typical size is around 10 pc. A possible interpretation is that
large clouds are more likely to undergo external influence, like
turbulent motion. Also, they may be younger objects which kept
their initial filamentary shape. It may also be just a resolution
effect: clouds with a smaller number of cells tend to be more
spherical while a higher number of cells allows for more com-
plex shapes. However, Fig. 8b demonstrates that this is not the
case: expect for M51, the distribution of the sphericity does not
vary much when looking at objects with a different number of
cells. Although clouds with very high sphericity are found only
when the number of cells is small, they represent only a tiny frac-
tion of the sphericity distribution, which is well sampled for all
cell numbers. In the stratified box simulations, there is no signif-
icant decrease of the sphericity with an increase of the number
<figure><img src="image_10.png" /><figcaption>Fig. 9. Cloud velocity dispersion as a function of the size. The dashed line illustrates a scaling of $\sigma$ with size $R^{0.75}$ while the curved solid line illustrates the trend of the centroids of each distributions. The contour lines and symbols are the same as in Fig. 6.</figcaption></figure>

of cells. Such a trend is however visible for the less resolved full-galaxy simulations, especially in M51, and in a lesser extent in Ramses-F20.

### 3.3.4. Shapes in observations

In observations, clouds can be approximated by the 2D equivalent of our ellipsoids, ellipses, since they appear as 2D projections on the sky. In this case a minor and major axis is reported and their ratio incorporates the shape. Since there is no information about the size along the line of sight<fnref n="2" />, it is unclear whether this ratio is more comparable to our *a/c* or *b/c*. If we assume that clouds are randomly oriented, then the observed value would represent an average of *a/c*, *b/c* and *a/b*. Taking the geometric mean of the peak value of the distribution for each of these ratios results in $(0.35 \times 0.55 \times 0.65)^{1/3} = 0.5$. The SEDIGISM (Duarte-Cabral et al. 2021), FQS (Benedettini et al. 2020, 2021) and Miville-Deschênes et al. (2017) catalogues provide the two projected sizes for their objects. Their ratio is typically between 0.2 and 0.8 and the distribution peaks between 0.4 and 0.6. For SEDIGISM and FQS, which have extracted clouds in a similar way, the peak is closer to 0.45, in line with the averaged value we find. For the Miville-Deschênes et al. (2017) catalogue, where the sizes were derived differently, the peak is more towards 0.6, close to what we find for *b/c*. As a test, we can also project our 3D clouds on a 2D plane and calculate the projected size as a geometric average of the major and minor axis of an approximate ellipse, equivalently to our ellipsoid approximations. This results indeed in a broad distribution with a peak around 0.5. Thus,

while an accurate comparison for the shape seems rather impossible, we can confidently say that simulation and observations agree that clouds are typically elongated.

## 4. Analysis of internal properties of clouds

### 4.1. Size – velocity dispersion relation

The velocity dispersion is also correlated with the cloud size (Hennebelle & Falgarone 2012; Heyer & Dame 2015). This is the so-called Larson relation (Larson 1981), which is considered a signature of the turbulent cascade (Elmegreen & Scalo 2004). Indeed, for compressible supersonic turbulence a scaling of $\sigma \propto l^\beta$ with $\beta$ between 1/3 and 1/2 is expected (Kritsuk et al. 2007). Observationally, a large scatter is found in the line-width – size relation. Some studies agree with this theoretically predicted slope (Solomon et al. 1987; Kauffmann et al. 2013), but some find steeper slopes (Miville-Deschênes et al. 2017), which shows that the relation may not be as universal as generally believed.

The velocity dispersion of our extracted clouds has been calculated as a mass averaged quantity (see Eq. (A.8)). The size is derived from the ellipsoid approximation as mentioned in the previous section and detailed in Appendix A.

#### 4.1.1. Continuous relation across simulations

Figure 9 shows the cloud internal velocity dispersion as a function of their size for all simulations. The dashed line illustrates the power-law scaling with 0.75. In addition, we show the actual scaling with the solid line that passes through the innermost contour line of all sets of clouds. Despite the overall large scatter, the majority of the clouds (inner contours) form a continuous distribution across the set of simulations.

We note that the distribution is bounded towards the low-$\sigma$ end approximately half an order of magnitude below the dashed line. There are hardly any clouds with a lower velocity dispersion. This is in line with the turbulent motions in the gas and the gravitational collapse in massive clouds, which increases the internal velocity dispersion. In contrast, there is no clear cut at the high-$\sigma$ end of the distribution. In particular for the small clouds (mainly the SILCC simulations) and the clouds from galaxy Ramses-F20, there is a population of clouds with strong internal motions. Overall, for small clouds the internal velocity dispersion is significantly enhanced compared to the indicated power-law scaling – in particular the inner contours of the cloud distributions from the SILCC simulations.

#### 4.1.2. Observed scaling and break in the power law

Canonically, the size–velocity dispersion relation has been described by a power law with slope 0.5. However, the simulation results indicate that the relation is somewhat more complex, in a sense that it in fact does not follow a simple power law. In Fig. 10 we compare our results to the catalogues from Miville-Deschênes et al. (2017) and Benedettini et al. (2020). We keep in mind that the velocity dispersion is determined in a very different way in the simulations compared to the observations. Nonetheless, the observations of Benedettini et al. (2020) show a flattening of the relation below sizes of about 5 pc, with the velocity dispersion saturating at a value of about 1.4 km s⁻¹ (when converted from 1D to 3D dispersion assuming isotropy), in agreement with the global trend found in the simulations. The bulk of the SILCC clouds even seem to slightly underestimate

<footnote n="2">² This information is in principle retrievable from clouds detected in 3D dust. However, so far no study has reported sizes along three orthogonal axes.</footnote>
<figure><img src="image_10.png" /><figcaption>Fig. 9. Cloud velocity dispersion as a function of the size. The dashed line illustrates a scaling of $\sigma$ with size $R^{0.75}$ while the curved solid line illustrates the trend of the centroids of each distributions. The contour lines and symbols are the same as in Fig. 6.</figcaption></figure>

of cells. Such a trend is however visible for the less resolved full-galaxy simulations, especially in M51, and in a lesser extent in Ramses-F20.

### 3.3.4. Shapes in observations

In observations, clouds can be approximated by the 2D equivalent of our ellipsoids, ellipses, since they appear as 2D projections on the sky. In this case a minor and major axis is reported and their ratio incorporates the shape. Since there is no information about the size along the line of sight<fnref n="2" />, it is unclear whether this ratio is more comparable to our *a/c* or *b/c*. If we assume that clouds are randomly oriented, then the observed value would represent an average of *a/c*, *b/c* and *a/b*. Taking the geometric mean of the peak value of the distribution for each of these ratios results in $(0.35 \times 0.55 \times 0.65)^{1/3} = 0.5$. The SEDIGISM (Duarte-Cabral et al. 2021), FQS (Benedettini et al. 2020, 2021) and Miville-Deschênes et al. (2017) catalogues provide the two projected sizes for their objects. Their ratio is typically between 0.2 and 0.8 and the distribution peaks between 0.4 and 0.6. For SEDIGISM and FQS, which have extracted clouds in a similar way, the peak is closer to 0.45, in line with the averaged value we find. For the Miville-Deschênes et al. (2017) catalogue, where the sizes were derived differently, the peak is more towards 0.6, close to what we find for *b/c*. As a test, we can also project our 3D clouds on a 2D plane and calculate the projected size as a geometric average of the major and minor axis of an approximate ellipse, equivalently to our ellipsoid approximations. This results indeed in a broad distribution with a peak around 0.5. Thus,

while an accurate comparison for the shape seems rather impossible, we can confidently say that simulation and observations agree that clouds are typically elongated.

## 4. Analysis of internal properties of clouds

### 4.1. Size – velocity dispersion relation

The velocity dispersion is also correlated with the cloud size (Hennebelle & Falgarone 2012; Heyer & Dame 2015). This is the so-called Larson relation (Larson 1981), which is considered a signature of the turbulent cascade (Elmegreen & Scalo 2004). Indeed, for compressible supersonic turbulence a scaling of $\sigma \propto l^\beta$ with $\beta$ between 1/3 and 1/2 is expected (Kritsuk et al. 2007). Observationally, a large scatter is found in the line-width – size relation. Some studies agree with this theoretically predicted slope (Solomon et al. 1987; Kauffmann et al. 2013), but some find steeper slopes (Miville-Deschênes et al. 2017), which shows that the relation may not be as universal as generally believed.

The velocity dispersion of our extracted clouds has been calculated as a mass averaged quantity (see Eq. (A.8)). The size is derived from the ellipsoid approximation as mentioned in the previous section and detailed in Appendix A.

#### 4.1.1. Continuous relation across simulations

Figure 9 shows the cloud internal velocity dispersion as a function of their size for all simulations. The dashed line illustrates the power-law scaling with 0.75. In addition, we show the actual scaling with the solid line that passes through the innermost contour line of all sets of clouds. Despite the overall large scatter, the majority of the clouds (inner contours) form a continuous distribution across the set of simulations.

We note that the distribution is bounded towards the low-$\sigma$ end approximately half an order of magnitude below the dashed line. There are hardly any clouds with a lower velocity dispersion. This is in line with the turbulent motions in the gas and the gravitational collapse in massive clouds, which increases the internal velocity dispersion. In contrast, there is no clear cut at the high-$\sigma$ end of the distribution. In particular for the small clouds (mainly the SILCC simulations) and the clouds from galaxy Ramses-F20, there is a population of clouds with strong internal motions. Overall, for small clouds the internal velocity dispersion is significantly enhanced compared to the indicated power-law scaling – in particular the inner contours of the cloud distributions from the SILCC simulations.

#### 4.1.2. Observed scaling and break in the power law

Canonically, the size–velocity dispersion relation has been described by a power law with slope 0.5. However, the simulation results indicate that the relation is somewhat more complex, in a sense that it in fact does not follow a simple power law. In Fig. 10 we compare our results to the catalogues from Miville-Deschênes et al. (2017) and Benedettini et al. (2020). We keep in mind that the velocity dispersion is determined in a very different way in the simulations compared to the observations. Nonetheless, the observations of Benedettini et al. (2020) show a flattening of the relation below sizes of about 5 pc, with the velocity dispersion saturating at a value of about 1.4 km s⁻¹ (when converted from 1D to 3D dispersion assuming isotropy), in agreement with the global trend found in the simulations. The bulk of the SILCC clouds even seem to slightly underestimate

<footnote n="2">² This information is in principle retrievable from clouds detected in 3D dust. However, so far no study has reported sizes along three orthogonal axes.</footnote>
<figure><img src="image_11.png" /><figcaption>**Fig. 10.** Comparison of the 3D size-velocity dispersion relation found in the simulation with several cloud catalogues from observations. Coloured contours are as in Fig. 6. Observational data feature only the 1D line-of-sight velocity dispersion, which was here multiplied by $\sqrt{3}$ for comparison with the 3D velocity dispersion computed in the simulation.</figcaption></figure>

the small scale velocity dispersion. This is not surprising since the simulations do not take into account processes that inject turbulence at small scales, such as stellar winds and jets. Another possible interpretation is that the observations are biased towards gravitationally collapsing objects. The outermost contour of LS, simulations which preferentially resolve the collapsing objects within them as discussed in Sect. 5.1, indeed seems to follow the observations.

While we do recover a global power-law behaviour for large clouds, the slope from the simulations is somewhat steeper, that is to say about 0.75, compared to the already fairly steep 0.6 found by Miville-Deschênes et al. (2017). For large clouds, the extraction from the simulations measures on average a slightly larger velocity dispersion than the observations from Miville-Deschênes et al. (2017) for a certain size bin. The simulation results may overestimate the dispersion due to several reasons. We calculate the dispersion in a somewhat naive way which does not take into account organised motions within the cloud such as rotation, radial gravitational inflow or colliding flows. These are not forms of turbulence and should thus be accounted for when one wants to quantify the turbulence velocity dispersion. An alternative interpretation is a systematically different size estimate, such that the sizes for the observed clouds are larger than the ones from the simulations. This is in fact a plausible explanation given that clouds are overlapping along the line of sight as illustrated in the simulations in Fig. 2 and cloud blending can be seen especially at large distances.

What the simulations and observations do agree on is the large scatter. For a specific size, the estimated velocity dispersion can vary by an order of magnitude.

### 4.1.3. Origin of high velocity dispersion clouds

A high velocity dispersion can have several origins, among which externally driven turbulence and self-gravitating collapse

are two prominent examples. In the case of gravitational collapse the clouds are likely to be embedded in a dense environment with cold temperatures. In the case of externally driven turbulence, that propagates into the clouds, the gas temperature can be enhanced if the external turbulence is driven by SN feedback and penetrates into the clouds from the warm or hot phase of the ISM.

In order to identify the mechanism at work for the small clouds we investigate the locations of the clouds in Fig. 1a for simulation SILCC-0.5pc. The corresponding SILCC simulations with higher and lower resolution are very similar. The grey-scale colour-map shows the gas column density. The clouds identified by HOP are over-plotted with the colour indicating the internal velocity dispersion. There is no clear correlation between the position of the clouds in the global gas structure and their velocity dispersion. However, the lowest values of $\sigma$ are found for points that tend to be more embedded in large gas structures. On the other hand, clouds closer to the low-density bubbles that have formed as a result of SN feedback have intermediate to high velocity dispersions. We further illustrate this in Fig. 11, where we show slices of density, temperature and velocity at the positions of two example clouds with high velocity dispersion, both located at the interface between cold gas and a hot bubble.

A more systematic analysis is shown in Fig. 12 for the clouds in SILCC-0.5 pc (again, different resolutions do not differ) for which we investigate a spherical region around the clouds with a radius twice the average size obtained from the cloud algorithm. Both panels show scatter plots of $\sigma$ as a function of cloud mass. In the left-hand panel the clouds are colour coded by the mass-weighted temperature in the investigated volume. The right-hand counterpart shows the fraction of molecular gas in the volume. We note that the clouds with higher velocity dispersion clearly show enhanced temperatures and lower fractions of molecular hydrogen. This illustrates that the clouds which are surrounded by warm or hot gas have statistically higher velocity dispersions. In the SILCC setups analysed here, the only source of hot turbulent gas are SN explosions. We therefore conclude that the enhanced motions in the clouds are due to SN-driven turbulence. We note that these warm clouds with low $H_2$ fractions would also contain little CO and as a result could be missing from classical observational catalogues.

In the full galaxy simulations RAMSES-F20 and M51, there is also a population of clouds with high-velocity dispersion, but the mechanisms at play are different. Most of these clouds, which are also very large, are rotating disks, which explains the high velocity dispersion. The existence of such disks with radius of around 100 pc is a spurious effect of the lack of resolution in these regions, as we would expect these structures to collapse and fragment into smaller objects.

### 4.2. Virial parameter

In addition to the size-velocity dispersion relationship, we also compute the virial parameter $\alpha_{vir}$ for the clouds in our simulations, using the definition given in Eq. (1). The results are shown in Fig. 13. We can see immediately that $\alpha_{vir}$ has a clear dependence on mass: $\alpha_{vir} \propto M^{-n}$, with $n = 0.4$. In other words, lower mass clouds on average have much larger values of $\alpha_{vir}$ than higher mass clouds. Therefore, although there is generally good agreement between the simulations for the mass ranges where their cloud populations overlap, the mean value of $\alpha_{vir}$ differs significantly between different simulations. This demonstrates that care must be taken when using simulations to draw conclusions about the gravitational boundedness of interstellar
<figure><img src="image_11.png" /><figcaption>**Fig. 10.** Comparison of the 3D size-velocity dispersion relation found in the simulation with several cloud catalogues from observations. Coloured contours are as in Fig. 6. Observational data feature only the 1D line-of-sight velocity dispersion, which was here multiplied by $\sqrt{3}$ for comparison with the 3D velocity dispersion computed in the simulation.</figcaption></figure>

the small scale velocity dispersion. This is not surprising since the simulations do not take into account processes that inject turbulence at small scales, such as stellar winds and jets. Another possible interpretation is that the observations are biased towards gravitationally collapsing objects. The outermost contour of LS, simulations which preferentially resolve the collapsing objects within them as discussed in Sect. 5.1, indeed seems to follow the observations.

While we do recover a global power-law behaviour for large clouds, the slope from the simulations is somewhat steeper, that is to say about 0.75, compared to the already fairly steep 0.6 found by Miville-Deschênes et al. (2017). For large clouds, the extraction from the simulations measures on average a slightly larger velocity dispersion than the observations from Miville-Deschênes et al. (2017) for a certain size bin. The simulation results may overestimate the dispersion due to several reasons. We calculate the dispersion in a somewhat naive way which does not take into account organised motions within the cloud such as rotation, radial gravitational inflow or colliding flows. These are not forms of turbulence and should thus be accounted for when one wants to quantify the turbulence velocity dispersion. An alternative interpretation is a systematically different size estimate, such that the sizes for the observed clouds are larger than the ones from the simulations. This is in fact a plausible explanation given that clouds are overlapping along the line of sight as illustrated in the simulations in Fig. 2 and cloud blending can be seen especially at large distances.

What the simulations and observations do agree on is the large scatter. For a specific size, the estimated velocity dispersion can vary by an order of magnitude.

### 4.1.3. Origin of high velocity dispersion clouds

A high velocity dispersion can have several origins, among which externally driven turbulence and self-gravitating collapse

are two prominent examples. In the case of gravitational collapse the clouds are likely to be embedded in a dense environment with cold temperatures. In the case of externally driven turbulence, that propagates into the clouds, the gas temperature can be enhanced if the external turbulence is driven by SN feedback and penetrates into the clouds from the warm or hot phase of the ISM.

In order to identify the mechanism at work for the small clouds we investigate the locations of the clouds in Fig. 1a for simulation SILCC-0.5pc. The corresponding SILCC simulations with higher and lower resolution are very similar. The grey-scale colour-map shows the gas column density. The clouds identified by HOP are over-plotted with the colour indicating the internal velocity dispersion. There is no clear correlation between the position of the clouds in the global gas structure and their velocity dispersion. However, the lowest values of $\sigma$ are found for points that tend to be more embedded in large gas structures. On the other hand, clouds closer to the low-density bubbles that have formed as a result of SN feedback have intermediate to high velocity dispersions. We further illustrate this in Fig. 11, where we show slices of density, temperature and velocity at the positions of two example clouds with high velocity dispersion, both located at the interface between cold gas and a hot bubble.

A more systematic analysis is shown in Fig. 12 for the clouds in SILCC-0.5 pc (again, different resolutions do not differ) for which we investigate a spherical region around the clouds with a radius twice the average size obtained from the cloud algorithm. Both panels show scatter plots of $\sigma$ as a function of cloud mass. In the left-hand panel the clouds are colour coded by the mass-weighted temperature in the investigated volume. The right-hand counterpart shows the fraction of molecular gas in the volume. We note that the clouds with higher velocity dispersion clearly show enhanced temperatures and lower fractions of molecular hydrogen. This illustrates that the clouds which are surrounded by warm or hot gas have statistically higher velocity dispersions. In the SILCC setups analysed here, the only source of hot turbulent gas are SN explosions. We therefore conclude that the enhanced motions in the clouds are due to SN-driven turbulence. We note that these warm clouds with low $H_2$ fractions would also contain little CO and as a result could be missing from classical observational catalogues.

In the full galaxy simulations RAMSES-F20 and M51, there is also a population of clouds with high-velocity dispersion, but the mechanisms at play are different. Most of these clouds, which are also very large, are rotating disks, which explains the high velocity dispersion. The existence of such disks with radius of around 100 pc is a spurious effect of the lack of resolution in these regions, as we would expect these structures to collapse and fragment into smaller objects.

### 4.2. Virial parameter

In addition to the size-velocity dispersion relationship, we also compute the virial parameter $\alpha_{vir}$ for the clouds in our simulations, using the definition given in Eq. (1). The results are shown in Fig. 13. We can see immediately that $\alpha_{vir}$ has a clear dependence on mass: $\alpha_{vir} \propto M^{-n}$, with $n = 0.4$. In other words, lower mass clouds on average have much larger values of $\alpha_{vir}$ than higher mass clouds. Therefore, although there is generally good agreement between the simulations for the mass ranges where their cloud populations overlap, the mean value of $\alpha_{vir}$ differs significantly between different simulations. This demonstrates that care must be taken when using simulations to draw conclusions about the gravitational boundedness of interstellar
<figure><img src="image_12.png" /><figcaption>**Fig. 11.** Two examples showing the positions of the clouds with their respective density (left), temperature (middle), and velocity environment (right). In both cases the clouds are located at the edge between cold gas and a hot SN-driven bubble. The resulting high velocity dispersion is thus due to SN-driven turbulence that mixes hot gas into the clouds rather than by gravitational collapse.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Fig. 12.** Scatter plot of the velocity distribution as a function of cloud mass for simulation SILCC-0.5pc. Left: colours indicate the mass-weighted temperature in a sphere around the centre of the HOP cloud with a radius corresponding to twice the average size of the HOP cloud. Right: colours indicate the molecular fraction in the surrounding sphere.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>**Fig. 11.** Two examples showing the positions of the clouds with their respective density (left), temperature (middle), and velocity environment (right). In both cases the clouds are located at the edge between cold gas and a hot SN-driven bubble. The resulting high velocity dispersion is thus due to SN-driven turbulence that mixes hot gas into the clouds rather than by gravitational collapse.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Fig. 12.** Scatter plot of the velocity distribution as a function of cloud mass for simulation SILCC-0.5pc. Left: colours indicate the mass-weighted temperature in a sphere around the centre of the HOP cloud with a radius corresponding to twice the average size of the HOP cloud. Right: colours indicate the molecular fraction in the surrounding sphere.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>**Fig. 13.** Virial parameter as a function of cloud mass for all simulations. Coloured contour lines are the same as in Fig. 6. Over-plotted are the slope of the scaling (thick grey line), as well as the observational data compiled for the PPVII review (Chevance et al. 2023).</figcaption></figure>

clouds<fnref n="3" />. Different simulations may lead to very different conclusions about whether or not clouds are bound if they differ significantly in terms of resolution or box size, and hence do not probe the same portion of the cloud mass hierarchy.

We also see that there is substantial scatter in the value of $\alpha_{vir}$ that we recover for any given cloud mass. For the majority of clouds, the value of $\alpha_{vir}$ varies by about one order of magnitude for a specific mass range, but in a small fraction of clouds the value of $\alpha_{vir}$ can be as much as 100 times the typical value for a given mass. At fixed mass, $\alpha_{vir}$ depends only on the cloud size $R$ and velocity dispersion $\sigma$. However, we have already seen that at fixed $M$, the cloud size varies by only around a factor of two in most cases (see Fig. 6), and so the majority of this scatter must be driven by variations in the velocity dispersion. This conclusion is consistent with what we have already found for the velocity dispersion-size relation: for most clouds, $\sigma$ varies by around a factor of three at fixed cloud size (see Fig. 9), corresponding to an order of magnitude variation in $\alpha_{vir}$ if we keep $M$ and $R$ fixed, but for a small number of clouds we find much larger values for $\sigma$ and hence for $\alpha_{vir}$. It hence seems likely that the factors responsible for the outliers in the $\alpha_{vir}-M$ distribution are the same as those responsible for the outliers in the $\sigma-R$ distribution, namely stellar feedback at low cloud masses and lack of resolution leading to artificial rotational support at high cloud masses (see Sect. 4.1).

Finally, it is interesting to compare the $\alpha_{vir}-M$ relation that we recover from the simulations with the one measured for real interstellar clouds. As is clear from the compilation of data in Chevance et al. (2023), different galactic and extra-galactic cloud surveys recover a similar scaling $\alpha_{vir} \propto M^{-0.5}$ at low cloud masses, but with a normalisation that varies between surveys by

up to an order of magnitude. These results are also in agreement with observations of the densest regions within molecular clouds, parsec-scales structures called clumps. Several surveys dedicated to the study of the dynamics at these scales using high density tracers (e.g. NH₃ or N₂H⁺) report a scaling $\alpha_{vir} \propto M^{-0.5}$ (Urquhart et al. 2018; Traficante et al. 2018a). Although the behaviour of this scaling relationship on the scale of individual clumps is debated (Kauffmann et al. 2013; Traficante et al. 2018b; Singh et al. 2021), it is remarkable to note that the average slope observed using high-density tracers is in overall agreement with CO cloud surveys and with the results of the simulations in this work. In fact, overall, the range of values found in observations are in good agreement with the location of most of our simulated clouds in the $\alpha_{vir}-M$ plane (see grey crosses in Fig. 13), but the small population of clouds we recover that have very high virial parameters do not appear to have observational counterparts in CO-based surveys. This may simply reflect the fact that many of these clouds have been strongly affected by stellar feedback and hence would in reality have little molecular gas associated with them. The slope that we recover for the $\alpha_{vir}-M$ relation also agrees well with the observed slope: our best fit value is slightly flatter than that found in the observations, but as we show in Fig. 13, a slope of -0.5 also provides a relatively good fit to the data from all of the simulations apart from Ramses-F20. However, we caution that this agreement may be strongly influenced by selection effects. Observationally, clouds with an integrated intensity close to the detection limit of a survey and a velocity dispersion close to the velocity resolution will all have approximately the same molecular gas surface density, implying that for these clouds, $R \propto M^{0.5}$. For fixed $\sigma$, this then yields $\alpha_{vir} \propto M^{-0.5}$. A similar combination of limited sensitivity and velocity resolution may also explain the similarity between the scaling relationships recovered for dense clumps and for molecular clouds. Our population of simulated clouds does not suffer from this selection effect, but the fact that, as previously discussed, HOP selects clouds with a relatively narrow range of mean densities also strongly biases us against finding clouds with low $\alpha_{vir}$ at low $M$.

### 4.3. Thermal state

In order to quantify the thermal state of the gas in the clouds identified in the different simulations, we have computed the mean temperature $T$ for each cloud. We compare the distributions of mass-weighted average cloud temperatures in Fig. 14, where we show results from the LS, M51 and Ramses-F20 simulations, plus a representative example from the SILCC simulations.

From the figure, we see that there are some significant differences in the average cloud temperatures in the simulations. All of the simulations show a clear peak at low temperatures, but they do not agree on the location of this peak: it is found at $T \sim 25$ K in the M51 simulation, $T \sim 100$ K in the LS and Ramses-F20 simulations, and at $T \sim 60$ K in the SILCC simulation. This appears to be primarily due to a difference in the effectiveness of photoelectric heating in the various simulations. In the LS and Ramses-F20 simulations, the photoelectric heating rate is computed assuming no local attenuation of the interstellar radiation field. Since much of the cloud mass is located at a density close to the HOP threshold density, we would therefore expect the temperature of this gas to be comparable to that of a patch of unshielded ISM with a density $n \sim 30 \text{ cm}^{-3}$. For an interstellar radiation field strength comparable to the solar neighbourhood value, this is $T \sim 100$ K (Wolfire et al. 2003), in good agreement with the low $T$ peak we find in these simulations. In

<footnote n="3">³ The question of whether the conventional definition of $\alpha_{vir}$ is an adequate diagnostic of whether or not clouds are gravitationally bound is an interesting one (see e.g. Ramírez-Galeano et al. 2022), but is out of the scope of this study.</footnote>
<figure><img src="image_14.png" /><figcaption>**Fig. 13.** Virial parameter as a function of cloud mass for all simulations. Coloured contour lines are the same as in Fig. 6. Over-plotted are the slope of the scaling (thick grey line), as well as the observational data compiled for the PPVII review (Chevance et al. 2023).</figcaption></figure>

clouds<fnref n="3" />. Different simulations may lead to very different conclusions about whether or not clouds are bound if they differ significantly in terms of resolution or box size, and hence do not probe the same portion of the cloud mass hierarchy.

We also see that there is substantial scatter in the value of $\alpha_{vir}$ that we recover for any given cloud mass. For the majority of clouds, the value of $\alpha_{vir}$ varies by about one order of magnitude for a specific mass range, but in a small fraction of clouds the value of $\alpha_{vir}$ can be as much as 100 times the typical value for a given mass. At fixed mass, $\alpha_{vir}$ depends only on the cloud size $R$ and velocity dispersion $\sigma$. However, we have already seen that at fixed $M$, the cloud size varies by only around a factor of two in most cases (see Fig. 6), and so the majority of this scatter must be driven by variations in the velocity dispersion. This conclusion is consistent with what we have already found for the velocity dispersion-size relation: for most clouds, $\sigma$ varies by around a factor of three at fixed cloud size (see Fig. 9), corresponding to an order of magnitude variation in $\alpha_{vir}$ if we keep $M$ and $R$ fixed, but for a small number of clouds we find much larger values for $\sigma$ and hence for $\alpha_{vir}$. It hence seems likely that the factors responsible for the outliers in the $\alpha_{vir}-M$ distribution are the same as those responsible for the outliers in the $\sigma-R$ distribution, namely stellar feedback at low cloud masses and lack of resolution leading to artificial rotational support at high cloud masses (see Sect. 4.1).

Finally, it is interesting to compare the $\alpha_{vir}-M$ relation that we recover from the simulations with the one measured for real interstellar clouds. As is clear from the compilation of data in Chevance et al. (2023), different galactic and extra-galactic cloud surveys recover a similar scaling $\alpha_{vir} \propto M^{-0.5}$ at low cloud masses, but with a normalisation that varies between surveys by

up to an order of magnitude. These results are also in agreement with observations of the densest regions within molecular clouds, parsec-scales structures called clumps. Several surveys dedicated to the study of the dynamics at these scales using high density tracers (e.g. NH$_3$ or N$_2$H$^+$) report a scaling $\alpha_{vir} \propto M^{-0.5}$ (Urquhart et al. 2018; Traficante et al. 2018a). Although the behaviour of this scaling relationship on the scale of individual clumps is debated (Kauffmann et al. 2013; Traficante et al. 2018b; Singh et al. 2021), it is remarkable to note that the average slope observed using high-density tracers is in overall agreement with CO cloud surveys and with the results of the simulations in this work. In fact, overall, the range of values found in observations are in good agreement with the location of most of our simulated clouds in the $\alpha_{vir}-M$ plane (see grey crosses in Fig. 13), but the small population of clouds we recover that have very high virial parameters do not appear to have observational counterparts in CO-based surveys. This may simply reflect the fact that many of these clouds have been strongly affected by stellar feedback and hence would in reality have little molecular gas associated with them. The slope that we recover for the $\alpha_{vir}-M$ relation also agrees well with the observed slope: our best fit value is slightly flatter than that found in the observations, but as we show in Fig. 13, a slope of -0.5 also provides a relatively good fit to the data from all of the simulations apart from Ramses-F20. However, we caution that this agreement may be strongly influenced by selection effects. Observationally, clouds with an integrated intensity close to the detection limit of a survey and a velocity dispersion close to the velocity resolution will all have approximately the same molecular gas surface density, implying that for these clouds, $R \propto M^{0.5}$. For fixed $\sigma$, this then yields $\alpha_{vir} \propto M^{-0.5}$. A similar combination of limited sensitivity and velocity resolution may also explain the similarity between the scaling relationships recovered for dense clumps and for molecular clouds. Our population of simulated clouds does not suffer from this selection effect, but the fact that, as previously discussed, HOP selects clouds with a relatively narrow range of mean densities also strongly biases us against finding clouds with low $\alpha_{vir}$ at low $M$.

### 4.3. Thermal state

In order to quantify the thermal state of the gas in the clouds identified in the different simulations, we have computed the mean temperature $T$ for each cloud. We compare the distributions of mass-weighted average cloud temperatures in Fig. 14, where we show results from the LS, M51 and Ramses-F20 simulations, plus a representative example from the SILCC simulations.

From the figure, we see that there are some significant differences in the average cloud temperatures in the simulations. All of the simulations show a clear peak at low temperatures, but they do not agree on the location of this peak: it is found at $T \sim 25$ K in the M51 simulation, $T \sim 100$ K in the LS and Ramses-F20 simulations, and at $T \sim 60$ K in the SILCC simulation. This appears to be primarily due to a difference in the effectiveness of photoelectric heating in the various simulations. In the LS and Ramses-F20 simulations, the photoelectric heating rate is computed assuming no local attenuation of the interstellar radiation field. Since much of the cloud mass is located at a density close to the HOP threshold density, we would therefore expect the temperature of this gas to be comparable to that of a patch of unshielded ISM with a density $n \sim 30 \text{ cm}^{-3}$. For an interstellar radiation field strength comparable to the solar neighbourhood value, this is $T \sim 100$ K (Wolfire et al. 2003), in good agreement with the low $T$ peak we find in these simulations. In

<footnote n="3">³ The question of whether the conventional definition of $\alpha_{vir}$ is an adequate diagnostic of whether or not clouds are gravitationally bound is an interesting one (see e.g. Ramírez-Galeano et al. 2022), but is out of the scope of this study.</footnote>
<figure><img src="image_15.png" /><figcaption>**Fig. 14.** Distribution of cloud-averaged temperatures.</figcaption></figure>

the SILCC and M51 simulations, we attempt to account approx-
imately for the effects of local attenuation using the TreeCol
algorithm (Clark et al. 2012; Wünsch et al. 2018). The clouds
formed in these simulations therefore feel an effective photoelec-
tric heating rate that is lower than in the Ramses simulations, and
hence have lower characteristic temperatures. Finally, the differ-
ence between the SILCC and M51 simulations stems from the
fact that the two populations of clouds have similar mean den-
sities but different sizes (see Fig. 6). The M51 clouds therefore
have mean surface densities that are a factor of a few larger than
the clouds in the SILCC 0.5 pc simulation, giving them higher
mean extinctions and hence lower photoelectric heating rates.

The temperature distributions shown in Fig. 14 also dif-
fer significantly at the high *T* end. In the SILCC and M51
simulations, very few clouds have average temperatures as
high as 1000 K, while in the LS and Ramses-F20 simulations,
there is a substantial population of clouds with high mean
temperatures, extending up to *T* ~ 10⁴ K. This behaviour is a
consequence of differences in the way that stellar feedback is
modelled in the simulations. In the LS simulations, the effects
of photo-ionisation feedback are included. As HOP identifies
clouds based purely on their density structure, without reference
to their chemical composition, it therefore extracts some clouds
that have been partially or totally photo-ionised, but that have
not yet expanded sufficiently to bring their mean density below
the cloud identification density threshold. The peak in the tem-
perature distribution at *T* ~ 10⁴ K corresponds to fully ionised
clouds, while clouds with temperatures between a few 100 K and
10⁴ K have only been partially ionised, and hence have a mean
temperature that is intermediate between the low temperature of
their neutral portion and the ~10⁴ K temperature of their ionised
portion. In the Ramses-F20 simulations, ionising radiation is
not explicitly followed with a radiative transfer code, but is
accounted for approximately by maintaining the temperature
of the gas in any grid cells containing massive stars at 10⁴ K
(Kretschmer & Teyssier 2020). Therefore, the mean cloud
temperature in this simulation depends on the number of cells
in the cloud that contain massive stars, with the highest mean

temperatures corresponding to clouds in which most or all of the
grid cells are affected by the local heating. Finally, we do not see
this population of warm clouds in the SILCC or M51 simulations
because these simulations do not include any treatment of photo-
ionisation feedback. In principle, SN feedback could produce a
similar effect if we were to happen to identify a cloud in which
a supernova had just exploded, but in practice, clouds in which
SNe have exploded are rapidly disrupted and so the chances of
this occurring are relatively small. That said, we do find exactly
one example of such a cloud in the M51 simulation, with mean
temperature of ~3 × 10⁵ K, visible on the far right in Fig. 14.

The fact that we see such clear differences between the
distributions of cloud-averaged temperatures in the different sim-
ulations indicates that this is a quantity that is very sensitive
to how we chose to model the thermal physics of the ISM and
the impact of stellar feedback in our simulations. However, the
fact that many of the other cloud properties that we have exam-
ined agree well between different simulations demonstrates that
even large differences in cloud temperatures have a relative small
effect on the other properties of the clouds, consistent with the
idea that these are shaped more strongly by turbulence than by
thermal pressure.

**5. Discussion**

We already discussed individual results extensively in their
corresponding sections. Here we address further the issue of
resolution and explore the dependence on the cloud extraction
algorithm.

5.1. *Resolution dependence*

The influence of the resolution of the simulation on the extracted
cloud population is most easily seen in the cloud mass spec-
trum. We already discussed in Sect. 3 that the peak of the
mass spectrum depends on the resolution. However, if we nor-
malise according to the number of clouds per unit area, we see
that all simulations follow the same distribution. The effect of
increasing the resolution is thus purely to be able to resolve an
additional population of smaller clouds; the absolute number of
large clouds remains the same. When the resolution is low, it
is possible that some clouds become blended. In this case one
would expect to find more massive objects. The universality of
the obtained mass spectrum slope for simulations with different
resolutions shows that this is not an issue, at least when structures
are extracted with the algorithm used in this work. Some small
clouds may be counted as being part of a larger cloud, but since
their mass and size are negligible compared to the properties of
the larger cloud, this does not change the statistics.

We also saw that the shape of the low-mass part of the mass
spectrum depends on the simulation code. The LS and Ramses-
F20 simulations show a gradual decline towards smaller and
smaller masses, while the SILCC runs show a sharp cut-off at
the resolution limit. The differences between the grid code sim-
ulations run with RAMSES (LS and Ramses-F20) and FLASH
(SILCC) are related to the different refinement strategies used
in the two codes. While both codes use refinement criteria that
trigger grid refinement in dense regions, the SILCC simulations
have block refinement, such that the entire rectangular block of
8³ cells around the dense region is also refined in order to keep
the grid structure simple. This results in a significant fraction
of cells having high refinement levels, as illustrated in Table 3,
which in turn allows for well-resolved turbulence inside the
majority of the domain. By contrast, the RAMSES runs were
refined in a very localised fashion, and hence although they
<figure><img src="image_15.png" /><figcaption>**Fig. 14.** Distribution of cloud-averaged temperatures.</figcaption></figure>

the SILCC and M51 simulations, we attempt to account approximately for the effects of local attenuation using the TreeCol algorithm (Clark et al. 2012; Wünsch et al. 2018). The clouds formed in these simulations therefore feel an effective photoelectric heating rate that is lower than in the Ramses simulations, and hence have lower characteristic temperatures. Finally, the difference between the SILCC and M51 simulations stems from the fact that the two populations of clouds have similar mean densities but different sizes (see Fig. 6). The M51 clouds therefore have mean surface densities that are a factor of a few larger than the clouds in the SILCC 0.5 pc simulation, giving them higher mean extinctions and hence lower photoelectric heating rates.

The temperature distributions shown in Fig. 14 also differ significantly at the high *T* end. In the SILCC and M51 simulations, very few clouds have average temperatures as high as 1000 K, while in the LS and Ramses-F20 simulations, there is a substantial population of clouds with high mean temperatures, extending up to *T* ~ 10⁴ K. This behaviour is a consequence of differences in the way that stellar feedback is modelled in the simulations. In the LS simulations, the effects of photo-ionisation feedback are included. As HOP identifies clouds based purely on their density structure, without reference to their chemical composition, it therefore extracts some clouds that have been partially or totally photo-ionised, but that have not yet expanded sufficiently to bring their mean density below the cloud identification density threshold. The peak in the temperature distribution at *T* ~ 10⁴ K corresponds to fully ionised clouds, while clouds with temperatures between a few 100 K and 10⁴ K have only been partially ionised, and hence have a mean temperature that is intermediate between the low temperature of their neutral portion and the ~10⁴ K temperature of their ionised portion. In the Ramses-F20 simulations, ionising radiation is not explicitly followed with a radiative transfer code, but is accounted for approximately by maintaining the temperature of the gas in any grid cells containing massive stars at 10⁴ K (Kretschmer & Teyssier 2020). Therefore, the mean cloud temperature in this simulation depends on the number of cells in the cloud that contain massive stars, with the highest mean

temperatures corresponding to clouds in which most or all of the grid cells are affected by the local heating. Finally, we do not see this population of warm clouds in the SILCC or M51 simulations because these simulations do not include any treatment of photo-ionisation feedback. In principle, SN feedback could produce a similar effect if we were to happen to identify a cloud in which a supernova had just exploded, but in practice, clouds in which SNe have exploded are rapidly disrupted and so the chances of this occurring are relatively small. That said, we do find exactly one example of such a cloud in the M51 simulation, with mean temperature of ~3 × 10⁵ K, visible on the far right in Fig. 14.

The fact that we see such clear differences between the distributions of cloud-averaged temperatures in the different simulations indicates that this is a quantity that is very sensitive to how we chose to model the thermal physics of the ISM and the impact of stellar feedback in our simulations. However, the fact that many of the other cloud properties that we have examined agree well between different simulations demonstrates that even large differences in cloud temperatures have a relative small effect on the other properties of the clouds, consistent with the idea that these are shaped more strongly by turbulence than by thermal pressure.

# 5. Discussion

We already discussed individual results extensively in their corresponding sections. Here we address further the issue of resolution and explore the dependence on the cloud extraction algorithm.

## 5.1. Resolution dependence

The influence of the resolution of the simulation on the extracted cloud population is most easily seen in the cloud mass spectrum. We already discussed in Sect. 3 that the peak of the mass spectrum depends on the resolution. However, if we normalise according to the number of clouds per unit area, we see that all simulations follow the same distribution. The effect of increasing the resolution is thus purely to be able to resolve an additional population of smaller clouds; the absolute number of large clouds remains the same. When the resolution is low, it is possible that some clouds become blended. In this case one would expect to find more massive objects. The universality of the obtained mass spectrum slope for simulations with different resolutions shows that this is not an issue, at least when structures are extracted with the algorithm used in this work. Some small clouds may be counted as being part of a larger cloud, but since their mass and size are negligible compared to the properties of the larger cloud, this does not change the statistics.

We also saw that the shape of the low-mass part of the mass spectrum depends on the simulation code. The LS and Ramses-F20 simulations show a gradual decline towards smaller and smaller masses, while the SILCC runs show a sharp cut-off at the resolution limit. The differences between the grid code simulations run with RAMSES (LS and Ramses-F20) and FLASH (SILCC) are related to the different refinement strategies used in the two codes. While both codes use refinement criteria that trigger grid refinement in dense regions, the SILCC simulations have block refinement, such that the entire rectangular block of 8³ cells around the dense region is also refined in order to keep the grid structure simple. This results in a significant fraction of cells having high refinement levels, as illustrated in Table 3, which in turn allows for well-resolved turbulence inside the majority of the domain. By contrast, the RAMSES runs were refined in a very localised fashion, and hence although they