© The Authors 2025

1. Introduction

Observations as well as theoretical models suggest that galaxies are surrounded by a halo of diffuse gas extending out to approximately the virial radius. Often referred to as the circumgalactic medium (CGM), this reservoir is believed to play an important role in galaxy growth and evolution (see Donahue & Voit 2022 for a recent review of the CGM).

The phase-structure of the CGM is complex. While a large fraction of the volume is dominated by a virialized warm-hot component, it also contains denser, colder gas. The high-velocity clouds (HVCs) observed in the Milky Way (e.g., Muller et al. 1963; Wakker & van Woerden 1997) as well as analogs in external galaxies (e.g., Gim et al. 2021; Weng et al. 2022) are common examples. Although such gas clouds were first detected several decades ago, most key questions surrounding their formation and evolution remain open.

In particular, there is no consensus as to how these clouds form. While theory suggests that the condensation of hot halo gas triggered by density perturbations may form dense pockets of cold gas (Binney et al. 2009; Fraternali et al. 2015; Nelson et al. 2020; Faerman et al. 2025), ram pressure stripping of satellite gas as they infall toward the center of their host halo may also be a source of formation (Mayer et al. 2006; Rohr et al. 2023). The clouds may also be a result of galactic-scale outflows driven from galaxies by stellar and active galactic nuclei (AGN) feedback (Voit et al. 2015a; Thompson et al. 2016; Ramesh et al. 2023c). Cool CGM gas could also be inflowing as part of a galactic fountain cycle (Fraternali & Binney 2006) or due to filamentary cosmological accretion, particularly at a high redshift (Kereš & Hernquist 2009; Nelson et al. 2016; Mandelker et al. 2020).

The lifetimes of cool clouds are similarly unknown. Early “cloud crushing” simulations have suggested that the clouds are expected to be short-lived, either as a result of fragmentation (Mellema et al. 2002; Zhang et al. 2017) or complete destruction (Klein et al. 1994; Schneider & Robertson 2017). However, more recent results suggest that their survivability may be enhanced by certain physical processes. As an example, the Kelvin-Helmholtz (KH) instability, driven by velocity shear as clouds pass through their ambient media, creates an interface layer of warm gas. Rapid cooling of this interface can then increase the cold mass of the cloud, effectively improving its lifetime (Scannapieco & Brüggen 2015; Gronke & Oh 2018; Fielding et al. 2020).

Other non-thermal processes also play a role. For instance, magnetic fields may provide additional pressure support to clouds (Sharma et al. 2010; Nelson et al. 2020; Fielding et al. 2023), suppress fluid instabilities through their associated magnetic tension (Ji et al. 2018; Sparre et al. 2020), influence their shape (Kwak et al. 2009) and kinematics (Banda-Barragán et al. 2016; Brüggen et al. 2023), or boost the drag force the clouds experience in the presence of a draped layer (Ramesh et al. 2024), thereby decreasing the time needed for them to begin comoving with their surroundings (McCourt et al. 2015). Similar to magnetic fields, cosmic rays may also contribute pressure support (Butsky et al. 2020; Brüggen & Scannapieco 2020).

Theoretical studies have greatly improved our understanding of cloud growth and survival; however, most are based on idealized (i.e., wind tunnel-type) simulations. Such setups are typically inexpensive from a computational standpoint and thus allow one to realize a higher numerical resolution and/or explore varying physics, but they have key limitations. First, their outcomes are sensitive to the initial conditions, which are ad hoc. The region(s) of the parameter space probed may not correspond to true clouds, and generalization is difficult (Jennings et al. 2022). Second, these setups generally assume clouds to be uniform density spheres embedded in uniform background media¹, and hence they lack the diversity and complexity of a realistic CGM. Lastly, since a preexisting cloud is initialized at the start of the simulation, these setups cannot address the question of origin.

An alternative approach is fully cosmological simulations. These evolve halo gas over cosmic epochs and cosmic scales and thus offer a better chance at self-consistently reproducing a realistic CGM, albeit typically at a coarser numerical resolution with respect to their idealized counterparts. However, recent projects have attempted to take a step forward in narrowing this resolution gap in order to allow for the study of small-scale cold gas structures evolved in the full ΛCDM context (see, e.g., Fig. 1 of Pillepich et al. 2024 and Weinberger & Hernquist 2023; Butsky et al. 2024 (Hummels et al. 2024) for on the fly (postprocessing) sub-grid approaches to model unresolved cold gas). In particular, the TNG50 simulation has been shown to resolve cold clouds in the gaseous halos of massive ellipticals (Nelson et al. 2020) and of Milky Way-like galaxies (Ramesh et al. 2023c; Lucchini et al. 2024). However, the simulation resolution limit dictates that only clouds above a mass threshold of ≳10⁶ M_⊙ can be studied. Given that this is close to the upper limit of known masses of HVCs (Thom et al. 2008) and that clouds down to masses of ≲10⁵ M_⊙ are routinely observed in the Milky Way sky (Wakker 2001; Adams et al. 2013), higher levels of numerical resolution are certainly required to simulate a complete set of analogs of observed CGM clouds.

Cosmological hydrodynamical simulations that specifically aim to increase the numerical resolution in the CGM offer a possible path forward (e.g., Suresh et al. 2019; Hummels et al. 2019; Rey et al. 2024). In this work, we use Project GIBLE (Ramesh & Nelson 2024), a suite of such “CGM-refinement” simulations of (z = 0) Milky Way-like galaxies. The runs we utilize in this work achieve a mass resolution of ∼10³ M_⊙ in the CGM², i.e. roughly 100 times better than TNG50, allowing us to study smaller structures, while also better resolving clouds in the high-mass end. As we elaborate below, the GIBLE runs also contain a substantial number of Monte Carlo tracer particles, which combined with a snapshot cadence as good as ∼30 Myr close to z = 0, makes it possible to quantify the growth and origin of cold clouds to an unprecedented level compared to earlier efforts with cosmological simulations.

The rest of the paper is organized as follows. In Sect. 2, we provide an overview of the simulation suite, the cloud identification algorithm, Monte Carlo tracers, and the use of these tracers to link clouds between snapshots. Results are presented and discussed in Sect. 3 and summarized in Sect. 4.

2. Methods

2.1. Sample and galaxy formation model

In this paper we use Project GIBLE (Ramesh & Nelson 2024), a suite of cosmological magneto-hydrodynamical zoom-in simulations of Milky Way-like galaxies with added preferential resolution refinement in the CGM, defined as the region bounded between 0.15 R_200, c and R_200, c (virial radius). Our sample consists of eight halos, all drawn from the TNG50 Milky Way sample described in detail in Ramesh et al. (2023b) (see also Pillepich et al. 2024). We achieve CGM gas mass resolutions of ∼10³, 10⁴ and 10⁵ M_⊙ in each of these halos, with the galaxy, i.e. the region within the inner boundary of the CGM (0.15 R_200, c), always maintained at ∼8.5 × 10⁵ M_⊙ (TNG100 resolution). Throughout the paper, we refer to these three resolution levels as ‘RF512’, ‘RF64’ and ‘R8’, respectively. All our main results are derived from the RF512 runs, while the other two are used for resolution convergence tests. In all cases, the simulation output is saved at 30 distinct snapshots: five at relatively high redshift (z ∼ [0.50, 0.67, 1.00, 1.50, 2.00]), fifteen between z ∼ 0.02 and z ∼ 0.15, i.e. with an average snapshot spacing of ∼100 Myr, and ten between z ∼ 0.02 and z = 0.0, i.e. with an even finer time cadence of ∼30 Myr. The latter ten snapshots are particularly useful to assess short time-scale phenomena, while the next fifteen make it possible to explore the question of cloud growth and evolution over a significant period of time (∼2 Gyr). We note that this is typically not possible with cloud crushing setups, which are usually evolved for a few hundred Myr at best, since simulating further is not realistic owing to their idealized nature.

Run with the AREPO code (Springel 2010), Project GIBLE uses the IllustrisTNG model (Weinberger et al. 2017; Pillepich et al. 2018) to account for the physical processes that regulate galaxy formation and evolution. This includes the physics of metal line cooling, star formation and subsequent stellar evolution and enrichment, supermassive black hole (SMBH) formation, and feedback from stars (supernovae) and SMBHs. Included in the TNG model is also ideal magnetohydrodynamics (Pakmor et al. 2011, 2014) – a uniform primordial field of 10⁻¹⁴ comoving Gauss seeded in the initial conditions (z ∼ 127) is self-consistently amplified through a combination of structure formation, small scale dynamos and feedback processes (Pakmor et al. 2020; Arámburo-García et al. 2021). Predictions by the TNG model suggest a typical field strength of 1−10 μG within galaxies at z ∼ 0 (Marinacci et al. 2018; Nelson et al. 2018; Ramesh et al. 2023a, b), consistent with recent observational inferences using the Faraday rotation of polarized background sources and fast radio bursts (FRBs) (e.g., Mao et al. 2017; Prochaska et al. 2019; Heesen et al. 2023; Böckmann et al. 2023).

2.2. Cloud and interface identification algorithm

Natural neighbor cells of the unstructured Voronoi tessellation of space provide an ideal spatial definition for clouds. Following Nelson et al. (2020) and Ramesh et al. (2023c), we define and identify clouds as spatially contiguous sets of cold (T ≤ 10^4.5 K) Voronoi cells which contain at least 10 member resolutionelements³. Clouds composed of fewer than 10 cells are considered unresolved, and are also referred to as cold seeds in subsequent sections of the paper. Cells that are gravitationally bound to satellite galaxies, as identified by the substructure identification algorithm SUBFIND (Springel et al. 2001) are excised.

The cloud-background interface is defined as the layer of gas cells immediately surrounding the cloud, and is determined using the connectivity of the Voronoi tessellation. It is thus the first cocoon of gas cells around the cloud. Our cloud definition and identification algorithm allows for clouds and interfaces of arbitrary sizes and shapes (see, e.g., Fig. 1 of Ramesh et al. 2023c).

2.3. Monte Carlo tracers and cloud merger trees

Following an identification of clouds at discrete snapshots, we link clouds with their progenitors and/or descendants with the use of Monte Carlo tracers (Genel et al. 2013; Nelson et al. 2015). In short, these tracers ‘reside’ within resolution elements, and are exchanged between neighboring cells and/or particles with a probability that scales with the mass flux. As Voronoi cells refine and/or de-refine (Springel 2010), tracers are also re-distributed based on the relevant volume ratio. This approach has been shown to be more accurate than velocity field tracers at tracing the underlying mass flow, particularly in regions of strong shocks and/or turbulence (Genel et al. 2013).

At the start of a given RFx simulation, we inject 4x tracers in every gas cell, resulting in an average of 4 tracers in every CGM gas cell⁴ following the start of the CGM refinement⁵. Clouds above the limit of ten member cells thus typically have ≳40 tracers associated with them. For a given cloud in snapshot n, its progenitors in snapshot n − 1 are defined as clouds that contain at least 1 of its tracers, and the progenitor with the largest number of tracers is termed the main progenitor. A similar procedure is followed to identify descendants in snapshot n + 1, and an analogous definition for the main descendant. In effect, we construct the merger trees of clouds (see also Lucchini et al. 2024), similar to galaxy and subhalo merger tree algorithms (e.g., Springel et al. 2005; Behroozi et al. 2013; Rodriguez-Gomez et al. 2016). A key difference is that interaction events between clouds may have multiple descendants (during a fragmentation), while this is typically rare for galaxies (see also Bahé et al. 2019). We note that while cold objects composed of fewer than ten cells (i.e. unresolved clouds) are included in the tree, mergers and/or fragmentations with them were not counted when computing the associated event rates.

2.4. Cloud origins and mass flux rates

In addition to linking clouds across snapshots and creating their merger trees, we also use tracers to assign one of three origin channels to clouds. Starting from a cloud at a given snapshot (always z = 0 in this work), we followed the main progenitor backward to a snapshot n such that none of its tracers were present in any CGM cold gas cell in snapshot n − 1. If the progenitor (in snapshot n) is linked to the galaxy, i.e. the central (ISM) disk, or other clouds within 0.15 R_200, c (in snapshot n − 1), we label the cloud as cold gas arising from the central galaxy. If not, we cross-match all tracers of the progenitor (in snapshot n) with those from snapshot n − 1. The cloud is considered to arise from satellites if any⁶ were part of satellite gas. We note that this could refer to gas stripped due to ram-pressure (e.g., Gunn & Gott 1972; Kenney & Young 1989; Rohr et al. 2023), or to those driven out by feedback within the satellite (e.g., Peluso et al. 2023). The z = 0 clouds in this category are thus direct descendants of cold gas structures that were ejected from satellites. On the other hand, if majority were in the non-cold phase of the CGM, we refer to it as a condensate of the hot halo.

