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

HAL
open science

Cosmic evolution of the star formation efficiency in Milky Way-like galaxies

Álvaro Segovia Otero, Oscar Agertz, Florent Renaud, Katarina Kraljic, Alessandro B Romeo, Vadim A Semenov

► To cite this version:

Álvaro Segovia Otero, Oscar Agertz, Florent Renaud, Katarina Kraljic, Alessandro B Romeo, et al. Cosmic evolution of the star formation efficiency in Milky Way-like galaxies. Monthly Notices of the Royal Astronomical Society, 2025, 538 (4), pp.2646 - 2659. 10.1093/mnras/staf423. insu-05023696

HAL Id: insu-05023696

https://insu.hal.science/insu-05023696v1

Submitted on 7 Apr 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Cosmic evolution of the star formation efficiency in Milky Way-like galaxies

Álvaro Segovia Otero<sup>1,★</sup> Oscar Agertz<sup>1</sup>,<sup>1</sup> Florent Renaud<sup>1,2,3</sup> Katarina Kraljic,<sup>2</sup> Alessandro B. Romeo<sup>4</sup> and Vadim A. Semenov<sup>5</sup>

Accepted 2025 March 9. Received 2025 March 9; in original form 2024 October 11

## ABSTRACT

Current star formation models are based on the structure of the interstellar medium (ISM), yet the details on how local physics propagates to galactic-scale properties are still debated. To investigate this, we use VINTERGATAN, a high-resolution cosmological zoom-in simulation of a Milky Way-like galaxy. We study how the velocity dispersion and density structure of the cold neutral ISM on 50–100 pc scales evolve with redshift and quantify their impact on the star formation efficiency per free-fall time-scale, $\epsilon_{ff}$. During starbursts velocity dispersions can reach ~50 km s⁻<fnref n="1" />, especially throughout last major merger events (1.3 < z < 1.5). After a merger-dominated phase (1 < z < 5), VINTERGATAN transitions into evolving secularly, featuring velocity dispersion levels of ~10 km s⁻<fnref n="1" />. Despite strongly evolving density and turbulence distributions over cosmic time, $\epsilon_{ff}$ at the resolution limit is found to change by only a factor of a few: from median efficiencies of 0.8 per cent at z > 1 to 0.3 per cent at z < 1. The mass-weighted average shows a universal $\langle\epsilon_{ff}\rangle \approx 1$ per cent, caused by an almost invariant virial parameter distribution in star-forming clouds. Changes in their density and turbulence levels are coupled, so the kinetic-to-gravitational energy ratio remains close to constant. We show that a theoretically motivated $\epsilon_{ff}$ is intrinsically different from its observational estimates adopting tracers of star formation, e.g. H$\alpha$. Since the physics underlying star formation can be lost on short time-scales (~10 Myr), caution must be taken when constraining star formation models from observational estimates of $\epsilon_{ff}$.

**Key words:** methods: numerical – ISM: structure – galaxies: star formation.

## 1 INTRODUCTION

Current star formation theories strive to connect, and track across cosmic time, the local properties of star-forming regions ($\lesssim$100 pc), the global parameters characterizing galaxies (> tens of kpc), and their cosmological environment. Nearby star-forming spiral galaxies follow the canonical main-sequence (MS) and Kennicutt–Schmidt (KS) relations: at a given stellar mass, disc galaxies with larger molecular gas surface densities present greater star formation rate (SFR) surface densities, depleting their gas reservoirs on time-scales of ~ few Gyr (gas depletion times are estimated as $\tau_{dep} = M_{gas}/\text{SFR}$; Kennicutt 1998; Bolatto et al. 2008; Leroy et al. 2013; Speagle et al. 2014). The absence of major mergers combined with low gas fractions ($f_g \lesssim 10$ per cent for $M_\star \gtrsim 10^{10}$ M$_\odot$; e.g. Scholte et al. 2024) allow these isolated systems to develop stable rotationally supported thin discs (Brinchmann et al. 2004; Wuyts et al. 2011; Wisnioski et al. 2015), where star formation on galactic scales is regulated by gas accretion, disc dynamics and instabilities, and feedback processes.

In the Local Universe, it is mainly within giant molecular clouds (GMCs) that the ISM is shielded from background radiation, reaching the required density and temperature conditions to trigger star formation (for a review, see McKee & Ostriker 2007). CO observations of GMCs in the Milky Way have revealed cloud scaling relations that relate their masses, sizes, and velocity dispersions (Larson 1981; Solomon et al. 1987; Heyer et al. 2009; Miville-Deschênes, Murray & Lee 2017). These have typical molecular gas surface densities of $\Sigma_{gas} \approx 100 \text{ M}_\odot \text{ pc}^{-2}$, effective radii of $R \sim 10–100 \text{ pc}$, and velocity dispersions of $\sigma_{gas} \sim 1–10 \text{ km s}^{-1}$. These properties translate into cloud virial parameters of $\alpha_{vir} \approx 1–10$ (Leroy et al. 2016; Miville-Deschênes et al. 2017; Sun et al. 2022), where $\alpha_{vir}$ is defined as the ratio between the kinetic to gravitational energy such that $\alpha_{vir} = 2E_{\text{kin}}/|E_{\text{grav}}|$. A constant $\alpha_{vir}$ of $\gtrsim 1$ is often understood as GMCs reaching dynamical equilibrium, although external cloud confinement by thermal and turbulent pressure complicates this picture (e.g. Elmegreen 1989; Grisdale et al. 2018). Star formation under these conditions is inefficient, which is commonly quantified using the star formation efficiency per free-fall time parameter $\epsilon_{ff}$, i.e. the fraction of stars formed from the total gas mass of a GMC on a free-fall time-scale $t_{ff} = \sqrt{3\pi/32G\rho}$. Observations of star-forming GMCs show $\epsilon_{ff}$ distributions with a universal average of

© 2025 The Author(s).

<footnote n="1"><sup>1</sup>Lund Observatory, Division of Astrophysics, Department of Physics, Lund University, Box 43, SE-221 00 Lund, Sweden</footnote>

<footnote n="2"><sup>2</sup>Observatoire Astronomique de Strasbourg, Université de Strasbourg, CNRS UMR 7550, F-67000 Strasbourg, France</footnote>

<footnote n="3"><sup>3</sup>University of Strasbourg Institute for Advanced Study, 5 allée du Général Rouvillois, F-67083 Strasbourg, France</footnote>

<footnote n="4"><sup>4</sup>Department of Space, Earth and Environment, Chalmers University of Technology, SE-41296 Gothenburg, Sweden</footnote>

<footnote n="5"><sup>5</sup>Center for Astrophysics | Harvard & Smithsonian, 60 Garden St, Cambridge, MA 02138, USA</footnote>

<footnote>*E-mail: alvaro.segovia@astro.lu.se</footnote>
1 per cent, followed by a spread of several orders of magnitude (see review by Krumholz, McKee & Bland-Hawthorn 2019). Preliminary systematic changes in $\epsilon_{ff}$ and $\alpha_{vir}$ have been found as a function of galactic environment (bulge, bars, spiral arms, disc), but the scatter around these parameters is dominated by intrinsic cloud-to-cloud differences (Sun et al. 2020a, b; Rosolowsky et al. 2021). The way in which gas densities and velocity dispersions in the ISM are affected by environmental factors, and whether this has a significant impact on the local star formation properties, still remains unclear even at low redshift.

At high redshifts ($z > 1$), star-forming galaxies are predominantly rotating discs, and their high gas fractions make them gravitationally unstable to clump formation (for reviews see Förster Schreiber & Wuyts 2020; Saintonge & Catinella 2022). Such a clumpy ISM is subject to strong feedback from enhanced SFRs, along with more frequent galactic interactions (Duncan et al. 2019). All of this sets elevated velocity dispersions at high redshifts, as measured from kinematics of molecular, atomic, and ionized gas tracers (Übler et al. 2019, and references therein). Compared to spiral disc galaxies at present day, average gas depletion times in high-$z$ galaxies are found to be shorter ($\tau_{dep} \ge 500$ Myr; see Tacconi et al. 2018; Tacconi, Genzel & Sternberg 2020). In contrast to low-redshift discs, the global $\tau_{dep}$ in starburst galaxies is around an order of magnitude shorter (30–300 Myr; Daddi et al. 2010; Rodighiero et al. 2011), making them a natural testbed for star formation physics in more extreme ISM conditions. Starbursts can be driven by mergers, but can also be triggered by secular processes in isolated gas-rich turbulent discs (Ciesla et al. 2023).

The star-forming gas clumps in starbursts observed to date are denser and more turbulent than GMCs in Local spirals, nevertheless, they also show $\alpha_{vir} \approx 1 - 10$, albeit with a large scatter (Leroy et al. 2015; Rosolowsky et al. 2021; Dessauges-Zavadsky et al. 2023, see also Krahm et al. 2024 for GMCs in the overlap region between the Antennae galaxies). Unfortunately, observations of gas clumps in such objects are scarce and limited to bright mergers in the Local Universe or lensed high-redshift galaxies (Dessauges-Zavadsky et al. 2019). It is therefore uncertain whether ISM conditions at high redshift or in extreme environments such as starbursts imply efficient local star formation, i.e. a higher $\epsilon_{ff}$.

Mapping the density and turbulent structure of the ISM to its star formation properties has been the goal of theoretical and computationally driven studies in the last couple of decades. Pioneering analytical work base the star formation process on the theory of supersonic isothermal gas (Krumholz & McKee 2005; Hennebelle & Chabrier 2011; Padoan & Nordlund 2011; Federrath & Klessen 2012; Burkhart 2018). These models, further calibrated on controlled magnetohydrodynamic (MHD) simulations of GMCs (5–500 pc in size; Grudić et al. 2018, 2019) and highly resolved ISM boxes (~ few pc in length, Federrath et al. 2010; Federrath & Klessen 2013; Federrath 2015), can predict the commonly observed average $\langle\epsilon_{ff}\rangle \approx 1$ per cent even though they locally present large $\epsilon_{ff}$ variations (Krumholz, Dekel & McKee 2012). More realistic set-ups face the challenge of resolving the large dynamic range that connects the inner structure of GMCs and their galactic and extra-galactic environment in isolated disc or cosmological simulations respectively. To that aim, subgrid star formation recipes are either based on the aforementioned analytical prescriptions where $\epsilon_{ff}$ directly depends on the local gas density and velocity dispersion (Semenov, Kravtsov & Gnedin 2016; Trebitsch et al. 2017; Kretschmer & Teyssier 2020; Nuñez-Castiñeyra et al. 2021), or have a fixed $\epsilon_{ff}$ typically in the range of 1–100 per cent (Renaud et al. 2013; Grisdale et al. 2018, 2019; Hopkins et al. 2018). These simulations reveal a complex interplay

among the extragalactic and galactic environment, star formation models, and stellar feedback prescriptions (e.g. Agertz, Teyssier & Moore 2011; Agertz & Kravtsov 2015; Semenov, Kravtsov & Gnedin 2018). For this reason, different star formation models can yield significant differences in ISM and GMC properties (e.g. Grisdale et al. 2017, 2018), but be degenerate in terms of integrated quantities on galactic scales such as SFRs or global gas depletion times (Hopkins, Narayanan & Murray 2013a).

