We present an analysis of the cold gas phase in a lowmetallicity starburst generated in a highresolution hydrodynamical simulation of a gasrich dwarf galaxy merger as part of the GRIFFIN project. The simulations resolve 4 Minlineformulatexmath idTM0001 notationLaTeXodottexmathinlineformula gas phase mass resolution inlineformulatexmath idTM0002 notationLaTeXsimtexmathinlineformula0.1 pc spatial resolution the multiphase interstellar medium with a nonequilibrium chemical heatingcooling network at temperatures below 10inlineformulatexmath idTM0003 notationLaTeX4texmathinlineformula K. Massive stars are sampled individually and interact with the interstellar medium ISM through the formation of H II regions and supernova explosions. In the extended starburst phase the ISM is dominated by cold inlineformulatexmath idTM0004 notationLaTeXTmathrmgas lt 300texmathinlineformula K filamentary clouds with selfsimilar internal structures. The clouds have masses of inlineformulatexmath idTM0005 notationLaTeX102.6texmathinlineformulainlineformulatexmath idTM0006 notationLaTeX105.6texmathinlineformula Minlineformulatexmath idTM0007 notationLaTeXodottexmathinlineformula with a powerlaw mass function inlineformulatexmath idTM0008 notationLaTeXmathrm dNmathrm dM propto Malphatexmathinlineformula with inlineformulatexmath idTM0009 notationLaTeXalpha  1.78 pm 0.08texmathinlineformula. They also follow the Larson relations in good agreement with observations. We trace the lifecycle of the cold clouds and find that they follow an exponential lifetime distribution and an efolding time of inlineformulatexmath idTM0010 notationLaTeXsimtexmathinlineformula3.5 Myr. Clouds with peak masses below inlineformulatexmath idTM0011 notationLaTeX104texmathinlineformula Minlineformulatexmath idTM0012 notationLaTeXodottexmathinlineformula follow a powerlaw relation with their average lifetime inlineformulatexmath idTM0013 notationLaTeXtau mathrmlife propto M0.3mathrmmaxtexmathinlineformula which flattens out for higher cloud masses at inlineformulatexmath idTM0014 notationLaTeXlt 10texmathinlineformula Myr. A similar relation exists between cloud size at peak mass and lifetime. This simulation of the evolution of a realistic galactic cold cloud population supports the rapid formation and disruption of starforming clouds by stellar radiation and supernovae on a timescale less than 10 Myr.
['Astrophysics - Astrophysics of Galaxies']
The masses structure and lifetimes of cold clouds in a highresolution simulation of a lowmetallicity starburst
{'No Header': ', 1-16 (2024)', 'The masses, structure and lifetimes of cold clouds in a high-resolution simulation of a low metallicity starburst': 'Constantina M. Fotopoulou 1 ★ , Thorsten Naab 1 , Natalia Lahén 1 , Miha Cernetic 1 , Tim-Eric Rathjen 2 , Ulrich P. Steinwandel 3 , Jessica M. Hislop 4 , Stefanie Walch 2 ,', 'Peter H. Johansson 4': '1 Max Planck Institute for Astrophysics, Karl-Schwarzschild-Str. 1, D-85740, Garching, Germany  \n- 2  \nI. Physikalisches Institut, Universitat zu Köln, Zülpicher Str. 77, D-50937 Köln, Germany Center for Computational Astrophysics, Flatiron Institute, 162 5th Avenue, New York, NY 10010, USA  \n- 3  \n4 Department of Physics, University of Helsinki, Gustaf Hällströmin katu 2, FI-00014 Helsinki, Finland  \nAccepted XXX. Received YYY; in original form ZZZ', 'ABSTRACT': 'We present an analysis of the cold gas phase in a low metallicity starburst generated in a highresolution hydrodynamical simulation of a gas-rich dwarf galaxy merger as part of the /g.pc/r.pc/i.pc/f.pc/f.pc/i.pc/n.pc project. The simulations resolve (4 M ⊙ gas phase mass resolution, ∼ 0.1 pc spatial resolution) the multi-phase interstellar medium with a non-equilibrium chemical heating/cooling network at temperatures below 10 4 K. Massive stars are sampled individually and interact with the ISM through the formation of HII regions and supernova explosions. In the extended starburst phase, the ISM is dominated by cold ( 𝑇 gas < 300 K) filamentary clouds with self-similar internal structures. The clouds have masses of 10 2 . 6 - 10 5 . 6 M ⊙ with a power law mass function, 𝑑𝑁 / 𝑑𝑀 ∝ 𝑀 𝛼 with 𝛼 = -1 . 78 (± 0 . 08 ) . They also follow the Larson relations, in good agreement with observations. We trace the lifecycle of the cold clouds and find that they follow an exponential lifetime distribution and an e-folding time of ∼ 3.5 Myr. Clouds with peak masses below 10 4 M ⊙ follow a power law relation with their average lifetime 𝜏 life ∝ 𝑀 0 . 3 max which flattens out for higher cloud masses at < 10 Myr. A similar relation exists between cloud size at peak mass and lifetime. This simulation of the evolution of a realistic galactic cold cloud population supports the rapid formation and disruption of star-forming clouds by stellar radiation and supernovae on a timescale less than 10 Myr.  \nKey words: ISM: clouds - ISM: evolution - galaxies: dwarf - galaxies: starburst - methods: numerical', '1 INTRODUCTION': "There is clear observational evidence that galactic star formation is closely connected to the galactic gas content and in particular to molecular gas (see e. g. Kennicutt 1998; Bigiel et al. 2008a; Leroy et al. 2008; Kennicutt & Evans 2012). Most molecular gas is organised in giant clouds with sizes of a few to ∼ 100 pc and masses of 10 to ∼ 10 6 M ⊙ (see e.g. Heyer & Dame 2015) which are observed to follow several scaling relations. The cloud mass functions (CMFs) follow a power law with slope 𝛼 , i.e. 𝑑𝑁 / 𝑑𝑀 ∝ 𝑀 𝛼 . The slope varies between 𝛼 = -1 . 5 and 𝛼 = -2 . 1 (Rosolowsky 2005) with typical values of 𝛼 ∼ -1 . 8 (Heyer et al. 2001) for the Milky Way. The formation of dense cold clouds and star clusters therein can be triggered and enhanced in extreme environments (Elmegreen & Efremov 1997) such as starburst galaxies (Holtzman et al. 1992;  \nLeroy et al. 2018), interacting massive galaxies (Whitmore et al. 2010; Herrera et al. 2012; Johnson et al. 2015) or merging dwarf galaxies (Johnson et al. 2000; Kepley et al. 2016; Ochsendorf et al. 2017).  \nMolecular clouds also follow relations between their kinematic and structural properties (Larson 1981). The so-called 'Larson relations' connect the cloud velocity dispersion along the line of sight, 𝜎 v , (from 12 CO line-width observations) to cloud size, 𝑅 c as 𝜎 𝑣 [ kms -1 ] ∝ 𝑅 0 . 38 c [ pc ] , the velocity dispersion to cloud mass as 𝜎 v [ kms -1 ] ∝ 𝑀 0 . 2 c , and the cloud density 𝑛 c to cloud size as 𝑛 c [ cm -3 ] ∝ 𝑅 -1 . 1 c [ pc ] . More recent estimates of the slope of the first Larson relation indicate values close to ∼ 0 . 5 (see e.g. Solomon et al. 1987; Heyer et al. 2001; Heyer & Dame 2015). However, there seem to be considerable variations in the observed slope ranging between 0.15 - 0.91, depending on the respective study (e.g. Bolatto et al. 2008a; Imara et al. 2020; Duarte-Cabral et al. 2021). The third Larson relation implies approximately constant cloud surface", '2 C. M. Fotopoulou et al.': 'densities (see e.g. Lombardi et al. 2010). Several other studies indicate surface densities varying from ∼ 10 M ⊙ pc -2 to 1000 M ⊙ pc -2 (Heyer et al. 2009; Hughes et al. 2013; Duarte-Cabral et al. 2021). It seems, therefore, unclear how universal observed molecular/cold cloud scaling relations behave.  \nThe very existence of a physical scaling relation between the density and the size of molecular clouds has been challenged by simulations (Mac Low 1999; Ballesteros-Paredes & Mac Low 2002; MacLow&Klessen2004),attributing its existence to observational artefacts, namely the unavoidable existence of a minimum threshold in column density. On the other hand, the size-velocity dispersion relation is not disputed since it is reliably reproduced by both observations and simulations. The debate lies in whether its existence confirms the presence of virial equilibrium or even energy equipartition, as argued by previous studies (Ballesteros-Paredes & Mac Low 2002; Mac Low & Klessen 2004).  \nThe lifetimes of molecular clouds are even more difficult to determine observationally (see Heyer & Dame 2015, for a discussion). Already early studies have estimated either short lifetimes ∼ 10 Myr or long lifetimes ∼ 100 Myr (Scoville & Solomon 1975; Scoville & Hersh 1979). More recent work seems to favour maximum lifetimes of several 10 Myr (e.g. Blitz et al. 2007; Meidt et al. 2015; Kruĳssen et al. 2019a). Knowledge about the cloud lifetimes is of fundamental astrophysical interest as it can reveal information about the formation (gravity, turbulence, galactic gas thermodynamics) and destruction (galactic shear, star formation and stellar feedback) mechanisms of the cold and dense phase of the interstellar medium (ISM) which is the site of star formation in all galaxies.  \nAlarge number of numerical studies have investigated the formation and disruption of the cold cloud population in massive, Milky Way-like, galaxies. The different simulation studies have been carried out with all major available (magneto-)hydrodynamical galaxy evolution codes (see e.g. Naab & Ostriker 2017). The studies have used different mass and spatial resolutions as well as different sub-grid models for ISM heating and cooling processes, chemistry, star formation and stellar feedback when included. In some studies the stellar disk (with/without spiral structure) and/or the galactic halo potential is analytic (Tasker & Tan 2009; Dobbs et al. 2011a; Smith et al. 2020), in some the halo is represented by particles and can dynamically respond (e.g. Hopkins et al. 2012). A general conclusion from these studies is that the formation and disruption of cold clouds is complex and deeply connected to the larger scale galactic ISM properties (see e.g. Chevance et al. 2020a, for a review). Clouds simulated in the galactic environment can have highly irregular structure (e.g. Benincasa et al. 2013; Smith et al. 2020; Treß et al. 2021) and in particular massive clouds can rapidly grow by cloud collisions (Li et al. 2018; Skarbinski et al. 2023). Isolated cloud formation/disruption studies starting with more or less idealised initial conditions are well suited for understanding the cloud internal processes but cannot deliver the full picture.  \nMost simulations, even without stellar feedback, indicated that cold clouds form with a declining CMF above the resolution limit with power law slopes of 𝛼 ∼ -1 . 8 (e.g. Tasker & Tan 2009; Dobbs et al. 2011a; Hopkins et al. 2012; Grisdale et al. 2018; Li et al. 2020; Jeffreson et al. 2021). It has been noticed (Hopkins et al. 2012) that this slope is similar to analytical expectations from gravitational collapse (Press & Schechter 1974; Bond et al. 1991). However, other factors could affect the mass function slope, e.g. Williamson et al. (2014); Li et al. (2020) find a slope increase for models with stronger stellar feedback.  \nMany studies indicate that simulated cold and dense clouds are marginally gravitationally bound or unbound (e.g. Hopkins et al.  \n2012; Treß et al. 2021) and only the densest regions might collapse (Dobbs et al. 2011c). Galaxy interactions are reported to produce clouds which are less bound (Pettitt et al. 2018). Whether the simulated clouds "automatically" follow the Larson relations remains uncertain. In some high-resolution galaxy evolution simulations, the forming clouds seem to follow the linewidth-size relation (e.g. Hopkins et al. 2012; Jeffreson et al. 2021). Nickerson et al. (2019) use one of the most complete modelling approaches including radiation transfer and molecular chemistry. Their linewidth-size relation appears to be steeper than observed (see also Treß et al. 2021, for the same effect). Grisdale et al. (2018) argue that the slope of the linewidth-size relation depends on the implementation of feedback (steeper without feedback) and on whether the analysis is performed in projection (2D) or 3D.  \nAll galaxy scale simulation studies indicate short cloud lifetimes of several Myr not exceeding a few 10 Myr (e.g Tasker 2011; Hopkins et al. 2012; Dobbs & Pringle 2013; Benincasa et al. 2020; Jeffreson et al. 2021). For example, Jeffreson et al. (2021) studied the time evolution of molecular clouds in simulated Milky Way-like galaxies and found self-gravitating clouds with lifetimes that scale with their sizes as 𝜏 life ∝ ℓ -0 . 3 . In this context, it is worth noting that most if not all simulation based cold (molecular) cloud studies suffer from two major uncertainties - the cloud identification algorithm (e.g. Nickerson et al. 2019) and the (sub-grid) modelling of the ISM (e.g. Li et al. 2020), in particular the modelling of molecular gas. For cloud identification different methods have been used like group-finding algorithms (Dobbs 2008; Hopkins et al. 2012; Nickerson et al. 2019), coherent peaks of high-density gas (Tasker & Tan 2009) or dendrogram-based algorithms (Li et al. 2020; Treß et al. 2021; Jeffreson et al. 2021). To identify \'molecular clouds\' different studies have used density thresholds (e.g. Tasker & Tan 2009), various post-processing methods for molecular hydrogen (e.g. Tasker & Tan 2009), or direct chemical networks (e.g. Treß et al. 2021). The simulations presented here also use a chemical network to follow molecular gas directly. However, we use a clump finding base algorithm as we are simulating a low metallicity system where the main reservoir for star formation is the cold and dense neutral gas phase (Hu et al. 2016).  \nIn this work, we investigate the formation, the evolution and the properties of simulated cold clouds. For this purpose we analyse the star-forming cold dense gas clouds emerging after a massive starburst caused by the merger of two gas-rich and metal-poor dwarf galaxies in the smoothed particle hydrodynamics (SPH) simulation carried out by Lahén et al. (2019, 2020a,b). The simulation is part of the /g.pc/r.pc/i.pc/f.pc/f.pc/i.pc/n.pc project 1 (Galaxy Realizations Including Feedback From INdividual massive stars), which aims to resolve current challenges in galaxy formation (see e.g. review by Naab & Ostriker 2017) implementing a detailed model for the ISM (Hu et al. 2014, 2016, 2017; Lahén et al. 2019, 2020a,b; Steinwandel et al. 2020; Bisbas et al. 2022; Hislop et al. 2022; Szakacs et al. 2022). The high-resolution simulations have solar mass gas phase resolution and sub-parsec spatial resolution and include a non-equilibrium chemical network/heating and cooling, individually realised massive stars and stellar lifetimes as well as a model for HII regions and resolved supernova explosions. We focus especially on the evolution of the emerging cold clouds and the comparison to observed scaling relations.  \nThis paper is structured as follows. In Sec. 2, we introduce  \nthe underlying simulation with its numerical realisation and initial conditions, as well as the developed analysis tools. A description of the cold interstellar medium (ISM) in our numerical experiment follows in Sec. 3. We analyse the mass and size functions of the detected dense cold gas clouds in Sec. 5 and test them for the aforementioned scaling relations in Sec. 5. The life cycle of the cold clouds is studied in Sec. 6. We discuss our findings in context, as well as possible caveats of our methods, in Sec. 7 and conclude this work with a summary in Sec. 8.', '2 SIMULATION AND METHODS': 'The dwarf galaxy merger simulation analysed in this paper is presented in Lahén et al. (2019) as part of the G/r.pc/i.pc/f.pc/f.pc/i.pc/n.pc project 1 . It utilises the SPHG/a.pc/l.pc implementation (Hu et al. 2014, 2016, 2017) of the smoothed particle hydrodynamics (SPH) code /g.pc/a.pc/d.pc/g.pc/e.pc/t.pc-3 (Springel 2005) and accounts for artificial viscosity and thermal conduction as described in Hu et al. (2014). For full details about the simulation, we refer to Lahén et al. (2019, 2020a). In the following, we briefly describe the setup and star formation, stellar feedback, and the chemistry model.', '2.1 Initial conditions': 'The simulation computes the merger of two identical gas-rich lowmetallicity isolated dwarf galaxies (see e.g. Hu et al. 2016, 2017). Each dwarf galaxy has a dark matter halo following a Hernquist (1990) profile with a virial mass of 𝑀 vir = 2 × 10 10 M ⊙ , a virial radius of 𝑟 vir = 44 kpc, a spin parameter of 𝜆 = 0 . 03, and an equivalent Navarro et al. (1997) concentration parameter of 𝑐 = 10. The embedded rotationally supported gas disk has a mass of 𝑀 gas = 4 × 10 7 M ⊙ , a scale radius of 𝑟 gas = 1 . 46 kpc and an initial metallicity of 𝑍 gas ≈ 0 . 1 Z ⊙ , which is a typical value for low-mass dwarf galaxies. The stellar disk has a mass of 𝑀 ∗ = 2 × 10 7 M ⊙ , a scale radius of 𝑟 ∗ = 0 . 73 kpc and a scale height of ℎ ∗ = 0 . 35 kpc. The dwarf galaxies are set up on an initially parabolic orbit with a pericentric distance of 𝑑 peri = 1 . 46 kpc and separation of 𝑑 init = 5 kpc.  \nThe mass resolution is 𝑚 DM = 7 × 10 3 M ⊙ for the dark matter particles and 𝑚 baryonic = 4 M ⊙ for the baryonic matter (gas and stars) resulting in 5 × 10 6 dark matter, 10 7 gas and 4 × 10 6 stellar particles per dwarf galaxy. The gravitational softening length is 𝜖 DM = 62 pc and 𝜖 baryonic = 0 . 1 pc for the dark matter and baryonic matter, respectively.', '2.2 Star formation and stellar feedback': 'A gas particle 𝑖 becomes eligible for star formation once its associated Jeans mass,  \n𝑀 J , i = 𝜋 5 / 2 𝑐 3 s 6 𝐺 3 / 2 𝜌 1 / 2 i , (1)  \ncrosses a threshold density of 8 SPH kernel masses ( 𝑀 𝐽,𝑖 < 3200 M ⊙ ) (Hu et al. 2017), where 𝑐 s is the local sound speed, 𝐺 is the gravitational constant, and 𝜌 i is the local gas density. Each SPHkernel contains approximately 100 neighbour particles. For gas particles with Jeans masses between 0.5 and 8 SPH kernel masses, star formation is modelled following Schmidt (1959), such that  \n𝑑𝜌 ∗ 𝑑𝑡 = 𝜖 ff 𝜌 gas 𝑡 ff , (2)  \nwhere 𝜌 ∗ and 𝜌 gas are the stellar and gas densities, respectively, and 𝑡 ff is the free-fall time of the star-forming gas. For the efficiency per free-fall time, 𝜖 ff , we assume a typical value of 0.02 (see Hislop et al. 2022, for a discussion on the impact of different efficiencies). For Jeans masses below 0.5 kernel masses ( 𝑀 𝐽,𝑖 < 200 M ⊙ ), corresponding to number densities of 𝑛 H > 10 3 . 5 cm -3 at temperatures between 𝑇 = 10 -100 K, we enforce instantaneous star formation. Once a gas particle is labelled as star-forming, we sample stars from a Kroupa (2001) IMF in the mass range of 0.08 M ⊙ to 50 M ⊙ . Stars with masses exceeding 4 M ⊙ are represented individually as distinct star particles, while stars with lower masses are aggregated into stellar particles of approximately 4 M ⊙ . Each massive star ( 𝑀 > 4 M ⊙ ) injects stellar feedback and contributes via farultraviolet (FUV) emission to the interstellar radiation field (ISRF). The FUV from the stars is integrated from the B/a.pcS/e.pcL library (Westera et al. 2002) and propagated into the medium considering dust shielding calculated with T/r.pc/e.pc/e.pcC/o.pc/l.pc (Clark et al. 2012).  \nAdditionally, young massive stars ( > 8 M ⊙ ) also contribute with hydrogen ionising radiation and the resulting HII regions are modelled as Strömgren spheres. The gas within these regions is heated to 10 4 Kandthe size of the regions is iterated until ionisation equilibrium is reached. All details about the modelling of FUV and hydrogen ionising radiation are presented in Hu et al. (2017).  \nMassive stars explode as core-collapse supernovae using the lifetimes of Georgy et al. (2013). We realise Type II supernovae (SNII) by injecting thermal energy equivalent to 𝐸 SN = 10 51 erg into the surrounding gas particles alongside with ejecta of mass and metals obtained from Chieffi & Limongi (2004). The environmental densities at SN explosion sites are pre-processed by the ionising radiation and previous SN explosions which enables a stronger coupling of the injected SN energy to the gas dynamics (see e.g. Walch & Naab 2015; Hu et al. 2016, 2017; Lahén et al. 2020b; Rathjen et al. 2021; Smith et al. 2021; Hislop et al. 2022). We also account for stellar winds from asymptotic giant branch (AGB) stars with metal yields from Karakas (2010).', '2.3 Chemistry': 'The code models non-equilibrium chemistry with heating and cooling while explicitly following the species H, H + , H 2 , C, C + , O, CO and e -, based on the chemical network by Nelson & Langer (1997) and adapted by Glover & Mac Low (2007). We account for photoelectric heating on dust grains, cosmic ray ionisation and photo-dissociation of H 2 , as well as fine-structure and metal cooling. However, in the high-temperature regime ( 𝑇 > 3 × 10 4 K), we switch to metallicity-dependent equilibrium cooling tables by Wiersma et al. (2009) under the assumption of an optically thin ISM subjected to an ionising UV field from the galactic background (Haardt & Madau 1996).', '2.4 Cold cloud identification': 'In metal rich galaxies, stars typically form within molecular clouds (see e.g. Lada & Lada 2003) that consist mainly of H 2 . Observational data from nearby galaxies have revealed a correlation between the surface density of the SFR ( Σ SFR ) and the molecular gas surface density ( Σ H 2 ) (Bigiel et al. 2008b, 2011). Dwarf galaxies, however, are typically metal-poor with a shallow H 2 content (Tremonti et al. 2004; Hunter et al. 2012). For many star-forming low mass dwarf galaxies very little or no H 2 is detected (Schruba et al. 2012; Michałowski et al. 2015). Therefore, the larger scale reservoir for  \n2.0  \n1  \nFigure 1. Gas surface density maps (top panels) and mass-weighted gas temperature maps (bottom panels) of the central 3 kpc of the dwarf galaxy merger at three different stages: first encounter at 80 Myr (left), the second encounter and peak of star formation (middle panel), and at the new equilibrium configuration at 300 Myr (right panel). Dense (yellow) and cold (blue) gas is structured in clouds and filaments. The hot (red, bottom panel) gas is generated by supernova explosions.  \n<!-- image -->  \nstar formation in dwarf galaxies might be cold atomic gas, which is also supported by the results from numerical simulations (Glover & Clark 2012; Krumholz 2012; Hu et al. 2014, 2016). In our simulated system only small amounts of H 2 form, which is typical for metalpoor systems, such as dwarf galaxies (see e.g. Kirby et al. 2013). To be more specific, the formation of H 2 peaks around the time of the major starburst, i.e. between 150-180 Myr and even during this peak period, H 2 accounts for - on average - 0.04 per cent of the total gas mass of the system. Therefore, we do not use H 2 to identify the sites of star forming clouds, but instead we use gas density and temperature criteria, i.e. gas which is sufficiently dense and cold to form stars.  \nTo identify the cold clouds in each simulation snapshot we use the friends-of-friends algorithm (FoF) implemented in the Python based /p.pc/y.pc/g.pc/a.pc/d.pc analysis package (Röttgers et al. 2020). For the cold and dense ISM gas phase we only consider gas above a defined density and below a defined temperature threshold of 𝑛 gas > 100 cm -3 and 𝑇 gas < 300 K (e.g. Walch et al. 2015; Rathjen et al. 2021). We varied the FoF linking length over an order of magnitude and found that a linking length of 1 pc results in stable cloud masses (see Appendix A).  \nWe only consider clouds (FoF groups) more massive than 100 SPH gas particles (approximately 400 𝑀 ⊙ ). We track the individual clouds through the simulation snapshots (separated by 1 Myr) via backward-tracing. A cloud in any given simulation snapshot is considered to be the descendant of a cloud in a previous snapshot if it contains at least 20 per cent of the particles of the parental cloud.  \nWetested a wide range of thresholds between 10 and 90 per cent and found that 20 per cent gives the most stable results, mapping both the growth and the destruction phases of the cold cloud life-cycles.', '3 EVOLUTION OF THE COLD ISM': 'In Fig. 1, we show an overview of the gas distribution (top panels) and the density-weighted temperature distribution (bottom panels) of the central 3 kpc of the dwarf galaxy merger at the first encounter at 80 Myr, at the peak of star formation during the second encounter after 170 Myr and at 300 Myr, when the system is more relaxed. At the first encounter (left panels) the gas gets compressed and reaches surface densities in excess of Σ gas ≳ 1000 M ⊙ pc -2 (top panel), more than two orders of magnitude higher than for isolated systems (Hu et al. 2016, 2017; Hislop et al. 2022). At these high surface densities, cold clouds (the blue regions in the bottom panel) are already forming in the turbulent ISM. At the second encounter (middle panels), the merger has the highest SFR and the effect of supernova feedback is visible in the form of filaments and shells at high gas surface density and SN driven bubbles at low surface densities. These bubbles are filled by hot gas (red-coloured, middle bottom panel) generated by SN explosions. The central super-bubble originates from a massive star cluster as discussed in Lahén et al. (2019, 2020a). Gas surface densities here can reach extreme values of Σ gas ≳ 10 3 M ⊙ pc -2 when the system reaches its highest SFR (see e.g. Lahén et al. 2020a). The ISM now shows a clear multi-phase structure with clumpy cold (blue), volume-filling warm (green) and  \n0  \n1  \nFigure 2. Zoom sequence (top left, top right, bottom right) of the gas distribution colour-coded by the mass weighted temperature (first three panels) of the simulation at 300 Myr (right panel in Fig. 1) from 10 kpc down to 100 pc with one of the most massive cold clouds in the system (5 . 7 × 10 4 M ⊙ ). The black circles in the top right panel indicate the positions and half-mass radii of individual cold clouds identified with the FoF algorithm. We focus solely on the cold ISM gas phase targeting gas with 𝑛 gas > 100 cm -3 and 𝑇 gas < 300 K. In the bottom right panel, we show a zoom onto a massive cloud with its environment and the raw particle distribution identified by FoF is shown in the bottom left panel. The highly irregular shape is typical for the massive cold clouds in the simulation.  \n<!-- image -->  \nhot (red) gas covering a temperature range of 10 K ≲ 𝑇 ≲ 10 7 K. The characteristic filamentary multi-phase structure of the cold ISM in the post-merger state is depicted in the right panels of Fig. 1. After the peak starburst at ∼ 170 Myr, the remnant remains in a turbulent state, with mass-weighted velocity dispersion of the cold gas exceeding 10 km s -1 , for several dynamical times until the end of the simulation at ∼ 450 Myr.  \nWe identify coherent structures in the cold ISM with the FoF identification procedure described in Sec. 2.4. We term these structures (FoF groups) as "cold clouds" and we identify hundreds of cold clouds at every simulation snapshot. An example of the cloud identification process is shown in Fig. 2 for the merger system at 𝑡 sim = 300 Myr (e.g. the right panels in Fig. 1). The top left panel shows the large scale gas distribution colour-coded by temperature on a 10 kpc scale. The hot gas (red) is generated by SNII explosions. A zoom into the central kiloparsec makes the filamentary cold cloud structure visible (blue regions, top right panel). Each black circle represents an identified cold cloud and the maximum spatial extent of the clouds is indicated by the size of the respective circle. A zoom into the region of one of the most massive clouds in this snapshot (bottom right panel) on a 200 pc scale depicts the typical highly irregular structure of the cold gas (blue). In the bottom left panel, we show the underlying SPH gas particle distribution of one identified cloud with a mass of 𝑀 = 5 . 7 × 10 4 M ⊙ . The substructure within the clouds with several dense sub-clumps becomes visible in this representation. This zoom sequence from 10 kpc down to 50 pc showcases the large dynamic range covered by the simulation.  \nA gallery of the morphology of the eight most massive cold clouds identified at 𝑡 = 300 Myr is presented in Fig. 3. Each cloud  \nFigure 3. Visualisation of the FoF identified SPH gas particle distribution of the eight most massive cold clouds at 300 Myr. Each cloud is centred to its respective densest region in 3D projection. The massive cold clouds typically have filamentary and irregular shapes with density substructure. The colourmap indicates the surface density of each cold gas cloud, from low density areas in dark blue to high density areas in yellow.  \n<!-- image -->  \nis centred on its densest region. The visualisation of the SPH gas particle distribution using the gas surface density as a colour map highlights the variations in shape and substructure. At this snapshot the clouds can become as massive as 9 . 5 × 10 4 M ⊙ . The most massive cloud identified in the entire simulation has a mass of 3 . 9 × 10 5 M ⊙ at 156 Myr. The cloud shown in Fig. 2 is the third most massive cloud in this simulation snapshot (second row, left panel of Fig. 3). At the more massive end, the clouds are typically elongated with sizes up to several tens or hundred parsecs and show density substructure.  \nWith the results of the cloud analysis, we separate the dense and cold ( 𝑛 gas > 100 cm -3 and 𝑇 gas < 300 K) gas phase in two categories: a \'cloud component\', i.e. the gas associated in clouds which are identified using the FoF algorithm (see Sec. 2.4) and a \'diffuse component\', i.e. the gas that is not identified as part of clouds. We present the time evolution of the cold gas for the full simulation in Fig. 4. In the top panel, the cold gas (black) is divided into a cloud (blue) and a diffuse (no-cloud, red) component. The mass in cold gas is low ( ∼ 10 4 M ⊙ ) and dominated by diffuse gas  \nFigure 4. Top panel: Time evolution of the total mass in dense and cold gas ( 𝑛 gas > 100 cm -3 and 𝑇 gas < 300 K, black) split up into cold clouds (blue) and the component not bound in clouds, i.e. the diffuse gas (red). The time of the first encounter at 50 Myr, the second encounter at 170 Myr and the post-merger state at 300 Myr are indicated by dashed, dotted and dasheddotted lines, respectively. The mass in cold gas strongly increases around the peak starburst. Middle panel: The mass ratio of cold clouds and cold diffuse gas. Unity is indicated by the horizontal line. The starburst drives the cold medium, which is initially dominated by diffuse gas, into a new equilibrium where the mass of the cold gas bound in clouds is now about four times the mass of the diffuse gas. Bottom panel: Star formation rate (SFR) as a function of time. The interaction increases the SFR by about two orders of magnitude and drives the system into an extended starburst. The SFR evolution is qualitatively similar to the mass evolution of cold gas in clouds (see top and middle panel).  \n<!-- image -->  \nat the beginning of the simulation when the galaxies are entirely dominated by warm gas. At this time, cold gas constitutes only ∼ 0 . 01 per cent of the total ∼ 10 7 M ⊙ of gas in the two merging galaxies. Already after the first encounter (vertical dashed line) the gas is compressed and cools so that the cold gas mass increases by about one order of magnitude. Now the cold phase ISM becomes dominated by clouds as can also be seen in the middle panel of Fig. 4. The effect becomes even stronger after the second encounter (vertical dotted line). The cold gas mass reaches ∼ 10 6 M ⊙ and the cold ISM is now predominantly in a clumpy phase and stays in such a configuration until the end of the simulation. The SFR of the system closely follows the evolution of the cold gas mass as well as the mass fraction of gas in clouds (bottom panel in Fig. 4).  \nFigure 5. Mass function of cold ( 𝑇 < 300K) clouds at 170 Myr (top panel) and 300 Myr (bottom panel) and increasing density cuts from 𝑛 gas > 100 cm -3 to 𝑛 gas > 2000 cm -3 shown in blue, orange, green, and red, respectively. The CMFs follow a power law 𝑑𝑁 / 𝑑𝑀 ∝ 𝑀 -𝛼 for all density thresholds, indicating self-similar structure. The power law slope for these two snapshots ranges between -1 . 6 ≲ 𝛼 ≲ -1 . 8. A typical observational estimate of 𝑑𝑁 / 𝑑𝑀 ∝ 𝑀 -1 . 8 from CO is shown by the black line (Heyer et al. 2001). The error bars represent the errors in the x- and y- axis. As errors in the x-axis 𝛿𝑥 , we use the widths of the respective mass bins: 𝛿𝑥 = 𝑑𝑀 . As errors in the y-axis 𝛿𝑦 , we calculate the standard errors of the mass distribution: 𝛿𝑦 = √︁ 𝑑𝑁 / 𝑑𝑀 / 𝑑𝑀 .  \n<!-- image -->', '4 MASS AND SIZE FUNCTION': 'In this Section and Section 5 we compare the two-dimensional projected properties of the cold clouds identified from the simulation with a three-dimensional friends-of-friends finder as explained in Section 2.4. It has been indicated in the literature that twodimensional and three-dimensional identification methods from simulation might result in differences in scaling relations. We do not address this issue here as we would like to keep the cloud definition consistent for all sections of the paper. We refer the reader to Grisdale et al. (2018) for a detailed discussion on cloud identification.  \nIn Fig. 5, we present the cold CMF from the minimum mass of ∼ 400 M ⊙ to ∼ 4 × 10 5 M ⊙ at the peak of star formation at 170 Myr (top panel) and at a later phase of the system at 300 Myr (bottom panel, see Fig. 1 for the corresponding global gas distributions). We compare the CMFs for four increasing upper limits for the cold gas density. In addition to our fiducial value of 𝑛 gas = 100 cm -3 (blue) we also perform the FoF analysis for cold gas denser than 𝑛 gas = 500 cm -3 (orange), 𝑛 gas = 1000 cm -3 (green), and 𝑛 gas = 2000 cm -3 (red). In all cases, we study the cold gas, i.e. the temperature limit of 300 K remains constant. For all den-  \nFigure 6. Timeevolution of the cold CMF from pre-merger times (dark blue) towards the new post-merger equilibrium (yellow) plotted every 50 Myr. At all phases, the mass function follows a power-law distribution. For comparison, we indicate by the black line ( 𝑑𝑁 / 𝑑𝑀 ∝ 𝑀 -1 . 8 ) the observed slope from CO from (Heyer et al. 2001). The dotted line ( 𝑑𝑁 / 𝑑𝑀 ∝ 𝑀 -1 ) indicates the one cloud per mass bin threshold.  \n<!-- image -->  \nFigure 7. Time evolution of the power law slope 𝛼 of the cold CMFs shown in Fig. 6 starting at 150 Myr (i.e. shortly before the starburst). The mean value of the slope is indicated by the black horizontal dashed line ( 𝛼 mean = -1 . 78) and the scatter of 1 𝜎 (± 0 . 08 ) is indicated by the grey area. The value of the observational slope of the CO luminosity function ( Δ 𝑁 / Δ 𝐿 CO ∝ 𝐿 𝛼 CO ) is indicated by the red horizontal dashed line ( 𝛼 obs = -1 . 8) and the observational scatter ( ± 0 . 03) is indicated by the red area (Heyer et al. 2001). The mass function slope is similar to the observations at all stages of the simulation.  \n<!-- image -->  \ncold CMFs follow similar 𝑑𝑁 / 𝑑𝑀 ∝ 𝑀 𝛼 power law distributions with slopes in the range of -1 . 9 ≲ 𝛼 ≲ -1 . 5. As expected, the normalisation decreases with increasing density thresholds. The slopes of the cold CMFs are comparable to observational values 𝛼 ∼ -1 . 8 for molecular CMFs inferred from the CO luminosity observations (e.g. Heyer et al. 2001). The similar slopes for denser substructures indicate a self-similar density structure within the clouds on the spatial and density scales investigated here. Such self-similarity is typically attributed to the presence of turbulence (see e.g. Elmegreen & Falgarone 1996; Stutzki et al. 1998; Elmegreen & Scalo 2004)  \nThe time evolution of the cold CMF for the full simulation is shown in Fig. 6 in steps of 50 Myr. The time evolution of the  \n2  \nFigure 8. Differential size distribution of the cold clouds detected at 170 Myr (black) and at 300 Myr (red). The size of the clouds is represented by the projected half mass radius 𝑅 1 / 2 . The differential size distribution of the clouds is a power law distribution with index around -2.5, for both 170 Myr and 300 Myr. The vertical dotted line indicates the observational limit of Heyer et al. (2001) at 𝑅 1 / 2 ∼ 3 pc. The error bars in the x-axis represent the widths of the respective size bins: 𝛿𝑥 = 𝑑𝑅 1 / 2 . The error bars in the y-axis represent the standard errors of the size distribution: 𝛿𝑦 = √︁ 𝑑𝑁 / 𝑑𝑅 1 / 2 / 𝑑𝑅 1 / 2 .  \n<!-- image -->  \nbest-fitting slope is shown in Fig. 7. From the early pre-merger stages of the simulation (dark blue) towards the late post-merger equilibrium (yellow) the mass function approximately follows a power-law distribution 𝑑𝑁 / 𝑑𝑀 ∝ 𝑀 𝛼 at every stage. The limiting minimum slope 𝛼 = -1 for one cloud per log mass bin is indicated by the black dashed line in Fig. 6. The slope of the cold CMF is in good agreement with the observed CO luminosity function (Heyer et al. 2001) (solid black line) at all evolutionary stages (see Fig. 7).  \nFor more than 300 Myr, the best-fitting slope of the CMF varies little (see Fig. 7). The mean value is 𝛼 ∼ -1 . 8 with a scatter of ∼ 0 . 1 dex. To test the correlation between the mass function slope and normalisation and the SFR we use the Pearson correlation coefficient 𝜌 that is defined as the ratio between the covariance of two variables 𝑋, 𝑌 and the product of their standard deviations 𝜎 X , 𝜎 Y , i.e. 𝜌 ( 𝑋,𝑌 ) = cov ( 𝑋,𝑌 )/( 𝜎 X 𝜎 Y ) . The Pearson correlation coefficient is 𝜌 ( 𝛼, SFR ) = 0 . 3 suggesting that there is only a weak correlation between the slope of the CMF 𝛼 and the actual SFR of the system. The correlation coefficient for the normalisation of the CMF 𝛽 and the SFR is 𝜌 ( 𝛽, SFR ) = 0 . 2, also suggesting a weak correlation. Therefore, the star formation activity of the system does not appear to have a strong impact on cold cloud mass function. However, the mass of the most massive cloud identified in each snapshot between 150 Myr and 458 Myr correlates strongly with the SFR with 𝜌 ( 𝑀 max , SFR ) = 0 . 8. This is in agreement with a similar correlation between the most massive forming star cluster and the SFR in the same simulation (Lahén et al. 2020a).  \nIn Fig. 8, we plot the cloud size function of the cold clouds at 170 Myr (black) and 300 Myr (red). Similarly to the CMF, the cloud size function is defined as the number of clouds in a specific size range divided by that range, i.e. 𝑑𝑁 / 𝑑𝑅 1 / 2 , plotted against the size 𝑅 1 / 2 . As a measure for the cloud size, we use the projected half mass radius 𝑅 1 / 2 , which is calculated from the respective centre of each cloud. Each cloud is centred on its densest region (see Fig. 3). To determine the half-mass radius 𝑅 1 / 2 of each cloud, we first calculate the total mass of the cloud and the radial distances of each gas particle from the densest point of the cloud, on the 𝑥 -𝑦 plane.', '8 C. M. Fotopoulou et al.': 'Subsequently, we sort the array of radial distances from minimum to maximum and the corresponding array of gas particle masses following the sorting of the radial distances array. This ensures that each element of both arrays corresponds to properties of the same gas particle. By summing the gas particle masses along the mass array until the half of the total mass of cloud is reached (or at least not exceeded), we identify the index at which this occurs. The radial distance associated with this index in the sorted array represents the projected half-mass radius on the x-y plane 𝑅 1 / 2 . This process is repeated for all clouds and for the 𝑦 -𝑧 and 𝑧 -𝑥 planes, with minimal variation observed between the resulting half-mass radii in each projection. The sizes of the simulated cold clouds range from a few parsecs to ∼ 50 pc. The cloud size function follows a powerlaw distribution 𝑑𝑁 / 𝑑𝑅 1 / 2 ∝ 𝑅 𝛼 1 / 2 with index 𝛼 = -2 . 47 (± 0 . 16 ) at 170 Myr (black) and 𝛼 = -2 . 52 (± 0 . 20 ) at 300 Myr (red). These values are shallower than that of the observations of Heyer et al. (2001) who find a slope of 𝛼 obs = -3 . 2 (± 0 . 1 ) . To directly compare with the observations with a detection threshold of 3 pc (Heyer et al. 2001), we have only used clouds with 𝑅 1 / 2 > 3 pc in our analysis, but we display the full range of our sample in Fig. 8. The reason for the shallower size distribution slopes of the simulated clouds are not obvious and will have to be investigated in future studies. Here, the cold cloud structure identification, the tracer, the stellar feedback implementation, and the galactic environment could be of importance. With observations, molecular cloud sizes are also difficult to measure and not many observations exist. The observational tracer, the algorithm to determine the sizes, or the galactic environment can have an impact on the size distribution (see e.g. discussions in Heyer & Dame 2015; Dib et al. 2020). From now on we will use only clouds with 𝑅 1 / 2 > 3 pc in our analysis, unless stated otherwise and we will display a vertical dotted line at 𝑅 1 / 2 ∼ 3 pc in the respective figures to indicate the observational limit of Heyer et al. (2001).', '5 COLD CLOUD SCALING RELATIONS': "The observed cold cloud scaling relations introduced by Larson (1981) can be understood as an indication that cold clouds are selfgravitating and, possibly, in some form of virial equilibrium (see Heyer et al. 2001; Heyer et al. 2009; Heyer & Dame 2015). To determine if a system - like a cold cloud or clump - is in virial equilibrium we use the virial parameter 𝛼 vir , which is a dimensionless parameter that quantifies the ratio of thermal and kinetic energies to gravitational energy in the system. It is mathematically expressed as (Bertoldi & McKee 1992):  \n𝛼 vir ≡ 5 𝜎 2 v 𝑅 𝐺𝑀 vir , (3)  \nwhere 𝐺 is the gravitational constant, 𝑅 is the radius, 𝑀 vir is the corresponding virial mass and 𝜎 𝑣 is the average one dimensional velocity dispersion of the system. When the gravitational energy of a cloud becomes comparable to its kinetic energy, then 𝛼 vir ∼ 1. However, when the surface pressure confines the cloud and selfgravity is negligible, then the virial parameter is significantly greater than one: 𝛼 vir > 1. The virial mass 𝑀 vir estimate from Eq. 3 and the mass of a cold cloud 𝑀 cc estimated from a measured effective (projected half-mass) radius 𝑅 1 / 2 , an effective surface mass density Σ 1 / 2 (the mass within the effective radius divided by the area) are:  \n𝑀 vir = 5 𝜎 2 v 𝑅 1 / 2 𝐺 and 𝑀 cc = 𝜋 Σ 1 / 2 𝑅 2 1 / 2 , (4)  \n2  \nFigure 9. Line-of-sight velocity dispersion in the 𝑧 -direction 𝜎 z , as a function of the size of the identified cold clouds in terms of the projected halfmass-radius 𝑅 1 / 2 on the 𝑥 -𝑦 plane at 170 Myr (top panel) and 300 Myr (bottom panel). The vertical dashed line marks the 3 pc detection threshold of Heyer et al. (2001). The empty circles represent clouds with 𝑅 1 / 2 ≤ 3 pc, whereas the filled circles represent clouds with 𝑅 1 / 2 > 3 pc. The solid red line shows the best-fit to the Solomon et al. (1987) points and the shaded red area represents the scatter. The solid black line represents the best fit to the simulation data while the shaded black area denotes the standard deviation of the cloud distribution with sizes larger than 3 pc.  \n<!-- image -->  \nrespectively. For clouds in virial equilibrium 𝛼 vir = 𝑀 vir / 𝑀 cc = 1 and  \n𝜎 𝑣 = √︂ 𝜋𝐺 5 𝑅 0 . 5 1 / 2 Σ 0 . 5 1 / 2 . (5)  \nIn Fig. 9, we show the first Larson relation for the clouds identified in the simulations with our fiducial density and temperature thresholds ( 𝑛 gas > 100 cm -3 and 𝑇 gas < 300 K). The onedimensional line-of-sight velocity dispersion ( 𝜎 z ) for each cloud is computed from the mass-weighted gas particle velocities along the 𝑧 -axis of the computational domain. The half mass radius 𝑅 1 / 2 is the projected radius enclosing half of the cloud mass from the densest region of the cloud. The effective mass surface density Σ 1 / 2 is the half-mass of each cloud divided by the half-mass projected area of the cloud around its densest point (Eq. 4). The relationship between the line-of-sight velocity dispersion and the size of the identified cold clouds shown in this figure was investigated at two characteristic times: right after the starburst at 170 Myr (top panel) and at 300 Myr (bottom panel) when the system has reached post-merger equilibrium. We show the best power-law fit to our simulated clouds and the standard deviation. For comparison, we show the canonical relation derived in Solomon et al. (1987) and discussed in the review by Heyer & Dame (2015). It is represented by the solid red line and the shaded red area represents their scatter. The fitted  \nFigure 10. Line-of-sight velocity dispersion in the 𝑧 -direction 𝜎 z vs. the product of the size of the identified cold clouds in terms of the face-on projected half-mass-radius 𝑅 1 / 2 on the 𝑥 -𝑦 plane and the effective surface density Σ 1 / 2 at 170 Myr (top panel) and 300 Myr (bottom panel). The empty circles represent clouds with 𝑅 1 / 2 ≤ 3 pc, whereas the filled circles represent clouds with 𝑅 1 / 2 > 3 pc. The solid line represents Eq. 5 for virial parameter 𝛼 vir = 1. The dashed line is derived from Eq. 4 for 𝛼 vir = 3 and the dotted line for 𝛼 vir = 10.  \n<!-- image -->  \nslopes ( 𝛼 170 Myr = 0.33( ± 0.03) and 𝛼 300 Myr = 0.46( ± 0.04) are somewhat shallower than the 'canonical' slope of 0.5 reported in Solomon et al. (1987). The slopes are well within the range reported from various observations (0.15 - 0.91, see Sec. 1). This indicates that the cold clouds in our simulation model have properties which are consistent with observed cold clouds.  \nFig. 10 shows the correlation between the one-dimensional line-of-sight velocity dispersion 𝜎 z and the size of identified cold clouds in terms of their projected half-mass-radius 𝑅 1 / 2 and their effective surface density Σ 1 / 2 at 170 Myr (top panel) and 300 Myr (bottom panel). We show the virial parameters of 1, 3, and 10, respectively, assuming a virial constant of 5 (see Eq. 4 and 5). For the sake of completeness, we should note that when testing the fits on our whole sample of simulated cold clouds, the parameters of the fits (i.e. slope, normalisation of the power laws) did not change significantly for neither Fig. 9 nor Fig. 10. In addition, the results shown Fig. 9 and Fig. 10 are similar for line-of-sight projections along the x- and y-axis (i.e. 𝜎 x , 𝜎 y ).  \nIn Fig. 11 we show the average one dimensional velocity dispersion 𝜎 v of all identified clouds as a function of the product of cold cold clouds size, given by the projected half-mass-radius 𝑅 1 / 2 , and the effective surface density Σ 1 / 2 , for different gas number density thresholds: 𝑛 gas > 100 cm -3 (grey), 𝑛 gas > 500 cm -3 (blue), and 𝑛 gas > 1000 cm -3 (red). The solid line in the plot represents Eq.  \nFigure 11. Average one dimensional velocity dispersion 𝜎 v = √︃ 𝜎 2 x + 𝜎 2 y + 𝜎 2 z / √ 3 as a function of the product of the identified cold clouds size in terms of the projected half-mass-radius 𝑅 1 / 2 and the effective surface density Σ 1 / 2 at 300 Myr for different gas density thresholds: 𝑛 gas > 100 cm -3 (grey), 𝑛 gas > 500 cm -3 (blue) and 𝑛 gas > 1000 cm -3 (red). The solid line represents Eq. 5 for virial parameter 𝛼 vir = 1. A grey shaded area of 0.3 dex around the 𝛼 vir = 1 line is plotted to visually indicate the number of clouds in approximate virial equilibrium. As the total gas number density threshold increases from 𝑛 gas > 100 cm -3 to 𝑛 gas > 1000 cm -3 , the total number of identified clouds decreases, but the number of clouds in virial equilibrium (i.e. the clouds in the shaded area) increases.  \n<!-- image -->  \n5 for a virial parameter 𝛼 vir = 1 and the grey-shaded indicates the range of 0 . 5 < 𝛼 vir < 2 which might be considered close to virial equilibrium. Values of 𝛼 vir < 0 . 5 would clearly indicate collapse, while 𝛼 vir > 2 points towards expansion. As the gas density threshold increases from 𝑛 gas > 100 cm -3 to 𝑛 gas > 1000 cm -3 , the total number of identified clouds decreases. With the higher gas density threshold, the identified cold clouds are the very dense cores of the previously identified clouds at a lower threshold, without their less dense outer regions. Some clouds are already so disrupted by stellar feedback or not located in suitable environmental conditions to form such dense cores, which leads to a decrease of the total number of identified clouds. However, the relative number of clouds in virial equilibrium, which are represented by the clouds in the shaded area, increases with the increase in density threshold. The dense regions of clouds, when outer less dense regions are ignored, are closer to virial equilibrium. Using only the identified clouds with densities higher than 1000 cm -3 and in virial equilibrium, i.e. the cold clouds represented by red scatter points that fall into the shaded area in Fig. 11, we calculate the SFR. To do so, first, we calculate the free-fall time of the selected clouds using the formula (Evans et al. 2021):  \n𝑡 ff = GLYPH<18> 3 𝜋 32 𝐺𝜌 GLYPH<19> 0 . 5 , (6)  \nwhere 𝐺 is the gravitational constant and 𝜌 is the cloud's average volume density, weighted by the SPH-particles mass for each cloud. The SFR cc derived by the selected cold clouds can be calculated using the sum of the masses of the clouds divided by the free-fall times:  \nSFR cc = 𝑁 ∑︁ 𝑖 = 1 𝑀 cc 𝑡 ff , (7)  \nwhere 𝑁 is the number of the selected clouds. By performing these calculations for the selected dense clouds ( 𝑛 gas > 1000 cm -3 ) in  \nFigure 12. The scaling coefficient 𝜐 o = 𝜎 v / 𝑅 0 . 5 1 / 2 as a function of the effective surface density Σ 1 / 2 at 300 Myr for different gas density thresholds: 𝑛 gas > 100 cm -3 (grey), 𝑛 gas > 500 cm -3 (blue) and 𝑛 gas > 1000 cm -3 (red). The average line-of-sight velocity dispersion is 𝜎 v = √︃ 𝜎 2 x + 𝜎 2 y + 𝜎 2 z / √ 3 and 𝑅 1 / 2 is the projected half-mass-radius of the identified cold clouds. The solid line shows the loci of gravitationally bound clouds.  \n<!-- image -->  \nvirial equilibrium we find that SFR cc = 0.04 M ⊙ yr -1 . The actual SFR in our simulation at 300 Myr is indeed 0.04 M ⊙ yr -1 (see also the bottom panel of Fig. 4 for a quick estimate of the actual SFR at 𝑡 = 300 Myr). The fact that the SFR cc calculated from the freefall time of the selected clouds is equal to the actual SFR of the simulation indicates that the star formation in our simulation comes directly and exclusively from high-density clouds close to virial equilibrium. Alternatively, we identified all the star particles at 𝑡 = 300 Myr that were created during the last 1 Myr and identified their progenitor gas particles at the previous snapshot of our simulation, i.e. at 𝑡 = 299 Myr. In the same snapshot, we identified the cold dense clouds with 𝑛 gas > 1000 cm -3 . These are the dense regions of the cold clouds, as we explained above. This analysis showed that these dense cores include 57 per cent of the gas particles that will turn into stars at 𝑡 = 300 Myr. By performing the same analysis at 𝑡 = 299 Myr for the cold gas with the lower density threshold of 100 cm -3 , we find that around 94 per cent of the gas particles that will crate stars at the next 1 Myr belong to the identified clouds. It is evident that in our simulation star formation takes place in the cold clouds and more than half of the stars are located in the dense regions of these clouds. To gain a statistical overview, we performed the same analysis for five snapshots around the time of the merger (170 Myr) and for five snapshots around the more quiescent stage of the system (300 Myr). Our analysis showed that on average around 60 per cent of the gas particles that create stars within 1 Myr are located in the denser regions of the identified cold clouds.  \nWhile the original third Larson relation indicates that all clouds have similar surface densities, Σ 0 , this is not the case in many observed systems and in our simulations (see Sec. 1). From Eq. 5, a constant surface density would imply that the scaling coefficient 𝜐 o , G = 𝜎 v / 𝑅 0 . 5 1 / 2 = ( 𝜋𝐺 Σ 1 / 2 / 5 ) 1 / 2 , which relates the velocity dispersion and cloud size in virial equilibrium, is constant with a value of ( 𝜋𝐺 Σ 0 / 5 ) 1 / 2 . In Fig. 12 we show the scaling coefficient 𝜐 o = 𝜎 v / 𝑅 0 . 5 1 / 2 for the simulated clouds at 𝑡 = 300 Myr at the same density thresholds as Fig. 11. The simulated clouds cover more than one order of magnitude in surface density, similar to observational estimates (e.g. Heyer et al. 2009; Hughes et al. 2013; Duarte-Cabral  \nFigure 13. The mass evolution for 10 cold clouds with lifetimes larger than 10 Myr, randomly selected from our sample of traced clouds. In their growth phase the clouds gain enough mass for gravitational collapse and the beginning of star formation. In their destruction phase the energy and momentum input from newly formed stars lead to the dispersal of the clouds. Both the growth and destruction phases of the clouds are visible on the selected mass traces.  \n<!-- image -->  \net al. 2021) and they also do not seem to follow the equilibrium scaling in general. In particular low surface density clouds deviate significantly from the equilibrium expectation indicated in Fig. 12.", '6 THE COLD CLOUD LIFE CYCLE': 'We have used the method described in Sec. 2.4 to follow the tracks of 31222 individual clouds starting shortly before the peak starburst (150 Myr) until the end of our simulation (457 Myr). The cold cloud life cycle ends as soon as the cloud cannot be detected anymore. 61 per cent of the cold clouds (i.e 19123 traced cold clouds) have lifetimes longer than 1 Myr. 91 per cent of the cold clouds that are detected in two consecutive snapshots (1 Myr) have lifetimes shorter than 10 Myr and no cloud lives longer than ∼ 30 Myr. As an example, we show the mass evolution of 10 randomly chosen clouds with lifetimes longer than 10 Myr is shown in Fig. 13. The cloud mass at first identification can be relatively high as the gas has to cross the density threshold to be detected. Clouds typically grow in mass after they are first detected (e.g. the first cloud on the left in Fig. 13). Thereafter, the clouds lose mass due to star formation and stellar feedback. Individual cold clouds, however, can have more complex histories. Once the clouds are not detected any more their life-cycle is complete.  \nIn Fig. 14, we show the lifetime distribution for all traced cold clouds that survive for at least 3 Myr (top left panel) as well as the distribution in three different mass bins covering the mass range of identified clouds, i.e. in mass bins of 10 2 . 6 M ⊙ ≲ 𝑀 cl < 10 3 . 6 M ⊙ (top right panel), 10 3 . 6 M ⊙ ≲ 𝑀 cl < 10 4 . 6 M ⊙ (bottom left panel), and 10 4 . 6 M ⊙ ≲ 𝑀 cl < 10 5 . 6 M ⊙ (bottom right panel). In general, the cold clouds follow an exponential lifetime distribution with an e-folding timescale of ∼ 3 . 5 (± 0 . 4 ) Myr. In the mass range of 10 2 . 6 M ⊙ ≲ 𝑀 cl < 10 3 . 6 M ⊙ the e-folding timescale is ∼ 3 . 2 (± 0 . 4 ) Myr, in the mass range 10 3 . 6 M ⊙ ≲ 𝑀 cl < 10 4 . 6 M ⊙ the e-folding timescale is ∼ 3 . 6 (± 0 . 4 ) Myr and on the high mass end of 10 4 . 6 M ⊙ ≲ 𝑀 cl < 10 5 . 6 M ⊙ the e-folding timescale is ∼ 6 . 1 (± 0 . 9 ) Myr. There is a moderate trend for lower- and medium-  \nFigure 14. Top left : Histogram of the lifetimes of the total number of cold clouds that survived for at least 3 Myr traced in our simulation, i.e. in the mass range 10 2 . 6 ≤ 𝑀 cl [ M ⊙] < 10 5 . 6 . The e-folding timescale is 3 . 5 (± 0 . 4 ) Myr. Top right : Histogram of the cold clouds traced in the mass range 10 2 . 6 ≤ 𝑀 cl [ M ⊙] < 10 3 . 6 with e-folding timescale 3 . 2 (± 0 . 4 ) Myr. Bottom left : Histogram of the cold clouds traced in the mass range 10 3 . 6 ≤ 𝑀 cl [ M ⊙] < 10 4 . 6 with e-folding timescale 3 . 6 (± 0 . 4 ) Myr. Bottom right : Histogram of the cold clouds traced in the mass range 10 4 . 6 ≤ 𝑀 cl [ M ⊙] < 10 5 . 6 with e-folding timescale 6 . 1 (± 0 . 9 ) Myr.  \n<!-- image -->  \nsive clouds have slightly longer typical lifetimes.  \nTo quantify this and check if the trend holds, we make a finer binning of the cold cloud masses in our sample and we calculate the average lifetime of the clouds in each mass bin 𝜏 life . As a representative value of the masses, we choose the peak mass that each cloud reaches during its lifetime 𝑀 max 2 . The results of our analysis are shown in Fig. 15 and for the binning, we use 20 logarithmic bins. Approximating the slope of the power-law distributions that fit the low mass range in our sample (i.e. clouds with peak masses 𝑀 max < 10 4 M ⊙ ), we find that there is a power-law dependency between the peak mass 𝑀 max and the average lifetime 𝜏 life such that 𝜏 life ∝ 𝑀 0 . 39 (± 0 . 03 ) max . This dependency indicates that the cold  \nclouds at the lower mass end of our sample have longer lifetimes the more massive they are. For the higher mass range in our sample (i.e. clouds with peak sizes 𝑀 max > 10 4 M ⊙ ), we have that 𝜏 life ∝ 𝑀 0 . 06 (± 0 . 01 ) max 3 . This indicates that the more massive cold clouds in our sample only have slight differences between their average lifetimes which tend to stabilise around 8 Myr. This might be explained by a maximum lifetime for the most massive clouds which is regulated by the onset of SNe. We analyse this possibility further down. We also find no evidence for a change of the cloud lifetime distribution depending on the state of the merger (peak starburst vs. later evolution, see Fig. B1 in the Appendix). This indicates that the cloud lifetimes are not affected by the environment.  \nSimilar as in Fig. 15, we perform a logarithmic binning of the sizes of the traced cold clouds. As a measure of size, we use the maximum distance that a SPH gas particle belonging to the cloud has from the centre of mass of the cloud. As a representative size value for the binning we chose the size each cloud has when it reaches its peak mass 𝑅 M max and we calculate the average lifetime  \nFigure 15. Peak mass of the traced clouds during their lifetime 𝑀 max binned in 20 logarithmic bins vs. the average lifetime of the clouds 𝜏 life in each bin. Each scatter point represents each bin. The error bars in the 𝑥 -axis represent the size of each bin, while in the 𝑦 -axis they stand for the standard error of the average value of the lifetimes in each bin. The solid lines indicate the indexes of the power-law distributions that represent our sample. For clouds with peak masses 𝑀 max < 10 4 M ⊙ we have 𝜏 life ∝ 𝑀 0 . 39 (± 0 . 03 ) max . For clouds with peak sizes 𝑀 max > 10 4 M ⊙ we have 𝜏 life ∝ 𝑀 0 . 06 (± 0 . 01 ) max . The trend is that more massive clouds have longer lifetimes and for clouds with peak masses larger than 10 4 M ⊙ the lifetimes tend to stabilise around 8 Myr.  \n<!-- image -->  \nof the clouds in each size bin 𝜏 life (we work with 20 logarithmic bins as in Fig. 15). The results are shown in Fig. 16. Each scatter point represents each size bin. The error in the 𝑥 -axis is represented by the horizontal error bars which are drawn based on the width of each logarithmic bin. For the error in the 𝑦 -axis we use the standard error of the average value of the lifetimes 𝜏 life in each bin which is visualised by the vertical error bars. As in Fig. 15, we approximate the slope of the power-law distributions that fit the smaller size range in our sample (i.e. clouds with sizes at peak mass 𝑅 M max < 20 pc) and we find a power-law dependency between the size at peak mass 𝑅 M max and the average lifetime 𝜏 life , i.e. 𝜏 life ∝ 𝑅 0 . 42 (± 0 . 06 ) M max . This dependency indicates that in the small-size-end of our sample the larger the clouds the longer they live. For the larger cold clouds in our sample (that is clouds with sizes at peak mass 𝑅 M max > 20 pc) we find a shallower 𝜏 life ∝ 𝑅 0 . 13 (± 0 . 05 ) M .  \nmax In Fig. 17 we plot the peak mass, 𝑀 max , against the corresponding size, 𝑅 M max , for each traced cold cloud. The color indicates the number of clouds in each pixel. The cloud peak size scales very well with the square root of its mass, following the relation 𝑅 M max ∝ 𝑀 0 . 527 (± 0 . 002 ) max . There is no obvious change in the masssize relation at the mass and size scales ( ∼ 10 4 𝑀 ⊙ , ∼ 20 𝑝𝑐 ) at which the lifetime slopes are changing (see Figs. 15 and 16).  \nWeinspect how many SNe explode in each cloud and correlate this number with the clouds peak mass and lifetime in Fig. 18. We find a moderate correlation ( 𝜌 = 0 . 53) for the 𝑀 max -𝑁 SN relation with a sub-linear power-law slope 𝛼 = 0 . 58 (± 0 . 04 ) . The 𝜏 life -𝑁 SN relation is not correlated ( 𝜌 = 0 . 12) and follows a flat distribution with 𝛼 = 0 . 20 (± 0 . 08 ) . The majority of clouds withstand not more than 1-2 SN explosions before they disperse (see accumulation of 𝑁 clouds at the bottom of both panels in Fig. 18). The moderate correlation between 𝑁 SN and 𝑀 max with the sublinear power-law slope indicates that more massive clouds can host more SNe. This is understandable with the fact that the SNe tend  \nFigure 16. Size at peak mass 𝑅 Mmax of the traced clouds during their lifetime binned in 20 logarithmic bins vs. the average lifetime of the clouds 𝜏 life in each bin. Each scatter point represents each bin. The error bars in the 𝑥 -axis indicate the size of the bin and in the 𝑦 -axis they represent the standard error of the average value of the lifetimes in each bin. The solid lines indicate the indexes of the power-law distributions that represent our sample. For clouds with sizes at their peak mass 𝑅 Mmax < 20 pc we have 𝜏 life ∝ 𝑅 0 . 42 (± 0 . 06 ) Mmax . For clouds with sizes at their peak mass 𝑅 Mmax > 20 pc we have 𝜏 life ∝ 𝑅 0 . 13 (± 0 . 05 ) Mmax . The trend is that smaller clouds have shorter lifetimes whereas for clouds with sizes larger than 20 pc the lifetimes tend to stabilise around 8 Myr.  \n<!-- image -->  \nFigure 17. Peak mass 𝑀 max of each cold cloud during its lifetime vs. its size at that given point of its life 𝑅 Mmax . The colour indicates the number of clouds in the respective region. The best fitting mass size relation, 𝑅 Mmax ∝ 𝑀 0 . 527 (± 0 . 002 ) max , is indicated.  \n<!-- image -->  \nto be very closely clustered in space and time. The more massive clouds host a bigger stellar population. Those stars, however, form at approximately the same time and multiple SNe can therefore explode in very short succession while the cloud is being dispersed. This finding is supported by the fact that there is no correlation between 𝜏 life and 𝑁 SN . The average time between the first and the last SNe in a cloud for clouds with at least 2 SNe occurring in them is 𝑑𝑡 avg SN = 2 . 26 Myr, and the median time difference is 𝑑𝑡 med SN = 1 Myr. In conclusion, the lifetime of the clouds is regulated by the onset of SNe. Multiple SNe, closely clustered in space and time, can occur within a cloud. More massive clouds have a slight tendency to host more SNe, because they also host a bigger stellar population.  \nFigure 18. Number of SNe, 𝑁 SN , which explode within a cold cloud as a function of the cloud peak mass, 𝑀 max , and lifetime, 𝜏 life . We indicate the slope, 𝛼 , of power-law fits as well as the Pearson correlation coefficients, 𝜌 .  \n<!-- image -->', '7 DISCUSSION': 'As we mentioned in Sec. 1 the three Larson relations are not independent but algebraically bound, i.e. if any two are valid the third one is valid, too. In Sec. 5, using this bond and equating the mass of the clouds to the virial mass we derived the cohesive Eq. 5 which we use in Fig. 10 to check the virialisation of the identified cold clouds. Observed molecular clouds in the Milky Way disk typically fall in between a virial parameter 𝛼 vir of 1 and 3 (Heyer et al. 2009; Heyer & Dame 2015) and the trend is similar for molecular clouds in nearby galaxies (Bolatto et al. 2008b; Leroy et al. 2015). However, in the Milky Way centre the values of 𝛼 vir are between 3 and 10 (Oka et al. 2001). Our simulated cold clouds have a virial parameter of 𝛼 vir ∼ 10 (see Fig. 10) and they follow the first Larson relation, but not Eq. 5. This means that they still follow a size-velocity dispersion relationship, but one that is independent of the surface density introduced in Eq. 5. By examining also Figs. 11 and 12, which provide a succinct illustration of the Larson scaling relations, we can infer that the complete sample of our simulated clouds does not satisfy the condition of virial equilibrium. However, when we increase the gas density threshold to 1000 cm -3 (red scatter points in Figs. 11 and 12), we notice that only a fraction of dense substructures (13 per cent, see Fig. 11) are in approximate virial equilibrium, indicating that denser regions within our identified cold clouds are gravitationally bound. This finding suggests that the clouds contain dense virialised cores surrounded by diffuse, lower density cold gas that forms a sort of cocoon around the cores. Our result is in agreement with the observational results of Evans et al. (2021), who confirmed - among others - that the value of the virial parameter varies in general with the surface density and the mass of the bound structures. Calculating the free-fall time of the dense cores that are in approximate virial equilibrium (see Eq. 6) and from it the SFR  \n(see Eq. 7), we find that this SFR is equal to the actual SFR of our simulation. This suggests that the dense cores within the cold clouds are the sites of star formation in our system.  \nPrevious studies in the literature suggest that clouds similar to ours, i.e. clouds which are not gravitationally bound but follow the linewidth-size relation exist if they are bound by external pressure (Elmegreen 1989; Bertoldi & McKee 1992; McKee & Tan 2003; Field et al. 2011; McKee et al. 2010; Hennebelle & Falgarone 2012; Heyer & Dame 2015). In more detail, Bertoldi & McKee (1992) showed that in giant molecular clouds (GMCs), smaller substructures with high virial parameter values ( 𝛼 vir ≫ 1) are bound by the external pressure of the surrounding interstellar medium. Some of these clumps contain cores which are in gravitational virial equilibrium, but these cores contain only a small fraction of the mass (Myers 1985). Another idea explored in the literature is that such clouds with high virial parameters can still follow the size-velocity dispersion relation if they have short lifetimes compared to their dynamical timescales (Field et al. 2011; Heyer & Dame 2015). Most of our simulated clouds live for a couple of Myr with an average e-folding timescale of 3.5 Myr (see Fig. 14), so this theory could describe the regime of our clouds and the reason of the discrepancy between non-virialisation and adherence to the first Larson relation.  \nThe variance between 𝜐 o = 𝜎 v / 𝑅 0 . 5 1 / 2 with Σ 0 . 5 1 / 2 that we encounter in Fig. 12 agrees with the results of Heyer et al. (2009) and suggests a deviation from strict adherence to the Larson relations. Heyer et al. (2009) attributes the displacement 𝜐 o = 𝜎 v / 𝑅 0 . 5 1 / 2 with  \nΣ 0 . 5 1 / 2 to the systematic underestimation of the surface densities of the clouds due to observational constraints (more specifically they refer to their values as \'lower limits"). Dobbs et al. (2011b) propose that the displacement is because the clouds have too large 𝜎 v to satisfy virial equilibrium and they resort to arguments of external pressure confinement of the gravitationally unbound clouds. Ballesteros-Paredes et al. (2011) argue that the displacement of 𝜐 o = 𝜎 v / 𝑅 0 . 5 1 / 2 with Σ 0 . 5 1 / 2 is not due to turbulence, but is rather caused by internal motions of the gas in the clouds undergoing gravitational collapse as a whole, while experiencing individual gravitational collapse in their cores.  \nFollowing the identification of the cold clouds in our simulation and the detailed study of their physical properties, we proceeded to track their life cycles to construct a comprehensive picture of their lifespan. Cold clouds are destroyed by stellar feedback thus having lifetimes in the range of 10-30 Myr and within ∼ 5 Myr after the emergence of very massive stars (see e.g. Kruĳssen et al. 2019b; Chevance et al. 2020b, 2022). The simulations of Milky Way-like systems by Jeffreson et al. (2023) put the average lifetime of molecular clouds at ∼ 25 Myr independent from the typical chemical lifetime of H 2 , which they found to be only 4 Myr and independent from the lifetime of the host cloud. Our simulated clouds live up to 28 Myr with most clouds in our sample being shortlived and surviving for only a couple Myr and ∼ 8 Myr after the emergence of very massive stars. This difference could be attributed to the more massive galaxies that Jeffreson et al. (2023) simulate compared to our simulation of a dwarf merger system.  \nIn Jeffreson et al. (2021) it is supported that the distribution of the number of clouds has exponential form provided that the rate of cloud formation and the lifespan of molecular clouds are time-invariant (e.g. if the total SFR remains relatively stable over time and the population of clouds does not change significantly over time, then this hypothesis is valid). Indeed, the distribution of the number of clouds in our sample has exponential form, with an average e-folding timescale 𝜏 life =3.5 Myr (see Fig. 14). Calculating', '14': 'the e-folding timescales for different mass ranges, we find that on the larger mass end (10 4 . 6 ≤ 𝑀 cl [ M ⊙ ] < 10 5 . 6 ) the average efolding timescale rises to 6.1 Myr (bottom right panel of Fig. 14). Exploring this correlation between the mass of the cold clouds and their lifetimes further, we find a trend: more massive clouds have longer lifetimes (see Fig. 15). In more detail, the lifetime 𝜏 life of clouds with peak masses ( 𝑀 max ) below 10 4 M ⊙ scales with 𝑀 max as: 𝜏 life ∝ 𝑀 0 . 39 max , whereas the lifetime of clouds with 𝑀 max > 10 4 M ⊙ tends to stabilise around 8 Myr. This happens because in our simulation the delay times of SNe are at around 8 Myr, therefore this is the timescale that even the most massive clouds will be disrupted not being able to resist to the explosive feedback of the SNe.  \nWe perform a similar analysis for the cloud size 𝑅 M max , which is the maximum distance of a SPH gas particle of the cloud from the centre of mass of the cloud at the time that the cloud has reached its peak mass. We find that the average lifetime for clouds with sizes below 20 pc correlates with the size as: 𝜏 life ∝ 𝑅 0 . 42 M max , whereas for clouds with sizes larger than 20 pc the correlation is 𝜏 life ∝ 𝑅 0 . 13 M max and their lifetime stabilises around 8 Myr. Jeffreson et al. (2021), using simulations of Milky Way-like isolated disc galaxies, find an anti-correlation between the lifetime and the size ℓ of their traced molecular clouds below the disc scale-height: 𝜏 life ∝ ℓ -0 . 3 . The disagreement is possible to lie in the differences in the tracing methods we use: our method is three-dimensional and we cannot account for cloud mergers with our tracing algorithm, whereas the Jeffreson et al. (2021) method is two-dimensional and accounting for cloud mergers. Both studies use similar star formation efficiencies 𝜖 ff : Jeffreson et al. (2021) use 𝜖 ff of 1 per cent and in our simulation 𝜖 ff = 2 per cent is used. Furthermore, there is a difference in the simulated systems: Jeffreson et al. (2021) simulate more massive and isolated disc galaxies, whereas we simulate a dwarf starburst which creates cold clouds very efficiently after the major merger.', '8 CONCLUSIONS': 'In this work, we analysed the cold and dense interstellar medium of a low metallicity starburst generated through a high-resolution hydrodynamical simulation of a gas-rich dwarf galaxy merger. The simulations resolve the multi-phase interstellar medium at 4 M ⊙ gas phase mass resolution and ∼ 0.1 pc spatial resolution. The properties of the ISM are followed with a non-equilibrium chemical heating/cooling network at temperatures below 10 4 K. Massive stars are sampled individually and interact with the ISM through a model implementation for the formation of HII regions and mostly resolved supernova explosions. Based on the definitions and methods presented in the paper, we identified and traced the formation and disruption of cold clouds. We summarise our main findings below and emphasise the remaining uncertainties in a consistent comparison with observations:  \n- · The cold ISM of the simulation evolves from a sub-dominant component in the warm gas-dominated isolated galaxy ISM to a highly filamentary dominant and clumpy ISM component in the merger-triggered starburst. We identified hundreds of cold clouds with masses of 10 2 . 6 - 10 5 . 6 𝑀 ⊙ following a power law mass function, 𝑑𝑁 / 𝑑𝑀 ∝ 𝑀 𝛼 with 𝛼 = -1 . 78 (± 0 . 08 ) and realistic size distributions. We find no statistical evidence for a change in mass function slope with the post-merger SFR of the system. The mass of the most massive cloud formed at any given time, however, correlates with the SFR of the system.\n- · The identified cold clouds are highly structured. The slope of the CMF remains roughly constant independent of the assumed density threshold. This indicates a self-similar density substructure.\n- · The identified clouds have virial parameters 𝛼 vir > 1 and approximately follow the Larson relations. The effective surface densities of the clouds is not constant. It varies from a few to a few hundred solar masses per square parsec.\n- · All cold clouds are short-lived. Their lifetime distributions are approximately exponential with e-folding timescales in the range from 3.2 Myr for lower-mass clouds up to 6.1 Myr for highermass clouds. The cloud lifetimes do not seem to be affected by the environment, i.e. the actual state of the ongoing merger.\n- · The lifetime of the cold clouds and their peak mass, 𝑀 max , follows a power-law distribution 𝜏 life ∝ 𝑀 𝛼 max . Cold clouds with peak masses 𝑀 max up to 10 4 M ⊙ have a power-law index 𝛼 = 0 . 39, whereas clouds with peak masses higher than 10 4 M ⊙ have 𝛼 = 0 . 06. The more massive the clouds the longer they live up to 10 4 M ⊙ and then their lifetimes saturate at ∼ 8 Myr.\n- · The lifetime of the cold clouds and their size at peak mass, 𝑅 M max , follows a power-law distribution 𝜏 life ∝ 𝑅 𝛼 M max . Cold clouds with sizes at peak mass 𝑅 M max up to 20 pc have a power-law index 𝛼 = 0 . 48, whereas clouds with sizes at peak mass larger than 20 pc have 𝛼 = 0 . 13. The larger the clouds the longer they live up to 20 pc and then their lifetimes tend to stabilise around 8 Myr.\n- · The size of the cold clouds correlates approximately with the square root of their peak mass: 𝑅 M max ∝ 𝑀 1 / 2 max .  \nOur analysis indicates that the assumed interstellar medium model forms realistic cold clouds with short lifetimes embedded in a multi-phase interstellar medium. In combination with the realistic star cluster properties forming in these clouds the same simulation (Lahén et al. 2020a) it appears that the basic cycle of cold cloud formation and clustered star formation in the galactic ISM is well captured by the GRIFFIN model. The model supports the rapid formation and disruption of star-forming clouds by stellar radiation and supernovae on timescale less than 10 Myr in the low metallicity interstellar medium.', 'ACKNOWLEDGEMENTS': "We would like to thank Chia-Yu Hu for his help in the realisation of the present work. T. Naab acknowledges support from the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany's Excellence Strategy - EXC-2094 - 390783311 from the DFG Cluster of Excellence 'ORIGINS'. J.M. Hislop and P.H. Johansson acknowledge the support by the European Research Council via ERC Consolidator Grant KETJU (no. 818930) and the Academy of Finland through grant 339127. Part of the simulations were carried out at CSC-IT Center for Science Ltd. in Finland. We acknowledge the Gauss Centre for Supercomputing e.V. (www.gauss-centre.eu) for providing computing time on the GCS Supercomputer SUPERMUC-NG at Leibniz Supercomputing Centre (www.lrz.de) under project numbers pn49qi and pn72bu. The data analysis was performed on the FREYA cluster hosted by The Max Planck Computing and Data Facility (MPCDF) in Garching, Germany. T.-E. Rathjen acknowledges support by the project 'NRW-Cluster for data-intensive radio astronomy: Big Bang to Big Data (B3D)' funded through the programme 'Profilbildung 2020', an initiative of the Ministry of Culture and Science of the State of North Rhine-Westphalia. The sole responsibility for the content of this publication lies with the authors.", 'DATA AVAILABILIY': 'The data underlying this article will be shared on reasonable request to the corresponding author.', '16 C. M. Fotopoulou et al.': 'Rosolowsky E., 2005, PASP, 117, 1403 Röttgers B., Naab T., Cernetic M., Davé R., Kauffmann G., et al., 2020, MNRAS, 496, 152 Schmidt M., 1959, ApJ, 129, 243 Schruba A., et al., 2012, AJ, 143, 138 Scoville N. Z., Hersh K., 1979, ApJ, 229, 578 Scoville N. Z., Solomon P. M., 1975, ApJ, 199, L105 Skarbinski M., Jeffreson S. M. R., Goodman A. A., 2023, MNRAS, 519, 1887 Smith R. J., et al., 2020, MNRAS, 492, 1594 Smith M. C., Bryan G. L., Somerville R. S., Hu C.-Y., Teyssier R., Burkhart B., Hernquist L., 2021, MNRAS, 506, 3882 Solomon P. M., Rivolo A. R., Barrett J., Yahil A., 1987, ApJ, 319, 730 Springel V., 2005, MNRAS, 364, 1105 Steinwandel U. P., Moster B. P., Naab T., Hu C.-Y., Walch S., 2020, MNRAS, 495, 1035 Stutzki J., Bensch F., Heithausen A., Ossenkopf V., Zielinsky M., 1998, A&A, 336, 697 Szakacs R., Péroux C., Zwaan M. A., Nelson D., Schinnerer E., Lahén N., Weng S., Fresco A. Y., 2022, MNRAS, 512, 4736 Tasker E. J., 2011, ApJ, 730, 11 Tasker E. J., Tan J. C., 2009, ApJ, 700, 358 Tremonti C. A., et al., 2004, ApJ, 613, 898 Treß R. G., Sormani M. C., Smith R. J., Glover S. C. O., Klessen R. S., Mac Low M.-M., Clark P., Duarte-Cabral A., 2021, MNRAS, 505, 5438 Walch S., Naab T., 2015, Mon. Not. R. Astron. Soc., 451, 2757 Walch S., et al., 2015, MNRAS, 454, 238 Westera P., Lejeune T., Buser R., Cuisinier F., Bruzual G., 2002, A&A, 381, 524 Whitmore B. C., Chandar R., Schweizer F., Rothberg B., Leitherer C., et al., 2010, AJ, 140, 75 Wiersma R. P. C., Schaye J., Smith B. D., 2009, MNRAS, 393, 99 Williamson D. J., Thacker R. J., Wurster J., Gibson B. K., 2014, MNRAS, 442, 3674', 'APPENDIX A: LINKING LENGHT STUDY': 'In the paper we use a linking length of 1 pc to identify cold clouds. In Fig. A1 we show the histogram of cold cloud masses for a representative snapshot identified with linking lengths of 0.5 pc (blue), 1 pc (orange), 5 pc (green) and 10 pc (red). A linking length of 1 pc allows for the best identification of low mass clouds ( ≲ 10 3 𝑀 ⊙ ) while the total mass detected in cold clouds ( 𝑀 tot ) changes very little, mainly for larger linking lengths for which we loose the information about low mass clouds again. Therefore we consider a linking length if 1 pc to be the best choice.', 'APPENDIX B: CLOUD LIFETIME STUDY': 'In Fig. B1 we show a comparison of the cloud lifetime distributions at the peak starburst, when we would expect the strongest impact of the environment and 130 Myr later when the system is in a more relaxed state. We do not find evidence for a change in the lifetime distribution with respect to the environment.  \nThis paper has been typeset from a T E X/L A T E X file prepared by the author.  \nFigure A1. Histograms of the masses of the identified clouds for different linking lengths used in the FoF algorithm: ℓ = 0.5 pc (blue), 1 pc (orange), 5 pc (green) and 10 pc (red).  \n<!-- image -->  \nFigure B1. Histograms of the lifetimes of cold clouds forming during the merger: 160 Myr < 𝑡 form < 180 Myr (red) and forming during a more quiescent later stage of the evolution of the system: 290 Myr < 𝑡 form < 310 Myr (blue).  \n<!-- image -->'}