We note that while tracers may originate from the IGM, these never form a majority (≲0.3%). Similarly, tracers linked to the warm-hot phase of the galaxy (e.g., hot outflows) never form a majority, although it is possible that the warm-hot gas in the CGM in snapshot n − 1 is linked to an outflow from an earlier epoch (snapshot ≤ n − 2). Since the time cadence of the simulation output is coarse beyond z ≳ 0.15 (corresponding to a lookback time of ∼2 Gyr), we terminate the assignment of origin channels at this z, and label all still existing clouds as ‘long-lived’. Lastly, at z ∼ 0.15, the progenitors of a small fraction (≲3%) of z = 0 clouds are below the limit of ten member cells to be defined as a cloud. These are not assigned any of the above flags, and are not considered in any of the subsequent analyses.

Finally, the tracers also enable us to compute mass fluxes between clouds and their environments. ‘New’ tracers in a cloud in snapshot n, which were not present in any cold gas cells in snapshot n − 1, are assumed to trace a mass flux into the cloud from the warm-hot phase. An analogous definition is followed to compute the mass out-flux from a cloud to the surrounding gas. Given the limited time cadence between snapshots, this mass flux need not only represent the cycling of material between the cloud and its interface, but also the flow of material into/from larger scales.

3. Results

3.1. The origin of z = 0 clouds

We begin with a visualization⁷ highlighting the different origin channels defined above. Fig. 1 shows the distribution of clouds in the CGM of a random GIBLE halo at z = 0. The image stretches ±R_200, c from edge-to-edge, where R_200, c is the virial radius of the halo, also marked by the white circle. While the projection extends across the same range along the perpendicular direction, clouds within a galactocentric distance of 0.15 R_200, c, i.e. inside the inner boundary of the CGM, are not considered in this work, and are not shown here. The clouds are colored by the different origin channels assigned to them: those that arise from the central galaxy are shown in white. Condensates of the hot CGM phase are shown in blue, while those ejected from satellites are colored in orange. We note that the long-lived clouds in this halo are not shown in this image. This simple visualization suggests that the white and blue channels dominate, while clouds arising through stripping of satellite gas are less abundant.

Fig. 1. Visual impression of the three origin sources of z = 0 cold CGM gas clouds. The white circle marks the virial radius of this GIBLE halo (S105; M200c ∼ 1012.29 M⊙; see Table 2 of Ramesh & Nelson 2024 for more details). Of the many thousand clouds in this CGM, most came to existence in the recent past (≲2 Gyr). A majority of these can be traced back to cold gas that was previously present in the central galaxy (white). On average across the sample of eight halos, these form ∼40–60% of the z = 0 population. A small fraction of clouds (∼3–5%) arise from satellites as they spiral in toward the center of the halo (orange). Colored in blue, most of the remainder of the clouds are formed through the in-situ condensation of hot halo gas. These account for ∼15–50%, depending on cloud mass. This simple visualization suggests that the inner regions are dominated by white clouds, while their blue and orange counterparts are situated farther out in the halo. This radial distribution is typical, as we expand upon in subsequent figures.

In Fig. 2, we quantify this visual inference by showing the fraction of clouds originating from each channel, as a function of z = 0 cloud mass. The cloud fraction arising from the galaxy is shown in black, condensation of the hot halo in red, and ejected from satellites in orange. To increase statistics, we combine all clouds across all eight GIBLE halos.

Fig. 2. Main panel: Fraction of clouds arising from the three origin channels as a function of (z = 0) cloud mass (indicated by thick solid curves). A large fraction (∼40–60%) were previously present in the central galaxy and have been driven out into the CGM in the last ≲2 Gyr (black). The ejection of gas from satellites gives rise to a small portion of the population (∼3–5%; orange). A majority of the remaining clouds are formed through condensation of the hot phase of the CGM in the last ∼2 Gyr. The thin gray curve shows the fraction of clouds that are long-lived, i.e. those whose progenitors already exist at least ∼2 Gyr ago (∼2–10%). Top panel: Assessment of numerical convergence. The fraction of clouds arising through in situ condensation of the hot phase in the RF512 suite (solid curves) is contrasted against the lower resolution RF64 (dashed) and RF8 (dotted) runs. Results between resolution levels approach convergence toward the high mass end, i.e. well above the resolution limit.

As suggested by Fig. 1, clouds from the central galaxy and hot halo condensation dominate. On average, ∼40% of ∼10^4.2 M_⊙ clouds at z = 0 arise as a result of cold gas transition from the galaxy into the CGM in the last ∼2 Gyr. This number increases with cloud mass, to as high as ∼60% for their massive counterparts (∼10^6.6 M_⊙). While ∼67% of Milky Way-like galaxies in TNG50 are in the regime where kinetic kicks by the central super-massive black hole (SMBH) is the major mode of feedback at z ≲ 0.15 (Pillepich et al. 2021; Ramesh et al. 2023b), all eight halos re-simulated in Project GIBLE are in the kinetic mode by z = 0. In our sample, outflowing clouds from the galaxy into CGM over the last ∼2 Gyr are primarily driven by SMBHs (i.e., AGN feedback). We note, however, that some of these clouds may not be ejected due to feedback processes, but rather be dislodged from the extended or outer disk due to rapid rotation and dynamics alone. Furthermore, galactic-scale winds driven by supernovae feedback will also produce cool outflowing gas (Fraternali & Binney 2006; Smith et al. 2024), a phenomenon more readily apparent at higher redshifts and/or lower galaxy masses (see also Suresh et al. 2019).

At the low-mass end (∼10^4.2 M_⊙), clouds are more likely to arise through hot-phase condensation (∼50%). This origin channel portrays the inverse mass trend as above, with the fraction dropping to ∼25% as one transitions to the high mass end (∼10^6.6 M_⊙). While theory suggest that such precipitation can take place if the cooling time of gas is sufficiently short (McCourt et al. 2012; Voit et al. 2017), this condition is typically not satisfied (globally) across the halo (Nelson et al. 2020). However, strong density perturbations can lead to run away thermal instability and cooling (Sharma et al. 2012b; Choudhury et al. 2019; Dutta et al. 2022), the exploration of which we return to shortly.

Descendants of cold gas structures ejected from satellites are relatively rare. The fraction of such clouds is relatively flat at ∼2–3% in the mass range of ∼10^4.2–10^6.2 M_⊙, with a slightly larger fraction of ∼5% for clouds as massive as ∼10^6.6 M_⊙. While past studies have shown that stripping of satellite gas may substantially add to the cold mass budget of the host CGM (e.g., Olano 2008; Steinhauser et al. 2016; Yun et al. 2019), only four of the eight halos in our sample have had gaseous satellites in-fall through them in the last ∼2 Gyr, thus possibly leading to an under-estimate. Indeed, for a larger sample of 132 Milky Way halos in TNG50, clouds preferentially cluster along the past trajectories of satellites (Ramesh et al. 2023c; see Salem et al. 2015; Lucchini et al. 2021; Krishnarao et al. 2022; Grønnow et al. 2024 for analogous inferences by observational studies).

In addition, we note that the total mass of these satellites at first infall was ∼[10^10.7, 10^10.8, 10^11.4, 10^11.45] M_⊙. All four first crossed the halo virial radius between z ∼ 0.15 and 0.5, and owing to the coarse time resolution beyond z ∼ 0.15, it is likely that these numbers are lower limits. In either case, the latter two of these interactions resulted in major mergers (merger mass ratio greater than 1/4, defined as the ratio of stellar masses when the secondary has its maximum mass; Rodriguez-Gomez et al. 2016). It is thus possible that while the CGM of these satellites was stripped as they spiraled toward the central galaxy, only a fraction of their cold-phase ISM was lost (e.g., (Ghosh et al. 2024), see also discussion below). This selective stripping could reduce the fraction of cold gas clouds that arise from theseinteractions.

The thin gray curve shows the fraction of long-lived clouds, i.e. those whose progenitors were already assembled ∼2 Gyr ago. The fraction of such clouds is relatively small (∼2%) at the low mass end (∼10^4.2 M_⊙), and gradually increases to ∼10% for their more massive counterparts (∼10^6.6 M_⊙). Larger clouds are thus more likely to survive over a long time scale. Studies with idealized simulations have reached a similar conclusion, arguing that clouds survive above some critical minimum mass (Gronke et al. 2022). The evolution of long-lived clouds is a main focus of our current work, which we return to and explore in Sect. 3.2.

In the top panel, we assess numerical convergence. To do so, we compare the fraction of clouds with a condensation origin in the RF512 suite (solid line) with the lower resolution RF64 (dashed) and RF8 (dotted) runs. We note that the latter two curves only begin at ∼10⁵ and 10⁶ respectively, since clouds below this threshold are unresolved in these runs. While there is a clear offset between RF512 and RF64 at low masses (∼10⁵ M_⊙), they approach convergence toward the high mass end, i.e. for massive clouds well above the resolution limit. The RF8 run, i.e. TNG50-1 resolution, is far from convergence at these mass scales, but approaches convergence for higher mass scales, although such clouds are rare in the sample (Ramesh & Nelson 2024).

Before proceeding further, we contrast our findings against a recent study that likewise aimed to quantify the origin of cold circumgalactic gas (Decataldo et al. 2024). This study is based on a single Milky Way like halo ‘Eris2k’ at z ∼ 0.3, simulated with the SPH code GASOLINE (Wadsley et al. 2004) at an average baryonic mass resolution of ∼2 × 10⁴ M_⊙, i.e. a factor of ∼10 coarser than our RF512 runs. Their results suggest that condensation out of the hot phase is the primary source of cold gas (∼80%), with contributions by satellites (∼2%) and outflows from the galaxy (∼15%) being less prominent. Apart from the differences in redshift and resolution between Eris2k and GIBLE, we note that the underlying numerical methods and feedback recipes are also different. We speculate that a combination of these factors is the leading cause for differing results. Moreover, while we show results for clouds stacked across all halos in Fig. 2, we mention that the fractions of the two dominant channels typically vary by ∼ ± 15% across the set of halos. A larger sample size with the Eris2k model may change their average result.

In Fig. 3, we explore physical properties of z = 0 clouds, split based on the three origin channels discussed above. Similar to Fig. 2, properties of clouds arising from satellites are shown in black, condensation of the hot halo in red, and ejected from satellites in orange. Clouds across the sample are again stacked.

Fig. 3. PDFs of properties of z = 0 CGM clouds split based on their origin channel. Curves corresponding to those arising as from the central galaxy (within the last ∼2 Gyr) are colored in black, condensates of the hot halo gas in red and those ejected from satellites in orange, as shown by the legend in the top-left panel. Distributions of cloud radius, metallicity, radial velocity and galactocentric distance are shown in the top-left, top-right, lower-left and lower-right panels, respectively. In the right panels, we show individual distributions of the S98 and S105 halos (see Table 2 of Ramesh & Nelson 2024 for more details) through dotted and dashed curves, respectively (vertically offset for better visibility). Distributions are broad for all four properties shown here, with PDFs of clouds ejected from satellites typically different as compared to the other two cases.

In the top-left panel, we show distributions of cloud sizes, measured as the radius of the volume equivalent sphere. While both the black and red curves peak around ∼650 pc, the former is more broad, extending to sizes as large as ≳2 kpc. This is a result of clouds in the latter category being more spherical, as opposed to elongated, with smaller major-to-minor axis ratios on average (not shown explicitly, although visible in Fig. 1). Clouds ejected from satellites are typically smaller, peaking at ∼450 pc. We interpret that this is due to densities being larger in the ISM of satellites than in the cold-phase CGM, leading to smaller clouds for a given mass. The distribution, however, extends to as large as ≳2 kpc, similar to the red curve, suggesting that this may be a signature of cold CGM gas of satellites being stripped, as we explore below.

The top-right panel shows PDFs of cloud metallicity. The red and black curves are largely similar, peaking at ∼25% Z_⊙. The distributions are, however, broad, with clouds as enriched as ∼2 Z_⊙, and as pristine in their metal content as ≲0.1 Z_⊙. The likeness between the two curves suggests a substantial level of mixing between halo gas and metal enriched outflows driven outward from the galaxy by feedback processes (Fielding et al. 2017; Nelson et al. 2019, see also Appendix A). A similar variance of metallicity is observed for high velocity clouds (HVCs) in the Milky Way halo, with values ranging between ≲0.3 Z_⊙ (Richter et al. 2001; Collins et al. 2003) to ≳1.5 Z_⊙ (Zech et al. 2008; Yao et al. 2011).