In this paper, we use VINTERGATAN (Agertz et al. 2021; Renaud et al. 2021a, b), a cosmological zoom-in simulation of a Milky Way-like galaxy to investigate the evolution of the ISM and local star formation properties, i.e. $\epsilon_{ff}$ as a function of redshift. In our simulation, the Padoan, Haugbølle & Nordlund (2012) effective model connects the density and velocity dispersion of the ISM with $\epsilon_{ff}$, which we further compare to Local Universe GMC observations from the PHANGS collaboration (Physics at High Angular resolution in Nearby GalaxieS, Schinnerer & Leroy 2024). This work is a follow-up of Segovia Otero, Renaud & Agertz (2022) and Renaud, Segovia Otero & Agertz (2022), where we demonstrated the crucial role played by galaxy mergers and morphological transformations in setting the global gas depletion time. Here, we focus on investigating the physics on smaller ($\lesssim 50-100$ pc) scales and how it connects to the global properties of a galaxy across cosmic time. Section 2 briefly summarizes the simulation set-up, Section 3 highlights the connection between the cosmological environment and the local star formation properties of Milky-Way like galaxies, and Sections 4 and 5 conclude with the main takeaways of this work and contextualize them with respect to other available models for star formation.

## 2 METHODOLOGY

The VINTERGATAN cosmological zoom-in simulation was run using the hydrodynamics + N-body code RAMSES (Teyssier 2002), assuming $\Lambda$ cold dark matter ($\Lambda$CDM) cosmological parameters.¹ The initial conditions correspond to the ‘m12i’ halo of the AGORA project (Assembling Galaxies Of Resolved Anatomy, Kim et al. 2014, 2016, Roca-Fàbrega et al. 2021, see also Wetzel et al. 2023), generated with the MUSIC code (Hahn & Abel 2011). Within a periodic box of 85 Mpc and $512^3$ dark matter particles, the progenitor Lagrangian region is selected to be $3R_{200,m}$ around the halo at $z=0$, where the virial radius and mass with respect to the mean cosmic background density at $z=0$ are $R_{200,m} = 334$ kpc and $M_{200,m} = 1.3 \times 10^{12}$ M$_\odot$, respectively. Its resolution was further enhanced to reach a dark matter particle resolution of $3.5 \times 10^4$ M$_\odot$, and a mass and spatial resolution of $7070$ M$_\odot$ and $\lesssim 20$ pc, respectively for the gas. A more exhaustive explanation of the technical aspects of the simulation can be found in Agertz et al. (2021). Here, we briefly describe the most relevant models within VINTERGATAN.

Star formation is treated as a Poisson process on a cell-by-cell basis (see Agertz et al. 2013), where the star formation rate density $\dot{\rho}_*$ follows

$$ \dot{\rho}_* = \epsilon_{ff} \frac{\rho_g}{t_{ff}}, \quad \text{with } \rho_g > 100 \text{ cm}^{-3}, \text{ and } T_g < 100 \text{ K.} \quad (1) $$

Here, $t_{ff}$ is the free-fall time-scale and $\epsilon_{ff}$ is the star formation efficiency per free-fall time from Padoan et al. (2012) (henceforth, PN12). Star particles with an initial mass of $10^4$ M$_\odot$ form when gas is colder and denser than the density and temperature thresholds. These conditions are set to prevent thermally supported gas from forming

<footnote>$^1H_0 = 70.2\text{km s}^{-1}\text{Mpc}^{-1}$, $\Omega_m = 0.272$, $\Omega_\Lambda = 0.728$, $\Omega_b = 0.045$.</footnote>
stars, a choice that has a minor impact on the overall star formation history as the free-fall time-scales associated to such diffuse gas are long and do not contribute to the global SFR.

In PN12, they simulate star formation in supersonically turbulent, magnetized gas with a wide range of initial conditions and find that the star formation efficiency per free-fall time is well approximated by

$$ \epsilon_{ff} = \epsilon_w \exp(-1.38\sqrt{\alpha_{vir}}) = \epsilon_w \exp\left(-1.6\frac{t_{ff}}{t_{turb}}\right). \quad (2) $$

The normalization constant $\epsilon_w$, here set to 0.5, accounts for processes such as mass loss via protostellar jets. The virial parameter is defined as

$$ \alpha_{vir} = \frac{2E_{\text{kin}}}{|E_{\text{grav}}|} = \frac{15\sigma_{\text{ID}}^2}{\pi G \rho_g L^2} = 1.35 \left( \frac{t_{ff}}{t_{\text{turb}}} \right)^2, \quad (3) $$

where

$$ t_{\text{turb}} = \frac{L}{2\sqrt{3}\sigma_{\text{ID}}} \quad (4) $$

is commonly referred to as the turbulence crossing time, which is on the order of the time-scale for turbulent energy dissipation in the case of driven supersonic turbulence (Mac Low et al. 1998; Stone, Ostriker & Gammie 1998; Mac Low 1999). L is chosen to be the cell size (≈20 pc), and $\sigma_{\text{ID}}$ is the one-dimensional velocity dispersion calculated as the standard deviation of each component of the velocity vector independently for all six adjacent cells, i.e. effectively computed on ~50 pc scales, $\sigma_{\text{ID}}^2 = (\sigma_x^2 + \sigma_y^2 + \sigma_z^2)/3$.

The star formation model in PN12 is exclusively virial parameter dependent. Alternatives exist, e.g. with additional Mach number dependencies (for a comprehensive discussion see Federrath & Klessen 2012), but have not been directly implemented in this work. In Section 4, we comment on the potential implications of using different prescriptions for $\epsilon_{ff}$. It is also worth mentioning that even though sophisticated, these star formation laws do not directly account for the galactic scale environment surrounding the star-forming region, hence not accounting for mechanisms such as shear and large-scale compression.

Each star particle represents a single stellar population with a universal initial mass function (IMF; Chabrier 2003) and age-, mass-, and gas/stellar metallicity-dependent feedback processes such as radiation pressure, stellar winds, core collapse, and type Ia supernovae calibrated on the STARBURST99 code (Raiteri, Villata & Navarro 1996; Leitherer et al. 1999; Kim & Ostriker 2015). Gas metallicity is initialized from a floor of $Z = 10^{-3}$ Z$_\odot$ in the zoomed-in region<fnref n="2" /> and subsequently enriched with iron (Fe) and oxygen (O) by supernovae (hereafter SNe), where we adopt the yields from (Woosley & Heger 2007). This allows for metallicity-dependent cooling for $T < 10^4$ K (Rosen & Bregman 1995) and $10^4 < T < 10^8.5$ K gas (Sutherland & Dopita 1993). Gas heating from an ultraviolet (UV) radiation background was added assuming reionization starting at $z = 8.5$ (Haardt & Madau 1996; Courty & Alimi 2004; Aubert & Teyssier 2010).

To extract reliable information from our data, we use robust statistics (Müller 2000; Romeo, Agertz & Renaud 2023). Namely, the median and median absolute deviation (MAD) provide reliable information on the ‘central value’ and ‘width’ of the sample even when a large fraction of them are outliers. The robust $1\sigma$ scatter is

derived as MAD/0.6745. We further emphasize that using the median rather than the mean is convenient under logarithmic transformations given that $\text{med}[\log(X)] = \log[\text{med}(X)]$. The robust scatter in logarithmic scale is approximated to $\sigma_{\log(X)} \approx \sigma_X / [X \ln(10)]$, which is more representative in the limit where $\sigma_X / X \ll 1$.

## 3 RESULTS

VINTERGATAN proposes a formation channel general to nearby star-forming spirals that reproduces Milky Way features at $z \approx 0$ (Agertz et al. 2021; Renaud et al. 2021a, b). Segovia Otero et al. (2022) demonstrate that Milky Way-like galaxies experience epochs of starburst or more quiescent star formation activity depending on the global gas depletion time of the cold ISM. They identify three evolutionary stages with order of magnitude changes in $\tau_{\text{dep}}$: an early stage at high redshift ($5 < z < 9$) with $\tau_{\text{dep}} \sim 1$ Gyr; a starburst phase ($1 < z < 5$) when $\tau_{\text{dep}}$ drops by an order of magnitude driven by mergers impacting on an assembled galactic disc; a secular stage ($z < 1$) when $\tau_{\text{dep}}$ increases to a few Gyr in the absence of mergers. In the following, we explore how the density and turbulent structure of the ISM evolve with time depending on these epochs, and assess how much this affects the local star formation properties in the galaxy. Note that the word *local* is used in this work to refer to ISM properties derived on GMC scales when speaking about observations and for parameters computed on scales comparable to the size of gas cells in simulations. Conversely, the word *global* refers to galactic scales.

### 3.1 ISM density and velocity dispersion across cosmic time

Fig. 1 shows the ISM in the aforementioned evolutionary stages of VINTERGATAN. An additional panel representing the last major merger (LMM) exemplifies the properties of the ISM from its first pericenter passage ($z \approx 1.5$) to its final coalescence with the main galaxy ($z \approx 1.3$). Maps showing the gas density of the main galaxy (face- and edge-on), its environment, and its main global properties at each phase are on the bottom. In the top panels, contours show the density and velocity dispersion levels of the ISM, computed by stacking simulation snapshots with an output resolution of 100 Myr within each redshift interval. The dashed red lines mark regions of constant virial parameters of $\alpha_{vir} = 1$ (thin) and $\alpha_{vir} = 10$ (thick) using equation (3).

For every output, we select a spherical volume with radius equal to three times the stellar half-mass radius of the galaxy<fnref n="3" /> and split it in cubes of 50 pc in size.<fnref n="4" /> Within each cube we divide the gas into two phases by applying density and temperature cuts:

(i) The cold + dense ISM (blue contours) corresponds to gas that fulfills the star formation criteria from equation (1), i.e. above $n > 100 \text{ cm}^{-3}$ and below $T < 100$ K. This phase traces molecular species (Klessen & Glover 2016). Densities are computed by dividing the corresponding mass of cold + dense gas by the volume of the 50 pc cube. Velocity dispersion values $\sigma$ represent the one-dimensional mass-weighted standard deviation of velocities of cold + dense gas cells within 50 pc cubes that contain more than one cold + dense cell.

(ii) The diffuse ISM (green contours) includes gas cells at temperatures of $T < 10^4$ K excluding those with cold + dense gas properties.

<footnote n="2">²This metallicity floor accounts for unresolved Population III stars (Wise et al. 2012; Agertz et al. 2020).</footnote>

<footnote n="3">³Galactic centres in each output are found by running a shrinking sphere algorithm, and half-mass radii are calculated from the mass of stars younger than 100 Myr within the virial radius of VINTERGATAN at each epoch.</footnote>

<footnote n="4">⁴All length units reported in this work correspond to physical coordinates.</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Top: response of the ISM in VINTERGATAN as a function of redshift to its extragalactic environment. For every redshift interval, each panel shows the distribution of the density *n* and velocity dispersion *σ* calculated from gas in (50 pc)³ cubes. Contours represent the mass-weighted two-dimensional histogram with probability densities of 0.6, 0.3, 0.1, 0.05, and 0.01 from thickest and darkest to thinnest and brightest contours, respectively. Two sets of contours are presented: cold + dense (blue) for gas *T* < 100 K and *n* > 100 cm⁻³ (star-forming gas according to equation 2); diffuse (green) for gas *T* < 10⁴ K (excluding cold + dense gas). Dashed red lines of constant αvir are included (equation 3), computed on 50 pc scales for αvir = 1 (thin) and αvir = 10 (thick). Bottom: snapshots of the evolution of VINTERGATAN. Each column represents the edge- and face-on views of the galaxy, colour-coded by its gas density at every cosmic epoch: Early (5 < *z* < 9), Starburst (1 < *z* < 5) including the LMM event from the first pericenter passage (*z* ≈ 1.5) to the final coalescence of the merger (*z* ≈ 1.3), and Secular (*z* < 1). Included in each panel is the SFR, computed from the mass of stars younger than 100 Myr, the total stellar mass M*, the gas fraction *f*<sub>g</sub>, and twice the half radius calculated accounting for the spatial distribution of stars younger than 100 Myr (and outlined in black dashed circles on the bottom row panels). Takeaway: density and velocity dispersion go hand-in-hand, specially for cold and dense gas, and respond to the cosmological environment: mergers stir the ISM towards high *σ* and *n*, when a more isolated disc has more lower turbulent levels and density values.</figcaption></figure>

