Clouds of Theseus: long-lived molecular clouds are composed of short-lived H2 molecules

We use passive gas tracer particles in an Arepo simulation of a dwarf spiral galaxy to relate the Lagrangian evolution of star-forming gas parcels and their H2 molecules to the evolution of their host giant molecular clouds. We find that the median chemical lifetime of H2 is just 4 Myr, independent of the lifetime of its host molecular cloud, which may vary from 1 to 90 Myr, with a substantial portion of all star formation in the galaxy occurring in relatively long-lived clouds. The rapid ejection of gas from around young massive stars by early stellar feedback is responsible for this short H2 survival time, driving down the density of the surrounding gas, so that its H2 molecules are dissociated by the interstellar radiation field. This ejection of gas from the H2-dominated state is balanced by the constant accretion of new gas from the galactic environment, constituting a"competition model"for molecular cloud evolution. Gas ejection occurs at a rate that is proportional to the molecular cloud mass, so that the cloud lifetime is determined by the accretion rate, which may be as high as 4 x 10^4 Msol/Myr in the longest-lived clouds. Our findings therefore resolve the conflict between observations of rapid gas ejection around young massive stars and observations of long-lived molecular clouds in galaxies, that often survive up to several tens of Myr. We show that the fastest-accreting, longest-lived, highest-mass clouds drive supernova clustering on sub-cloud scales, which in turn is a key driver of galactic outflows.