The orange curve, however, portrays a different behavior: it is bimodal with a primary peak at ∼0.45 Z_⊙ and a smaller secondary peak at a lower value of ∼0.15 Z_⊙. This is a result of superposing the distributions from the different halos with their distinct satellite populations. To dig deeper, we also show the individual distributions from two particular halos: the dashed curve corresponds to the halo with a ∼10^11.40 M_⊙ satellite discussed above, while the dotted line shows the halo with the ∼10^11.45 M_⊙ satellite, both vertically offset downward for visibility. Clouds from these two halos dominate the orange curve, so considering these alone is sufficient to explain its structure. The dashed line peaks at higher metallicity, suggesting that these clouds are signatures of ISM gas being stripped. On the contrary, the dotted curve peaks at lower metallicity, and we interpret that these clouds may instead primarily arise from the cold CGM of this satellite. At the same time, a small feature at higher Z hints that a portion of the ISM is also stripped, as we discuss more below. The low-Z tail of clouds arising from satellites is also more prominent compared to the black and red curves, suggesting that much of the CGM gas of these satellites is metal poor in comparison to the CGM of the central galaxy (Pasquali et al. 2012). We note that the gas phase metallicites of both satellites and centrals, particularly that of the CGM, is sensitive to feedback episodes (Péroux et al. 2020), which further adds to the diversity in distributions.

In the lower-left panel, we show distributions of radial velocities. The black curve peaks at ∼ − 40 km/s, i.e. most clouds that moved outward from the galaxy in the last ∼2 Gyr are now inflowing, tracing the second phase of the halo fountain cycle (Péroux & Howk 2020). The distribution, however, extends to large positive radial velocities, with clouds moving from the galaxy into the CGM at speeds as large as ∼200 km/s, corresponding to the first phase of a halo-scale fountain flow (Oppenheimer & Davé 2008).

The red curve peaks at a comparatively larger negative radial velocity of ∼ − 60 km/s, i.e. a larger fraction of clouds are inflowing as compared to the black curve, as they precipitate out of the hot halo and begin infalling under the influence of gravity (Voit et al. 2015b). Although outflowing clouds in this category are less common as compared to the black curve, a small fraction of clouds are moving away from the galaxy, with outward velocities as large as ∼200 km/s. These are clouds that likely precipitated out of warm-hot outflows from the galaxy, or fromgalactic cold gas that dissolved in the surrounding wind and thereafter condensed (e.g., Zhang et al. 2017, see also Fig. 5). Clouds stripped from satellites are also preferentially inflowing, with the distribution peaking at ∼ − 60 km/s. Some clouds are, however, outflowing, either as a result of their rotational vectors aligning with the direction of ram pressure (Sparre et al. 2024), or possibly due to feedback processes within the satellite driving gas along the stripped tail (Peluso et al. 2023).

Lastly, in the lower-right panel, we focus on the distributions of galactocentric distance. The black curve is relatively flat between the inner radius of the CGM (0.15 R_200, c) and ∼0.3 R_200, c, and begins to decline toward larger distances, i.e. these clouds are centrally concentrated. The red curve peaks at ∼0.4 R_200, c: these clouds dominate the inner halo. However, the distribution is broad, with clouds as far as ∼0.9 R_200, c, but also close to the inner boundary of the CGM. Similar to the top-right panel, clouds ejected from satellites show a bimodal behavior, with peaks at ∼0.4 and ∼0.8 R_200, c, once again as a result of superposing distributions across halos. The dashed curve (the ∼10^11.4 M_⊙ satellite) peaks at relatively large galactocentric distances, suggesting that this satellite was stripped soon after infall. The dotted curve (the ∼10^11.45 M_⊙ satellite) also shows a primary peak at large galactocentric distances, and this corresponds to the low-metallicty CGM gas. A second less prominent peak in the inner region of the halo is seen, signifying the onset of stripping of ISM gas as ambient densities grow larger (Ayromlou et al. 2019). Furthermore, these two satellites with roughly the same total mass at infall seem to loose their gas in varying ways, suggesting that orbital diversity and large-scale environment play important roles in such hydrodynamic interactions (e.g., Roediger & Brüggen 2007; Sales et al. 2015).

3.1.1. Precipitation of the hot phase

As discussed above, cold clouds of gas may be able to condense out of the hot phase of the CGM (e.g., Fraternali & Binney 2006). Indeed, this is the origin channel for ∼25–50% of clouds in our sample (Fig. 2). In this sub-section, we take a deeper look into their growth.

We begin with evolutionary tracks of mass growth in Fig. 4. We group clouds into four bins: those that first came into existence as cold CGM gas cells⁸ ∼250, 750, 1300 and 1850 Myr ago, shown in black, red, orange and blue, respectively. These numbers are thus also upper limits on the formation times of their progenitors. Solid curves show the median, while the shaded bands correspond to the 16th-84th percentile regions. As before, we stack clouds across all halos. The gray band shows the mass range within which gas cells are maintained by AREPO through refinement and de-refinement, i.e. corresponding to the mass of one gas cell, while the red dashed line shows ten times the average mass resolution.

Fig. 4. Evolutionary tracks of clouds that condense out of the hot halo. Solid curves correspond to the median while the shaded bands show the 16th-84th percentile regions. Colors correspond to different epochs where the cloud first came into existence. In all cases, the origin of the cloud can be traced back to a ‘cold seed’, typically comprised of only a few gas cells. Following that, clouds accrete matter to grow to a median mass of ∼105.2 M⊙ by z = 0. The percentile regions, however, are broad, implying that evolutionary tracks of individual clouds can be diverse.

Going backward in time starting from z = 0, all evolutionary tracks lead back to a small cold gas object, typically only comprised of a few cells. The fact that it is not always one cell at the start is likely due to limited time resolution between snapshots⁹. In all cases, starting from this cold over-density, curves grow steadily to reach a median mass of ∼10^5.2 M_⊙ by z = 0. However, the percentile regions are broad, implying that evolutionary tracks may be diverse. Furthermore, while the median curves are smooth and roughly monotonous, the growth of individual objects is typically not (see upper panels of Fig. 8 for analogous growth tracks).

Precipitation of clouds out of the hot phase thus begins with the formation of a ‘cold seed’. Such an idea was discussed earlier by Nelson et al. (2020), where it was suggested that thermal instability and runaway cooling could thereafter lead to accretion of cold material (e.g., Field 1965; Sharma et al. 2012a; Dutta et al. 2022), ultimately resulting in the growth of the seed into a bigger cloud. However, it is presently unclear as to how such a central cold over-density is formed in the first place.

We aim to address this question in Fig. 5. We focus on seeds, i.e. contiguous cold objects composed of fewer than ten cells (Sect. 2.2) that first appeared at z = 0, and we note that the result is similar for seeds that appeared at higher z. To boost statistics, we stack tracers across all seeds, and also across all halos. The y-axis shows the fractional contribution to the mass of these seeds as a function of lookback time, all the way back to z ∼ 2 (t_lookback ∼ 10.5 Gyr). We consider five distinct sources: (a) gas from the central galaxy (all phases of gas included) shown in black, (b) the cold CGM, i.e. clouds that existed in the past but subsequently dissolved into the hot phase, in red, (c) the warm-hot CGM in orange, (d) gas from satellites i.e. other galaxies in general¹⁰ (all phases included) in blue, and (e) the IGM¹¹ in gray.

Fig. 5. Dissecting the different sources of gas that seed the condensation of cold clouds out of the hot halo at z ∼ 0. On short timescales (≲500 Myr), the dissolution of past clouds into the hot phase, and the subsequent mixing, is the primary source, and it is also potentially stirred by outflows from the central galaxy into the CGM. Transitioning out into previous cosmic epochs as early as z ∼ 2 (tlookback ∼ 10.5 Gyr), the mass giving rise to these seeds was predominantly smoothly accreted from the IGM into the halo (∼85%), while clumpy accretion through satellites contributes the remaining ∼15%.

At a lookback time of ∼30 Myr, i.e. just before the seed was formed, ∼95% of the mass was in the warm-hot phase of the CGM. The remainder (∼5%) arose from the central galaxy, which by definition corresponds to warm-hot gas being ejected outward. While the contribution of other cold CGM clouds is 0% at this first snapshot, it rises to as large ∼25% at earlier times. The dissolution of clouds into the background halo, and the eventual mixing between phases, thus plays an important role in seeding the next generation of clouds (see also Nelson et al. 2020). In addition, the mass contribution due to outflows¹² from the galaxy rises to as high as ∼10%. Turbulence driven by the interaction of these outflows with the CGM gas may further contribute to precipitation (Mohapatra & Sharma 2019).

Zooming out to previous cosmic epochs as early as z ∼ 2, the mass contribution of IGM gas rises to ∼85%. A similar increase in the contribution of gas from other galaxies is also seen, to ∼15% at t_lookback ≳ 2 Gyr. Almost all of the gas seeding condensation at z = 0 (and in general low-z) is thus extragalactic at z ∼ 2, after which it accretes, either smoothly or in a clumpy fashion (Kereš et al. 2005; Nelson et al. 2013). At least some of this gas ultimately reaches the galaxy (Appendix A), and is recycled back into the CGM, possibly multiple times (Anglés-Alcázar et al. 2017; Suresh et al. 2019). This halo fountain, as well as the motion of satellite galaxies through the CGM, likely produces additional density perturbations in the halo (Gauthier 2013; O’Sullivan et al. 2018), giving rise to clouds that thereafter dissolve in the hot phase, eventually leading to a new population of clouds, and so on. Overall, the short time-scale (≲500 Myr) dissolution of pre-existing cold structures into the hot phase primarily drives the formation of cold seeds. Simultaneously, the source of their baryonic mass is ultimately of cosmological origin (see also Appendix A).

3.1.2. Predicting the origin of Milky Way clouds

We conclude the origin-segment of this paper with a simple prediction for the origin of high velocity clouds (HVCs), and in general cold clouds, in the Milky Way halo. Following Ramesh et al. (2023c), we place a hypothetical observer at a random point along the galactic plane, at a galactocentric distance of 8.34 kpc. Further, the observer is considered to be moving along a perfectly circular orbit around the galactic center, at a velocity of 240 km/s. This observer is consistent with the solar location and motion in the Milky Way (Reid et al. 2014). We repeat this randomization 100 times for each halo, resulting in a total of 800 random realizations for the eight GIBLE halos stacked together.

In Fig. 6, we show the fraction of cloud origin as a function of the line-of-sight velocity to this hypothetical observer. Clouds arising from the galaxy are shown in black, condensation of the hot halo in red, and ejected or stripped from satellites in orange. Solid lines show the median, and the shaded bands the 16th-84th percentile regions across the random realizations. As the ultra-refined region in our simulation is bounded by 0.15 R_200, c (∼ 30–40 kpc), only clouds beyond this distance are considered. Known distances to HVCs are typically ∼4–15 kpc (Wakker et al. 2008; Peek et al. 2016), although these are estimated using absorption-line spectra along the line-of-sight of background stars, and thus the lack of known clouds at larger distance estimates is likely a direct result of the dearth of distant halo stars, i.e. an observational bias.

Fig. 6. Probabilities of the three origin channels as a function of line-of-sight velocities of clouds in the frame of the local standard of rest. The vertical dashed lines delineate commonly used definitions to classify clouds in the Milky Way halo. Outflowing clouds predominantly trace cold gas arising from the central in the recent past (≲2 Gyr), while their inflowing counterparts typically condense out of the hot phase. Clouds ejected from satellites are relatively rare (≲5%), but roughly equally likely at all (line-of-sight) velocities.

The different vertical dashed lines and labels correspond to commonly used definitions to classify clouds based on their kinematics. The low velocity clouds (LVCs) are those with (absolute) line-of-sight velocities in the range ≲40 km/s, the intermediate velocity clouds (IVCs) are in the range of 40–90 km/s, and HVCs are ≳90 km/s (e.g., Wakker 2001).

Figure 6 suggests that LVCs and IVCs are roughly equally likely to arise as a result of cold gas moving outward from the galaxy, as well as through condensation of the hot phase. At larger velocities, the two median curves begin to separate out. Specifically, at line-of-sight velocities ≲ − 100 km/s, the red curve overtakes the black, i.e. inflowing HVCs predominantly trace clouds that condense out of the hot halo. On the contrary, at velocities ≳100 km/s, the inverse trend is observed, i.e. outflowing HVCs are more likely descendants of cold gas structures that were ejected from the galaxy. We note, however, that the two percentile regions overlap out to velocities as large as ∼ ± 200 km/s, particularly in the regime of positive velocities. Unambiguous identification of the origin of clouds in the Milky Way halo is thus not possible with kinematics alone.