We include gas denser than 100 cm⁻³ and warmer than 100 K, yet its relative contribution is negligible as implied by the contour lines. These density and temperature ranges are typical of the warm neutral or ionized medium traced by neutral hydrogen (H1) and recombination lines such as Hα. Its densities and velocity dispersions have been computed in the same way as for the cold + dense gas, minding the new temperature and density constraints.

The velocity dispersion and density distributions for both *diffuse* and *cold + dense* gas phases evolve with cosmic epoch. This evolution happens at a constant virial parameter (see Section 3.2), with αvir ≥ 10 for *diffuse* gas and at 1 < αvir < 10 for *cold + dense* gas. In fact, *cold + dense* gas shows a positive correlation between the velocity dispersion and density with σ ∝ n^0.5, in agreement with previous studies of disc galaxies (Semenov et al. 2016). Such a scaling is compatible with the Larson scaling relations (Larson

1981; Grisdale et al. 2018) when a *fixed spatial scale* is considered. *Diffuse* gas has less inertia and is easier to accelerate, making it more susceptible to stellar feedback (Ejdetjärn et al. 2022) and the extragalactic environment. Thus, it deviates from the *σ − n* relation, with more scatter at *z* > 1.

During the *early* stage (5 < *z* < 9), most of the gas is found at velocity dispersion values of *σ* ≤ 30 km s⁻¹, which is in line with C II observations of galaxies in the 10⁸ ≤ M* ≤ 10¹⁰ M☉ mass range at these redshifts (left panel in Fig. 1, see also Pope et al. 2023; Posses et al. 2023). Although gas-rich and immersed in a merger-dominated environment, SFRs in VINTERGATAN are the lowest in its star formation history (Segovia Otero et al. 2022, for observations of gas fractions see Heintz et al. 2022; Aravena et al. 2023, and Fudamoto et al. 2022; Leethochawalit et al. 2023; Trussler et al. 2023 for SFRs of coeval systems). It is not until
the galaxy enters the *starburst* stage ($1 < z < 5$) that the highest velocity dispersion levels are reached, especially in the case of the *LMM* event for both gas phases (central panels in Fig. 1). This is in agreement with observed starbursting discs slightly more massive than VINTERGATAN at similar redshifts ($\sigma \approx 20-80\,\text{km s}^{-1}$ also from C II measurements; Herrera-Camus et al. 2022; Parlanti et al. 2023; Rizzo et al. 2023; Roman-Oliveira, Fraternali & Rizzo 2023). With the assembly of a galactic disc at $z \approx 4.8$, mergers tidally compress gas, triggering star formation on short depletion times due to an excess of gas in the *cold + dense* phase (Renaud et al. 2022; Segovia Otero et al. 2022). Nevertheless, even if an excess of gas at high densities leads to high SFRs (which inevitably goes hand-in-hand with increased stellar feedback), disentangling the roles of feedback- and gravity/compression-driven turbulence is challenging in cosmological zoom-in simulations (but see Renaud et al. 2014, who tell both effects apart using high time-resolution simulations of idealized mergers). In the absence of major mergers (*secular* stage, $z \le 1$), and with the depletion of gas reservoirs, SFRs and velocity dispersion values decrease most notably in the *diffuse* gas ($\sigma < 10\,\text{km s}^{-1}$ for H I observations in Local spirals; Tamburro et al. 2009; Eibensteiner et al. 2023). Consequently, its scatter also decreases substantially.

## 3.2 Properties of the cold and dense ISM

Densities reached in VINTERGATAN imply free-fall time-scales at least a factor of 10 shorter than the output frequency of the simulation (~100 Myr). Thus, to get insight into the actual ISM properties upon star formation we use the *star formation events* (hereafter *SF-events*) gas. This data set is different to the previous two gas phases in that they have not been extracted directly from the low-cadence outputs of the simulation. *SF-events* properties are logged on the fly from gas cells that have undergone star formation: densities are directly read from the star-forming cell, and velocity dispersion values are calculated using the velocity vectors of the six adjacent cells. Stars are spawned from gas cells with the highest resolution, which is approximately kept at ~20 pc throughout the simulation, with adjacent cells at a similar level of refinement. Hence, the velocity dispersion is computed on scales that roughly correspond to sizes of 50–100 pc. *SF-events* gas is represented in solid black in Fig. 2, revealing very turbulent (≥ 50 km s⁻¹), dense (≥ 10⁴ cm⁻³), and cold (≈10 K) properties that resemble those that make molecules such as HCN, HCO⁺, NH₃, CH, CN, and CS (Klessen & Glover 2016; Gallagher et al. 2018b; Wilson et al. 2023).

From this point on, we focus our analysis on the *cold + dense* and *SF-events* gas in order to connect the ISM structure with its star formation properties at each cosmic epoch. Fig. 2 shows the evolution of $\alpha_{vir}$ (equation 3), Mach number,<fnref n="5" /> $t_{ff}$, $t_{turb}$ (equation 4), density, and velocity dispersion for *SF-events* gas in black and *cold + dense* gas in blue. *Cold + dense* gas is further split in two: gas within 50 pc cubes containing stars younger than 10 Myr in solid blue; gas cubes without star particles younger than 10 Myr in dotted blue. It is worth noting that the velocity dispersion values in both *cold + dense* and *SF-events* gas are computed at approximately the same scales, but other

properties such as $t_{ff}$ (or $c_s$) are derived from the density (or mass-weighted average temperature) on 50 pc cubes compared to reading off the star-forming cell properties on ≲20 pc before spawning a star particle. We emphasize that our goal here is to: (i) analyse the properties of the ISM at the resolution limit of the simulation where the PN12 star formation model actually comes into play; (ii) compare the distributions of *SF-events* and *cold + dense* properties; and (iii) compare VINTERGATAN to GMC scaling relations from the PHANGS collaboration, who provide GMC properties on a variety of scales (45–150 pc; Rosolowsky et al. 2021; Sun et al. 2023; Schinnerer & Leroy 2024).

The most striking result in Fig. 2 is the difference between the *SF-events* and the *cold + dense* gas distributions. Note that the distributions are normalized, so *SF-events* gas actually conforms a smaller sample of cells representing the more extreme regions of the $\sigma-n$ space. An order of magnitude shift to shorter $t_{ff}$ and $t_{turb}$ dynamical time-scales is seen, highlighting the importance of capturing these short-lived events, and implying that *SF-events* gas is indeed a biased sub-set of denser and more turbulent gas intimately connected to stellar feedback. To explore this idea, we compare *cold + dense* gas with and without recent star formation (solid and dotted blue in Fig. 2). We find that gas embedding stars younger than 10 Myr reach denser and more turbulent states, i.e. stronger resemblance to *SF-events* distributions. The same trends are recovered when defining recent star formation as stars younger than 5 Myr and 100 Myr (time-scales probed by Hr and UV observations, respectively), with stronger shifts between recent and past star-forming ISM the older the stellar age cut is.

Looking at the *cold + dense* distributions, the shortest $t_{ff}$ and $t_{turb}$ are achieved during the *LMM* event compared to the other epochs (third and fourth rows in Fig. 2), expected from the evolution of the corresponding density and velocity dispersion distributions. Nevertheless, both dynamical time-scales conspire to render an $\alpha_{vir}$ probability density function (PDF) that is roughly constant in all four cosmic epochs (top row in Fig. 2), with median values of $\alpha_{vir} \approx 4.4$ and robust scatter of 3.3. Approximately stationary $\alpha_{vir}$ distributions with time is also noticed for *SF-events* gas, which overlap substantially with *cold + dense* gas ones despite the stark differences between them. This is particularly the case at $z > 1$, with both *SF-events* and *cold + dense* $\alpha_{vir}$ distributions showing marginally larger median values. In merger-driven environments at these redshifts, higher $\alpha_{vir}$ values are observed in recent simulations by He et al. (2023), yet CO observations of GMCs in mergers and lensed turbulent discs at $z > 1$ cannot conclusively estimate a systematic change. Some studies point towards $\alpha_{vir} \approx 1$ in both mergers and isolated galaxies (Dessauges-Zavadsky et al. 2019, 2023), whereas others report increased $\alpha_{vir}$ values in starbursts, low-density regions, galactic centre environments, and even mergers (Leroy et al. 2015). For the particular case of the Antennae galaxies, the order of magnitude increase in physical resolution to ~10 pc from Krahm et al. (2024) compared to Wei, Keto & Ho (2012) renders an average of $\langle \alpha_{vir} \rangle = 13.80 \pm 7.84$.

In the absence of mergers during the *secular* phase, the median $\alpha_{vir}$ is around a factor of 3 larger in *SF-events* gas with respect to *cold + dense*. Such moderate shift to lower virial parameters in *cold + dense* gas is due to lower velocity dispersion values increasing $t_{turb}$, and is in good agreement with CO-based observations of GMCs in nearby galaxies<fnref n="6" /> (highlighted in orange boxes on the top right

<footnote n="5">⁵The expression for the Mach number is $M = \sqrt{3\sigma_{1D}/c_s}$, where $c_s$ is the thermal sound speed determined by the temperature of gas in the cell T, the Boltzmann constant k_B, a mean molecular weight of $\mu = 1.3$, and the mass of hydrogen $m_H$ so that $c_s = \sqrt{k_B T / \mu m_H}$. Note that the sound speed for both *cold + dense* and *SF-events* gas is $c_s \lesssim 1\,\text{km s}^{-1}$, making turbulence supersonic for all measured values of the velocity dispersion.</footnote>

<footnote n="6">⁶Our equation for $\alpha_{vir}$ is almost identical to that of Sun et al. (2018), (2020b), except that in their work they use the projected gas surface density, $\Sigma_g = \pi R^2$, where R is the depth of the line of sight. They also add a geometrical</footnote>
<figure><img src="image_4.png" /></figure>

**Figure 2.** Gas mass-weighted PDFs of $\alpha_{vir}$, $M$, $t_{ff}$, the turbulent time-scale $t_{turb}$, and velocity dispersion at each cosmic epoch in VINTERGATAN. PDFs plotted in solid black correspond to SF-events while those in blue represent cold + dense gas. Curves in solid blue show ISM properties evaluated on 50 pc cubes with star particles younger than 10 Myr. Those in dotted blue PDFs characterize the ISM with no recent star formation. Additionally, vertical lines correspond to the median values of each distribution with same colour and line style. The orange rectangles in the top right panel represent the scatter (robust standard deviation, see last paragraph in Section 2) in $\alpha_{vir}$ values evaluated on 80 pc in Sun et al. (2018) and 90 pc in Sun et al. (2020b). Black arrows are just included to clarify the extent of the scatter. **Takeaway:** even if the density and turbulence of the ISM on 50 pc scales changes with cosmic epoch, the virial parameter shows a roughly constant distribution as a function of redshift.
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Cosmic evolution of the $\epsilon_{ff}$ gas mass-weighted PDFs (**left**) and local star formation time-scales (**right**) in 15 redshift bins. The on-the-fly $\epsilon_{ff}$ is calculated with equation (2) (PN12), using SF-events properties as input. Small line markers at the bottom the left and right most panels indicate the evolution of the median $\epsilon_{ff}$ and $t_{sf}$, respectively. Two additional panels display the re-calculated PDFs for only cells with $\epsilon_{ff} > 5$ per cent (**centre left**). For clarity, four colour-coded PDFs show the averaged distributions for each cosmic epoch. **Takeaway:** $\langle\epsilon_{ff}\rangle \approx 1$ per cent is constant throughout cosmic time, and therefore the local time-scales for star formation are mostly dependent on the changes in the gas density PDFs, through $t_{ff}$.</figcaption></figure>

panel in Fig. 2; Sun et al. 2018, 2020b; Rosolowsky et al. 2021). On the other hand, SF-events experiences an increase in $\alpha_{vir}$ due to the fact that less gas is compressed to dense states, creating longer $t_{ff}$ on average. As we mentioned previously, if we take SF-events to resemble dense gas tracers, its $\alpha_{vir}$ values overestimate the HCN-based virial parameters (Gallagher et al. 2018a; Jiménez-Donaire et al. 2019; Neumann et al. 2023), which agree better with CO-based observations.