INTRODUCTION
In recent years, a new, dynamical picture of star formation has emerged.Numerical simulations of disk galaxies have demonstrated that Lagrangian parcels of gas in the interstellar medium undergo constant and rapid cycles of compression into, and disruption out of, a high-density 'star-forming state' (Semenov et al. 2017(Semenov et al. , 2018;;Shin et al. 2022).The average time spent by gas in the star-forming state (∼ 5-15 Myr) is much shorter than the time spent outside the star-forming state ( > ∼ 100 Myr).This picture can explain several of the key observed attributes of star formation in disc galaxies.The fact that galaxies convert only a small fraction of their gas to stars per galactic rotation (e.g.Kennicutt 1998;Wyder et al. 2009;Daddi et al. 2010) is consistent with the small fraction of time spent by gas in the star-forming state.The observed spread of star-formation rates (SFRs) at a given molecular gas surface density found when galaxies are observed at high resolution (Onodera et al. 2010;Schruba et al. 2011), and the relatively poor correlation between molecular masses and star formation rate indicators for clouds in the Milky Way (Mooney & Solomon 1988;Lee et al. 2016), reflect the fact that any particular parcel of H 2 rich gas may have just started forming stars (in which case its luminosity per unit gas mass in SFR tracers such as H or infrared emission will be low) or may have just completed a star formation cycle and left the star-forming state (in which case the luminosity per unit mass will be high).The spatial decorrelation of molecular gas and young stars on small scales ( 100 pc) relative to galactic scales (∼ 1 kpc, Schruba et al. 2010;Kruĳssen et al. 2019a) and the rapidity with which star clusters become optically revealed (Hollyhead et al. 2015;Sokal et al. 2016) are consistent with the disruption of dense, star-forming gas on short (< 5 Myr) time-scales by the radiation and thermal pressure from young, massive stars (presupernova stellar feedback).Finally, the near-proportionality of the star formation rate surface density Σ SFR and the molecular gas surface density Σ H 2 observed on kpc scales (e.g., Wong & Blitz 2002;Bigiel et al. 2008;Leroy et al. 2013) can also be explained by a selfregulation process, in which the density distribution of molecular gas is shaped by the stars it forms (Semenov et al. 2019).
However, the connection between this Lagrangian picture of star formation and the properties of observable star-forming regions or 'giant molecular clouds' remains unclear.Molecular clouds are observed to have a large range of masses, sizes and densities, spanning over two orders of magnitude in Milky Way-like disc galaxies (e.g.Sun et al. 2018;Colombo et al. 2019).These star-forming regions are not universally destroyed on time-scales of 5-15 Myr.Both observations (e.g.Scoville & Hersh 1979;Engargiola et al. 2003;Blitz et al. 2007;Murray 2011;Corbelli et al. 2017) and simulations (e.g.Dobbs & Pringle 2013;Grisdale et al. 2019;Jeffreson et al. 2021a) show a large range of survival times for giant molecular clouds, ranging over two orders of magnitude, from 1 to > 100 Myr.While short-lived (∼ 10 Myr), low-mass (∼ 10 4 M ) molecular clouds dominate by number, a substantial fraction of molecular mass (and thus galactic star formation) resides in the most massive ( > ∼ 10 6 M ) molecular clouds (e.g.Murray 2011;Miville-Deschênes et al. 2017;Faesi et al. 2018), which tend to live for longer periods of time (e.g.Jeffreson et al. 2021a).As such, the mass-weighted molecular cloud lifetime is typically up to a factor of 10 longer than the star-forming lifetime of a Lagrangian gas parcel.We may therefore ask: how do we reconcile the short Lagrangian star formation time-scale with the long lifetimes of observable molecular clouds?
This question is closely-related to another unanswered question in the field of star formation: what is the chemical lifetime of molecular hydrogen and how does it relate to the molecular cloud lifetime?Based on observations of massive molecular clouds in the low-density inter-arm regions of spiral galaxies, Scoville & Hersh (1979) and Koda et al. (2009Koda et al. ( , 2016) ) argue that the survival time of hydrogen molecules must be comparable to the inter-arm crossing time.Given that the long time-scale for the formation of new hydrogen molecules from the low-density, inter-arm gas prohibits the formation of new massive clouds in these regions, these authors argue that high-mass molecular agglomerations must form within the spiral arms, then later fragment into smaller (but still massive) molecular clouds as they transit into the inter-arm regions.However, this hypothesis is at odds with the short Lagrangian lifetimes of starforming gas (Semenov et al. 2017): when cold, dense, star-forming gas parcels are ejected into a warm, diffuse non-star-forming state by stellar feedback on time-scales of ∼ 5-15 Myr, their hydrogen molecules should also be dissociated by the interstellar radiation field.
The chemical lifetime of H 2 has important consequences for interpreting the observed spatial distribution of young star clusters and dense molecular gas, and in particular their decorrelation on sub-kiloparsec scales (e.g.Blitz et al. 2007;Kawamura et al. 2009;Miura et al. 2012;Corbelli et al. 2017;Schruba et al. 2010;Kruĳssen et al. 2019a).If the H 2 survival time-scale is short, then the observed spatial decorrelation of young stars and H 2 on sub-kiloparsec scales implies that stellar feedback rapidly dissociates molecular gas and destroys the star-forming gas in its vicinity, providing evidence for short lifetimes of star-forming regions.However, if the H 2 survival time-scale is long, the same observations imply that the molecular gas is simply pushed away from young stars, where it can continue to form new stars.
In this work, we seek to reconcile the apparent contradiction between the Lagrangian and Eulerian views of molecular cloud lifetime, and to illuminate the relationship of cloud to chemical lifetimes.We use passive Lagrangian tracer particles in a chemodynamic simulation of an isolated disc galaxy in the moving-mesh code A to directly compare the star-forming time-scale of Lagrangian gas parcels, the chemical lifetime of molecular hydrogen, and the lifetimes of giant molecular clouds.We describe our numerical prescription in Section 2, and examine the distribution of molecular mass and star formation among short-and long-lived molecular clouds in Section 3. We compare our molecular cloud lifetimes to the Lagrangian star formation time-scale and the molecular hydrogen survival time in Section 4, and use this result to develop a simple picture of molecular cloud evolution in Section 4.3.In Section 4.4, we also show that the massive, long-lived molecular clouds formed in this picture account for the clustering of supernova feedback in our simulation, and are therefore the key drivers of galactic winds.We discuss our results in the context of the existing literature in Section 5, and present our conclusions in Section 6.

SIMULATION OF A DWARF SPIRAL GALAXY
Here we describe our simulation setup and methods.

Initial conditions
We simulate a dwarf flocculent spiral galaxy that is analagous in its gas and stellar mass distribution to the bulgeless nearby galaxy NGC300.The initial condition includes a dark matter halo at a mass resolution of 1.254 × 10 7 M , a stellar disc at a mass resolution of 3437 M , and a gas disc at a mass resolution of 859 M .The dark matter halo follows the profile of Navarro et al. (1997), with a concentration parameter of  = 15.4,a spin parameter of  = 0.04, a mass of 8.3 × 10 10 M and a circular velocity of  200 = 76 km s −1 at the virial radius.The stellar disc is of exponential form, with a mass of 1 × 10 9 M , a scale-length of 1.39 kpc, and a scale-height of 0.28 kpc.The corresponding exponential gas disc extends beyond the stellar disc, with a mass of 2.2 × 10 9 M (giving a gas fraction of 68 per cent) and a scale-length of 3.44 kpc.

Hydrodynamics and chemistry
The initial condition is evolved using the moving-mesh hydrodynamics code A (Springel 2010).Within A , the gas dynamics is modelled on the unstructured moving mesh defined by the Voronoi tesselation about a discrete set of points, which move according to the local gas velocity.The gravitational acceleration vectors of the Voronoi gas cells, stellar particles and dark matter particles are computed using a hybrid TreePM gravity solver.
The temperature and chemical composition of the gas in our simulation is modelled using the simplified network of hydrogen, carbon and oxygen chemistry described in Nelson & Langer (1997) and in Glover & Mac Low (2007a,b).The fractional abundances of the chemical species H, H 2 , H + , He, C + , CO, O and e − are computed and tracked for each gas cell, and self-consistently coupled to the heating and cooling of the interstellar medium via the atomic and molecular cooling function of Glover et al. (2010).The gas equilibrates to a state of thermal balance between line-emission cooling and heating due to the photo-electric emission from polycyclic aromatic hydrocarbons and dust grains, as they interact with the interstellar radiation field (ISRF) and with cosmic rays.We assign a value of 1.7 Habing fields to the UV component of the ISRF (Mathis et al. 1983), a value of 3 × 10 −17 s −1 to the cosmic ray ionisation rate (van der Tak & van Dishoeck 2000), and we assume the solar value for the dust-to-gas ratio.We note that this results in a higher metallicity than the sub-solar value that is observed for NGC300 (Bresolin et al. 2009), but we demonstrate in Section 3.1 that our simulation nevertheless produces a realistic interstellar medium structure, with a realistic distribution of molecular cloud properties.
We use the T C algorithm introduced by Clark et al. (2012) to model the dust-and self-shielding of molecular hydrogen from dissociation by the ISRF.This allows us to accurately model the nonequilibrium abundance of molecular hydrogen during the run-time of the simulation, and so to compute its value for each gas cell as a function of time.

Tracer particles
We introduce passive tracer particles to the simulation following the Monte Carlo prescription of Genel et al. (2013), which allows us to track the Lagrangian mass flow and molecular fraction of gas as it moves through simulated GMCs, despite the fact that A 5 kpc is not a Lagrangian code.Via this prescription, tracer particles are moved along with the gas cells in the simulation, and are exchanged between gas cells according to a probability set by the mass flux between them.When a gas cell is converted to a star particle, the tracer particles associated with that gas cell are moved to the star particle with a probability set by the ratio of the stellar mass to the original gas cell mass.Similarly, tracer particles attached to star particles are transferred back to the gas reservoir according to the masses ejected in stellar winds and supernovae.As such, the mass of tracers in the gas and stellar reservoirs remains equal to the masses of these reservoirs throughout the simulation.
In this work, we assign one tracer particle to each gas cell in the initial condition, which sets an initial mass distribution for the tracer particles, equal to the initial mass distribution of gas cells.Since the cells are not of exactly equal mass initially, this means that the tracers are not uniform in mass either.However, we show in Appendix A that the distribution of effective tracer masses converges to a uniform value after after < 100 Myr.We analyse the simulated disc only after it is in a state of dynamical equilibrium, between simulation times of 500 and 800 Myr.The effective tracer mass during this period is stable at a value of ∼ 450 M , corresponding to 1.9 tracer particles per gas cell.In what follows, we also make sure to analyse only large populations of tracer particles to mitigate the effects of the Poisson noise: our results are derived for samples containing 10 2 to 10 3 distinct molecular clouds, the smallest of which sees the passage of ∼ 100 tracer particles over its lifetime.

Star formation and feedback
The star formation efficiency  ff per free-fall time of the gas in our simulations follows the parametrisation of Padoan et al. (2017).These authors conduct high-resolution simulations of turbulent fragmentation and find that  ff depends on the local virial parameter  vir of the gas, according to  ff = 0.4 exp (−1.6 0.5 vir ). (1) To use such a star formation recipe, simulations need a model for unresolved turbulent velocity dispersion, , to calculate  vir ∝  2 /( 2 ).To this end, existing methods in galaxy formation simulations include explicit dynamic models for sub-grid turbulence (e.g., Braun & Schmidt 2015;Semenov et al. 2016;Kretschmer & Teyssier 2020), estimation of  from the gradients of resolved gas velocity (e.g., Hopkins et al. 2013), or calculation of  over some resolved scale (e.g., Gensior et al. 2020).As our simulations do not include a sub-grid model for turbulence, we follow the prescription of Gensior et al. (2020).In brief, the scale over which  vir is calculated is computed using Sobolev approximation:  = |  g / ∇ g |, where ∇ g is the cubic spline kernel-weighted average of the gas volume density gradient, with respect to the radial distance from the central gas cell.The smoothing length of the cubic spline kernel is chosen to enclose the 32 nearest-neighbour cells.We refer the reader to Gensior et al. (2020) for a more detailed explanation.The star formation rate volume density of each gas cell in the simulation is therefore given by d * , d where  ff, = √︁ 3/(32   ) is the local free-fall time-scale for the gas cell  with a mass volume density of   , and  ff is given by Equation (1).
We set a lower limit of  thresh / H  = 100 cm −3 on the volume density of hydrogen atoms above which star formation is allowed to occur, as well as an upper limit of  thresh = 100 K on the temperature; here  ≈ 1.4 is the mean mass per H atom.The value of  thresh is the density of Jeans-unstable (collapsing) gas at our mass resolution of 859 M and at the temperature of ∼ 30 K reached by the molecular gas in our simulation.For a spherical gas cell, the radius associated with this density threshold is 3 pc, and so we employ the adaptive gravitational softening scheme in A with a minimum softening length of 6 pc and a gradation of 1.5 times the Voronoi gas cell diameter.We set the softening length of the stellar particles to the same Our simulated dwarf spiral galaxy reproduces the spatial decorrelation observed for NGC300 between young stars and dense gas on the scales of molecular clouds-the so-called "tuning fork diagram" (Kruĳssen et al. 2019a, thin lines with markers).The branches of the tuning fork show the average bias of gas depletion times measured in the apertures of a given size (-axis) placed on the peaks of dense gas (top branch) or young stars (lower branch).To quantify the snapshot-to-snapshot variation of the tuning fork in our simulation, the thick solid lines show medians and the shaded regions show 16-84 and 2.5-97.5 interpercentile ranges.The diagram is computed over the range of galactocentric radii of  = 2-6 kpc (see the text for details).
value, and choose a softening length of 280 pc for the dark matter particles, according to the convergence tests presented in Power et al. (2003).Because our simulations resolve the gas disc scale-height and the Toomre mass at all scales, the adaptive gravitational softening avoids the majority of artificial fragmentation at scales larger than the Jeans length (Nelson 2006).The stellar feedback in our simulation is modelled via injection of energy and momentum by supernovae and pre-supernova HII regions.To compute the number of supernovae generated by each star particle formed during the simulation, we assign a stellar population drawn stochastically from a Chabrier (2003) initial stellar mass function (IMF), using the Stochastically Lighting Up Galaxies (SLUG) stellar population synthesis model (da Silva et al. 2012(da Silva et al. , 2014;;Krumholz et al. 2015).Within SLUG, the resulting stellar populations are evolved along Padova solar metallicity tracks (Fagotto et al. 1994a,b;Vázquez & Leitherer 2005) during run-time, using S 99-like spectral synthesis (Leitherer et al. 1999).This modelling provides the number of supernovae  * ,SN generated by each star particle during each simulation time-step, as well as the ionising luminosity of the cluster and the mass Δ * it has ejected.
We use the values of  * ,SN and Δ * for each star particle to compute the momentum and thermal energy injected by supernova explosions at each time-step.In the case of  * ,SN = 0, we assume that all mass loss results from stellar winds.In the case of  * ,SN > 0, we assume that all mass loss results from supernova explosions, and we model the corresponding kinetic and thermal energy injection due to the expanding blast-wave.At our mass resolution of 859 M per gas cell, the energy-conserving/momentum-generating phase of supernova blast-wave expansion is unresolved, and so we follow the prescription introduced by Kimm & Cen (2014): we explicitly inject the terminal momentum of the blast-wave into the set of gas cells  that share faces with the nearest-neighbour cell to the star particle.We use the unclustered parametrisation of the terminal momentum derived from the high-resolution simulations of Gentry et al. ( 2017 which is given by with an upper limit imposed by the condition of kinetic energy conservation, as the shell sweeps through the gas cells .The momentum is distributed among the facing cells as described in Jeffreson et al. (2021b).
The pre-supernova feedback from HII regions is implemented according to the model of Jeffreson et al. (2021b).This model accounts for the momentum injected by both radiation pressure and the 'rocket effect', whereby warm ionised gas is ejected from cold molecular clouds, following the analytic work of Matzner (2002) and Krumholz & Matzner (2009).Momentum is injected for groups of star particles with overlapping ionisation-front radii, which are identified using a Friends-of-Friends grouping prescription.The momentum is received by the gas cell closest to the luminosity-weighted centre of the Friends-of-Friends group, and is distributed to the set of adjoining neighbour cells.The gas cells inside the Strömgren radii of each Friends-of-Friends group are heated to a temperature of 7000 K, and are held above this temperature floor for as long as they receive ionising photons from the group.In contrast to Jeffreson et al. (2021b), we also explicitly and fully ionise the gas inside these Strömgren radii, rather than relying on the chemical network to do so.

PROPERTIES OF THE SIMULATED GALAXY
Here we compute the key observable properties of our simulated galaxy and its interstellar medium, and demonstrate that their values are close to those obtained via direct observations in comparable galactic environments.

Gas morphology and global interstellar medium structure
Figure 1 shows the spatial distribution of the total (left), atomic (centre-left), total molecular (centre-right) and CO-luminous molecular (right) gas reservoirs at face-on and edge-on viewing angles, at a simulation time of 800 Myr.The atomic and total molecular gas abundances are computed within our simplified chemical network during simulation run-time, as described in Section 2. To partition the molecular gas into its CO-luminous and CO-dark components, we post-process our simulation outputs using the D astrochemistry and radiative transfer model (Krumholz 2013(Krumholz , 2014)), as described in Appendix B, and delineate CO-dark from CO-bright material using a threshold that we detail below.
The qualitative structure of the CO-bright interstellar medium is similar to that observed by Faesi et al. (2018) in the nearby dwarf spiral galaxy NGC300, displayed in Figure 1 of Kruĳssen et al. (2019a).
In Figure 2, we quantify the agreement between the structure of our simulated interstellar medium and the observed interstellar medium in NGC300, showing that the simulation reproduces the observed spatial decorrelation between regions of recent star formation (traced by young stars) and regions of high molecular gas surface density (traced by CO emission).To this end, we use the so-called "tuning fork diagram" that quantifies the relative bias of depletion time measured in the apertures centered either on dense gas or young stars (e.g., Schruba et al. 2010;Kruĳssen et al. 2019a).The relative excess of dense gas or young stars in apertures of a given size (-axis) results in the upper (blue) and lower (green) branches, respectively.The values calculated from observations by Kruĳssen et al. (2019a) are given by the thin lines and open data points.
We calculate the dense gas and young stellar branches for our simulation as outlined in Semenov et al. (2021).In particular, we use molecular gas with a surface density of Σ H 2 ,CO > Σ H 2 ,CO,min = 13  pc −2 and a line-of-sight velocity dispersion of < 1 km s −1 × Σ H 2 ,CO /Σ H 2 ,CO,min as a proxy for the dense CO-bright gas; here Σ H 2 ,CO is the H 2 surface density that would be inferred based on the computed CO luminosity, rather than the true H 2 surface density.We define this quantity precisely in Appendix B. The line-of-sight velocity dispersion is computed as the CO-weighted mean velocity dispersion of the Voronoi gas cells along rays perpendicular to the galactic mid-plane.We assume that the H signal in young starforming regions is produced by star particles with ages 2-5 Myr.The shaded regions show the 1 and 2 spread of values for our simulated interstellar medium between simulation times of 500 and 800 Myr.
The convergence, on large scales, of the upper and lower branches in Figure 2 is a general property of star-forming galaxies (e.g.Kim et al. 2022), indicating the tight correlation between dense gas and young stars that is manifested in the near-linear relation between the molecular gas density and the star formation density observed in normal (non-starburst) star-forming galaxies on >kiloparsec scales (e.g., Wong & Blitz 2002;Bigiel et al. 2008;Leroy et al. 2013).On small scales, the branches diverge, as dense, CO-emitting gas and ionised, H-emitting gas become spatially decorrelated.
The precise shape of the branches in Figure 2 is sensitive to the modelling of star formation and stellar feedback used in the numer- ical simulation (Fujimoto et al. 2019;Semenov et al. 2021).The agreement of this shape with the relation computed from observational data by Kruĳssen et al. (2019a) indicates that our models produce a reasonable interstellar medium structure for galactocentric radii of  = 2-6 kpc.The inner 2 kpc of the simulation displays a stronger correlation (narrower branch opening) than is computed by Kruĳssen et al. (2019a).Although this is likely due to the existence of several high-mass molecular clouds in this region, which are removed ('masked') in their analysis, we exercise caution and exclude the region  < 2 kpc from our subsequent analysis, as we cannot verify its agreement with the observational data.
For the remaining galactocentric radii  > 2 kpc, we demonstrate in Figure 3 that our simulated dwarf spiral galaxy is similar in its atomic and molecular depletion times (upper panel) and its gas and stellar surface densities (lower panel) to NGC300.The simulated values are given by the bold lines, while the corresponding observed values from Westmeier et al. (2011) and Kruĳssen et al. (2019a) are given by the thin lines.

Identification
In this work, we identify CO-bright, observable clouds using isocontours of value log (Σ H 2 ,CO /M pc −2 ) = −1.5 on the surface density Σ H 2 ,CO of CO-bright molecular hydrogen perpendicular to the galactic mid-plane, shown in dark pink in Figure 4.This threshold corresponds to the natural break in the distribution of Σ H 2 ,CO produced by our chemical post-processing, which is described in detail in Appendix B. At surface densities higher than the thresh-old, gas cells contain at least some shielded, CO-dominated gas.At lower densities, CO exists only as a uniformly-mixed, unshielded, low-abundance component.
The yellow contours in Figure 4 enclose regions with a projected H 2 abundance of Σ H 2 /Σ gas > 0.3, accounting for 70 per cent of the total H 2 in the simulation, including CO-dark H 2 .We find that the majority of star formation (78 per cent) occurs in the CO-luminous H 2 , despite the fact that this gas reservoir accounts for only 26 per cent of the galactic molecular hydrogen mass.In what follows, we will therefore focus on the CO-luminous giant molecular clouds.We note that 99 per cent of the gas tracer particles occupying the CO-luminous state inside the pink contours of Figure 4 have previously resided in the CO-dark, H 2 -rich state denoted by the yellow contours.That is, the gas in our simulation almost always enters CO-luminous clouds by way of the CO-dark envelope.

Observational checks
An important check for the applicability of our simulation to the real Universe is whether the simulated population of CO-bright molecular clouds has similar properties to the observed population in a dwarf spiral galaxy.In Figure 5, we compare the instantaneous mass distribution (left panel) and size distribution (centre panel) of our simulated molecular clouds across all simulation times to the population of resolved clouds in the inner disc ( < ∼ 3 kpc) of NGC300 (Faesi et al. 2018) and in the outer regions ( > ∼ 2 kpc) of M33 (Gratier et al. 2012).We divide the simulated cloud population into four mass bins, to highlight the correspondence between the cloud mass, surface density and velocity dispersion.
We see that the slope of the simulated cloud mass distribution agrees well with the observed slope of  = 2 ± 0.1 in the outer regions of M33, but is significantly steeper than the corresponding slope in NGC300.It does not display a truncation at high masses, as seen for NGC300 (not shown in the figure).Given that our simulated cloud population is made up of clouds at  > 2 kpc, we would expect that our results are more closely-comparable to the M33 sample.Note however that such a comparison should be interpreted with caution, given the sensitivity of the mass function slope to the method of cloud identification used (Pineda et al. 2009;Hughes et al. 2013;Rosolowsky et al. 2021).Nevertheless, we note that the distribution of masses for the molecular clouds in our simulation is in reasonable agreement with observed values.
Comparing the coloured contours to the coloured lines in the righthand panel, we see that the cloud population generally follows a line of constant virial parameter, as seen in the observations of resolved clouds in NGC300 (black data points, Faesi et al. 2018).Higher-mass clouds are closer to virial equilibrium on average, with higher levels of gravitational boundedness.

Tracking and evolution
We track the temporal evolution of the simulated molecular cloud population via the algorithm described in Jeffreson et al. (2021a).In brief, the position of the two-dimensional isocontour enclosing a star-forming region at simulation time  is projected forward by the time-step Δ = 1 Myr of our simulation output, using the velocities of the gas cells in the region.A pair of star-forming regions is temporally linked as parent and child if there exists any overlap between the projected contour of the parent at time  and the contour outlining the child at time  + Δ.
Via this tracking procedure, we produce the cloud evolution network for the entire simulation, composed of ∼ 8000 complete, independent segments between simulation times of  = 500 and 800 Myr, and between galactocentric radii of  = 2 and 6 kpc.These segments correspond to time-evolving molecular clouds, which we use in the analysis presented in Section 4. The lifetimes of these clouds are calculated as the end-to-end time between the formation of the first parent cloud and the destruction of the last child, accounting for all cloud mergers and splits, as shown in Figure 6.

RESULTS
In this section, we combine our analysis of the molecular cloud network with results derived from the tracer particles in our simulation to study the formation and destruction of H 2 molecules within CObright molecular clouds, as a function of the cloud mass and cloud lifetime.We begin with the distribution of masses and lifetimes derived from the tracking of Eulerian clouds in Section 4.1, turn to the results from the tracer particles in Section 4.2, and combine the two viewpoints in Section 4.3.In Section 4.4, we discuss the impact of the Eulerian cloud mass distribution on the clustering of supernovae.

The Eulerian view: long-lived clouds
In the upper panel of Figure 7, we show the fraction of the simulated galactic H 2 reservoir (orange) and star formation rate (purple) that is accounted for by CO-bright molecular clouds of different lifetimes.We see that half of the star formation occurring in CO-bright molecular clouds (40 per cent of the galactic total) is accounted for by molecular clouds that live longer than 25 Myr, despite the fact that the most-common lifetime for molecular clouds, by number, is < 10 Myr (centre panel of Figure 7).Clouds with lifetimes > 25 Myr also contain about half the total CO-luminous H 2 mass in the galaxy.
In the lower panel of Figure 7, we demonstrate that molecular cloud lifetime is tightly correlated and monotonically increasing with the peak cloud mass achieved throughout this lifetime.We will discuss this proportionality in detail in Section 4.3, but for now we simply note that 'long-lived' molecular clouds are synonymous with 'massive' clouds.Therefore, a substantial fraction of the CO-bright molecular gas in our dwarf spiral galaxy simulation is contained in clouds that are both massive and long-lived.

The Lagrangian view: short-lived H 2 molecules
We now examine the chemical survival time of molecular hydrogen in CO-bright molecular clouds, by tracking the H 2 abundance of Lagrangian gas parcels (tracer particles) as they transit through these clouds.We consider a gas parcel to have passed through a cloud if it moves within the isocontour defining its edge at any time.We show the time evolution of the H 2 mass fraction  H 2 =  H 2 / gas for the collection of gas parcels transiting through a single example cloud in Figure 8.We omit tracer particles that remain at molecular hydrogen abundances below  H 2 = 10 −5 throughout the simulation run-time from our analysis; this cut excludes those tracer particles that are far from the galactic mid-plane but appear inside the isocontour in projection.
As the plot shows, individual fluid elements experience a wide range of variations in  H 2 , with many experiencing rapid fluctuations in the H 2 abundance, particularly as the cloud is disrupted.To quantify this behaviour, for each tracer particle that passes through a CO-bright molecular cloud, we identify all the contiguous time periods for which it has  H 2 < 0.5 and > 0.5; for each such period, we record its duration, the median (in time) value of  H 2 during it, and the lifetime of the CO-bright cloud through which it passed. 1 In the left and right panels of Figure 9, we plot the median  H 2 versus duration, with points colour-coded by molecular cloud lifetimes; the left panel shows intervals of  H 2 > 0.5, while the right shows  H 2 < 0.5.The plot shows that gas rapidly cycles in and out of an H 2dominated state ( H 2 > 0.5), in qualitative agreement with Semenov et al. (2017).The time spent by gas parcels in the H 2 -rich state is short, ranging from 1 to ∼ 20 Myr with a median value of 4 Myr.The median abundance of molecular hydrogen over the duration of each cycle is 0.63.Crucially, the distributions of H 2 molecule lifetimes and  H 2 values are almost identical between the samples of molecular clouds with different lifetimes (different colours).Thus we conclude that independent of the lifetime of the host molecular cloud, the chemical survival time of of H 2 molecules is about 4 Myr.
The right-hand panel of Figure 9 shows that the time spent by gas parcels in the H 2 -poor state ( H 2 < 0.5) is correlated with the lifetime of the host molecular cloud, albeit weakly.Gas parcels ejected from higher-mass, longer-lived clouds (yellow) return to the H 2 -dominated state after shorter periods of time, and retain higher median H 2 abundances, than those associated with shorter-lived, lower-mass clouds (purple).That is, H 2 is not as efficiently or quickly dissociated in higher-mass molecular clouds.Gas parcels ejected from these clouds are more likely to become 'trapped', and cycle back to an H 2 -dominated state multiple times: the typical number of cycles for tracer particles transiting through clouds of lifetime < 25 Myr is 1, but for regions of lifetime > 25 Myr, the number of cycles is between 2 and 4.
In Figure 10, we demonstrate that the ejection of gas from the H 2 -dominated state is driven by early stellar feedback.To construct 1 For particles that have multiple episodes of  H 2 > 0.5 with a period of  H 2 < 0.5 between them, we take the molecular cloud lifetime for the H 2 -poor phase to be the lifetime of the molecular cloud through which the particle passed during the preceding H 2 -rich phase.Essentially no tracer particles pass through two different CO-bright molecular clouds during a single cycle of having  H 2 > 0.5 and < 0.5.
Figure 10.The behaviour of gas parcels (tracer particles) around all young star particles in our simulation.The colorbar corresponds to the lifetime of the host molecular cloud.Thin lines show the median values for the gas tracer particles that are located within 10 pc of a star particle at the instant of its formation, while thick lines show the medians across this population, computed in bins of molecular cloud lifetime.The left column shows the time evolution starting at the instant of star formation, while the right column shows the time evolution starting from the moment when the star particle produces its first supernova.We see that early, pre-supernova feedback (left-hand column) ionises and drives the gas away from young stars on time-scales between 2 and 5 Myr across all molecular clouds, with the ejection being slightly slower in the longer-lived, higher-mass regions.Supernovae (right-hand column) generally occur when the gas has already been expelled from an H 2 -rich state.
this figure, we tag every Lagrangian tracer particle that lies within 10 pc of a star particle at the moment of its formation.For each tagged tracer, we also identify the lifetime of the GMC within which that tracer resides.On the left-hand side of Figure 10, we plot the time evolution of the H 2 fraction, total gas density, and distance from the star for each of these tracers, terminating each line at the time of the star particle's first supernova.The thin lines in the figure correspond to the histories of individual tracers, and the thick lines to median values we obtain by binning the tracers by molecular cloud lifetime.Although the variation of individual tracks is large, we see that on average gas is efficiently ejected from the H 2 -dominated state by pre-supernova feedback on time-scales between 2 and 5 Myr.The gas is pushed out to a distance of ∼ 80 pc and down to a density of between 1 and 10 cm −3 (close to the bulk mean density of hydrogen in the interstellar medium), where its H 2 is dissociated by the interstellar radiation field.However, the efficiency of this ejection is a strong function of molecular cloud lifetime: gas parcels in highermass, longer-lived molecular clouds start at higher densities and H 2 fractions, and after star formation they retain higher density and H 2 abundance out to larger distances from the newly-formed star particle.On the right-hand side of Figure 10, we show the further time evolution of these same Lagrangian tracers, but starting at the instant of the first supernova.The plot demonstrates that the role of supernova feedback in destroying H 2 is secondary to that of pre-supernova feedback.By the time that stars explode as supernovae, the surrounding gas is already relatively H 2 -poor ( H 2 < 0.5) and close to the mean density of the bulk interstellar medium (∼ 1 cm −3 ).The main (partial) exception to this is is in the most massive and longest-lived clouds, where supernovae do further decrease the H 2 fraction.However, even for these clouds the mean H 2 fraction is below 50% by the time the first supernovae occur.
The increased difficulty of destroying H 2 molecules in highermass/long-lived molecular clouds by pre-supernova feedback results in an increase in the integrated star formation efficiency with cloud mass and lifetime.In Figure 11, we calculate the efficiency as the fraction of tracer particles that are transferred from gas to star particles as they transit through molecular clouds, and find an increase from ∼ 1 to ∼ 4 per cent for an increase in cloud lifetime from ∼ 10 to ∼ 90 Myr.The horizontal dashed line represents the median integrated star formation efficiency for the entire CO-luminous gas reservoir.This corresponds to the star formation efficiency associated with the most common cloud lifetime of ∼ 10 Myr, as reported in Figure 7.
Interestingly, Figure 11 shows that the lifetimes of molecular regions,  GMC , correlate near-linearly with their integrated SFE,  int .Such a correlation might explain the near-linear correlation between Σ H 2 and Σ ★ observed on kiloparsec scales in normal (non-starburst) galaxies.The depletion time of molecular gas can be expressed as  H 2 ≡ Σ H 2 / Σ ★ ∼  GMC / int , and therefore, if  int scales only with  GMC then  H 2 will be constant, i.e., Σ H 2 ∝ Σ ★ independent of kiloparsec-scale environment (see Semenov et al. 2019, for more detail).Indeed, we find that in the region where we apply our analysis,  H 2 is close to constant (see Figure 3).

Reconciling the views: molecular cloud evolution driven by the competition between accretion and ejection
We have shown in Section 4.1 that CO-bright molecular clouds are relatively long-lived.By contrast, in Section 4.2 we have shown that their constituent H 2 molecules are very short-lived.In Figure 12, we show how these two results can be reconciled.
In panel (), we explicitly show the mass evolution of all molecular clouds in our simulation, as a function of time.To aid visualisation, we divide the clouds into five bins of lifetime, indicated by colour.The median values at each time are given by the solid lines and the shaded regions indicate the corresponding interquartile ranges.In agreement with the trend shown in Figure 7, longer-lived molecular clouds achieve higher peak masses.The mass evolution of the molecular clouds in the simulation is remarkably symmetrical.All clouds begin and finish at the same mass, which is determined by the onset of effective H 2 self-and dust-shielding at our mass resolution of 859 M .The clouds increase in mass for the first half of their lives, and decrease in mass for the second half.
In panel () of Figure 12, we demonstrate that the mass evolution of the molecular clouds can be described by the competing accretion and ejection of gas into and out of a state of high H 2 abundance.The solid lines show the median mass of Lagrangian gas parcels (tracer particles) entering the H 2 -dominated state per unit time per molecular cloud; the dashed lines show the total mass exiting this state per unit time.The masses of the molecular clouds (upper panel) switch from increasing to decreasing when the rate of mass ejection crosses above the rate of mass accretion.From mass conservation, the evolution of cloud masses can be described by the simple formula where  () is the instantaneous mass of the cloud, and  accr and  ej are the instantaneous rates of H 2 mass accretion and ejection, respectively.The total lifetime of the cloud is therefore given by We note that this finding is broadly consistent with the analytic models proposed by Goldbaum et al. (2011), Inutsuka et al. (2015), Burkert (2017) and Kobayashi et al. (2017).These authors argue that molecular cloud lifetimes are determined primarily by external accretion rates, and that the galactic distributions of molecular cloud mass and star formation rate results from a competition between gas accretion driven by the compression of gas at the interfaces between large-scale converging flows and gas ejection driven by galactic shear and stellar feedback.Panels () and () of Figure 12 demonstrate that the rate of molecular gas ejection from clouds is correlated in time with the instantaneous star formation rate and, closely related, the momentum injection rate from HII regions.The median rate of star formation per cloud is indicated by the dashed lines in panel (), and the corresponding radial momentum injected by the HII regions is indicated by the dashed lines in panel ().The star formation rate tracks the cloud mass, and therefore the rate of molecular gas ejection by stellar feedback also tracks the cloud mass.It begins at the same value for all clouds, and diverges only after it crosses the mass accretion rate and so begins to destroy the cloud.
We demonstrate the proportionality between the mass ejection rate  ej and the rate of momentum injection   by HII regions more explicitly in Figure 13.The fact that  ej ∝   can be understood in terms of the momentum required to eject material from a (spherical) shows rate of radial momentum injection by HII regions; both of these rates are computed using a 5 Myr averaging window to suppress fluctuations due to the stochasticity of star formation.Note that the total mass follows a symmetrical pattern of mass evolution from formation to destruction, and that the mass loss rate from the H 2 -rich phase, star formation rate, and momentum injection rate are all strongly correlated.The schematic at the top shows a possible lifecycle for a short-lived (green) and a long-lived (yellow) molecular cloud, including the effects of accretion and feedback-induced gas ejection.bounding surface at radius  shield around each young star particle, where  shield is the radius at which  H 2 drops to 0.5, set by the extent of self-and dust-shielding.In this case, where  esc is the escape speed at  shield , which we write in terms of the mean volume density  of the gas inside  shield .The variation in the values of  and  shield is not very large, resulting in the proportionality  ej ∝   shown in Figure 13.
We therefore find that the rate of mass ejection  ej across all molecular clouds in our simulation can be straight-forwardly parametrised in terms of the rate of momentum injection   into these regions by early feedback.

The clustering of supernovae in massive, long-lived molecular clouds
The concentration of star formation in the most massive and longlived molecular clouds (Section 4.1), fed by the fast accretion of new molecular gas (Section 4.3), drives the spatial and temporal clustering of supernovae.We demonstrate this in Figure 14, which shows the two-point correlation function  (Δ) of supernova explosions occurring in more and less massive clouds.The black dashed line corresponds to the time-averaged median value of the two-point correlation function for all supernova explosions occurring across the simulated galaxy, within all time slices of 1 Myr over the 300 Myr simulation interval.This line indicates that on scales Δ < ∼ 80 pc, the supernovae are more clustered ( > 1) than would be expected for a Poisson (uniform) distribution of objects across the galactic midplane, and on scales Δ > 80 pc they are less clustered ( < 1).The degree of clustering of supernovae rises steeply on small scales of < ∼ 30 pc, and Figure 14 demonstrates that this clustering is accounted for almost exclusively by supernovae in clouds with masses for supernova explosions as a function of their separation Δ over time intervals of 1 Myr, averaged over all times (solid lines).The shaded regions give the interquartile ranges over these times.The black dashed line is calculated for all supernovae across the entire simulation, while the green and purple lines represent the supernovae associated with star particles born in long-lived, high-mass and short-lived, low-mass molecular clouds, respectively.Massive and long-lived clouds account for supernova clustering on scales < ∼ 30 pc.
> 10 5 M (green line).By contrast, the supernovae in lower-mass clouds (purple line) approach a random distribution at small scales.
The clustering of supernovae on such small scales, approaching the softening length of 6 pc used in our simulation, has been shown to enhance the momentum injected per supernova to the interstellar medium by a factor of at least 4 (e.g.Gentry et al. 2017).This, in turn, increases the mass-loading of outflows (Fielding et al. 2018) as well as the burstiness of star formation across galaxies (Smith et al. 2021).
Figure 15 helps to visualise this intuitive result, comparing the time-evolving H 2 column density of a low-mass cloud in our simulation (∼ 5 × 10 4 M , upper panels) to that of a high-mass cloud (∼ 10 6 M , lower panels).White filled stars represent young star particles of age < 5 Myr, which are emitting significant quantities of ionising radiation, and empty stars denote supernova explosions.The yellow arrows show the velocities of gas tracer particles.In the low-mass cloud, star formation within a radius of ∼ 100 pc is halted as the feedback from the young star particle destroys the entire molecular gas reservoir.By contrast, the high-mass cloud continues to accrete more molecular mass from the lower right-hand corner (inward-pointing yellow arrows), even as an entire cluster of supernovae are occurring in the top left-hand corner.Pre-supernova feedback destroys only small sections of the cloud, while new stars continue to form in other nearby sections.When massive stars eventually explode as supernovae, they are therefore surrounded by other nearby supernova explosions: a cluster of supernovae.

DISCUSSION
Here we discuss both the implications of our results in the context of the existing literature and the caveats associated with our calculations.3. High-mass star-forming regions cannot be entirely destroyed by distributed pre-supernova feedback, and so are responsible for the clustering of supernovae on small scales in the simulation.

Previous results for molecular cloud lifetimes
Our demonstration that accretion-ejection balance is the key process in molecular cloud evolution helps resolve a number of problems in the literature regarding molecular cloud lifetimes.First, it helps explain why molecular cloud lifetimes computed in simulations of isolated GMCs tend to be too short compared to observations.In such isolated cloud simulations, molecular cloud lifetimes for observationally-motivated initial densities are typically 5 − 10 Myr (e.g., Grudić et al. 2018Grudić et al. , 2021;;He et al. 2019;Fukushima et al. 2020;Kim et al. 2021), a factor of 2 − 5 shorter than the ≈ 20 − 30 Myr mean generally inferred from observations of nearby galaxies (e.g.Kawamura et al. 2009;Miura et al. 2012;Kruĳssen et al. 2019b;Chevance et al. 2020).Our result provides a natural explanation for this discrepancy: simulations of isolated clouds are necessarily missing the accretion part of the accretion-ejection balance.With no fresh mass supplies, these simulations produce lifetimes closer to the Lagrangian lifetime of individual molecules than to the lifetimes of real molecular clouds.
Our results also help to explain the seeming contradiction between the presence of molecular clouds in the inter-arm regions of grand design spiral galaxies and a typical observed cloud lifetimes of 20 − 30 Myr, much less than the inter-arm transit time.Our simulated dwarf spiral galaxy has an average gas density (1 − 5 cm −3 ) and velocity dispersion (∼ 10 km s −1 ), similar to the values found in the interarm regions of grand design spiral galaxies such as M51.Despite these conditions, and despite the short median survival time of an H 2 molecule in our simulation (∼ 4 Myr, see Figures 9 and 10), we nonetheless find that massive molecular clouds in our simulation can survive for up to 90 Myr.These long lifetimes are sustained by a high accretion rate of new molecular gas throughout the cloud lifetime: up to 4 × 10 4 M Myr −1 for the longest-lived clouds.This accretion rate is an order of magnitude faster than predicted by Koda et al. (2009Koda et al. ( , 2016) ) for the comparable inter-arm environment in M51, and therefore provides an explanation for the existence of clouds in this inter-arm region that is not in contradiction with a much shorter typical molecule survival time (∼ 4 Myr) and a much shorter typical cloud survival time (≈ 20 − 30 Myr, see Figure 7).Possible processes driving this fast accretion are converging flows due to supernova-driven bubbles or converging flows due to large-scale galactic-dynamical processes such as galactic shear.Pinpointing and describing this rate of accretion in our simulation will be the topic of a future paper.

Caveats
We conclude this discussion by pointing out two limitations of calculations, and considering their possible impact.

Unresolved shielding of molecular gas
The resolution of our simulation (859 M , probing down to ∼ 2 pc inside molecular clouds) presents a potential caveat to the result that the chemical survival time of H 2 molecules is 4 Myr in a dwarf spiral galaxy.Specifically, Joshi et al. (2019) show that the H 2 abundance derived using our chemical network and shielding prescription is converged only at a resolution of 0.2 pc, due to the existence of sub-pc dense clumps of molecular gas.
In such clumps both the ratio of the H 2 formation and destruction rates and the amount of self-shielding would be enhanced compared to those we measure at our ≈ 2 pc resolution, and as a result the value of 4 Myr that we have derived for this lifetime may be an underestimate.This would also mean that the rate  ej of gas ejection from molecular clouds would be lowered relative to the value we have calculated.The rate  accr of gas accretion would be unaffected.However, there is circumstantial evidence that this cannot be a major effect.We have shown in Figure 3 that the average HI and H 2 gas densities and depletion times in our simulated galaxy are in good agreement with observations, as are the masses and sizes of its giant molecular clouds (Figure 5).If the H 2 survival time were significantly longer, agreement with observations would worsen.This suggests that the rate limiting step in H 2 destruction is not so the timescale for chemical reactions, which we do not resolve, but the timescale for molecular material to be dynamically ejected from well-shielded clouds, which we do.Thus our failure to capture the details of the chemistry does minimal harm.Given this discussion, however, we do caution that the H 2 survival time of 4 Myr depends on a number of factors, including the resolution of our simulation, the properties of the large-scale galactic environment, and the metallicity of the simulated gas reservoir (see next section).We have shown that this survival time is short, but great significance should not be attached to its exact value.

Galactic-scale metallicity variations
We have assumed a globally-invariant, Solar value for the gas-todust ratio in our simulation.However, the observed metallicity of NGC300 in fact varies between 0.8 and 0.2 times the Solar value, decreasing monotonically with galactocentric radius (Bresolin et al. 2009).According to our competition model for the evolution of molecular clouds, lowering the gas-phase metallicity in the outer regions of the galactic disc may have several different (competing) effects.Reducing the metallicity will reduce the shielding of H 2 from the ISRF, and so reduce the distance from young stars at which H 2 is ionised as it is ejected by early stellar feedback.This will cause an increase in the rate  ej of H 2 ejection from molecular clouds.However, lowered metallicity will also reduce the degree of cooling and star formation in the dense gas, causing a decrease in  ej .The effect of lowered metallicity on the mass accretion rate  accr is more straight-forward: (1) the dust-catalysed formation of H 2 from atomic hydrogen will be slower, and (2) gas parcels will need to reach higher densities before becoming shielded from the ISRF.Both of these effects will act to decrease the rate of H 2 accretion onto molecular clouds.
As such, it is difficult to predict exactly how the inclusion of live metal-enrichment will affect our results, except to note that the environmentally-dependent metallicity will likely be a variable that plays a significant role in setting  ej .Further studies that vary the metallicity and control for galaxy morphology are necessary to determine the functional form of this dependence.

CONCLUSIONS
In Section 1 we posed two questions: (1) what is the chemical lifetime of H 2 in the interstellar medium of galaxies and how does it relate to the molecular cloud lifetime?And (2) how can we reconcile the long observed lifetimes of giant molecular clouds (e.g.Koda et al. 2009) and the short Lagrangian residence times of gas in the star-forming state (e.g.Semenov et al. 2017)?In this work, we have used passive gas tracer particles in an A simulation of a dwarf spiral galaxy to arrive at the following answers: (i) The typical survival time of H 2 molecules is 4 Myr in our simulation, independent of the lifetime of the host molecular cloud, which has a much longer median value of ∼ 25 Myr.The rapid destruction of H 2 molecules is driven by early stellar feedback from HII regions.
(ii) The reason that molecular clouds are able to survive much longer than their constituent molecules is that the molecular cloud lifetime is determined by the competition between the feedbackdriven ejection of molecular gas, and the accretion of new molecular gas from the large-scale galactic environment.This competition means that long-lived and massive molecular clouds can be sustained given a sufficiently-fast accretion rate.In our simulation, a value of  accr ∼ 4 × 10 4 M Myr −1 can feed a molecular cloud for 90 Myr, and up to a peak mass of 10 6 M .
Our answers have important implications for the dynamics of gas on both cloud and galactic scales.In particular, we find that (i) The time spent by gas in the H 2 -dominated state is much shorter than the time-spent in the H 2 -poor state, in qualitative agreement with the results of Semenov et al. (2017).A given Lagrangian parcel of gas spends most of its life transiting between molecular clouds, with only brief passages through clouds.
(ii) A partial exception to this is that gas can sometimes become 'trapped' around higher-mass molecular clouds, cycling back to a dense, H 2 -dominated state faster and more frequently than occurs for gas parcels passing through low-mass clouds.This 'trapping' raises the integrated star formation efficiency in massive, long-lived molecular clouds up to four times higher than that in low-mass, short-lived clouds.
(iii) Despite their small number, massive and long-lived molecular clouds account for a not insignificant amount of the total star formation budget of a galaxy.Moreover, they drive the clustering of supernova feedback on sub-cloud scales, which is in turn a key driver of galactic outflows.
In an upcoming paper, we will investigate the galactic-scale physics driving the fast accretion of new gas onto massive, long-lived molecular clouds, and so driving their higher star formation efficiencies and levels of supernova clustering.the native resolution of the simulation (6 pc).The black solid line has slope unity, and shows the expected relationship for a uniform distribution of effective tracer particle masses, while the dashed lines denote the Poisson error expected due to the Monte Carlo exchange of a small number of tracer particles between gas cells.The plot shows that, after a period of 100 Myr (green points and error bars), both the number of tracer particles per pixel, and the error in this number, conform to the expectation for a uniform distribution of effective tracer particle masses.

Figure 1 .
Figure1.Column density maps of the total (Σ gas , left), atomic (Σ HI , centre-left), total molecular (Σ H 2 , centre-right), and CO-luminous molecular (Σ H 2 ,CO right) gas distribution for the simulated dwarf spiral galaxy, viewed perpendicular to (top panels) and across (lower panels) the galactic mid-plane, at a simulation time of 800 Myr.
Figure 2. Our simulated dwarf spiral galaxy reproduces the spatial decorrelation observed for NGC300 between young stars and dense gas on the scales of molecular clouds-the so-called "tuning fork diagram"(Kruĳssen et al. 2019a, thin lines with markers).The branches of the tuning fork show the average bias of gas depletion times measured in the apertures of a given size (-axis) placed on the peaks of dense gas (top branch) or young stars (lower branch).To quantify the snapshot-to-snapshot variation of the tuning fork in our simulation, the thick solid lines show medians and the shaded regions show 16-84 and 2.5-97.5 interpercentile ranges.The diagram is computed over the range of galactocentric radii of  = 2-6 kpc (see the text for details).

Figure 3 .
Figure3.Our simulated dwarf spiral galaxy is similar to NGC300 in its gas and stellar surface densities and gas depletion times.Upper panel: Molecular and atomic gas depletion times as a function of galactocentric radius (solid lines), measured within annuli of width 500 pc.Thin lines show the values measured for NGC300 byKruĳssen et al. (2019a).Lower panel: Stellar, molecular and atomic gas surface densities as a function of galactocentric radius, measured within the same annuli.Thin lines show the values measured for NGC300 byKruĳssen et al. (2019a) and byWestmeier et al. (2011).

Figure 4 .
Figure 4. Thresholds for the identification of CO-luminous giant molecular clouds (pink isocontours) and their CO-dark envelopes down to a fractional projected abundance of 0.3 (yellow isocontours).The upper panel shows the CO-luminous molecular hydrogen column density (pink) and the lower panel shows the total molecular hydrogen column density (blue-yellow).

Figure 5 .
Figure 5. Left: Mass distribution of CO-bright molecular clouds in our simulation, divided into four mass bins.The solid black line gives the observed powerlaw slope d /d ∝  − ,  = 1.76 ± 0.07 found by Faesi et al. (2018) in NGC300.The dotted black line gives the powerlaw slope  = 2 ± 0.1 found by Gratier et al. (2012) in M33.Centre: Size distribution of CO-bright molecular clouds in our simulation, divided into the same four mass bins.Right: Molecular gas line-of-sight velocity dispersion as a function of the molecular gas surface density for the CO-luminous molecular clouds.The solid lines indicate virial parameters of  vir = 2 for spherical beam-filling clouds at the mean region size for each mass bin.Black data points and errorbars represent the resolved GMC sample of Faesi et al. (2018) in NGC300.

Figure 6 .
Figure 6.The lifecycle of a single simulated molecular cloud with a lifetime of 10 Myr.The actual evolution of the H 2 surface density as it moves across the galaxy is shown in the upper panel.Both the CO-bright and CO-dark H 2 are shown.This region undergoes a split and re-merger 4 Myr after its birth.The corresponding section of the cloud evolution network is shown in the lower half of the figure.

Figure 7 .
Figure 7. Upper panel: Time-averaged cumulative fraction of the total galactic molecular mass (orange) and star formation rate (SFR; purple) accounted for by CO-bright clouds, as a function of their lifetime.The solid lines denote the time-averaged median values, while the shaded regions denote the corresponding interquartile ranges.The normalisation of the CDFs reflects the fact that CO-dark gas accounts for 22 per cent of star formation and 74 per cent of galactic H 2 (see Section 3.2.1).The orange and purple dashed vertical lines represent the median cloud lifetimes, weighted by the H 2 mass and by galactic SFR, respectively.Centre panel: Probability distribution of the molecular cloud lifetime by number.Lower panel: Median and interquartile range of the molecular cloud peak mass, as a function of the cloud lifetime.Grey circles represent the values for individual GMCs.

Figure 8 .
Figure8.The evolution of the molecular hydrogen abundance  H 2 of individual Lagrangian gas parcels (thin black lines), as they transit through a CO-bright molecular cloud of lifetime 10 Myr.This is the same molecular cloud shown in Figure6.The black vertical lines denote the times at which the cloud appears and disappears, based on our cloud-tracking algorithm.

Figure 9 .
Figure 9.The time spent by gas tracer particles in the H 2 -rich ( H 2 ≥ 0.5; left) and H 2 -poor ( H 2 < 0.5; right) phases, for tracer particles that transit through CO-bright molecular clouds.The transparent data points give the values for 1/5000th of the cycles undergone by tracer particles, while the solid data points and lines give the median values and interquartile ranges along each axis, computed in bins of molecular cloud lifetime as indicated by the colours of the solid points.All points are coloured according to the lifetime of the cloud through which the tracer particle transits.The time-scales spent in the H 2 -poor state are strictly lower limits, as 60 per cent of the values abut the beginning or end of the simulation window at  = 500 or  = 800 Myr, and so are truncated.

Figure 11 .
Figure 11.Integrated star formation efficiency over the lifetime of a CObright molecular cloud, as a function of its lifetime.The black data points and error bars correspond to the median values and interquartile ranges in each lifetime interval.The horizontal dashed line gives the total stellar mass produced in all clouds, divided by the total mass of all clouds, which is equal to the total SFE within the simulation.The transparent data points correspond to the values for individual molecular clouds, coloured by the peak mass they attain during their lifetimes.

Figure 12 .
Figure 12.Cloud properties as a function of time since cloud formation, computed as medians over bins of total cloud lifetime, from the shortest-lived clouds (purple) to the longest-lived (yellow).Panel (a) shows instantaneous cloud mass, with the solid lines showing the medians and the shaded bands showing the interquartile range.Panel (b) shows mass flux into (solid) and out of (dashed) the H 2 -dominated phase,  H 2 > 0.5.Panel (c) shows star formation rate, and (d)shows rate of radial momentum injection by HII regions; both of these rates are computed using a 5 Myr averaging window to suppress fluctuations due to the stochasticity of star formation.Note that the total mass follows a symmetrical pattern of mass evolution from formation to destruction, and that the mass loss rate from the H 2 -rich phase, star formation rate, and momentum injection rate are all strongly correlated.The schematic at the top shows a possible lifecycle for a short-lived (green) and a long-lived (yellow) molecular cloud, including the effects of accretion and feedback-induced gas ejection.

Figure 13 .
Figure13.The rate of mass ejection  ej as a function of the radial momentum injection rate   , for all molecular clouds in the simulation that are receiving feedback momentum.The scatter plot shows 1/5th of the total measurements for the cloud population, coloured according to the lifetime of the cloud.The black points and errorbars give the median and interquartile range values in 10 different mass bins.To guide the eye, the black solid line corresponds to the case of direct proportionality between mass ejection and momentum injection (not a fit).The black dashed horizontal line shows the resolution limit of the simulation, which corresponds to ejecting one resolution element of mass Δ = 859 M over the time interval Δ = 1 Myr that we use in our numerical computation of  ej .

Figure 14 .
Figure14.The two-point correlation function  (Δ) for supernova explosions as a function of their separation Δ over time intervals of 1 Myr, averaged over all times (solid lines).The shaded regions give the interquartile ranges over these times.The black dashed line is calculated for all supernovae across the entire simulation, while the green and purple lines represent the supernovae associated with star particles born in long-lived, high-mass and short-lived, low-mass molecular clouds, respectively.Massive and long-lived clouds account for supernova clustering on scales < ∼ 30 pc.

Figure 15 .
Figure 15.Examples of stellar feedback occurring in a low-mass, short-lived star-forming region (∼ 5 × 10 4 M , upper panels) and in a high-mass, long-lived star-forming region (∼ 10 6 M , lower panels).Time runs from left to right across the page, with the region ages given in white.Star particles of age < 5 Myr (emitting pre-supernova feedback) are indicated by filled white stars, while supernova explosions are indicated by open white stars.The velocities of the tracer particles that transit through the clouds are indicated by yellow arrows.The colorbar indicates the total molecular gas surface density, and white contours in the upper panel enclose regions of projected molecular hydrogen abundance Σ H 2 /Σ g ≥ 0.3.High-mass star-forming regions cannot be entirely destroyed by distributed pre-supernova feedback, and so are responsible for the clustering of supernovae on small scales in the simulation.

Figure A1 .
Figure A1.Mean and standard deviation of the tracer number as a function of gas surface density in 2D columns of area (6 pc) 2 through the entire midplane of the simulated galaxy disc, at simulation times of 0 Myr, 10 Myr and 100 Myr.The pixel size corresponds to the native resolution of the simulation, at which molecular clouds are identified and analysed in two dimensions.The thick black line corresponds to one tracer particle per 450 M , while the dashed black lines show a range of ± √ no.tracers around this, the error expected from Poisson sampling for no.tracers 1.