Clouds ejected from satellites are relatively rare (≲5%), as previously discussed, but more or less equally likely at all line-of-sight velocities. We note, however, that these fractions are for clouds stacked across all masses and distances, and numbers may change if the selection of clouds is restricted based on a given metric. For instance, the fraction of clouds arising from the cold gas supply of the galaxy is a bit larger in the inner halo, as discussed earlier in the lower-right panel of Fig. 3.

Considering additional physical properties beyond kinematics could help break degeneracies and more clearly identify cloud origins. Metallicity and metal abundance ratios are one example. Mock absorption spectra through clouds and their interfaces, as a function of origin, will establish mappings between observable and unobservable cloud properties (e.g., Eisert et al. 2024; Weng et al. 2024) in future work (Guo et al. in prep.).

3.2. Evolution of long-lived clouds

Throughout the rest of this paper, we focus exclusively on the long-lived clouds in our sample. In addition, to limit the stochastic noise associated with Monte Carlo tracers while estimating quantities such as mass flux, we only consider clouds above a mass threshold of 10^4.8 M_⊙. In effect, we restrict the analysis to relatively large clouds, typically resolved by ≳50 cells and containing ≳200 tracers. We examine their evolution over the last ∼2 Gyr, a span of time accessible due to the fully cosmological nature of our simulations.

In Fig. 7, we begin with a visualization of cold clouds in a GIBLE halo at z = 0. The background image shows a ray-traced volume rendering of CGM gas clouds, with the virial radius (R_200c ∼ 270 kpc) indicated by the white circle (using a visualization method similar to Nelson et al. 2016). The various colored surfaces highlight isothermal contours sampled by five narrow Gaussian transfer functions at log₁₀(T[K]) = [4.05, 4.15, 4.25, 4.35, 4.45] with a standard deviation of σ = 0.006 in log₁₀(K), corresponding to different layers of cold gas within clouds. A temperature gradient into the center of clouds is visible in most cases, with the innermost region typically being the coldest (bluish colors), as opposed to outer layers of clouds that are relatively warmer (reddish colors). Although not seen here, this temperature gradient also extends outward, i.e. these clouds are surrounded by a warm-interface layer that is sandwiched between the cloud and the hot background, as previously shown by Nelson et al. (2020) and Ramesh et al. (2023c) for clouds in TNG50. This structure is similar to idealized simulations run at much higher numerical resolution (e.g., Gronke & Oh 2018; Fielding et al. 2020), and is believed to play a role in cloud growth and survival by facilitating rapid condensation onto the cloud through radiative cooling (Dutta et al. 2022).

Fig. 7. Visualization of cold clouds in a GIBLE halo (S98; see Table 2 of Ramesh & Nelson 2024 for more details) at z = 0 using a ray-traced volume rendering sampled by five narrow Gaussian transfer functions in temperature space, as shown by the color bar (background image). The white circle denotes the outer boundary of the CGM refinement region, i.e. the virial radius of the halo (∼270 kpc). In the foreground, we show the (past) trajectory of a ∼106.1 M⊙ cloud identified at z = 0, with the various translucent stars highlighting the positions of the main progenitor of this cloud at 25 distinct snapshots back to z ∼ 0.15. The insets at the corners of the image show a volume rendering of (the main progenitor of) this cloud at four distinct redshifts, going back in time in clockwise direction starting from the bottom left, with the cloud oriented such that its mean velocity vector points to the right. In the ∼2 Gyr between z ∼ 0.15 and 0, this cloud clearly evolves in terms of shape, size (characterized by the varying scale bars in the insets) and mass (shown by the different labels; in units of M⊙). Also visible in the top-right inset is an instance of an ongoing merger between clouds, a phenomenon explored and quantified later in this paper.

In the foreground, we show the (past) trajectory of a ∼10^6.1 M_⊙ cloud identified at z = 0, with the various translucent stars highlighting the positions of the main progenitor of this cloud at 25 distinct snapshots up to z ∼ 0.15, using the merger tree constructed as detailed in Sect. 2.3. About ∼2 Gyr ago (i.e. z ∼ 0.15), this cloud is situated close to the virial radius of the halo, and soon begins infalling toward the center. During this time, the cloud changes significantly in terms of shape, size and mass as it interacts with its surroundings as well as other clouds. This is highlighted in the insets at the corners of the image, which show volume renderings of gas that belongs to the cloud at four distinct snapshots, going back in time (clockwise from the bottom left), with the cloud oriented such that its mean velocity vector points to the right in each case. The non-monotonic evolution of the mass of the cloud is shown by the labels in the corresponding insets, in units of M_⊙. Also visible in the top-right inset is an example of an ongoing merger between two clouds, a phenomenon that is explored and quantified in subsequent parts of this paper. We also highlight the elongated, filamentary, and in general non-spherical geometry of this cloud throughout its evolution, a common feature of cool clouds in GIBLE halos.

Generalizing beyond this single example, Fig. 8 shows the average evolutionary history of cloud masses over the last ∼2 Gyr for a population of clouds. As before, we stack clouds across all eight GIBLE halos. In the main panel, we segregate clouds into four bins of width 0.4 dex centered on M_cl ∼ [10^5.0, 10^5.5, 10^6.0, 10^6.5] M_⊙, shown in black, red, orange and blue, respectively. Solid curves show the medians, while the shaded bands correspond to the 16th-84th percentile regions. While the four different median curves portray different behaviors up to a lookback time of around a Gyr, they begin to converge at earlier times, i.e. on average, the progenitors of clouds from all four mass bins had roughly the same mass at some point in the past. Specifically, on average, all progenitors were roughly ∼10^6.25 M_⊙ in mass about 2 Gyr ago. As we elaborate further in subsequent sections, this mass threshold is important with respect to merger and fragmentation events of clouds and thus to their growth (see also discussion below for resolution dependence). In all cases, the percentile regions are relatively broad: starting from a spread of ∼0.3 dex at z = 0, the percentile regions are as large as ≳2.0 dex at earlier times, suggesting that evolutionary tracks of individual clouds are diverse.

Fig. 8. Main panel: Evolution of masses of z = 0 clouds whose main progenitors already existed ∼2 Gyr ago (z ∼ 0.15). Solid curves correspond to medians, shaded bands to the 16th-84th percentile regions, and the colors to four different mass bins. On average, progenitors of clouds from all four bins were massive (∼106.25 M⊙) 2 Gyr ago, albeit with relatively broad percentile regions in all cases. In the top panel, we show individual tracks for a subset of clouds split into three categories, chosen randomly: those whose masses have overall increased (left), remained overall constant (center) and overall decreased (right) in this time frame. For exactly one curve in each case, we show all instances of mergers and fragmentations as identified in our cloud merger trees: the former are shown in blue dots, the latter in red, and overlapping dots of both colors correspond to instances where both these events take place between neighboring snapshots. As expected, sudden jumps in curves are associated with these events. Finally, in the inset of the main panel, we assess numerical convergence by contrasting results of RF512 (solid) with the lower resolution RF64 (dashed) and RF8 (dotted) runs for the case of (z = 0) ∼106 and 106.5 M⊙ clouds. While the three runs are qualitatively similar, curves converge at a higher mass for the lower resolution runs (see text for discussion). We note that the axes labels of the inset and the top panels are identical to that of the main panel.

In the inset of the main panel, we assess the impact of numerical resolution. Results from the RF512 suite (solid curves) are contrasted against the lower resolution RF64 (dashed) and RF8 (dotted) runs. Curves in the inset are limited to clouds of (z = 0) masses ∼10⁶ and 10^6.5 M_⊙. Both mass bins converge well at lower resolution. However, the mass at which curves converge is higher for lower resolution runs. Results are thus qualitatively similar across resolution levels, albeit quantitatively different, likely due to the absence of ‘small’ clouds at coarser resolution, thus effectively dragging the median to larger values.

To explore the possible diversity in evolutionary tracks of clouds further, in the top panel, we show the individual histories of a subset of (z = 0) 10^6.0 M_⊙ clouds, split into three categories: those whose masses are overall increasing (left), overall constant (center) and overall decreasing (right) over the last 2 Gyr. We define the first (third) category to be clouds whose progenitors always have a mass lower (greater) than their z = 0 value. For clouds in the second category, we consider those whose progenitors have a mass within ± 0.75 dex of the z = 0 value. The red dashed line toward the bottom of each of these panels signifies ten times the average baryonic resolution, i.e. the lower limit for the mass of clouds. Thus by definition, individual tracks cannot extend below this red line.

Evolutionary tracks are clearly unique in each case: clouds evolve through distinct pathways over time, with impulsive responses to merger events similar to galaxies (e.g., Rodriguez-Gomez et al. 2016). While smooth mass changes occur during certain intervals of time, sudden jumps in mass are also visible. To understand this better, for one curve in each panel, we show instances of mergers and fragmentations as identified during the construction of our cloud merger trees. The former are shown in blue dots, the latter in red, and overlapping dots of both colors correspond to instances where both these events take place between neighboring snapshots. As expected intuitively, sudden jumps in curves are associated with one of these two events, or in some cases, both. The mass growth or loss of clouds thus cannot be studied in isolation, if it is dominated (or significantly affected) by interactions at the cloud population level, as we explore further below. We note, however, that due to a large scatter in clouds, the trends in all three panels are overall rather weak.

3.2.1. Merger and fragmentation events

Episodes of mergers and fragmentations clearly seem to play an important role in how clouds evolve over time. Idealized cloud-crushing simulations can assess how a cloud fragments into smaller cloudlets over time (Klein et al. 1994), and how these fragments may thereafter coagulate to form a larger mass (Banda-Barragán et al. 2020, 2021). However, mergers between halo-scale distinct clouds are not captured in idealized setups. In this sub-section, we quantify the merger and fragmentation rates expected when a cosmological population of clouds interact with each other, and assess the physical conditions that drive these events.

In Fig. 9, we begin with a visualization of a merger (top row) and fragmentation (bottom row) event. In both cases, neighboring panels are separated in time by ∼30 Myr. The scale bars of all panels in a given row are the same, as shown in the bottom-right corner of each row. All panels are oriented with the mean velocity vector of the cloud (system) pointed toward the right. Panels with two colors show times when two distinct clouds exist, while panels with only one color contain a single cloud.

Fig. 9. Visualization of two clouds merging into one (top row) and of a cloud fragmenting into two smaller objects (bottom row). In both rows, neighboring panels are separated in time by ∼30 Myr. The scale bars of all panels of a given row are the same, shown in the lower right of each. All panels are oriented such that the mean velocity vector of the cloud (system) is pointed toward the right. Panels with two colors show times when two distinct clouds are present, while panels with only one color contain a single cloud.

In the top row, we see a case of two clouds moving toward each other, merging into a larger mass, and comoving as a single object thereafter. On the contrary, the bottom panel highlights an example of a single cloud spontaneously breaking into two smaller objects, which drift apart from each other soon after. We note that for simplicity and clarity, we here pick examples where the cloud system is relatively isolated, i.e. these clouds do not interact with any other clouds in this time frame. However, in general, between two snapshots, a cloud may fragment into multiple objects, or multiple clouds may merge together, or both a fragmentation and merger could take place (the latter already seen in Fig. 8). Such cases are probably due to limited time resolution between snapshots, and it is likely that cloud interactions take place in a series of binary steps (Fig. B.1), much like the case shown in Fig. 9.

In the main panel of Fig. 10, we quantify the rate at which merger (black) and fragmentation (orange) events occur. Curves show the mean values, while shaded bands correspond to the ± 1σ variation. As these events are short time-scale phenomena and their rates are sensitive to time cadence between snapshots (Fig. B.2), we here only derive rates between z = 0 and ∼0.02, i.e. where we have a time resolution of ∼30 Myr, and compute mean rates averaged over these snapshots. However, we find no evidence for a trend of events rates with redshift, and thus results from this analysis can likely also be extended out to z ∼ 0.15.

Fig. 10. Main panel: Merger and fragmentation rates as a function of (z = 0) cloud mass. The former are shown in black, and the latter are in orange. A clear trend of event rates with cloud mass is visible, with lower rates for smaller clouds. Top panel: Difference in rates between merger and fragmentation events. On average, massive clouds merge more often than they fragment, while the opposite is true for their less massive counterparts.