Additionally, we have checked that our results are robust to small changes to the spatial scale by repeating the analysis using 100 pc cubes, with no change to the main conclusions drawn here.

## 3.3 Implications for the local $\epsilon_{ff}$

We next turn to the resulting $\epsilon_{ff}$ distributions for the SF-events gas. The left most panel in Fig. 3 presents the evolution of the $\epsilon_{ff}$ and local depletion time PDFs with redshift. Given that the PN12 model is $\alpha_{vir}$-dependent, a weakly increasing virial parameter with time (solid black in Fig. 2) leads to a weak decrease of $\epsilon_{ff}$. The mass-weighted mean of $\epsilon_{ff}$ is approximately 1 per cent, in agreement with the universally observed and predicted values (Krumholz et al. 2012; Lee, Miville-Deschênes & Murray 2016; Semenov et al. 2016; Utomo et al. 2018; Grisdale et al. 2019; Semenov, Conroy & Hernquist 2024b, see also Polzin et al. 2024 where $\langle\epsilon_{ff}\rangle \approx 1$ per cent for a range of metallicities). We can rewrite the local star formation law in equation (1) as $\rho_g/\dot{\rho}_* = t_{ff}/\epsilon_{ff}$, where $t_{sf} = t_{ff}/\epsilon_{ff}$ is the local gas depletion time. Then, in our case where the gas upon star formation in our simulation features a nearly constant $\epsilon_{ff}$, $t_{sf} \propto t_{ff}$. This means that changes to the local rate of star formation are mainly driven by changes in the density PDF (as $t_{ff} \propto \rho_g^{-0.5}$). The right most panel of Fig. 3 shows the evolution of $t_{sf}$. During the starburst and LMM epochs, both $t_{sf}$ and $t_{ff}$ (third row of Fig. 2) shifts to shorter time-scales, which is coeval with the drop of global depletion time values seen in Segovia Otero et al. (2022), caused by tidal compression and shocks during galaxy interactions (Renaud et al. 2022). In this

figure, the role of increased levels of ISM turbulence is therefore not to impact the distribution of $\epsilon_{ff}$ on small scales, but rather to drive the ISM density PDF towards higher densities, which in turn leads to lower global gas depletion times (see also Kraljic et al. 2014, 2024). We note that the median $t_{sf}$ is an order of magnitude shorter ($\sim 100$ Myr) than the global depletion time $\tau_{dep}$ ($\sim 1$ Gyr) due to the fact that most gas is not star forming (Semenov et al. 2017; Polzin et al. 2024).

The shape of the distributions, along with the median values (markers on the bottom of Fig. 3), reveal a slight evolution with time: red distributions showing symbolically higher $\epsilon_{ff}$ medians at $z > 1$ (0.8 per cent) and a family of blue distributions with a factor of $\sim 3$ lower $\epsilon_{ff}$ after the LMM event (0.3 per cent). Given the span of the $\epsilon_{ff}$ PDFs (0.001 per cent–30 per cent) and the more drastic evolution of $t_{ff}$ and $t_{turb}$, we consider $\epsilon_{ff}$ to be roughly stationary. However, this change in efficiency is enough for parts of the ISM in the starburst and LMM periods to feature $\epsilon_{ff}$ above $\sim 20$ per cent. These star-forming regions contain the largest amount of dense turbulent gas, but as can be seen from the PDFs, they make out a small fraction of the total gas mass of the galaxy during that epoch ($\ge 2$ per cent during the starburst phase compared to $< 0.5$ per cent for the secular and early epochs). They also correspond to the shortest time-scales ($< 1$ Myr), potentially making them ideal locations for massive star cluster formation (Li, Gnedin & Gnedin 2018; Li et al. 2020, 2022).

In Section 4, we discuss how alternatives to PN12 could affect our conclusions.

## 3.4 Inferring $\epsilon_{ff}$ from observations

A large number of studies have been undertaken to constrain $\epsilon_{ff}$ using Milky Way GMC- and cloud-scale observations of Local Universe galaxies (for reviews, see Krumholz et al. 2019; Schinnerer & Leroy 2024) as well as in gravitationally lensed high-z galaxies (e.g. Dessauges-Zavadsky et al. 2019, 2023). An apples-to-apples comparison with theory and simulations is not trivial (see Grudić et al. 2018, 2019, for high resolution GMC simulations, Grisdale et al. 2019 for simulations of isolated disc galaxies, Nuñez-Castiñeyra et al. 2021 for cosmological simulations, and Bemis & Wilson 2023 for a more observational perspective). Observationally, $\epsilon_{ff}$ quantifies the ratio between recent star formation and total mass budget of

factor *f* that corrects for an assumed radial density profile of CO inside their spherical GMCs so that $f = (1 - \gamma/3)(1 - 2\gamma/5)$, where $\gamma$ is the slope of the power law of the profile (Bertoldi & McKee 1992; Rosolowsky & Leroy 2006).
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Theoretical (top row) and observational (bottom row) estimates of the star formation efficiencies per free-fall time as a function of Σ, σ, and αvir. Only cold + dense gas properties in outputs at z < 1 are used. The running median of the two-dimensional histograms analysed on a grid of cubes of 100 pc in size is shown in blue, with its robust standard deviation as a blue shaded band (see Section 3.2). In solid green, we have the PHANGS-like running median also with its associated robust standard deviation band evaluated over on a grid of 1 kpc cubes (see Section 3.4 for a more thorough explanation on how to calculate PHANGS-like properties). In addition, the bottom row contains, in orange, a linear fit of the form y = m(x - x₀) + b where the coefficients take the values m = 0.00, 0.01, 0.10, b = -2.44, -2.44, -2.45, x₀ = 1.66, 0.74, 0.52, respectively from left to right (L25, Leroy et al. 2025). The orange band indicates the scatter around the fit with values of 1σ = 0.28, 0.28, 0.27. The bottom left panel shows a hatched region where VINTERGATAN is limited by the mass resolution of star particles when sampling the star formation model. **Takeaway:** star formation models based on the theory of isothermal gravoturbulent gas relaxes the apparent contradiction seen between models of star formation and star formation scaling relations observed in GMCs.</figcaption></figure>

GMCs, normalized by $t_{ff}$, and probed by molecular gas and SFR tracers. Therefore, a time-lag exists between the actual trigger of star formation in a specific cloud, and the eventual signal of recent star formation (YSOs or Hα, where the latter has an associated time-lag on the order of ~10 Myr). In contrast, in the context of theories of star formation in gravoturbulent isothermal gas, $\epsilon_{ff}$ measures the amount of gas mass that is turned into stars on a free-fall time-scale by integrating a log-normal density PDF above a certain density threshold. Such an instantaneous description of $\epsilon_{ff}$ is derived from the current density and turbulent structure of a piece of the ISM. This is not necessarily equal to an $\epsilon_{ff}$ inferred from gas tracers which encodes star formation imprinted at an earlier instance of a gas cloud.

In Fig. 4, we use VINTERGATAN outputs in the *secular* phase of evolution ($z < 1$) to compute theoretically and observationally motivated estimates of $\epsilon_{ff}$. The top row shows two-dimensional histograms with their running medians and robust standard deviations in blue corresponding to the *theoretical* star formation efficiency per free-fall time, $\epsilon_{ff,100}$, against the gas surface density, velocity dispersion, and virial parameter of *cold + dense* gas. This time, we calculate $\epsilon_{ff,100}$ on a grid of 100 pc-sized cubes using equation (2) from the velocity dispersion on these scales, and the density from the mass of *cold + dense* gas within each 100 pc cube. Note that these measurements are not necessarily the same as the on-the-fly data in Fig. 3 due to the larger sizes of the cubes with respect to the maximally resolved cells, and the inclusion of more diffuse and less turbulent gas (Fig. 3). We shift our analysis from 50 to 100 pc because we want to compare the predictions of the PN12 model

applied on 100 pc scales with the observational $\epsilon_{ff}$ estimates provided by the PHANGS collaboration, most of which are also on physical sizes of ~100 pc (Rosolowsky et al. 2021; Sun et al. 2022). The analysis behind Fig. 4 was also performed on 50 pc scales, but no major differences were observed.

The bottom row in Fig. 4 shows a similar analysis, only applying an *observational* estimator for the star formation efficiency per free-fall time (Leroy et al. 2017; Schinnerer & Leroy 2024), hereafter $\epsilon_{obs}$:

$$ \epsilon_{obs} = \frac{t_{ff}}{\tau_{dep}}. \qquad (5) $$

Here, $t_{ff}$ and $\tau_{dep}$ are the local free-fall and gas depletion time-scales computed from the *cold + dense* gas embedding stars younger than 10 Myr. Alternatives to equation (5) have been put forward by, e.g. Lee et al. (2016),⁷ but our findings and main trends are not affected by the change of $\epsilon_{obs}$ prescription.

It is clear from Fig. 4 that $\epsilon_{ff,100}$ and $\epsilon_{obs}$ are dramatically different measurements. Starting with the $\epsilon_{ff}-\Sigma$ relations, $\epsilon_{ff,100}$ increases with increasing $\Sigma$, as expected from a decreasing $\alpha_{vir}$ with increasing local density. The opposite trend is found for $\epsilon_{obs}$, which stems from dense gas often probing early stages of GMCs, when only low levels of star formation have occurred and before feedback has had time to disrupt them. This yields a low-mass ratio between young