Both merger and fragmentation rates show a clear trend with cloud mass. For instance, a merger rate of ∼20 Gyr⁻¹ for clouds of mass 10^5.0 M_⊙ increases gradually to ∼54 Gyr⁻¹ for their more massive counterparts with ∼10^6.5 M_⊙. A similar trend is observed for fragmentation events: a rate of ∼26 Gyr⁻¹ for low mass clouds (10^5.0 M_⊙) increases to ∼46 Gyr⁻¹ for clouds of mass 10^6.5 M_⊙. However, as mentioned above, event rates are sensitive to time resolution: these numbers are likely an underestimate, and values may be larger by a factor of ∼4 at an improved snapshot cadence of ≲6.0 Myr (Fig. B.2). Regardless, the qualitative mass trends shown here remain unchanged at even higher time resolution.

In the top panel, we show the difference between merger and fragmentation rates (Δ_mer, frag; y-axis) as a function of z = 0 cloud mass (x-axis). The dashed horizontal line at a value of zero demarcates the regime where mergers are more frequent than fragmentations (above the line) versus the opposite. A net rate of ∼ − 7 Gyr⁻¹ for clouds of mass 10^5.0 M_⊙ (i.e. 7 fragmentations more than mergers every Gyr) increases to ∼8 Gyr⁻¹ for more massive clouds at 10^6.5 M_⊙ (i.e. 8 mergers more than fragmentations every Gyr), with the transition between negative and positive values occurring roughly at a mass of ∼10^6.2 M_⊙. We note that this is close to the mass above which clouds in the RF512 run experience a net growth in mass over the last ∼2 Gyr (Fig. 8), suggesting that merger and fragmentation events indeed play an important role in cloud mass growth over time¹³. In general, massive clouds (well above the resolution limit) experience a greater number of mergers than fragmentations, and vice versa for their smaller counterparts. The bigger clouds are thus more likely to survive for longer periods of time. We note that similar conclusions have been reached by studies with idealized setups, postulating that clouds above a critical threshold mass¹⁴ have greater chances of increased survivability (e.g., Gronke et al. 2022).

In Fig. 11, we correlate these event rates with two important physical quantities: local pressure, and a local measure of environmental cloud density i.e. neighbor number. Similar to Fig. 10, we show the rates of mergers (black) and fragmentations (orange) in the main panels, while the difference between the two (Δ_mer, frag) is shown in top panels. Solid curves show the average values, while the shaded regions correspond to the ±σ variation at fixed x-axis value. Unlike Fig. 10 where the event rates were averaged over the last ∼300 Myr (i.e. the last 10 snapshots), the rates here correspond to those between neighboring snapshots, and properties are computed at the snapshot preceding the event. This enabled us to answer the question of whether a given property of the cloud at the present time influences the chances of a merger and/or fragmentation in the near future. While we show results stacked over clouds of all masses, we mention that the following trends are qualitatively similar at roughly fixed cloud mass, although typically steeper for more massive clouds.

Fig. 11. Main panels: Merger and fragmentation rates as a function of pressure contrast between the cloud and the interface (left) and Δ10 (right). Top panels: Difference between merger and fragmentation rates. A strong trend with pressure contrast can be seen, with under-pressurized clouds being more likely to undergo a fragmentation event as opposed to a merger. The trend with Δ10 is weaker – clouds with a very small (large) number of neighbors are likely to fragment (merge), while the net rate is largely flat for intermediate values of Δ10.

We began by exploring the rate of interactions as a function of the total (i.e. thermal plus magnetic) pressure contrast of the cloud with respect to the interface layer Δ_cl, in. Positive (negative) values, i.e. those that fall to the right (left) of the vertical line marked at a value of zero, correspond to clouds that are over (under)-pressurized with respect to their local interfaces.

An average merger rate of ∼8 Gyr⁻¹ for clouds that are heavily under-pressurized increases to ∼48 Gyr⁻¹ for those that are over-pressurized with respect to their interface. While the fragmentation rate also portrays the same qualitative trend, the rise is less steep: a rate of ∼19 Gyr⁻¹ increases to ∼42 Gyr⁻¹ between the two extremities of this plot. The difference between the two, shown in the top panel, shows a clear trend: the net difference between merger and fragmentation events is negative (positive) when clouds are under (over)-pressurized, i.e, these clouds are more (less) likely to fragment that to merge, with the cross-over between the two occurring roughly at a pressure contrast of zero. The pressure contrast is thus an important quantity that has an impact of the occurrence of such interaction events and thus on the overall growth of clouds (Li et al. 2020). Earlier work with TNG50 has shown that clouds are typically severely thermally under-pressurized with respect to their ambient media, and the addition of magnetic pressure brings them closer to pressure balance with their surroundings (Ramesh et al. 2023c), highlighting the importance of magnetic fields.

In the right panel, we show trends with respect to Δ₁₀, a proxy for the clustering of clouds around other clouds introduced earlier in Ramesh et al. (2023c), where it was shown that clouds are not uniformly distributed throughout the halo, but are rather typically clustered around their analogs. It is defined as the number of neighboring clouds that lie within a sphere of radius 10kpc centered on a given cloud. A value of zero thus implies that there are no neighboring clouds within this sphere, a value of one corresponds to one neighbor, and so on. We note that, in this approach, we treat clouds as point-like objects located at their center of mass, which is a decent approximation since extended clouds are not as frequent as their smaller counterparts (Nelson et al. 2020; Ramesh et al. 2023c; Ramesh & Nelson 2024). While we only show trends with respect to Δ₁₀ here, we mention that results are qualitatively similar for Δ₂₀ and Δ₃₀.

While both the merger and fragmentation rates increase toward larger values of Δ₁₀, the trend is weak. An average merger rate of ∼26 Gyr⁻¹ at Δ₁₀ ∼ 3 increases gradually to ∼40 Gyr⁻¹ at Δ₁₀ ∼ 40. The rate of fragmentation events is rather similar, with slightly elevated (decreased) rates at low (high) Δ₁₀, as seen in the top panel: the difference between the two rates is negligible for most values of Δ₁₀ (i.e. 10 ≲ Δ₁₀ ≲ 35). However, at lower values, the net rate is negative, signifying that fragmentation events are more common than mergers when clouds are surrounded by very few counterparts, as intuitively expected. On the other extreme, when clouds are strongly clustered and have a larger number of neighbors, the merger rate overtakes that of fragmentations. In this case clouds have a net positive mass growth. The spatial clustering of clouds may thus play more than one important role in cloud evolution: regulating the net merger rate and thus mass growth while also increasing cloud lifetimes through the process of drafting, as seen in wind tunnel setups (Williams & Shelton 2022).

3.2.2. Mass flux between clouds and their surroundings

Results from various idealized simulations suggest that the mass flux into and out of clouds, and their net balance, is an important factor that governs cloud growth and survivability. Fluid instabilities such as the Kelvin-Helmholtz instability trigger mixing of the cloud with the ambient gas, thereby driving an outflux from the cloud (Agertz et al. 2007). This mixing gives rise to an intermediate layer of warm gas which, if radiative cooling is substantial, may condense down onto the cloud resulting in an influx of cold gas (e.g., Melioli et al. 2005; Marinacci et al. 2010; Armillotta et al. 2016; Gronke & Oh 2018). In this sub-section, we quantify the in- and out-flux rates of CGM clouds.

In Fig. 12, we show the flux rates as a function of cloud mass at z = 0 for the same sample of long-lived clouds as above. Similar to merger rates, we find that flux rates are sensitive to the time cadence between snapshots and hence restrict the computation of fluxes to the final ten snapshots (0.0 ≲ z ≲ 0.02) which have a time resolution of ∼30 Myr. As before, fluxes are averaged over these 10 snapshots (∼300 Myr). However, we note that we find no strong trend of flux rates with redshift, and the results of this analysis can likely be extended out to z ∼ 0.15. Solid curves show the median, while the shaded bands correspond to the 16th–84th percentile regions.

Fig. 12. Main panel: Mass in- and out-flux as a function of (z = 0) cloud mass. The former are shown in black, while the latter are in orange. Both rates correlate with cloud mass and are higher for more massive clouds. In the inset, we show the difference between the two rates: the net flux, i.e. in-flux minus out-flux. On average, clouds of all masses experience a net inflow of cold gas into them, with a higher net positive rate for more massive clouds.

The main panel shows the mass flux rates into (black) and out of (orange) clouds. While both are roughly ∼1.5 × 10³ M_⊙ Myr⁻¹ for clouds of mass ∼10⁵ M_⊙, the in-flux increasing sharply by a factor of around ten to ∼15 × 10³ M_⊙ Myr⁻¹ for clouds of mass ∼10^6.5 M_⊙. While the out-flux also scales with mass, the rise is less steep, with an average out-flux of ∼11 × 10³ M_⊙ Myr⁻¹ for clouds of mass ∼10^6.5 M_⊙.

The difference between the two, i.e. the net flux, is shown in the inset. Interestingly, on average, the net rate is positive for all mass clouds we study: a flux of ∼0.15 × 10³ M_⊙ Myr⁻¹ for low mass clouds (∼10⁵ M_⊙) increases to ∼4.5 × 10³ M_⊙ Myr⁻¹ for their more massive counterparts (∼10^6.5 M_⊙). While the net flux averaged over ∼300 Myr⁻¹ is positive, note that there may be clouds that experience a net out-flow of mass between snapshots, as we explore shortly. In conjunction with Fig. 8, this result suggests that, for small clouds, the in-flow of cold material has a sub-dominant impact on mass growth in comparison to cloud-cloud interactions: despite the (small) net in-flow of gas, their overall mass decreases as a result of repeated fragmentations.

In Fig. 13, we explore four important physical properties that impact the mass flux into and out of clouds. The main panels shows the in-flux (black) and out-flux (orange) rates. Solid curves correspond to the median, and the top panels show the difference between the two, i.e. the net mass flux. The horizontal dashed gray line demarcates the regime of net in-flow into the cloud (above the line, positive values) from that of a net out-flow. Similar to Fig. 11, the rates here correspond to those between neighboring snapshots, and properties are computed at the snapshot preceding the mass exchange. As before, while we show results stacked over clouds of all masses, we mention that the following trends are qualitatively similar at roughly fixed cloud mass, although typically steeper for more massive clouds.

Fig. 13. Main panels: In-flux (black) and out-flux (orange) rates as a function of velocity contrast between the cloud and the interface (vrel; left), cloud crushing time (tcc; center left), vorticity in the interface (ω; center right), and the dot product between the magnetic tension vector and density gradient in the interface, i.e. the effective tension (fT, eff = fT ⋅ $\mathrm{\nabla }\widehat{\rho }$; right). The top panel shows the difference between the two, i.e. the net mass flow rate. While clouds with relatively low values of vrel and ω experience a net in-flux of material, larger values of these properties trigger a net loss of cold gas from the cloud. The inverse trend holds for tcc, as higher values correspond to larger shear times, thereby leading to a lower fraction of mass loss from the cloud. The strength and topology of magnetic fields may further be important, as a net magnetic tension force acting against the density gradient, i.e. away from the cloud and into the background, results in a net flux of mass into the cloud, likely by inhibiting the mixing of the background into the cloud.

The left panel shows flux rates as a function of the velocity contrast between the cloud and the interface (v_rel). An average in-flux of ∼2.5 × 10³ M_⊙ Myr⁻¹ for clouds at similar velocities as their interface (v_rel ≲ 5 km s⁻¹) increases to ∼3.25 × 10³ M_⊙ Myr⁻¹ for clouds moving at larger relative velocities (v_rel ∼ 80 km s⁻¹). The growth of the out-flux rate is more steep, increasing from ∼1.8 × 10³ to ∼4.0 × 10³ M_⊙ Myr⁻¹ between the two extremities of this plot.

On average, the net flux is positive for relatively small velocity contrasts (v_rel ≲ 40 km/s), signifying a net inflow of cold gas into the cloud. For larger contrasts, on average, clouds experience a net loss of mass into the ambient gas. These results are qualitatively consistent with those from higher resolution idealized studies which have argued that a large velocity shear may lead to an overall loss of cloud mass as a result of fluid instabilities such as the KH-instability (e.g., Scannapieco & Brüggen 2015), or by resulting in an increase of the turbulent velocity that may further have an impact on cloud growth over time (e.g., Abruzzo et al. 2024).

In addition to v_rel, the size (R_cl) and the overdensity of the cloud (χ) influence the time-scale over which the cloud is impacted by these fluid instabilities. Commonly referred to as the cloud crushing time-scale, this is defined as t_cc = ${\chi }^{1/2}{R}_{\mathrm{cl}}{v}_{\mathrm{rel}}^{-1}$ (Klein et al. 1994). In the center-left panel, we explore the dependence of the flux rates on t_cc. While both the in- and out-fluxes anticorrelate with t_cc, the net flux increases toward larger t_cc and eventually flattens out. We construe such a trend to arise since clouds with larger t_cc have larger shearing time scales, thereby leading to a lower fraction of mass loss from the cloud due to mixing with the background (Scannapieco & Brüggen 2015). It is nonetheless important to note that the cooling time-scales of the interface may vary across the sample of clouds, which we explore in Sect. 3.3.