$$ 7\epsilon_{obs} = \frac{t_{ff}}{t_{y,\star} M_{y,\star}} \frac{M_{y,\star}}{M_{g} + M_{y,\star}}, \text{ with } M_{y,\star} \text{ being the mass of stars younger than } t_{y,\star}. $$
stars and the total mass of clouds, hence a low $\epsilon_{\text{obs}}$. The high-$\epsilon_{\text{obs}}$, low-$\Sigma$ regions of the distribution tend to represent gas-poor clouds in their final stages of evolution, where large amounts of young stars drive large estimates of $\epsilon_{\text{obs}}$ (Feldmann & Gnedin 2011; Grudić et al. 2018; Grisdale et al. 2019; Grudić et al. 2019). In contrast, such diffuse environment makes the instantaneous $\epsilon_{\text{ff,100}}$ predict very low star formation efficiencies. We note however, that the picture is much more complex than this; for instance, as seen from the $\epsilon_{\text{obs}} - \alpha_{\text{vir}}$ relation, regions with apparent high $\epsilon_{\text{obs}}$ are in fact found to be both the most bound ($\alpha_{\text{vir}} < 1$) as well as unbound ($\alpha_{\text{vir}} > 10$), which indicates that not all regions are necessarily undergoing late stages of star formation and feedback disruption (for an example of diverse cloud evolutionary paths, see fig. 3 in Grisdale et al. 2019).

The $\epsilon_{\text{ff}}-\sigma$ relations predict a decreasing trend for both theoretical and observational efficiencies, yet additional efforts are required to precisely quantify whether such correlations are driven by the coupling between a multiphase ISM and different stellar feedback mechanisms (see Discussion section, but also Federrath 2015; Andersson et al. 2024) or whether more large-scale sources of turbulence ought to be considered, e.g. shear.<fnref n="8" />

To ascertain whether the star formation properties in VINTERGATAN at low redshifts are compatible with observations, we compare our results to those from the PHANGS collaboration. In particular, we follow the analysis outlined in Leroy et al. (2025) (hereafter L25). Cloud properties such as $\Sigma$, $\sigma$, and $\alpha_{\text{vir}}$, are measured on pixels of 50–150 pc in size further averaged over 0.5–1.5 kpc domains and weighted by the intensity of the CO (2-1) emission line (for more information on how to compute these parameters see Sun et al. 2022, 2023; Schinnerer & Leroy 2024). L25 evaluate $\epsilon_{\text{obs}}$ using equation (5), but $t_{\text{ff}}$ is now the CO luminosity-weighted average over the kpc domains and $\tau_{\text{dep}}$ is the ratio between the CO gas and SFR surface densities from the image at the native resolution convolved with a Gaussian beam to the desired kpc resolution (see appendix A in Sun et al. 2023; for improvements on the ‘aperture-scale averaged’ data products from Sun et al. 2022). Finally, the empirical relation seen in orange for all bottom panels in Fig. 4 with its $1\sigma$ scatter is the result of fitting a first-order polynomial between $\epsilon_{\text{obs}}$ and $\Sigma$, $\sigma$, $\alpha_{\text{vir}}$, respectively, and additionally bootstrapping the parent sample N times.

To closely mimic such analysis, each of the VINTERGATAN outputs at $z < 1$ are now split in 1 kpc-sized regions. We then find which of the 100 pc cubes are located within each kpc region and re-calculate the gas mass-weighted average over the kpc domain of $\Sigma$, $\sigma$, $\alpha_{\text{vir}}$, and $t_{\text{ff}}$. It is worth emphasizing that we do not evaluate mass-weighted averages of cell values on kpc scales, but rather use those calculated from the 100 pc-sized cubes. $\epsilon_{\text{obs}}$ is still computed using equation (5), but from a gas mass-weighted free-fall time, $\langle t_{\text{ff}} \rangle$ divided it by the gas depletion times of each kpc domain.

Solid green lines in Fig. 4 show the running median and associated scatter through these kpc domains. We see that VINTERGATAN and PHANGS agree well mostly at low $\epsilon_{\text{obs}}$ values, namely, in the high-$\Sigma$, high-$\sigma$ regime of the parameter space. The most striking result is appreciated in the last column of Fig. 4. While $\epsilon_{\text{ff,100}}$ naturally follows the PN12 star formation law with an exponential shape, the $\epsilon_{\text{obs}}$ estimator recovers an almost constant star formation efficiency

for virial parameters between $1 < \alpha_{\text{vir}} < 10$ in accordance with PHANGS analysis. While the match between $\epsilon_{\text{obs}}$ in VINTERGATAN and PHANGS is interesting, it remains to be seen if observed $\epsilon_{\text{obs}}$ trends can be used to constrain star formation models, or whether any signature of an underlying star formation law becomes washed out due to complex overlap of diverse cloud life cycles in a single ISM patch (see Grisdale et al. 2019). For instance, we note that the scatter in the simulation data is greater than that of the PHANGS fit. More analysis of simulations adopting a variety of star formation and feedback models is necessary to ascertain to what degree the scatter carries information of actual cloud diversity and ISM environment, as opposed to numerical aspects. Making use of outliers to cloud distribution (e.g. cloud mass) has the potential for providing more robust insights into the physics of star-forming clouds (Renaud, Agertz & Romeo 2024).

The match between VINTERGATAN and PHANGS at $\Sigma < 100 \text{ M}_\odot \text{ pc}^{-2}$ and $\sigma < 100 \text{ km s}^{-1}$ is poor. Part of the reason for this mismatch is inherent to VINTERGATAN and other cosmological simulations when sampling star formation using massive star particles, here with initial masses of $10^4 \text{ M}_\odot$. Such high masses cannot adequately sample low levels of star formation in low density gas, hence introducing a bias. To illustrate this we added a sampling limit<fnref n="9" /> with functional form $\epsilon_{\text{obs}} \propto 1/\Sigma^{1.5}$. This region of the parameter space is unreachable by our cosmological simulation and seen as a hatched grey region limited by a black dashed line in the bottom left panel of Fig. 4. Higher resolution simulations are necessary to better probe star formation in the low density parts of the ISM, which we leave for future work.

# 4 DISCUSSION

In this paper, we have studied how the density and velocity dispersion of the ISM in a Milky Way-mass galaxy evolve with cosmic time. Gravity and turbulence in gas clouds evolve conjointly, leading to an $\alpha_{\text{vir}}$ distribution that is almost constant in time. Applying the PN12 recipe for star formation, we get that the local theoretical efficiency is $\langle \epsilon_{\text{ff}} \rangle \approx 1$ per cent computed on 20 pc scales. This is in good agreement with observational efficiency estimators, e.g. in the PHANGS collaboration. In the following section, we discuss in-depth how mapping between the theoretical and observational efficiencies would be affected by the life cycle of clouds, merger-driven environments at high redshift, and alternative implementations of star formation.

## 4.1 Mapping between theoretical and observational efficiencies

Fig. 5 explores the mapping between $\epsilon_{\text{obs}}$ and $\epsilon_{\text{ff,100}}$ over time by presenting $\epsilon_{\text{obs}}-\epsilon_{\text{ff,100}}$ contours, and their corresponding one-dimensional histograms for the LMM event (left-hand side, $1.3 < z < 1.5$) and the secular epoch (right-hand side, $z < 1$). The distributions of counts appear as solid curves, gas mass-weighted distributions are illustrated in dotted lines, and stellar mass-weighted ones in dashed lines. Both main panels suggest a tentative anticorrelation between $\epsilon_{\text{obs}}$ and $\epsilon_{\text{ff,100}}$ (see also Khullar et al. 2024).

<footnote n="9">⁹The sampling limit is computed from equation (5), re-writing it as $\epsilon_{\text{obs}} = t_{\text{ff}}/\tau_{\text{dep}} = (t_{\text{ff}}/t_{y,\star}) (\Sigma_{y,\star}/\Sigma)$. Re-formulating the free-fall time as $t_{\text{ff}} \propto (\Sigma/L)^{-0.5}$ we get that $\epsilon_{\text{obs}} \propto (L^{0.5}/t_{y,\star}) (\Sigma_{y,\star}/\Sigma^{1.5})$. Here, $\Sigma_{y,\star}$ corresponds to a maximally resolved star particle ~ $10^4 \text{ M}_\odot$ inside a circle of radius equal to $L/2$ where $L = 100 \text{ pc}$, and $t_{y,\star}$ is the maximum age of young stars set to 10 Myr (Fig. 2).</footnote>

<footnote n="8">⁸Cold + dense gas in the secular phase of VINTERGATAN is often located within the inner half-mass radius of the galaxy, where galactic shear has been shown to play a crucial role in shaping the local star formation properties in the ISM of these regions ('the Brick'; Petkova et al. 2023).</footnote>
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Comparison between the theoretically (PN12) versus the observationally motivated efficiency (equation 5) for the LMM and the secular phases of VINTERGATAN. Both quantities are evaluated on 100 pc scales, with colour-coded contours representing the two-dimensional histogram between $\epsilon_{\text{obs}}$ and $\epsilon_{\text{ff,100}}$ indicating probability densities of 0.6, 0.3, 0.2, 0.1, 0.03, and 0.01 from most yellow to most purple, respectively. The corresponding one-dimensional histograms are placed on the top and sides of the main panels, where the solid curve representing the counts is weighted by stellar mass from particles younger than 10 Myr inside the 100 pc cubes (dashed) and by the enclosed gas mass (dotted). The diagonal dashed line crossing the two main panels just shows the one-to-one correlation between $\epsilon_{\text{obs}}$ and $\epsilon_{\text{ff,100}}$. For comparison, the robust scatter of $\epsilon_{\text{obs}}$ from Utomo et al. (2018) is highlighted as a green rectangle. Takeaway: $\epsilon_{\text{obs}}$ and $\epsilon_{\text{ff,100}}$ do not correlate with each other.</figcaption></figure>

Right-hand side panels in Fig. 5 show the $\epsilon_{\text{obs}}-\epsilon_{\text{ff,100}}$ relation for VINTERGATAN at $z < 1$. Surprisingly, an order of magnitude difference exists between a median $\epsilon_{\text{ff,100}}$ of ~4 per cent (right most panel) and a median $\epsilon_{\text{obs}}$ of ~0.4 per cent (top right panel) regardless of the distributions being stellar or gas mass-weighted. This suggests that: (i) from a PN12 point of view the ISM is considered to be more efficiently star-forming on 100 pc scales and free-fall time-scales of ~3.5 Myr than what $\epsilon_{\text{obs}}$ implies; (ii) the observational estimator $\epsilon_{\text{obs}}$ describes a less efficiently star-forming ISM with gas depletion times of ~300 Myr. From Fig. 4, we concluded that star-forming clouds in their initial stages of evolution are dense and turbulent, with gas masses larger than that of stars younger than 10 Myr, appearing inefficient when traced by observations yet efficient as predicted by theory. Contour lines in Fig. 5 do show a tendency towards the low-$\epsilon_{\text{ff,100}}$, high-$\epsilon_{\text{obs}}$ corner of the panel indicating such cloud evolution. Nevertheless, some of these early clouds are found in the central regions of VINTERGATAN and have the highest SFRs, so even after the onset of the first SNe, it is not clear whether feedback will be able to couple to such large cold and dense gas reservoirs to impact $\alpha_{\text{vir}}$ (Querejeta et al. 2019; see Hopkins et al. 2013b for feedback-regulated star formation, where SNe feedback weakly couples to cold and dense gas clouds with high HCN/CO compared to pre-SNe feedback; Ohlin, Renaud & Agertz 2019 for tests of SNe explosions in turbulent ISM boxes; and Zakardjian et al. 2023 for spatial offsets between HCN and CO peaks).

During the LMM phase, VINTERGATAN is undergoing a star-burst. Panels on the left-hand side of Fig. 5 illustrate how the

$\epsilon_{\text{obs}}-\epsilon_{\text{ff,100}}$ trend remains visible in a high-$z$ major merger, almost orthogonal to the one-to-one relation (black dashed line). If we look at the distribution of counts or gas mass-weighted $\epsilon_{\text{ff}}$ values (solid and dotted, respectively), there is a minor tendency towards lower $\epsilon_{\text{ff,100}}$ and higher $\epsilon_{\text{obs}}$ compared to the secular phase, with median $\epsilon_{\text{ff,100}} \sim \epsilon_{\text{obs}} \sim 1$ per cent (left most and top left panels). In turn, the stellar mass-weighted $\epsilon_{\text{obs}}$ distributions are significantly shifted towards high efficiencies of ~10 per cent (dashed lines). This is because shorter local gas depletion times of ~100 Myr (Fig. 2), imply that there is an enhancement of star formation on shorter time-scales when VINTERGATAN is observed as a starburst galaxy rather than a main-sequence galaxy. On 100 pc scales, this translates into larger stellar-to-gas mass ratios, which skews the stellar-mass weighted $\epsilon_{\text{obs}}$ distributions. Such efficiency values have been reported from YSOs (Young Stellar Objects) in dense clumps in the Milky Way (Heyer et al. 2016), dense gas tracers like HCN and CS from extragalactic sources (Wu et al. 2010), and from dense star-forming regions in starbursting high-redshift lensed galaxies (Dessauges-Zavadsky et al. 2019, 2023).

Interestingly, the $\epsilon_{\text{ff,100}}$ distribution indicates less efficient star formation during the LMM (left most panel) than in the secular phase (right most panel) due to enhanced $\alpha_{\text{vir}}$. This is at odds with Fig. 3, which characterizes VINTERGATAN as a marginally less efficient galaxy after the LMM as gas is no longer in high density states. Even though both $\epsilon_{\text{ff,100}}$ and $\epsilon_{\text{ff}}$ are computed using equation (2), differences arise due to the different scales in which these parameters
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Comparison among three $\epsilon_{ff}$ models and the evolution of their gas mass-weighted PDFs with redshift. While the solid lines outline the on-the-fly $\epsilon_{ff}$ PDFs in VINTERGATAN, the same $n$ and $\sigma$ from the SF-events gas is used as input in the other two models. No additional simulations have been run with different $\epsilon_{ff}$ models. From left to right we have: $\epsilon_{ff,PN12}$ solid lines represent the PN12 prescription, SFF<sub>KM05</sub> dashed distributions stand for the single free-fall model adopted by Federrath & Klessen (2012) from Krumholz & McKee (2005), and the MFF<sub>KM05</sub> in dotted PDFs correspond to the multi free-fall adaptation of the SFF<sub>KM05</sub> based on Hennebelle & Chabrier (2011). *Takeaway*: the same on-the-fly density and velocity dispersion values are used as inputs for all models, so it is just a prove of concept that, a priori, different models result in different $\epsilon_{ff}$ distributions.</figcaption></figure>

are computed, but most importantly due to the fact that SF-events and cold + dense describe inherently different states of gas. Namely, SF-events only include gas that has eventually formed a star, with moderately larger velocity dispersion values yet substantially denser gas. This makes up for an $\epsilon_{ff}$ evolution that follows changes in the gas density PDF.

Our results indicate that the multiscale aspect of deriving star formation efficiencies as well as the use of different SFR and molecular gas tracers ought to be considered. These will weigh $\epsilon_{ff}$ distributions differently and imply various degrees of feedback coupling to the surrounding ISM that must ultimately be taken into account when comparing star formation models to observational estimators (Utomo et al. 2018; Krumholz et al. 2019; Schinnerer & Leroy 2024, and references therein).

## 4.2 Alternative star formation efficiency models

The ISM in VINTERGATAN is kept in equilibrium, at a roughly constant $\alpha_{vir}$ distribution. Clouds with too low $\alpha_{vir}$ will be rapidly and efficiently turned into stars due to an exponentially increasing $\epsilon_{ff}$, while those with too large $\alpha_{vir}$ are likely at the end of their life cycle. In this section, we discuss how our findings could be affected when changing sub-grid recipes for star formation, and how can that modify the structure of the ISM. In order to illustrate the sensitivity of $\epsilon_{ff}$ to different star formation recipes, Fig. 6 shows three $\epsilon_{ff}$ models given a fixed gas distribution, i.e. for the same cell density and velocity dispersion input (a more rigorous comparison would imply rerunning VINTERGATAN with the mentioned models):

(i) $\epsilon_{ff,PN12}$: star formation prescription used in the VINTERGATAN simulation, presented in equations (1) and (2) (solid lines in Fig. 6). This model has been calibrated on high resolution simulations of magnetized ISM boxes (Padoan et al. 2012) that cover a range in virial parameters of $0.3 \ge \alpha_{vir} \ge 13$. The main take away from the PN12 model is that $\epsilon_{ff}$ exponentially decreases with increasing

$t_{ff}/t_{turb}$, and changes by less than a factor of two with magnetic field strength for typical star-forming regions.

(ii) **SFF<sub>KM05</sub>**: single free-fall models (SFF) define a family equations that explicitly integrate the log-normal distribution of gas densities to calculate the fraction of dense gas above a critical density that undergoes star formation on a fixed free-fall time-scale.¹⁰ Dashed curves in Fig. 6 outline the star formation model based on equation (20) in Krumholz & McKee (2005), and further adapted by Federrath & Klessen (2012) (see their table 1). In this model, the free-fall time is fixed to that of the average density of the cloud.

(iii) **MFF<sub>KM05</sub>**: multi-free-fall models (MFF) as dotted curves in Fig. 6 describe a family of star formation equation similar to SFF ones only that they allow the critical density threshold to vary given that gravitationally bound structures of different densities will collapse at different $t_{ff}$ (see Hennebelle & Chabrier 2011, and table 1 in Federrath & Klessen 2012).

In the future, we will evolve self-consistently these models using VINTERGATAN physics for strict comparison, but the main takeaway from Fig. 6 is that the choice of $\epsilon_{ff}$ prescription shifts the emerging $\langle\epsilon_{ff}\rangle$. From the $\alpha_{vir}$ distributions in Fig. 2 we see that a large fraction of the star-forming ISM has virial parameters that exceed those of the MHD ISM boxes used for calibrating $\epsilon_{ff,PN12}$, with Mach numbers approaching ~100. These are extreme conditions where $\alpha_{vir}$ is extrapolated beyond its limits underestimating $\epsilon_{ff,PN12}$ in comparison to the other two models. Both SFF<sub>KM05</sub> and MFF<sub>KM05</sub> are $\alpha_{vir}$ and $M$ dependent, meaning that for the same cell density input, additional dependencies on its velocity dispersion shift distributions towards efficiencies of 10 per cent and 100 per cent, respectively. That is particularly the case for the MFF<sub>KM05</sub> model, with an extra factor $\exp[(3/8)\sigma_s^2]$, where $\sigma_s = \ln(1+b^2M^2)$. In a recent study, Brucy, Hennebelle & Colman (2024) and Brucy et al. (2024) conclude that the density PDF at high Mach numbers can even depart from the often assumed log-normal PDF shapes (Castaing 1996; Hopkins 2013). Such high-resolution isothermal-ISM box simulations hold the promise of improving $\epsilon_{ff}$ predictions in high $M$ gas such as that found in VINTERGATAN (second row in Fig. 2).

The role of turbulence is critical in regulating star formation (Renaud, Kraljic & Bournaud 2012; Kraljic et al. 2024), but the details still remain unsettled. A variety of sub-grid recipes have been implemented to compensate for the effect of small-scale eddies unresolved in most galaxy simulations. Semenov et al. (2016) implemented a sub-grid scale turbulent model (SGS; Schmidt & Federrath 2011; Schmidt et al. 2014) where the production and dissipation of unresolved turbulent energy is followed more explicitly. Their isolated galaxy simulation well reproduces observed $\epsilon_{ff}$ distributions and the KS relation, and has been further extended to include a more realistic cosmological environment also recovering a near-universal $\epsilon_{ff}$ of ~1 per cent (Semenov et al. 2024a, b). Nevertheless, determining the exact role that more physically informed turbulence models play is challenging. Despite VINTERGATAN not including

<footnote>$^{10}\epsilon_{ff} = \frac{\epsilon_w}{\phi_t} \int_{s_{crit}}^{\infty} \frac{t_{ff}(\rho_0) \rho}{t_{ff}(\rho) \rho_0} p(s) \, ds$, with $p(s)$ being the log-normal distribution of $s = \ln(\rho/\rho_0)$, $s_{crit}$ the critical density equal to $(\pi^2/5)\phi_x^2\alpha_{vir}M^2$ for the hydrodynamics-only case in both the SFF<sub>KM05</sub> and MFF<sub>KM05</sub> models (Federrath & Klessen 2012), and $t_{ff}(\rho)$ is the free-fall time-scale set to $t_{ff}(\rho_0)$ in SFF<sub>KM05</sub> but kept inside the integral for MFF<sub>KM05</sub>. Both $\phi_t$ and $\phi_x$ are calibration factors set to $1/\phi_t = 3$ and $\phi_x = 0.12$ for SFF<sub>KM05</sub> and $1/\phi_t = 0.49$ and $\phi_x = 0.19$ for MFF<sub>KM05</sub> in Fig. 6. The turbulent forcing parameter is $b_{turb} = 0.38$, representing a mix of solenoidal and compressive modes of turbulence.</footnote>
an SGS model, average ISM densities and velocity dispersions in Semenov et al. (2024b), are consistent with our results, even when the sourcing for unresolved turbulence is changed (see Semenov 2024, for turbulence based on local numerical dissipation). The ISM can even develop order of magnitude differences in velocity dispersion values, only to self-regulate and maintain the same $\epsilon_{ff}$ distribution (see Kretschmer & Teyssier 2020, for shear- versus SNe-induced turbulence).

Evolving the ISM with the alternative $\epsilon_{ff}$ models such as those in Fig. 6 must go hand-in-hand with a tailored feedback recipe. In the FIRE (Feedback in Realistic Environments) collaboration $\epsilon_{ff}$ is assumed to be 100 per cent, making feedback the main agent regulating star formation (Faucher-Giguère, Quataert & Hopkins 2013; Hopkins et al. 2014; Orr et al. 2018; Khullar et al. 2024). At sufficiently high resolution gas collapses on exceedingly short $t_{ff}$, forms stars on local time-scales of $t_{sf} = t_{ff}/\epsilon_{ff}$, with clustered feedback quickly and efficiently dispersing the clouds (see discussion in Hopkins, Quataert & Murray 2012). Hopkins, Quataert & Murray (2011) argue that a lower value of $\epsilon_{ff}$ simply allows for gravitational collapse to proceed for longer. The actual choice of $\epsilon_{ff}$ is deemed unimportant at such high densities as it is just a matter of reaching the necessary feedback budget to disperse the cloud. When feedback is this efficient, changes to $\epsilon_{ff}$ or other star formation criteria lead to little-to-no impact on galactic-scale integrated quantities like stellar masses or SFRs (e.g. Hopkins et al. 2018). When feedback is less efficient, it is $\epsilon_{ff}$ that regulates star formation (Semenov et al. 2017) causing differences in global properties such as SFRs (Benincasa et al. 2016), the normalization of the KS relation (Agertz et al. 2013), bulge-to-disc mass ratios (Gensior, Kruijssen & Keller 2020), or morphology of H I discs (Gensior et al. 2024).