Fluid instabilities drive the mixing of material across the interface, effectively cycling matter between the cloud and the background (e.g., Fielding et al. 2020). In the center-right panel, we aim to quantify the impact of this on flux rates, through measurements of vorticity (ω = |∇ × v|) in the interface. Trends are similar to the left panel, as expected since v_rel and ω correlate strongly with each other (Ramesh et al. 2024): the in-flux rate increases from ∼1.9 × 10³ to ∼3.9 × 10³ M_⊙ Myr⁻¹ from low (∼10 km s⁻¹ kpc⁻¹) to high values of vorticity (∼95 km s⁻¹ kpc⁻¹). The out-flux rate increases more steeply, from ∼1.2 × 10³ to ∼4.2 × 10³ M_⊙ Myr⁻¹ over the same range of ω.

Overall, the net flow rate is positive for small values of ω, corresponding to a net flow of cold gas into the cloud. For large values of ω the net flow becomes negative, signifying a transfer of mass from the cloud into the ambient gas. We interpret these two limiting cases as arising from a tug-of-war between two key physical processes. As long as the rate of mass loss from the cloud due to mixing is smaller than the rate at which gas from the ambient gas cools and condenses onto the cloud, the cloud mass increases over time (e.g., Fielding & Bryan 2022).

In addition to radiative cooling, it is believed that non-thermal processes, such as magnetic fields, may be important. For instance, as clouds pass through the halo, they may sweep up ambient field lines to create a draped layer around them (e.g., Dursi & Pfrommer 2008; Pfrommer & Dursi 2010), which may thereafter suppress fluid instabilities along their interfaces (e.g., Sparre et al. 2020; Cottle et al. 2020). We next assess the possible impact that such a phenomenon may have on mass flux rates. We do so by first computing the magnetic tension, f_T = |B|²κ, where κ = (b ⋅ ∇)b is the magnetic curvature vector (Shen et al. 2003; Boozer 2005). Here, b = B/B is the unit vector in the direction of the magnetic field B. We have previously explored the amount of magnetic draping and its impact on CGM gas clouds (Ramesh et al. 2024).

Here, we take the dot product of the tension vector with the unit vector of the density gradient $\mathbf{\nabla }\widehat{\mathit{\rho }}$ to yield an effective tension force f_T, eff = f_T ⋅ $\mathbf{\nabla }\widehat{\mathit{\rho }}$. Since the gas density decreases from the cloud center toward its outskirts, and thereafter into the background (Nelson et al. 2020; Ramesh et al. 2023c), $\mathbf{\nabla }\widehat{\mathit{\rho }}$ points into the cloud from the background. Positive values of f_T, eff thus correspond to cases where the tension is directed into the cloud, while the opposite holds for negative values.

The right panel of Fig. 13 shows the different flux rates as a function of the effective tension, f_T, eff. On average, both the in- and out-flux rates are largely similar for positive values of f_T, eff, i.e. the net mass flux rate is ∼0 in this regime. However, for negative values of f_T, eff, we note a net in-flow of mass into clouds. We interpret this as the result of the tension force inhibiting a flow of gas from the background into the cloud (i.e. along the density gradient), thereby reducing mixing and vorticity, and eventually leading to a net in-flux of material. While draped configurations of magnetic fields are believed to possibly also suppress the impact of thermal conduction (Ettori & Fabian 2000; Brüggen et al. 2023), which might otherwise contribute to cloud evaporation and destruction (Marcolini et al. 2005; Vieser & Hensler 2007), the current GIBLE simulations do not include this physical mechanism. Future extensions that include conduction can explore its impact on cloud evolution in a cosmological context.

3.3. Comparison with theoretical expectations

So far we have qualitatively discussed our results in the context of previous work from the literature. We now provide a more quantitative comparison with analytic models and results from idealized simulation studies. In particular, we focus on the theoretical expectations for mass flux rates in to and out of clouds.

In Fig. 14, we show a number of quantities of interest as a function of cloud radius R_cl. The various solid colored curves correspond to medians for clouds segregated into different mass bins, scatter points to individual clouds, and the dashed gray lines to relevant scaling relations. Since the measurement of R_cl is sensitive to the galactocentric distance of the cloud, as this impacts the mean density of gas, we focus on clouds in the inner halo (∼0.3 R_200c) where statistics are the best (Ramesh et al. 2023c), and mention that these results and discussion hold true for clouds in other regions of the halo as well. As before, we restrict the selection of clouds to the final ten snapshots. Flux rates correspond to those between neighboring snapshots, and all properties of clouds are computed at the snapshot preceding the mass exchange.

Fig. 14. Mass in-flux (top panel), cloud surface area (center), and mean interface density (lower) as a function of radius for clouds in the inner halo (galactocentric distance ∼0.3 R200c). The various solid colored curves correspond to median relations for clouds segregated into different mass bins, scatter points to individual clouds, and the dashed gray lines to relevant scaling relations. See main text for more discussion.

In the top panel, we begin with the mass in-flux into the cloud, Ṁ_in. In all four mass bins, Ṁ_in increases gradually toward larger R_cl, and is well represented by a ∝R_cl¹ scaling in all cases. In a recent study by Tan & Fielding 2024, it was suggested that Ṁ_in ∝ ${R_{\mathrm{cl}}}^{11/4}$ for clouds in their tall box setup. We note, however, that this was the cumulative scaling for all clouds, i.e. not split by mass. In our simulations, we find that M_cl ∝ ${R_{\mathrm{cl}}}^{3/2}$ (not shown), implying a Ṁ_in ∝ ${R_{\mathrm{cl}}}^{5/2}$ scaling if clouds of all masses are stacked together, yielding a largely similar power law slope as reported by Tan & Fielding (2024).

To understand the Ṁ_in ∝ R_cl¹ scaling further, we considered a model for cloud growth where the in-flux into the cloud (Ṁ_in) is a function of the surface area, A_cl; density of the interface, ρ_int; and the inflow velocity of gas into the cloud v_in (see Gronke & Oh 2020; Tan & Fielding 2024):

$$\begin{array}{c}\hfill \dot{m_{\mathrm{in}}}\sim A_{\mathrm{cl}}·{\rho }_{\mathrm{int}}·v_{\mathrm{in}}.\end{array}$$(1)

The center panel of Fig. 14 shows A_cl as a function of R_cl. All curves follow a ∝ ${R_{\mathrm{cl}}}^{5/2}$ scaling, with a tight correlation of the scatter points around their respective medians. This is in remarkable agreement with results from idealized setups run at much higher resolution that have shown how the surface areas of clouds that continuously experience a shear force scale as A_cl ∝ ${R_{\mathrm{cl}}}^{5/2}$ (Gronke & Oh 2023; Tan et al. 2023), corresponding to a fractal dimension D = 2.5 in the mixing layer (Fielding et al. 2020).

The variation of ρ_int with R_cl, shown in the bottom panel, is non-trivial. While it is typically assumed that ρ_int is independent of R_cl, we find a reasonably strong ρ_int ∝ ${R_{\mathrm{cl}}}^{-3/2}$ scaling in all four mass bins. While not shown here, we mention that ρ_cloud, the mean density of the cloud, also is $\propto {R_{\mathrm{cl}}}^{-3/2}$, and ρ_int inherits this scaling as clouds aim to maintain pressure balance with their surroundings (Ramesh et al. 2023a). We find that these trends are driven by the inhomogeneities of the CGM at fixed galactocentric distance and fixed cloud mass, wherein clouds with different mean densities exist, which would then naturally lead todifferent interface densities. Such large pressure fluctuations imply that clouds embedded in low pressure environments have lower densities and thus larger radii at fixed mass (see also Stern et al. 2021).

Radiative turbulent mixing layer simulations suggest that the final ingredient of the model, v_in, can be written as (Fielding et al. 2020)

$$\begin{array}{c}\hfill v_{\mathrm{in}}\sim t_{\mathrm{cool},\min }^{-(1-\frac{1}{2(1-k)}))}·L_{\mathrm{turb}}^{1-\frac{1}{2(1-k)}}·v_{\mathrm{turb},\mathrm{L}}^{\frac{1}{2(1-k)}},\end{array}$$(2)

where t_cool, min is the minimum cooling time in the interface layer, L_turb is the outer scale of turbulence, and v_turb, L is the turbulent velocity at large scales, i.e. the typical length scale along of turbulent driving. Lastly, k is the power-law slope of the velocity structure function (VSF), i.e. VSF ∝ l^k, where l is the separation scale. To compute v_in, one thus first needs to estimate k.

In Fig. 15, we show VSFs for a subset of clouds from the M_cl ∼ 10^5.2 M_⊙ bin. In particular, to probe possible trends as a function of R_cl, we pick the extreme 10^th percentiles of clouds based on size. VSFs of the smaller (larger) set of clouds are shown on the left (right). Thin colored curves correspond to VSFs of individual clouds, while the thick black curve to the median. We compute the VSF with the first three layers of Voronoi cells around clouds, i.e. the interface plus two further background layers, although we have verified that the following results are not sensitive to this exact choice. The vertical dashed gray curves correspond to the average size of the cells selected for the analysis.

Fig. 15. Velocity structure functions for a subset of Mcl ∼ 105.2±0.2 M⊙ clouds in the inner halo (galactocentric distance ∼0.3 R200c). The extreme tenth percentile of clouds were picked based on cloud radius, with the VSFs of smaller (larger) clouds shown on the left (right). Thin colored curves correspond to individual clouds, thick solid curves to median, and the vertical dashed gray line to the average cell size. Both medians and a vast majority of the individual curves follow a ∝l1/2 scaling, suggesting that the power law slope of the VSF is k ∼ 1/2.

To guide the eye, the red dashed curves show a ∝ l^1/2 scaling. Medians in both panels, as well as most individual curves, seem to closely trace this scaling, suggesting that the power-law index k of the VSFs in the regions around these clouds is ∼1/2. Although not shown, VSFs of other clouds, i.e. other masses and galactocentric radii, show similar behavior. The VSFs in both panels also appear to flatten out at roughly the same separation scales, implying that L_turb is largely independent of R_cl. This is in agreement with Tan & Fielding (2024), where the same scaling was suggested for their sample of clouds.

We note, however, that resolution is limited on these scales, and it is possible that the inferred VSFs are influenced by numerical dissipation (e.g., Bauer & Springel 2012; Zier & Springel 2022). Comparison with higher resolution cosmological simulations run in the future will help assess and verify the accuracy of these inferences. For the present work, adopting that k = 1/2, Equation (2) reduces to $v_{\mathrm{in}}\sim v_{\mathrm{turb},\mathrm{L}}^1$, i.e. v_in in this case depends only on v_turb, L, and is independent of t_cool, min and L_turb. To estimate the dependence of v_in on R_cl, we need to assess the analogous scaling of v_turb, L.

Measuring v_turb, L directly from the simulation output requires knowledge of the associated turbulence mixing length scale, L_turb. While L_turb can be derived from the VSF (see, e.g., Abruzzo et al. 2024), we instead adopt the simpler approach of comparing the normalizations of the VSFs as a function of R_cl. As seen in the two panels of Fig. 15, and many other cases that we have tested, this is largely independent of R_cl, implying that v_turb, L ∝ ${R_{\mathrm{cl}}}^0$.

Plugging in these different scaling relations (A_cl ∝ ${R_{\mathrm{cl}}}^{5/2}$, ρ_int ∝ ${R_{\mathrm{cl}}}^{-3/2}$ and v_in ∼ v_turb, L ∝ ${R_{\mathrm{cl}}}^0$) into Equation (1) yields $\dot{m_{\mathrm{in}}}\sim {R_{\mathrm{cl}}}^{5/2}·{R_{\mathrm{cl}}}^{-3/2}·{R_{\mathrm{cl}}}^0$ ⟹ $\dot{m_{\mathrm{in}}}\sim {R_{\mathrm{cl}}}^1$, in excellent agreement with the top panel of Fig. 14. These results are consistent with k = 1/2 and the corresponding limit of the cloud growth model. We note once again, however, that resolution on these scales is limited in our simulations, and a convergence study toward even better numerical resolution will clarify this picture.

As a final point of comparison, we assess the relation between the mass in-flux into the cloud and the minimum cooling time in the interface. Turbulent mixing layer simulations predict the flux to scale as ∝ $t_{\mathrm{cool},\min }^{-1/4}$ in the fast cooling regime, with a larger negative power law slope (∝ $t_{\mathrm{cool},\min }^{-1/2}$) in the regime of slow cooling (Ji et al. 2018; Fielding et al. 2020; Tan et al. 2021).