In both feedback efficiency regimes the emerging structure of the ISM varies significantly as seen in density and temperature distributions (e.g. Hopkins et al. 2013a), nebular emission line ratios (see Katz et al. 2024, for radiative transfer cosmological simulations with on the fly non-equilibrium chemistry), or H$_2$ depletion time-scales on GMC scales (Semenov et al. 2018). Nevertheless, explicitly evolving the star formation-feedback cycle with various $\epsilon_{ff}$ and feedback recipes can render resembling ISM structures and GMC properties. For instance, enhancing the SNe energy output (Katz et al. 2024), switching $\epsilon_{ff}$ strategies from a constant value to a multi-freefall one (Kretschmer & Teyssier 2020), strengthening protostellar feedback, modifying feedback models (see Nuñez-Castiñeyra et al. 2021, for delayed cooling versus mechanical feedback prescriptions) or star formation schemes (see Kang et al. 2025, for a sink particle-versus turbulence-based star formation criteria) all produce a clumpy ISM with young and bright stellar populations forming from denser gas. We ultimately want to understand what determines how much gas mass and on which time-scales it piles up at the threshold of star formation. Existing and new observations of GMCs in the space defined by Fig. 4 along with additional diagnostics like those mentioned previously stand as stringent probes to calibrate the feedback-star formation balance in simulations.

## 5 CONCLUSIONS

In this paper, we analyse the density and velocity dispersion of the ISM and its implications on the local star formation efficiency per free-fall time in VINTERGATAN, a high resolution cosmological zoom-in simulation of a Milky Way-type galaxy (Agertz et al. 2021). We follow the evolution of their PDFs in the context of JWST observations (Heintz et al. 2022; Aravena et al. 2023), lensed galaxies and mergers around the cosmic peak of star formation (Dessauges-

Zavadsky et al. 2019, 2023; Rizzo et al. 2023; Roman-Oliveira et al. 2023), and GMC observables at $z \approx 0$ from the PHANGS collaboration (Sun et al. 2018, 2020b; Rosolowsky et al. 2021; Schinnerer & Leroy 2024; L25). Our findings can be summarized in four points:

(i) **The density and turbulent structure of the star forming ISM on 50 pc scales reacts to its extragalactic environment** (Fig. 1). Three cosmic epochs are identified: (i) *early on* ($5 < z < 9$), where the gas-rich protogalaxy is being assembled mostly from *in-situ* star formation and minor contributions from neighbouring galaxies that sustain moderate velocity dispersions levels ($\sigma < 30$ km s$^{-1}$; Pope et al. 2023; Posses et al. 2023); (ii) a *starburst* epoch ($1 < z < 5$) that ends with the *last major merger* event ($1.3 < z < 1.5$), where the galaxy builds a rotationally supported disc and experiences merger-driven starbursts that compress the ISM to $\sigma \approx 100$ km s$^{-1}$ (Herrera-Camus et al. 2022; Parlanti et al. 2023) with global gas depletion times of $\sim 0.1$ Gyr due to an excess of gas at high densities; and (iii) a final *secular* phase where VINTERGATAN settles into an extended main-sequence galaxy where $\sigma \approx 10$ km s$^{-1}$ (Tamburro et al. 2009).

(ii) **The density and velocity dispersion distributions on 50 pc scales evolve with cosmic time, but conspire to yield virial parameter distributions that are nearly constant** (Fig. 2). The resulting ISM at low redshift has virial parameter distributions within $1 < \alpha_{vir} < 10$, values that are in agreement with observations (Leroy et al. 2015; Miville-Deschênes et al. 2017; Dessauges-Zavadsky et al. 2019; Rosolowsky et al. 2021). The more exposed the ISM is to recent star formation, and therefore affected by the subsequent stellar feedback, the denser and more turbulent it is. The presence of mergers, or absence thereof, pushes the ISM to such extreme conditions, nevertheless, their exact role in shaping $\alpha_{vir}$ distributions is still unclear.

(iii) **The emerging $\epsilon_{ff}$ distribution spans several orders of magnitude, from 0.001 per cent up to 30 per cent, but averages to a roughly constant ($\epsilon_{ff} \approx 1$ per cent across cosmic time)** (Fig. 3) in line with observations (Lee et al. 2016; Utomo et al. 2018; Krumholz et al. 2019) and recent computationally driven studies (Semenov et al. 2016; Kretschmer & Teyssier 2020; Nuñez-Castiñeyra et al. 2021; Polzin et al. 2024). The model adopted in VINTERGATAN is based on PN12 and is exclusively $\alpha_{vir}$-dependent. Thus, for a constant $\epsilon_{ff}$ the local gas depletion time-scales as $t_{sf} \sim t_{ff}/\epsilon_{ff} \sim t_{ff} \propto \rho_g^{-0.5}$. As such, changes in the gas density PDFs become the main controlling factor in setting the local rate of star formation. The role of turbulence is to instead establish the density PDF in the first place (Kraljic et al. 2014, 2024; Renaud et al. 2014; Segovia Otero et al. 2022), rather than inducing large changes to $\epsilon_{ff}$.

(iv) **Theoretical, more instantaneous star formation efficiencies are inherently different to their observational estimators as they represent GMCs in different physical states** (Fig. 4). Different star-forming gas tracers (CO, HCN, CS, etc.) probe a range of densities and temperatures. Tracers of recent star formation (YSOs, H$\alpha$, UV, etc.) also provide information about the star-forming cloud, but it must be *integrated* over the lifetime of the tracer employed (Grudić et al. 2019). This often hinders us from accessing the local instantaneous conditions for star formation. Fig. 5 hints towards a possible anticorrelation between the theoretical and observational $\epsilon_{ff}$, proposing that clouds in their initial stages of evolution begin with high $\epsilon_{ff,100}$ and low $\epsilon_{obs}$, and end their lives with low $\epsilon_{ff,100}$ and high $\epsilon_{obs}$. Our conclusions cast light on to seemingly contradictory results between star formation models and observations of star-forming GMC (Schruba, Kruijssen & Leroy 2019). This behaviour
is consistent with PHANGS observations (Leroy et al. 2017, L25; Schinnerer & Leroy 2024).

Linking local star formation parameters from ISM properties is highly non-linear and complex and has to rely on well calibrated sub-grid models (Nobels et al. 2024). These are rapidly improving with time (see Girma & Teyssier 2024, who include magnetic fields and Polzin et al. 2024 who add a metallicity dependency), but still do not consider the effects of the larger-scale environment (e.g. shear, compression, etc.).

In summary, star formation models, while degenerate, can be singled out when comparing to observed ISM properties. Moreover, information can be obtained from the scatter in GMC scaling relations, but how exactly it can help constrain star formation models is yet to be pinned down. In our work, we stress the need for understanding the coupling among star formation, feedback, and a multiphase ISM. Thus, a change in star formation models must go hand-in-hand with a change in the feedback recipe. Further observational data will be soon available in the advent of upcoming JWST and ALMA (Atacama Large Millimetre Array) observations, and the further development of dedicated surveys like PHANGS and MaNGA (Mapping Nearby Galaxies at Apache Point Observatory), generating multiscale multiwavelength data products of cloud properties.

## ACKNOWLEDGEMENTS

We thank the anonymous referee for a helpful and constructive report. We are also deeply grateful to Adam K. Leroy for providing the empirical data from the PHANGS Collaboration and to Marta Reina Campos for insightful discussions. OA and ASO acknowledge support from the Knut and Alice Wallenberg Foundation, the Swedish Research Council (grant 2019–04659), and the Swedish National Space Agency (SNSA Dnr 2023–00164). FR acknowledges support provided by the University of Strasbourg Institute for Advanced Study (USIAS), within the French national program Investment for the Future (Excellence Initiative) IdEx-Unistra. Support for VS was provided by Harvard University through the Institute for Theory and Computation Fellowship.

## DATA AVAILABILITY

The data underlying this article will be shared on reasonable request to the corresponding author.

## REFERENCES

* Agertz O., Kravtsov A. V., 2015, ApJ, 804, 18

* Agertz O., Teyssier R., Moore B., 2011, MNRAS, 410, 1391

* Agertz O., Kravtsov A. V., Leitner S. N., Gnedin N. Y., 2013, ApJ, 770, 25

* Agertz O. et al., 2020, MNRAS, 491, 1656

* Agertz O. et al., 2021, MNRAS, 503, 5826

* Andersson E. P., Mac Low M.-M., Agertz O., Renaud F., Li H., 2024, A&A, 681, A28

* Aravena M. et al., 2023, A&A, 682, A24

* Aubert D., Teyssier R., 2010, ApJ, 724, 244

* Bemis A. R., Wilson C. D., 2023, ApJ, 945, 42

* Benincasa S. M., Wadsley J., Couchman H. M. P., Keller B. W., 2016, MNRAS, 462, 3053

* Bertoldi F., McKee C. F., 1992, ApJ, 395, 140

* Bolatto A. D., Leroy A. K., Rosolowsky E., Walter F., Blitz L., 2008, ApJ, 686, 948

* Brinchmann J., Charlot S., White S. D. M., Tremonti C., Kauffmann G., Heckman T., Brinkmann J., 2004, MNRAS, 351, 1151

* Brucy N., Hennebelle P., Colman T., Klessen R. S., Le Yhuelic C., 2024, A&A, 690, A44

* Hennebelle P., Brucy N., Colman T., 2024, A&A, 690, A43

* Burkhart B., 2018, ApJ, 863, 118

* Castaing B., 1996, J. Phys. II, 6, 105

* Chabrier G., 2003, PASP, 115, 763

* Ciesla L. et al., 2023, A&A, 672, A191

* Courty S., Alimi J. M., 2004, A&A, 416, 875

* Daddi E. et al., 2010, ApJ, 714, L118

* Dessauges-Zavadsky M. et al., 2019, Nat. Astron., 3, 1115

* Dessauges-Zavadsky M. et al., 2023, MNRAS, 519, 6222

* Duncan K. et al., 2019, ApJ, 876, 110

* Eibensteiner C. et al., 2023, A&A, 675, A37

* Ejdetjärn T., Agertz O., Östlin G., Renaud F., Romeo A. B., 2022, MNRAS, 514, 480

* Elmegreen B. G., 1989, ApJ, 338, 178

* Faucher-Giguère C.-A., Quataert E., Hopkins P. F., 2013, MNRAS, 433, 1970

* Federrath C., 2015, MNRAS, 450, 4035

* Federrath C., Klessen R. S., 2012, ApJ, 761, 156

* Federrath C., Klessen R. S., 2013, ApJ, 763, 51

* Federrath C., Roman-Duval J., Klessen R. S., Schmidt W., Mac Low M. M., 2010, A&A, 512, A81

* Feldmann R., Gnedin N. Y., 2011, ApJ, 727, L12

* Förster Schreiber N. M., Wuyts S., 2020, ARA&A, 58, 661

* Fudamoto Y. et al., 2022, ApJ, 934, 144

* Gallagher M. J. et al., 2018a, ApJ, 858, 90

* Gallagher M. J. et al., 2018b, ApJ, 868, L38

* Gensior J., Kruijssen J. M. D., Keller B. W., 2020, MNRAS, 495, 199

* Gensior J. et al., 2024, MNRAS, 531, 1158