In Fig. 16, we examine these trends for clouds in our simulation. As before, we select clouds in the inner halo (galactocentric distance ∼0.3 R_200c) to excise the radial trend of t_cool, min, and segregate clouds into bins based on their mass. The corresponding medians are shown through colored solid curves, while the dotted gray lines show scalings of $\propto t_{\mathrm{cool},\min }^{-1/4}$.

Fig. 16. Mass in-flux as a function of tcool, min for clouds in the inner halo (galactocentric distance ∼0.3 R200c). The colored curves correspond the medians of different mass bins, and the gray curves to relevant scaling relations. The mass in-flux anticorrelates with the cooling time, with a typical power law slope of ∼ − 1/4, suggesting that these clouds are in the fast cooling regime.

The mass in-flux in all cases anticorrelates with the minimum cooling time, i.e. clouds with faster cooling interfaces experience a larger in-flux of material. The power law-slope in all cases is ∼ − 1/4, suggesting that these clouds are in the regime of fast cooling. Although not shown explicitly, we find that the net flux also scales as $t_{\mathrm{cool},\min }^{-1/4}$ at low-t_cool, min, while the power law slope tends to more negative values at larger t_cool, min. A more quantitative comparison with the mass flux rates expected from theoretical models will be enabled by future simulations.

4. Summary

In this paper we have quantified the origin and evolution of cold gas clouds in Project GIBLE, a suite of cosmological magneto-hydrodynamical zoom-in simulations with super-Lagrangian refinement in the CGM down to a mass resolution of ∼10³ M_⊙. Focusing on halos of Milky Way-like galaxies, we used Monte-Carlo tracers to identify progenitors and descendants of clouds across time. We constructed cloud merger (interaction) trees and used them to identify cloud origins and explore the mass growth histories of clouds. Our main findings are as follows:

1. The population of z = 0 clouds primarily arises from two roughly balanced origins: recent (≲2 Gyr) outflows of cold gas driven from the galaxy into the CGM (∼40–60%) and precipitation of hot halo gas (∼25–50%). More massive clouds are more likely to arise from outflows, while low-mass clouds are more likely to originate from hot halo cooling. Only ∼5% of clouds are direct descendants of cold gas-stripped or otherwise removed from satellites (Figs. 1 and 2).

2. Properties of z = 0 clouds are diverse and depend of their origin. Clouds stripped from satellites are typically smaller in size and situated at larger galactocentric distances. The distribution of cloud metallicity is bimodal, with a primary peak at ∼0.45 Z_⊙ and a smaller secondary peak at ∼0.15 Z_⊙. These features reflect specific accretion histories of massive satellites. In contrast, clouds originating from the central galaxy or hot halo condensation have intermediate metallicities of ∼0.25 Z_⊙ (Fig. 3).

3. Clouds that precipitate out of the hot halo begin from a cold overdense seed. Rapid cooling and accretion leads to its growth over time, and such cold seeds eventually reach a mass of ∼10⁵ M_⊙ by z = 0. This exact value decreases as numerical resolution improves. Cold seeds are primarily created as a result of short timescale (≲500 Myr) dissolution of previous cold gas into the hot phase, which produces density fluctuations sufficient for the next generation of clouds to form (Fig. 4).

4. The baryonic material that builds cool CGM clouds, as with all baryonic matter in halos, is ultimately of cosmological origin. Roughly ∼85% of their mass comes from smooth accretion from the IGM at earlier cosmic epochs (Fig. 5).

5. We make predictions for the origin of cold clouds in the Milky Way halo. While LVCs and IVCs are roughly equally likely to arise as a result of cold galactic outflows and hot-phase precipitation, outflowing (inflowing) HVCs are predominantly signatures of the former (latter). Clouds arising from the cold gas of satellites are rare but equally likely at all line-of-sight velocities (Fig. 6).

6. Roughly 10% of the z = 0 cloud population is long-lived, with their progenitors having already assembled ∼2 Gyr ago. These clouds have unique and diverse evolutionary tracks as they pass through the halo, interact with other clouds, and exchange mass with their surroundings. Clouds portray a temperature gradient into their surroundings with a warm interface layer separating them from the hot background, which is believed to improve chances of their survivability (Figs. 7 and 8).

7. Long-lived clouds merge and fragment at an average rate of ≳20 Gyr⁻¹. Above some critical mass threshold, more massive clouds undergo more mergers than fragmentations (Fig. 10). This results in net mass growth over time. The frequency of these interactions is modulated by the pressure imbalance between clouds and their interfaces, with under-pressurized clouds being more likely to fragment. The clustering of clouds around other clouds plays a secondary role, and those with many neighbors are more likely to merge (Fig. 11).

8. Clouds also exchange matter with their surrounding gas, at a rate of ≳10³ M_⊙ Myr⁻¹, and massive clouds experience a larger net flux (Fig. 12). Small-scale physical properties have an impact on mass exchange: the net flux into the cloud anticorrelates (correlates) with the velocity contrast as well as interface-layer vorticity (cloud crushing timescale). The topology and magnitude of magnetic fields may further play a role, as we find that a net magnetic tension force acting against the density gradient can help suppress mixing of the background into the cloud (Fig. 13). In addition, the flux rates in and out of clouds may also depend on the surface area as well as the cooling time of gas in the interface (Figs. 14 and 16).

This work provides an initial set of cosmological inputs to questions related to the origin and evolution of cool clouds in the CGM of galaxies. Clear avenues for future extensions exist. What is the fate of these clouds? Do they fuel the cold gas supply of the central galaxy to sustain star formation or perish as they accrete toward the center? If they do increase the galactic cold mass budget, does this impact future outflows and the production of new clouds? We reserve such explorations to futurestudies.

Acknowledgments

RR and DN acknowledge funding from the Deutsche Forschungsgemeinschaft (DFG) through an Emmy Noether Research Group (grant number NE 2441/1-1). RR is a Fellow of the International Max Planck Research School for Astronomy and Cosmic Physics at the University of Heidelberg (IMPRS-HD). RR and DN thank Ralf Klessen and Annalisa Pillepich for valuable inputs. MB acknowledges funding by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy – EXC 2121 “Quantum Universe” – 390833306. This work has made use of the VERA supercomputer of the Max Planck Institute for Astronomy (MPIA) operated by the Max Planck Computational Data Facility (MPCDF), and of NASA’s Astrophysics Data System Bibliographic Services.

References

Appendix A: Tracing Cosmic Gas Accretion

In this section, we extend the analysis of Fig. 5, dissecting the various sources of gas that populate all clouds at z = 0, not just seeds as previously explored. We consider five distinct sources: (a) gas from the central galaxy (all phases of gas included) shown in black, (b) the cold CGM, i.e. the cold phase of the main progenitor or that of those clouds that merged, in red, (c) the warm-hot CGM in orange, (d) gas from satellite galaxies (all phases included) in blue, and (e) the IGM in gray. We note that the black curve has been multiplied by a factor of 25 to improve visibility.

Clouds that arise from the central galaxy are shown in the left panel of Fig. A.1. A very small fraction of tracers were in the galaxy in the past (≲ 2 %¹⁵). Although counter-intuitive, this simply reflects the rapid cycling of matter in and out of clouds, and that the present gas content of a cloud usually carries little to no information about its past (see also Nelson et al. 2020). Owing to poor time resolution at snapshots beyond z ∼ 0.15, the fraction of mass from the galaxy may be underestimated, since strong recycling could rapidly remove accreted gas from the central, biasing this contribution to instead be counted in the satellite and IGM categories. In either case, this suggests that the recycling of gas from the galaxy into the CGM predominantly took place ≳ 2 Gyr ago (Ford et al. 2016).

Fig. A.1. Similar to Fig. 5, we dissect the various sources of gas that populate z = 0 clouds, split based on their origin source: those arising from the central galaxy in the left panel, through condensation of hot halo gas in the center, and from satellites in the right. The five different colors correspond to the distinct pathways that we consider. Similar to cold seeds, ≳ 80 % of the z = 0 mass was in the IGM at z ∼ 2.

On short time scales (≲ 500 Myr), a large fraction of mass was either in the cold CGM, i.e. progenitor clouds, or in the warm-hot phase of the CGM. At t_lookback ≳ 1 Gyr, most of the mass is extragalactic, after which it is smoothly accreted from the IGM (∼ 85 %), or in a clumpy fashion (∼ 15 %). We note that ∼ 14 % of mass in the z = 0 central galaxy is in satellites at z ∼ 2, while the remaining ∼ 86 % is in the IGM. For the z = 0 warm-hot CGM, these numbers are ∼ 13 and 87 % respectively. Thus, (a) satellites contribute a slightly larger fraction to cold seeds (Fig. 5) and a fraction of z = 0 clouds (left and center panels of Fig. A.1) than to the central galaxy and the warm-hot CGM, (b) the IGM contributes a smaller fraction to the central than to the warm-hot CGM, suggesting that only a fraction of material from the IGM eventually collapses into the central, while the rest remains in the CGM due to pressure support, and (c) satellites contribute a slightly larger fraction to the central than to the warm-hot gas, implying that not all satellite gas is stripped during the passage through the CGM.

In the center panel, we examine clouds that arise through the precipitation of the hot halo. Numbers are largely similar to the left-panel, with the main exception that a greater fraction of their z = 0 mass is in the warm-hot CGM ∼ 500 Myr ago. The mass fraction from the central galaxy is ≲ 1 %, reflecting the fraction that the cold seed acquires from the galaxy (Fig. 5).

The right-panel focuses on clouds ejected from satellites. While a large fraction of gas is in other CGM clouds in the recent past (≲ 200 Myr), i.e. clouds with which merger events took place, ≳ 10 % of the current gas can be traced back to satellites. The inner core of these clouds is thus retained, while gas in the outer regions of the cloud is recycled. This is in contrast to clouds arising from the galaxy, in which case a greater amount of gas from the core is cycled out. We speculate that this is due to clouds ejected from satellites typically being denser (top-left panel of Fig. 3), thereby allowing them to retain their cores for a longer time. Lastly, ∼ 80 % of mass of clouds in this category is in the IGM at z ∼ 2, consistent with the scenario that IGM gas is first accreted into these satellites when they were centrals of their own halos, prior to their infall.

Appendix B: Impact of Limited Time Resolution

In this section, we assess the impact of (limited) time spacing between snapshots on a subset of results related to merger events discussed above in the main text. For this, we make use of a GIBLE halo for which we have saved a hundred snapshots between z = 0 and z ∼ 0.01 (average snapshot spacing of ∼ 1.4 Myr) and two hundred between z ∼ 0.01 and z ∼ 0.1 (average spacing of ∼ 6.0 Myr). These snapshots do not contain Monte Carlo tracers that we utilize in the main work, and hence we instead used the ParticleIDs of gas cells to construct cloud merger trees. We verified that this is a good alternative when the time cadence is as good as ≲ 6.0 Myr, particularly for those clouds that contain ≳ 50 gas cells, by comparing against the merger tree constructed using tracers at the snapshots used in the main body of this work. Since one can neither compute mass fluxes in and out of clouds without Monte Carlo tracers, nor unambiguously assign origin channels, we restrict the following analysis to merger and fragmentation events.

As mentioned in Section 3.2.1, between two snapshots, a cloud may fragment into multiple objects, or multiple clouds may merge together, or both a fragmentation and merger could take place. We therefore checked the prevalence of this issue in limited time cadence snapshots. To do so, we select z = 0 long-lived clouds, i.e. those z = 0 clouds that were already assembled ∼ 2.0 Gyr ago, with masses in the range [10^5.0, 10^5.4] M_⊙. Fig. B.1 shows the number of clouds that interact during a merger event between two neighboring snapshots, as a function of snapshot time cadence. Crosses show the mean, while error bars correspond to the 1σ variation. While the average is quite high (∼ 5) when the snapshot cadence is as coarse as ∼ 100 Myr, at the best available time resolution (∼ 1.4 Myr), the average number is close to 2, i.e. events are pairwise (binary) merger between clouds. Higher numbers at worse time resolution are indicative of multiple binary mergers occuring between two neighboring snapshots. Though not shown explicitly, a similar trend holds for fragmentation events as well.

Fig. B.1. Number of clouds that interact during a merger event between two neighboring snapshots (y-axis) as a function of time cadence (x-axis). Crosses show the mean, while error bars correspond to the 1σ variation. While the average is ∼ 5 when the snapshot cadence is as coarse as ∼ 100 Myr, it is close to 2 at out best available snapshot spacing, suggesting that all pairwise merger events are correctly identified.