* Girma E., Teyssier R., 2024, MNRAS, 527, 6779

* Grisdale K., Agertz O., Romeo A. B., Renaud F., Read J. I., 2017, MNRAS, 466, 1093

* Grisdale K., Agertz O., Renaud F., Romeo A. B., 2018, MNRAS, 479, 3167

* Grisdale K., Agertz O., Renaud F., Romeo A. B., Devriendt J., Slyz A., 2019, MNRAS, 486, 5482

* Grudić M. Y., Hopkins P. F., Faucher-Giguère C.-A., Quataert E., Murray N., Kereš D., 2018, MNRAS, 475, 3511

* Grudić M. Y., Hopkins P. F., Lee E. J., Murray N., Faucher-Giguère C.-A., Johnson L. C., 2019, MNRAS, 488, 1501

* Haardt F., Madau P., 1996, ApJ, 461, 20

* Hahn O., Abel T., 2011, MNRAS, 415, 2101

* He H., Bottrell C., Wilson C., Moreno J., Burkhart B., Hayward C. C., Hernquist L., Twum A., 2023, ApJ, 950, 56

* Heintz K. E. et al., 2022, ApJ, 934, L27

* Hennebelle P., Chabrier G., 2011, ApJ, 743, L29

* Herrera-Camus R. et al., 2022, A&A, 665, L8

* Heyer M., Krawczyk C., Duval J., Jackson J. M., 2009, ApJ, 699, 1092

* Heyer M., Guttermuth R., Urquhart J. S., Csengeri T., Wienen M., Leurini S., Menten K., Wyrowski F., 2016, A&A, 588, A29

* Hopkins P. F., 2013, MNRAS, 430, 1880

* Hopkins P. F., Quataert E., Murray N., 2011, MNRAS, 417, 950

* Hopkins P. F., Quataert E., Murray N., 2012, MNRAS, 421, 3488

* Hopkins P. F., Narayanan D., Murray N., 2013a, MNRAS, 432, 2647

* Hopkins P. F., Narayanan D., Murray N., Quataert E., 2013b, MNRAS, 433, 69

* Hopkins P. F., Kereš D., Oñorbe J., Faucher-Giguère C.-A., Quataert E., Murray N., Bullock J. S., 2014, MNRAS, 445, 581

* Hopkins P. F. et al., 2018, MNRAS, 480, 800

* Jiménez-Donaire M. J. et al., 2019, ApJ, 880, 127

* Kang C., Kimm T., Han D., Katz H., Devriendt J., Slyz A., Teyssier R., 2025, A&A, 693, A149

* Katz H., Rey M. P., Cadiou C., Kimm T., Agertz O., 2024, preprint (arXiv:2411.07282)

* Kennicutt R. C. J., 1998, ApJ, 498, 541

* Khullar S., Matzner C. D., Murray N., Grudić M. Y., Guszejnov D., Wetzel A., Hopkins P. F., 2024, ApJ, 973, 40

* Kim C.-G., Ostriker E. C., 2015, ApJ, 802, 99

* Kim J.-h. et al., 2014, ApJS, 210, 14
Kim J.-h. et al., 2016, ApJ, 833, 202

Klessen R. S., Glover S. C. O., 2016, *Saas-Fee Advanced Course*, 43, 85

Krahm G., Finn M. K., Indebetouw R., Johnson K. E., Kamenetzky J., Bemis A., 2024, ApJ, 964, 166

Kraljic K., Renaud F., Bournaud F., Combes F., Elmegreen B., Emsellem E., Teyssier R., 2014, ApJ, 784, 112

Kraljic K. et al., 2024, A&A, 682, A50

Kretschmer M., Teyssier R., 2020, MNRAS, 492, 1385

Krumholz M. R., McKee C. F., 2005, ApJ, 630, 250

Krumholz M. R., Dekel A., McKee C. F., 2012, ApJ, 745, 69

Krumholz M. R., McKee C. F., Bland-Hawthorn J., 2019, ARA&A, 57, 227

Larson R. B., 1981, MNRAS, 194, 809

Lee E. J., Miville-Deschênes M.-A., Murray N. W., 2016, ApJ, 833, 229

Leethochawalit N. et al., 2023, ApJ, 942, L26

Leitherer C. et al., 1999, ApJS, 123, 3

Leroy A. K. et al., 2013, AJ, 146, 19

Leroy A. K. et al., 2015, ApJ, 801, 25

Leroy A. K. et al., 2016, ApJ, 831, 16

Leroy A. K. et al., 2017, ApJ, 846, 71

Leroy A. K. et al., 2025, preprint (arXiv:2502.04481) (L25)

Li H., Gnedin O. Y., Gnedin N. Y., 2018, ApJ, 861, 107

Li H., Vogelsberger M., Marinacci F., Sales L. V., Torrey P., 2020, MNRAS, 499, 5862

Li H., Vogelsberger M., Bryan G. L., Marinacci F., Sales L. V., Torrey P., 2022, MNRAS, 514, 265

Mac Low M.-M., 1999, ApJ, 524, 169

Mac Low M.-M., Klessen R. S., Burkert A., Smith M. D., 1998, Phys. Rev. Lett., 80, 2754

McKee C. F., Ostriker E. C., 2007, ARA&A, 45, 565

Miville-Deschênes M.-A., Murray N., Lee E. J., 2017, ApJ, 834, 57

Müller J., 2000, J. Res. Natl. Inst. Stand. Technol., 105, 551

Neumann L. et al., 2023, MNRAS, 521, 3348

Nobels F. S. J., Schaye J., Schaller M., Ploeckinger S., Chaikin E., Richings A. J., 2024, MNRAS, 532, 3299

Nuñez-Castiñeyra A., Nezri E., Devriendt J., Teyssier R., 2021, MNRAS, 501, 62

Ohlin L., Renaud F., Agertz O., 2019, MNRAS, 485, 3887

Orr M. E. et al., 2018, MNRAS, 478, 3653

Padoan P., Nordlund Å., 2011, ApJ, 730, 40

Padoan P., Haugbølle T., Nordlund Å., 2012, ApJ, 759, L27 PN12

Parlanti E., Carniani S., Pallottini A., Cignoni M., Cresci G., Kohandel M., Mannucci F., Marconi A., 2023, A&A, 673, A153

Petkova M. A. et al., 2023, MNRAS, 525, 962

Polzin A., Kravtsov A. V., Semenov V. A., Gnedin N. Y., 2024, OJAp 7, 114

Pope A. et al., 2023, ApJ, 951, L46

Posses A. C. et al., 2023, A&A, 669, A46

Querejeta M. et al., 2019, A&A, 625, A19

Raiteri C. M., Villata M., Navarro J. F., 1996, A&A, 315, 105

Renaud F., Kraljic K., Bournaud F., 2012, ApJ, 760, L16

Renaud F. et al., 2013, MNRAS, 436, 1836

Renaud F., Bournaud F., Kraljic K., Duc P. A., 2014, MNRAS, 442, L33

Renaud F., Agertz O., Read J. I., Ryde N., Andersson E. P., Bensby T., Rey M. P., Feuillet D. K., 2021a, MNRAS, 503, 5846

Renaud F., Agertz O., Andersson E. P., Read J. I., Ryde N., Bensby T., Rey M. P., Feuillet D. K., 2021b, MNRAS, 503, 5868

Renaud F., Segovia Otero Å., Agertz O., 2022, MNRAS, 516, 4922

Renaud F., Agertz O., Romeo A. B., 2024, A&A, 687, A91

Rizzo F. et al., 2023, A&A, 679, A129

Roca-Fàbrega S. et al., 2021, ApJ, 917, 64

Rodighiero G. et al., 2011, ApJ, 739, L40

Roman-Oliveira F., Fraternali F., Rizzo F., 2023, MNRAS, 521, 1045

Romeo A. B., Agertz O., Renaud F., 2023, MNRAS, 518, 1002

Rosen A., Bregman J. N., 1995, ApJ, 440, 634

Rosolowsky E., Leroy A., 2006, PASP, 118, 590

Rosolowsky E. et al., 2021, MNRAS, 502, 1218

Saintonge A., Catinella B., 2022, ARA&A, 60, 319

Schinnerer E., Leroy A. K., 2024, ARA&A, 62, 9

Schmidt W., Federrath C., 2011, A&A, 528, A106

Schmidt W. et al., 2014, MNRAS, 440, 3051

Scholte D. et al., 2024, MNRAS, 535, 3

Schruba A., Kruijssen J. M. D., Leroy A. K., 2019, ApJ, 883, 2

Segovia Otero Å., Renaud F., Agertz O., 2022, MNRAS, 516, 2272

Semenov V. A., 2024, preprint (arXiv:2410.23339)

Semenov V. A., Kravtsov A. V., Gnedin N. Y., 2016, ApJ, 826, 200

Semenov V. A., Kravtsov A. V., Gnedin N. Y., 2017, ApJ, 845, 133

Semenov V. A., Kravtsov A. V., Gnedin N. Y., 2018, ApJ, 861, 4

Semenov V. A., Conroy C., Smith A., Puchwein E., Hernquist L., 2024a, preprint (arXiv:2409.18173)

Semenov V. A., Conroy C., Hernquist L., 2024b, preprint (arXiv:2410.09205)

Solomon P. M., Rivolo A. R., Barrett J., Yahil A., 1987, ApJ, 319, 730

Speagle J. S., Steinhardt C. L., Capak P. L., Silverman J. D., 2014, ApJS, 214, 15

Stone J. M., Ostriker E. C., Gammie C. F., 1998, ApJ, 508, L99

Sun J. et al., 2018, ApJ, 860, 172

Sun J. et al., 2020a, ApJ, 892, 148

Sun J. et al., 2020b, ApJ, 901, L8

Sun J. et al., 2022, AJ, 164, 43

Sun J. et al., 2023, ApJ, 945, L19

Sutherland R. S., Dopita M. A., 1993, ApJS, 88, 253

Tacconi L. J. et al., 2018, ApJ, 853, 179

Tacconi L. J., Genzel R., Sternberg A., 2020, ARA&A, 58, 157

Tamburro D., Rix H. W., Leroy A. K., Mac Low M. M., Walter F., Kennicutt R. C., Brinks E., de Blok W. J. G., 2009, AJ, 137, 4424

Teyssier R., 2002, A&A, 385, 337

Trebitsch M., Blaizot J., Rosdahl J., Devriendt J., Slyz A., 2017, MNRAS, 470, 224

Trussler J. A. A. et al., 2023, MNRAS, 523, 3423

Übler H. et al., 2019, ApJ, 880, 48

Utomo D. et al., 2018, ApJ, 861, L18

Wei L. H., Keto E., Ho L. C., 2012, ApJ, 750, 136

Wetzel A. et al., 2023, ApJS, 265, 44

Wilson C. D., Bemis A., Ledger B., Klimi O., 2023, MNRAS, 521, 717

Wise J. H., Turk M. J., Norman M. L., Abel T., 2012, ApJ, 745, 50

Wisnioski E. et al., 2015, ApJ, 799, 209

Woosley S. E., Heger A., 2007, Phys. Rep., 442, 269

Wu J., Evans Neal J. I., Shirley Y. L., Knez C., 2010, ApJS, 188, 313

Wuyts S. et al., 2011, ApJ, 742, 96

Zakardjian A. et al., 2023, A&A, 678, A171

This paper has been typeset from a *TEX/LATEX* file prepared by the author.