Similarly, limited snapshot time cadence also has an impact on the computed merger rates between clouds. In Fig. B.2, we show the merger rate as a function of the average time resolution, for the same selection of clouds as above. Crosses again show the mean, and error bars correspond to the 1σ variation. An average rate of ∼ 14 Gyr⁻¹ (11 Gyr⁻¹) at a snapshot spacing of ∼ 30 Myr (100 Myr) is clearly an underestimate, and average merger rates as high as ∼ 40 Gyr⁻¹ can be seen at improved time resolutions. Interestingly, the average merger rate is roughly consistent for snapshot spacings of ∼ 1.0 and 6.0 Myr, suggesting that these numbers reflect a converged value, although the 1σ variation continues to increase. Resolving collective cloud dynamics in time thus requires a snapshot cadence at least as good as a few Myr.

Fig. B.2. Merger rate of clouds (y-axis) as a function of snapshot time resolution (x-axis). Crosses correspond to the mean, while error bars show the 1σ variation. Merger rates at coarse time resolution are underestimated, and approach a converged value at snapshot spacings of ≲ 6.0 Myr.

All Figures

Fig. 1. Visual impression of the three origin sources of z = 0 cold CGM gas clouds. The white circle marks the virial radius of this GIBLE halo (S105; M200c ∼ 1012.29 M⊙; see Table 2 of Ramesh & Nelson 2024 for more details). Of the many thousand clouds in this CGM, most came to existence in the recent past (≲2 Gyr). A majority of these can be traced back to cold gas that was previously present in the central galaxy (white). On average across the sample of eight halos, these form ∼40–60% of the z = 0 population. A small fraction of clouds (∼3–5%) arise from satellites as they spiral in toward the center of the halo (orange). Colored in blue, most of the remainder of the clouds are formed through the in-situ condensation of hot halo gas. These account for ∼15–50%, depending on cloud mass. This simple visualization suggests that the inner regions are dominated by white clouds, while their blue and orange counterparts are situated farther out in the halo. This radial distribution is typical, as we expand upon in subsequent figures.

In the text

Fig. 2. Main panel: Fraction of clouds arising from the three origin channels as a function of (z = 0) cloud mass (indicated by thick solid curves). A large fraction (∼40–60%) were previously present in the central galaxy and have been driven out into the CGM in the last ≲2 Gyr (black). The ejection of gas from satellites gives rise to a small portion of the population (∼3–5%; orange). A majority of the remaining clouds are formed through condensation of the hot phase of the CGM in the last ∼2 Gyr. The thin gray curve shows the fraction of clouds that are long-lived, i.e. those whose progenitors already exist at least ∼2 Gyr ago (∼2–10%). Top panel: Assessment of numerical convergence. The fraction of clouds arising through in situ condensation of the hot phase in the RF512 suite (solid curves) is contrasted against the lower resolution RF64 (dashed) and RF8 (dotted) runs. Results between resolution levels approach convergence toward the high mass end, i.e. well above the resolution limit.

In the text

Fig. 3. PDFs of properties of z = 0 CGM clouds split based on their origin channel. Curves corresponding to those arising as from the central galaxy (within the last ∼2 Gyr) are colored in black, condensates of the hot halo gas in red and those ejected from satellites in orange, as shown by the legend in the top-left panel. Distributions of cloud radius, metallicity, radial velocity and galactocentric distance are shown in the top-left, top-right, lower-left and lower-right panels, respectively. In the right panels, we show individual distributions of the S98 and S105 halos (see Table 2 of Ramesh & Nelson 2024 for more details) through dotted and dashed curves, respectively (vertically offset for better visibility). Distributions are broad for all four properties shown here, with PDFs of clouds ejected from satellites typically different as compared to the other two cases.

In the text

Fig. 4. Evolutionary tracks of clouds that condense out of the hot halo. Solid curves correspond to the median while the shaded bands show the 16th-84th percentile regions. Colors correspond to different epochs where the cloud first came into existence. In all cases, the origin of the cloud can be traced back to a ‘cold seed’, typically comprised of only a few gas cells. Following that, clouds accrete matter to grow to a median mass of ∼105.2 M⊙ by z = 0. The percentile regions, however, are broad, implying that evolutionary tracks of individual clouds can be diverse.

In the text

Fig. 5. Dissecting the different sources of gas that seed the condensation of cold clouds out of the hot halo at z ∼ 0. On short timescales (≲500 Myr), the dissolution of past clouds into the hot phase, and the subsequent mixing, is the primary source, and it is also potentially stirred by outflows from the central galaxy into the CGM. Transitioning out into previous cosmic epochs as early as z ∼ 2 (tlookback ∼ 10.5 Gyr), the mass giving rise to these seeds was predominantly smoothly accreted from the IGM into the halo (∼85%), while clumpy accretion through satellites contributes the remaining ∼15%.

In the text

Fig. 6. Probabilities of the three origin channels as a function of line-of-sight velocities of clouds in the frame of the local standard of rest. The vertical dashed lines delineate commonly used definitions to classify clouds in the Milky Way halo. Outflowing clouds predominantly trace cold gas arising from the central in the recent past (≲2 Gyr), while their inflowing counterparts typically condense out of the hot phase. Clouds ejected from satellites are relatively rare (≲5%), but roughly equally likely at all (line-of-sight) velocities.

In the text

Fig. 7. Visualization of cold clouds in a GIBLE halo (S98; see Table 2 of Ramesh & Nelson 2024 for more details) at z = 0 using a ray-traced volume rendering sampled by five narrow Gaussian transfer functions in temperature space, as shown by the color bar (background image). The white circle denotes the outer boundary of the CGM refinement region, i.e. the virial radius of the halo (∼270 kpc). In the foreground, we show the (past) trajectory of a ∼106.1 M⊙ cloud identified at z = 0, with the various translucent stars highlighting the positions of the main progenitor of this cloud at 25 distinct snapshots back to z ∼ 0.15. The insets at the corners of the image show a volume rendering of (the main progenitor of) this cloud at four distinct redshifts, going back in time in clockwise direction starting from the bottom left, with the cloud oriented such that its mean velocity vector points to the right. In the ∼2 Gyr between z ∼ 0.15 and 0, this cloud clearly evolves in terms of shape, size (characterized by the varying scale bars in the insets) and mass (shown by the different labels; in units of M⊙). Also visible in the top-right inset is an instance of an ongoing merger between clouds, a phenomenon explored and quantified later in this paper.

In the text

Fig. 8. Main panel: Evolution of masses of z = 0 clouds whose main progenitors already existed ∼2 Gyr ago (z ∼ 0.15). Solid curves correspond to medians, shaded bands to the 16th-84th percentile regions, and the colors to four different mass bins. On average, progenitors of clouds from all four bins were massive (∼106.25 M⊙) 2 Gyr ago, albeit with relatively broad percentile regions in all cases. In the top panel, we show individual tracks for a subset of clouds split into three categories, chosen randomly: those whose masses have overall increased (left), remained overall constant (center) and overall decreased (right) in this time frame. For exactly one curve in each case, we show all instances of mergers and fragmentations as identified in our cloud merger trees: the former are shown in blue dots, the latter in red, and overlapping dots of both colors correspond to instances where both these events take place between neighboring snapshots. As expected, sudden jumps in curves are associated with these events. Finally, in the inset of the main panel, we assess numerical convergence by contrasting results of RF512 (solid) with the lower resolution RF64 (dashed) and RF8 (dotted) runs for the case of (z = 0) ∼106 and 106.5 M⊙ clouds. While the three runs are qualitatively similar, curves converge at a higher mass for the lower resolution runs (see text for discussion). We note that the axes labels of the inset and the top panels are identical to that of the main panel.

In the text

Fig. 9. Visualization of two clouds merging into one (top row) and of a cloud fragmenting into two smaller objects (bottom row). In both rows, neighboring panels are separated in time by ∼30 Myr. The scale bars of all panels of a given row are the same, shown in the lower right of each. All panels are oriented such that the mean velocity vector of the cloud (system) is pointed toward the right. Panels with two colors show times when two distinct clouds are present, while panels with only one color contain a single cloud.

In the text

Fig. 10. Main panel: Merger and fragmentation rates as a function of (z = 0) cloud mass. The former are shown in black, and the latter are in orange. A clear trend of event rates with cloud mass is visible, with lower rates for smaller clouds. Top panel: Difference in rates between merger and fragmentation events. On average, massive clouds merge more often than they fragment, while the opposite is true for their less massive counterparts.

In the text

Fig. 11. Main panels: Merger and fragmentation rates as a function of pressure contrast between the cloud and the interface (left) and Δ10 (right). Top panels: Difference between merger and fragmentation rates. A strong trend with pressure contrast can be seen, with under-pressurized clouds being more likely to undergo a fragmentation event as opposed to a merger. The trend with Δ10 is weaker – clouds with a very small (large) number of neighbors are likely to fragment (merge), while the net rate is largely flat for intermediate values of Δ10.

In the text

Fig. 12. Main panel: Mass in- and out-flux as a function of (z = 0) cloud mass. The former are shown in black, while the latter are in orange. Both rates correlate with cloud mass and are higher for more massive clouds. In the inset, we show the difference between the two rates: the net flux, i.e. in-flux minus out-flux. On average, clouds of all masses experience a net inflow of cold gas into them, with a higher net positive rate for more massive clouds.

In the text

Fig. 13. Main panels: In-flux (black) and out-flux (orange) rates as a function of velocity contrast between the cloud and the interface (vrel; left), cloud crushing time (tcc; center left), vorticity in the interface (ω; center right), and the dot product between the magnetic tension vector and density gradient in the interface, i.e. the effective tension (fT, eff = fT ⋅ $\mathrm{\nabla }\widehat{\rho }$; right). The top panel shows the difference between the two, i.e. the net mass flow rate. While clouds with relatively low values of vrel and ω experience a net in-flux of material, larger values of these properties trigger a net loss of cold gas from the cloud. The inverse trend holds for tcc, as higher values correspond to larger shear times, thereby leading to a lower fraction of mass loss from the cloud. The strength and topology of magnetic fields may further be important, as a net magnetic tension force acting against the density gradient, i.e. away from the cloud and into the background, results in a net flux of mass into the cloud, likely by inhibiting the mixing of the background into the cloud.

In the text

Fig. 14. Mass in-flux (top panel), cloud surface area (center), and mean interface density (lower) as a function of radius for clouds in the inner halo (galactocentric distance ∼0.3 R200c). The various solid colored curves correspond to median relations for clouds segregated into different mass bins, scatter points to individual clouds, and the dashed gray lines to relevant scaling relations. See main text for more discussion.

In the text

Fig. 15. Velocity structure functions for a subset of Mcl ∼ 105.2±0.2 M⊙ clouds in the inner halo (galactocentric distance ∼0.3 R200c). The extreme tenth percentile of clouds were picked based on cloud radius, with the VSFs of smaller (larger) clouds shown on the left (right). Thin colored curves correspond to individual clouds, thick solid curves to median, and the vertical dashed gray line to the average cell size. Both medians and a vast majority of the individual curves follow a ∝l1/2 scaling, suggesting that the power law slope of the VSF is k ∼ 1/2.

In the text

	Fig. 16. Mass in-flux as a function of tcool, min for clouds in the inner halo (galactocentric distance ∼0.3 R200c). The colored curves correspond the medians of different mass bins, and the gray curves to relevant scaling relations. The mass in-flux anticorrelates with the cooling time, with a typical power law slope of ∼ − 1/4, suggesting that these clouds are in the fast cooling regime.
In the text

Fig. A.1. Similar to Fig. 5, we dissect the various sources of gas that populate z = 0 clouds, split based on their origin source: those arising from the central galaxy in the left panel, through condensation of hot halo gas in the center, and from satellites in the right. The five different colors correspond to the distinct pathways that we consider. Similar to cold seeds, ≳ 80 % of the z = 0 mass was in the IGM at z ∼ 2.

In the text

Fig. B.1. Number of clouds that interact during a merger event between two neighboring snapshots (y-axis) as a function of time cadence (x-axis). Crosses show the mean, while error bars correspond to the 1σ variation. While the average is ∼ 5 when the snapshot cadence is as coarse as ∼ 100 Myr, it is close to 2 at out best available snapshot spacing, suggesting that all pairwise merger events are correctly identified.

In the text

	Fig. B.2. Merger rate of clouds (y-axis) as a function of snapshot time resolution (x-axis). Crosses correspond to the mean, while error bars show the 1σ variation. Merger rates at coarse time resolution are underestimated, and approach a converged value at snapshot spacings of ≲ 6.0 Myr.
In the text