Emergence and cosmic evolution of the Kennicutt-Schmidt relation driven by interstellar turbulence

The scaling relations between the gas content and star formation rate of galaxies provide useful insights into processes governing their formation and evolution. We investigate the emergence and the physical drivers of the global Kennicutt-Schmidt (KS) relation at $0.25 \leq z \leq 4$ in the cosmological hydrodynamic simulation NewHorizon capturing the evolution of a few hundred galaxies with a resolution of $\sim$ 40 pc. The details of this relation vary strongly with the stellar mass of galaxies and the redshift. A power-law relation $\Sigma_{\rm SFR} \propto \Sigma_{\rm gas}^{a}$ with $a \approx 1.4$, like that found empirically, emerges at $z \approx 2 - 3$ for the most massive half of the galaxy population. However, no such convergence is found in the lower-mass galaxies, for which the relation gets shallower with decreasing redshift. At the galactic scale, the star formation activity correlates with the level of turbulence of the interstellar medium, quantified by the Mach number, rather than with the gas fraction (neutral or molecular), confirming previous works. With decreasing redshift, the number of outliers with short depletion times diminishes, reducing the scatter of the KS relation, while the overall population of galaxies shifts toward low densities. Using pc-scale star formation models calibrated with local Universe physics, our results demonstrate that the cosmological evolution of the environmental and intrinsic conditions conspire to converge towards a significant and detectable imprint in galactic-scale observables, in their scaling relations, and in their reduced scatter.


Introduction
Decades of observational works have highlighted scaling relations between the gas content of galaxies and their star formation rate (SFR), as a key ingredient of galaxy formation.The pioneering work of Schmidt (1959) revealed a tight relation between the densities of gas and SFR.It has later been complemented by Kennicutt (1989), which proposed an empirical relation between the surface densities of neutral gas and SFR, of the form Σ SFR ∝ Σ a gas with an index a ≈ 1.4 for local star-forming galaxies.Since, a number of studies have extended the range of this Kennicutt-Schmidt (KS) relation by considering a more diverse population of galaxies like local starbursts (e.g.Kennicutt 1998), high redshift discs (e.g.Tacconi et al. 2010), submillimeter galaxies (e.g.Bouché et al. 2007), and sub-galactic scales in local galaxies (e.g.Bigiel et al. 2008) to name a few (see Kennicutt & Evans 2012 for a review).
The inferred slope of ∼ 1.4 for nearby normal spirals (Kennicutt 1989), recently confirmed by de los Reyes & Kennicutt (2019) in their revisited analysis, is consistent with measurements at higher redshifts (z = 1.5, Daddi et al. 2010).The slope of 1.4 − 1.5 has also been found for the combined sample of normal and starbursting local galaxies (Kennicutt & De Los Reyes 2021).On the other hand, dwarf galaxies yield slopes closer to unity (e.g.Filho et al. 2016;Roychowdhury et al. 2017), such that including them in the samples lowers the slope to ∼ 1.3 and increases the scatter of the KS relation (de los Reyes & Kennicutt 2019).Starburst galaxies alone appear to result in different values of a slope for different samples.For instance, Kennicutt & De Los Reyes (2021) suggests values of 1 − 1.2, shallower than previous findings (≈ 1.3 − 1.4, see Kennicutt 1989;Daddi et al. 2010).However, there seems to be a consistency in the findings of evidence for a bimodal (or even multimodal) relation for star-bursts and non-starbursting galaxies, yet with significant overlap (e.g.Daddi et al. 2010;Genzel et al. 2010;Kennicutt & De Los Reyes 2021).
As the neutral gas phase also includes diffuse atomic gas, that is yet to collapse, the SFR correlates more strongly with the molecular gas contents alone (e.g.Bigiel et al. 2011).This motivated the introduction of another KS-like relation, molecular one, with a slope empirically found to be close to unity in nearby star-forming galaxies both on galactic (e.g.Kennicutt 1998;Liu et al. 2015;de los Reyes & Kennicutt 2019) and sub-kpc scales (e.g.Bigiel et al. 2008;Leroy et al. 2008Leroy et al. , 2013;;Onodera et al. 2010;Schruba et al. 2011;Sun et al. 2023), in local (e.g.Liu et al. 2015) and high-redshift starbursts (e.g.Sharon et al. 2013;Rawle et al. 2014), and high-redshift galaxies both in galaxyaveraged (e.g.Genzel et al. 2010;Tacconi et al. 2013;Freundlich et al. 2019) and spatially resolved studies (e.g.Freundlich et al. 2013;Genzel et al. 2013).Considering the molecular gas alone reduces the variations in the measured slopes across these families of galaxies.
This wealth of observational studies comes with a vast diversity of resolutions, scales, tracers, conversion factors, and fitting methods, which make comparisons and compilations delicate (see e.g. de los Reyes & Kennicutt 2019).For instance, Sun et al. (2023) found that systematic uncertainties in the estimation of slopes due to different choices of SFR calibrations (see also e.g.Genzel et al. 2013, on the impact of the adopted extinction model) may be of about 10% to 15%, while CO-to-H 2 conversion factor may produce an additional 20% to 25% (in qualitative agreement with e.g.Liu et al. 2015;de los Reyes & Kennicutt 2019).Despite uncertainties on the values to adopt for such conversion factors (Bolatto et al. 2013), both observations and simulations report non-negligible variations across galactic disks (Teng et al. 2023), and from galaxy to galaxy (Narayanan et al. 2011), caused by the underlying range of the physical conditions.This is particularly important in starbursting galaxies which yield a significantly lower CO-to-H 2 conversion factor than the standard Milky Way value (Renaud et al. 2019b).This adds to uncertainties on the slope and scatter of the KS relation of heterogeneous samples.
Similarly, the choice of the fitting method was also found to have a significant effect on the derived slopes (e.g.Shetty et al. 2013;Kennicutt & De Los Reyes 2021).de los Reyes & Kennicutt (2019) recently revisited the KS relation for nonstarbursting galaxies comparing three widely used fitting techniques, ordinary linear regression, bivariate regression, and hierarchical Bayesian linmix model, finding changes for the inferred slope of up to ∼ 30%.
Nowadays, this variety of results seems to be acknowledged as being mostly due to systematics related to the abovementioned methodological choices.Yet, there is still a poor understanding of the physics behind the intrinsic scatter of the KS relation.Ongoing and future missions are opening new windows on the physics of the earlier Universe, in particular on the star formation activity of galaxies during their first few Gyr thanks to the James Webb Space Telescope.To accompany these efforts, cosmological simulations can provide insights into the behaviors of the current models in these high redshift conditions (e.g.Kravtsov 2003;Feldmann et al. 2012;Semenov et al. 2019).Models and sub-grid prescriptions for star formation and feedback are calibrated using detailed observations in the local Universe, and even mainly from the Solar neighborhood.It is thus important to understand how they behave when applied to different environments.In particular, identifying at which cosmic epoch a given scaling relation emerges is a crucial step in the interpretation of observations at high redshift, and the further confrontation with existing models.
This first paper of a series, intended to complement our understanding of the evolution of the physics of star formation across cosmic time, focuses on the question of when the KS relation emerges, how it evolves, and what physical parameters are primarily driving it.To address these questions, we use the large-scale zoom-in hydrodynamic simulation NewHorizon (Dubois et al. 2021), and perform an analysis of the KS relation on the galactic scale and at different cosmic epochs, from redshift 4 down to 0.25.
The remainder of this paper is structured as follows.Section 2 describes the simulated data set and methods used in the analysis.Section 3 presents the results on the emergence of the KS relation and its dependence on different physical properties of galaxies.These results are discussed in Section 4 and finally, Section 5 concludes.

The NewHorizon simulation
This work makes use of the NewHorizon 1 simulation (Dubois et al. 2021), a large-scale zoom-in simulation of a subvolume extracted from the large-scale cosmological simulation Horizon-AGN (Dubois et al. 2014;Kaviraj et al. 2017).Combining a relatively large volume with a resolution typical of standard zoom-in simulations, NewHorizon captures the structure of the cold interstellar medium of several hundreds of galaxies.This allows us to fully resolve the wider cosmic environment as well as emergently produce a realistic distribution of galaxy properties.Therefore, we are able to perform statistical studies on many galaxy properties, at an unprecedented resolution over such volumes.
NewHorizon reproduces reasonably well many observables (see Dubois et al. 2021), such as the galaxy stellar mass function, the cosmic SFR density, the stellar density, the stellar mass-star formation rate main sequence, galaxy gas fractions, the specific SFR-mass relation, the size-mass relation, the mass-metallicity relation, and the Tully-Fisher relation (see also e.g.Volonteri et al. 2020;Jackson et al. 2021b,a;Martin et al. 2021;Park et al. 2021;Grisdale et al. 2022).The details of the simulation can be found in Dubois et al. (2021), and we describe here only the features of interest for the analysis of the KS relation.
The NewHorizon simulation is run with the adaptive mesh refinement code Ramses (Teyssier 2002), with ΛCDM cosmology compatible with the WMAP-7 data (Komatsu et al. 2011).The mass resolution is 1.2 × 10 6 M ⊙ for the dark matter and 1.3 × 10 4 M ⊙ for the stars.The refinement strategy allows to reach the spatial resolution of up to 34 pc.
NewHorizon includes heating of the gas from a uniform UV background following Haardt & Madau (1996) and models the self-shielding of the ultraviolet background in optically thick regions following Rosdahl & Blaizot (2012).Gas cooling down to ≈ 10 4 K is allowed through collisional ionization, excitation, recombination, Bremsstrahlung, and Compton cooling.Further cooling of metal-enriched gas down to 0.1 K follows tabulated rates from Dalgarno & McCray (1972) and Sutherland & Dopita (1993).
Gas above the density threshold of 10 H cm −3 is converted into stars following the Schmidt relation ρ⋆ = ϵ ⋆ ρ/t ff , where ρ⋆ is the star formation rate density, ρ the gas mass density,  et al. 2017;Trebitsch et al. 2017Trebitsch et al. , 2021)).ϵ ⋆ is a function of the local turbulence Mach number M, and the virial parameter α vir = 2E kin /E grav (E kin and E grav are respectively the turbulent and gravitational energies): where s crit (α vir , M) is the critical logarithmic density contrast of the gas density probability distribution function with variance σ 2 s (M) (see Dubois et al. 2021, for details).The parameter ϕ t is set to the best-fit value between the theory and the numerical experiments (Federrath & Klessen 2012) and ϵ, set to 0.5, mimics proto-stellar feedback effects to regulate the amount of gas eligible to form stars (Matzner & McKee 2000;Alves et al. 2007;André et al. 2010).In short, this prescription favors the rapid formation of stars in dense, gravitationally collapsing medium with compressible turbulence.
A&A proofs: manuscript no.KS_aa NewHorizon includes feedback from type II supernovae (Thornton et al. 1998) following the mechanical supernova feedback scheme of Kimm & Cen (2014, Kimm et al. 2015) to ensure a correct amount of radial momentum transfer.NewHorizon also follows the formation, growth, and dynamics of massive black holes and the associated feedback from active galactic nuclei, following two different modes depending on the Eddington rate (Dubois et al. 2012).At low accretion rates, the massive black hole powers jets releasing mass, momentum, and total energy into the gas (the so-called radio mode feedback, Dubois et al. 2010), while at high rates, it releases only thermal energy (the so-called quasar mode, Teyssier et al. 2011).

Postprocessing and sample selection
Galaxies are identified with the AdaptaHOP halo finder (Aubert et al. 2004) run on the stellar particle distribution (see Dubois et al. 2021 for details).This work employs the 100% purity sample, i.e. halos and embedded galaxies devoid of low-resolution DM particles. 2ollowing the convention adopted in Dubois et al. (2021), we identify the neutral gas component (atomic and molecular), noted H i + H 2 , as denser than 0.1 H cm −3 and colder than 2 × 10 4 K, and the H 2 molecular component denser than 10 H cm −3 and colder than 2 × 10 4 K. Reproducing the ionisation and molecular states of the gas would require a detailed treatment of radiative transfer and molecular chemistry, out of the scope of this paper.The surface densities of neutral (Σ H i+H 2 ) and molecular gas (Σ H 2 ), and of the star formation rate (Σ SFR ) are computed within the (three-dimensional) effective radius R eff of each galaxy, defined as the geometric mean of the half-mass radius of the projected stellar densities along each of the Cartesian axes (see Dubois et al. 2021, for more details).
We do not de-project the galaxies when computing surface densities, and note that this is not expected to have a strong impact on statistical distributions of surface densities (e.g.Appleby et al. 2020).The star formation rate is estimated by considering only the stars younger than 10 Myr, consistent with the timescale probed by the Hα-based SFR indicator.This choice is a compromise, as longer time scales would tend to include the effects of stellar feedback on the properties of the interstellar medium and increase the intrinsic scatter of the Σ gas −Σ SFR relation (e.g.Feldmann et al. 2012;Andersson et al. 2021), in particular at higher redshifts.
The turbulence Mach number M of each gas cell is computed as M = σ g /( √ 3c s ), where σ g and c s are its velocity dispersion sound speed, respectively (see Kraljic et al. 2014, for more details).Then, the Mach number of the galaxy is given by the mass-weighted average of M of every neutral gas (H i + H 2 ) cell 3 within the effective radius R eff .
In this paper, we analyse the population of galaxies at the redshifts 4, 3, 2, 1, and 0.25.We consider only the galaxies with stellar masses M ⋆ above 10 7 M ⊙ , and having hosted star formation in the last 10 Myr.At z = 4, this corresponds to ∼ 90% of the entire sample of galaxies with M ⋆ ≥ 10 7 M ⊙ , while with decreasing redshift this fraction decreases to ∼ 40% at z = 0.25.This is essentially due to the lack of star formation activity dur-Table 1. Number of galaxies and median of their stellar mass (in log M ⊙ ) used in our analysis, at each redshift.Note that our sample is limited to galaxies of stellar mass ≥ 10 7 M ⊙ , with SFR ≥ 10 −3 M ⊙ yr −1 and containing neutral gas.ing the last 10 Myr and is limited to galaxies with M ⋆ ≲ 10 8.5 M ⊙ at z ≥ 2, while below z ∼ 2, more and more massive galaxies are concerned.Only ≲ 2% of galaxies at z = 4 − 1 and ∼ 7% at z = 0.25 do not host any neutral gas within their effective radius and these are limited to low mass range (≲ 10 8 M ⊙ ) at all redshifts.The resulting numbers of galaxies at each redshift are provided in Table 1.Examples of representative galaxies from the various stellar mass bins used in the analysis and at different redshifts are shown in Fig. 1.The stellar mass bins adopted throughout the paper are defined using the quartiles of the mass distribution at each redshift, and thus yield evolving ranges as the overall population grows.

Fitting method
To quantify the correlation between the surface densities of gas and SFR, we fit the distributions with the relation log Σ SFR = a(log Σ gas ) + b, with the best-fit values for the slope a and intercept b.In this paper, we do not attempt to provide a thorough study of the impact of different fitting methods used in the literature on the estimated values for the obtained parameters (we refer the readers to e.g.Hogg et al. 2010, for a discussion on fitting methods used in science).Nevertheless, we compare three different fitting methods: the ordinary least square (OLS) technique, the OLS bisector technique (Isobe et al. 1990), and the Bayesian linear regression.The results of the Bayesian regression are shown throughout the paper, as it provides a more robust treatment of errors and is thus particularly adapted to observational measures.In Appendix B, we report the results of the OLS bisector, together with a more detailed comparison of different methods.In short, all three fitting methods provide a qualitatively similar trend for the slopes and dispersions around the best fit as a function of redshift and stellar mass.We note however that quantitatively, the values of the slope differ: they are systematically higher for the bisector OLS method, and the dispersion around the best fit is also systematically higher.Overall, depending on the population of galaxies and the gas tracer under consideration, the choice of the fitting method can produce changes of 10% to 30% for the slope, in agreement with similar estimates of the impact of fitting algorithms on recent observational data (see de los Reyes & Kennicutt 2019; Kennicutt & De Los Reyes 2021).These differences should be kept in mind when comparing the values reported in the literature.

Distributions of galaxies in the KS plane
We start by investigating the diversity of star-forming galaxies and its evolution with cosmic time, by analysing the distributions of galaxies in the Kennicutt-Schmidt (KS) plane in different stellar mass bins, and as a function of redshift.distribution of galaxies of the NewHorizon simulation in the Σ H 2 -Σ SFR plane, in different stellar mass bins (columns) at different redshifts (rows), colour-coded by their specific star formation rate (sSFR = SFR/M ⋆ ), computed for the entire galaxy with the same timescale as Σ SFR .
Although the distributions vary quite significantly between panels, at a given redshift, there is a substantial overlap in the range of values of both Σ H 2 and Σ SFR within the KS plane.The lowest and highest tails of these distributions at each redshift are typically dominated by galaxies within the lowest and highest stellar mass bins, respectively.Overall, the distributions of galaxies of a given stellar mass quartile move within the KS pa-rameter space towards lower values of Σ H 2 and Σ SFR with decreasing redshift: fewer and fewer galaxies are found far above the canonical KS relation4 (solid line), at all stellar masses.This is accentuated after cosmic noon (z < 2) where only a handful of galaxies reach the sequence of starbursts (dotted line).
As expected, the sSFR of galaxies decreases with decreasing redshift, in particular at z ≤ 2. It also decreases with increasing stellar mass at each redshift.The sSFR of galaxies strongly correlates with Σ SFR at each redshift and in each stellar mass The correlation between sSFR and Σ SFR seen at in different stellar mass bins (Fig/ 2) is still apparent when stacking all galaxies.At all redshifts, the slope (a) and the dispersion (σ) around the best-fit relation (dash-dotted line) are larger for the neutral gas than for the molecular gas.At all redshifts, the correlation is stronger for molecular gas than for neutral gas.
bin, essentially because the sSFR is computed using stars with the same age as Σ SFR (< 10 Myr).This strong correlation vanishes when considering older stars (e.g.< 100 Myr).As a consequence, at fixed Σ H 2 , galaxies with shorter depletion times, have a higher sSFR than those with longer depletion times.We stress that this behavior, although not surprising, is not obvious.Starbursting systems have short depletion times, i.e. the normalisation of the SFR by the gas mass5 , while the sSFR is the SFR normalized by the stellar mass.Galaxies with a given sSFR but different gas fractions could then have significantly different depletion times.As a matter of fact, the systematic qualification of starburst galaxies as outliers above the main sequence of star formation is being questioned by observations (Gómez-Guijarro et al. 2022;Ciesla et al. 2023, see also e.g.Tacconi et al. 2018) and simulations (Renaud et al. 2022).
The trends and correlations from Fig. 2 persist when the neutral gas (Σ H i+H 2 ) is considered instead of the molecular gas alone (Σ H 2 , see Fig . A.1), but with a steepening of the slopes, weakening of correlations, and increased scatter, in agreement with observations, at all redshifts and in all stellar mass bins.
Figure 3 shows the distributions in the Σ H 2 -Σ SFR and Σ H i+H 2 -Σ SFR planes, but now for the entire population of galaxies at different redshifts, by stacking all galaxies from Figs. 2  and A.1, respectively, where 2D histograms are computed by averaging the colour-coded quantity in each bin.The gradients in sSFR seen in individual stellar mass bins are still apparent when stacking all stellar masses.Similarly, the entire galaxy population shows a stronger correlation, smaller dispersion, and shal-lower slope between the SFR density and the molecular gas, than with the neutral gas.
The shallower slope of the correlation with H 2 is consistent with studies of nearby galaxies both on galactic (e.g.Kennicutt 1998;Liu et al. 2015; de los Reyes & Kennicutt 2019; Kennicutt & De Los Reyes 2021) and sub-kpc scales (e.g.Kennicutt et al. 2007;Bigiel et al. 2008;Leroy et al. 2008).Although the value of the slope depends on the employed fitting method and types of galaxies under consideration, it is found to be approximately linear.The relation between SFR and total gas surface densities for a combined sample of normal and starburst galaxies is found to be superlinear with slopes 1.4−1.5 (Kennicutt 1998;Kennicutt & De Los Reyes 2021).A similar slope is found for a sample of nearby normal spiral galaxies (de los Reyes & Kennicutt 2019) and at higher redshifts (z ∼ 1.5; Daddi et al. 2010), while the inclusion of dwarf galaxies tends to produce a shallower slope of ∼ 1.3 (de los Reyes & Kennicutt 2019).When fitted separately, starburst galaxies appear to follow a relation with slope 1−1.2, as recently revealed by Kennicutt & De Los Reyes (2021), which is shallower compared to previous studies finding slopes of 1.3−1.4(e.g.Kennicutt 1989;Daddi et al. 2010), but confirms bimodal (or possibly multimodal) relation for the global star formation (e.g.Daddi et al. 2010;Genzel et al. 2010).
We will now explore when these observed relations emerge.

Emergence of the KS relation
The redshift dependence of the trends highlighted in the previous section suggests that the KS relation evolves with cosmic time.In this section, we explore its emergence and overall evolution.Fig. 4 shows the evolution of the slope a of the relations (Σ SFR ∝ Σ a gas ) fitted with the Bayesian linear regression method,  and A.1, i.e. using the Bayesian fitting method, for the four mass bins considered, from the low mass bin in a light color to the most massive one in dark color.The red points show the slope of the entire galaxy population at a given redshift (i.e.without accounting for their mass, as in Fig. 3).The emergence of the KS relation is shown by the convergence of the slope of the massive galaxies (from the two most massive bins) near the observed relation at z ≈ 2 − 3. Low mass galaxies do not show signs of convergence toward a fixed slope: their KS relation gets continuously shallower after z ≈ 2 − 3. and Fig. 5 displays the dispersion of the data around these fits (see Appendix B.3 for the equivalent plots using another fitting method).
The least massive galaxies (the lower half of the stellar mass distribution) clearly differ from the most massive cases (the highest stellar mass bin): at all redshifts, their KS relations are shallower and more dispersed.An examination of the distributions (Fig. 2) reveals that this originates from the presence of galaxies with short depletion times at low Σ H 2 .Such cases of rapid star formation in (relatively) diffuse gas possibly due to environmental triggers like mergers, or strong (or fast) gas outflows, are found at low mass at all redshifts, but only at high redshift for the massive quartiles.At cosmic noon (z ∼ 2), this regime disappears at all masses but reappears at z ≲ 2 at low masses.This explains the overall "bell" shape at low mass in Fig. 4.This effect is significantly more pronounced in the neutral gas (H i + H 2 ), which indicates that the fraction of molecular over neutral gas plays a role in star formation in diffuse gas.It is likely that at the lowest masses, the galaxies comprise only one active starforming region at a given instant, i.e. a molecular cloud with an extended atomic envelope (recall Fig. 1), which would favour the star formation regime noted here.
The slope of the KS relation stabilizes below z ≈ 2 for the overall population (red thick line on Fig. 4) and the most massive galaxies (M ⋆ ≳ 10 8 M ⊙ at this redshift).This is also the epoch when the dispersion around the relation reaches its final, minimum plateau (Fig. 5).Therefore, the present-day KS relation emerges at cosmic noon (z ≈ 2) in the most massive galaxies.Our results predict that populating the KS plane with observational data from the top 50% most massive galaxies at redshifts ≳ 3 would result in a different and significantly more dispersed relation than the one currently established at low and intermediate redshifts (when considering local spirals, z = 1 − 2 discs, and BzK galaxies, Daddi et al. 2010;Tacconi et al. 2010;Salmi et al. 2012).
However, this is not the case for the least massive galaxies, for which no stabilization of the slope is seen, neither in molecular nor neutral gas.Interestingly, the dispersion around the best fit of these galaxies still yields a behaviour very similar, qualitatively and quantitatively, to that of the most massive ones, i.e. a decrease until z ≈ 2 − 3 followed by a relatively flat plateau.Hence, the relation for the low mass galaxies becomes simultaneously tighter and shallower below z ≈ 2−3.This could indicate that the extreme cases at high redshift either evolve to a more massive quartile via rapid growth or conversely get quenched and disappear from the star-forming sample.The remaining lowmass objects would then display a more homogeneous behavior.
These shallower KS relations of low-mass galaxies are in qualitative agreement with the observations of local dwarf galaxies, which report slopes around unity (e.g.Filho et al. 2016;Roychowdhury et al. 2017) 6 .The underlying reason is still debated and probably consists of an interplay between galaxy interactions and the low-metallicity contents of these dwarf galaxies which caps their efficiency at forming molecular gas (Cormier et al. 2014).Such hypotheses are in line with our measurements of star formation in diffuse gas that we interpret as star-forming regions with extended atomic envelopes.Confirming these ideas requires a resolved analysis of these galaxies, instead of the statistical approach we follow here.Thus, we will explore these hypotheses in a forthcoming paper.
When considering the relation between Σ SFR and the neutral gas surface density Σ H i+H 2 , we retrieve qualitatively the same evolution of the slope and the dispersion with redshift and stellar mass (for all mass quartiles), but with larger slopes.The reason for this steepening of the relations is the presence of sub-efficient star-forming regions in galaxies at low Σ H i+H 2 , often referred to as the "break" of the KS relation (see an illustration in Bigiel et al. 2008).The physical origin of the break has been shown analytically (Renaud et al. 2012) and numerically (Kraljic et al. 2014) to be caused by low levels of turbulence which do not efficiently promote the formation of dense gas, or in other words, by a low filling factor of star-forming gas in the volumes considered.In turn, the break becomes more apparent when including the atomic component in our analysis, by increasing the gas surface density without altering Σ SFR , which bends the distribution of galaxies in the KS plane below the canonical KS relation.
At low redshifts, the low mass galaxies of our sample correspond to dwarfs of which the low-metallicities (Dubois et al. 2021) could explain the inefficient formation of molecular gas (at small scales), and in turn slow down star formation, even at high Σ H 2 at galactic scales (but see the discussion of Roychowdhury et al. 2017 on the relatively small effect of the metallicity on the KS relation).Interestingly, Figure 13 of Dubois et al. (2021) indicates that the relation between the stellar mass and the metallicity in NewHorizon varies only very weakly with the redshift.We have checked that this remains true when selecting the star forming galaxies only.We confirm that the gas and stellar metallicities7 for the selection of galaxies within the KS plane increase with stellar mass and with decreasing redshift, as expected, but the redshift evolution of the mass-metallicity relation is only weak.Moreover, at all redshifts, at a given stellar mass, the metallicity does not show any gradient within the KS plane.This implies that low redshift dwarfs have a similar metallicity as the galaxies with the same stellar mass at z ≳ 2, but which are then in our upper mass bin, and already follow distributions close to the canonical KS relation.This demonstrates that the stellar mass and the metallicity are not key parameters in driving the emergence of the KS relation.The role of other physical quantities is explored in the next section.

Molecular and total gas content
We now investigate whether the relation between Σ H 2 and Σ SFR is driven by the gas content of galaxies.
Figure 6 (top panel) shows the molecular gas fraction, defined as the fraction of H 2 mass over the neutral gas, i.e.M H 2 /(M H i + M H 2 ), and its evolution within the KS plane with the cosmic time for the entire galaxy population.The fraction of molecular gas decreases with decreasing redshift.Moreover, it correlates with Σ H 2 resulting in vertical contours (within the KS plane) at all redshifts.The same trends are seen at a given stellar mass, although the fraction of molecular gas increases with stellar mass (Fig. A.2, see also Dubois et al. 2021, their figure 19).
Molecular gas content of galaxies may also be defined in terms of baryonic fraction, i.e.M H 2 /(M H i + M H 2 + M ⋆ ).The correlation between the baryonic molecular gas fraction and Σ H 2 is maintained at all redshifts (Fig. 6, middle panel), although it weakens at z ≲ 2 when it is only carried by the low-mass galaxies -the massive galaxies having very low baryonic molecular gas fraction, independently of Σ H 2 (see Fig. A.3).This lack of correlation at low redshift results from the baryonic fraction of molecular gas vs stellar mass relation getting shallower at these late times (see top right panel of figure 19 of Dubois et al. 2021).
We finally consider the neutral baryonic gas fraction, i.e.
As already reported by Dubois et al. (2021) for the NewHorizon galaxies, this fraction strongly anticorrelates with the stellar mass, but only mildly depends on the redshift at a given mass.This is confirmed in the KS plane for galaxy stacks (Fig. 6, bottom panel) and individual stellar mass bins (Fig. A.4).At high redshift (z = 4), the neutral gas fraction increases with Σ H 2 for all stellar masses, but this correlation disappears at later epochs, where this fraction varies weakly with Σ H 2 at fixed stellar mass but varies strongly with stellar mass.The combination of these relations between the neutral gas fraction and Σ H 2 at high z, and M ⋆ at all z, translates into a non-trivial evolution of the distributions of the neutral gas fraction in the KS plane when the stellar mass is marginalized out.The reversal of the trend between the neutral gas fraction and Σ H 2 between high and low redshift is thus a direct consequence of the dependencies highlighted above, and of how early the galaxies build up their stellar masses.
In conclusion, both the molecular and neutral gas fractions vary with one of the parameters of the KS plane (Σ H 2 ), but have close to no influence on the other (Σ SFR ), except in the very diffuse gas of the low mass galaxies, as noted above.As such, at galactic scales, the KS relation does not originate from the gas fractions of the galaxies, at any redshift, nor at any stellar mass.The quantity that correlates better within the KS plane is the fraction of cold gas that is in the dense phase, however, it does not fully capture the variation with both the surface density of gas and the star formation rate of galaxies.

Turbulence
Previous works (e.g.Vazquez-Semadeni 1994;Renaud et al. 2014) have pointed out the role of turbulence in setting the density distributions of gas, and thus the amount of star-forming gas in galaxies.Here, we test whether turbulence explains the emergence of the KS relation over cosmic time.
Figure 7 shows the distribution of galaxies in the KS plane, now color-coded with their turbulence Mach number M. At all

])
Fig. 6.Same as Fig. 3, but color-coded by the molecular over cold gas fraction (top), and baryonic molecular (middle) and neutral (bottom) gas fractions.Molecular gas fraction (top panel) correlates with Σ H 2 at all redshifts.Baryonic molecular fraction (middle panel) correlates with Σ H 2 at high redshifts, below z∼2 this correlation weakens and is essentially carried by low-mass galaxies.The correlation between the neutral gas fraction and Σ H 2 is apparent only at z = 4.At z ≲2 the trend reverses such that this fraction decreases with increasing Σ H 2 .
redshifts, M evolves monotonically along the best-fit relation in the KS parameter space, by increasing with Σ H 2 and Σ SFR .Furthermore, at fixed Σ H 2 , Σ SFR is positively correlated with M.
As shown in Fig. A.5, these trends are independent of the galaxy's stellar mass.At fixed stellar mass, galaxies at high redshifts are more turbulent than their low redshift counterparts, and at fixed redshift, more massive galaxies tend to be more turbulent than their lower mass counterparts, a direct consequence of increasing gas richness in galaxies with redshift (e.g.Bournaud et al. 2010;Renaud et al. 2012).
To investigate the physical origin of these trends, Fig. 8 shows the distributions of the Mach number M ∝ σ vel / √ T , and the underlying quantities which are the gas velocity dispersion (σ vel ) and the temperature (T )8 .
While the velocity dispersion decreases with redshift for all mass bins, only massive galaxies maintain high values at low z (∼ 10 km s −1 ), leading to significantly higher median values.This is confirmed by observations at low redshifts of lower dispersion in dwarf galaxies (∼ 1 km s −1 ) than in massive galaxies (∼ 10 km s −1 , e.g.Hunter et al. 2021).This general trend likely originates in parts from the overall lowering of the star formation activity with cosmic time, and possibly the less efficient coupling of feedback with the local interstellar medium (ISM), as opposed to intergalactic medium due to low escape velocity, in low-mass galaxies (i.e. with shallow potential wells).
The locations of the galaxies with high σ vel in the KS plane (Fig. A.6) reveal complex correlations with the indicators of star formation: while galaxies with short depletion times tend to have high σ vel , high-velocity dispersion are found across the entire KS plane.This indicates that stellar feedback is not the only factor in setting the velocity dispersion, and therefore the KS relation (see also Agertz et al. 2011;Agertz & Kravtsov 2015), at the galactic scale.Galactic dynamics and interaction-triggered stirring are likely important drivers of the velocity dispersion (see Renaud et al. 2014 for an illustration that increased velocity dispersion is a cause and not a consequence of starbursts in mergers).
The temperature of all mass quartiles decreases with decreasing redshift, but Fig. 8 (middle-column) reveals that the stellar mass only discriminates the distributions of T at late times (z ≲ 1).Contrary to the velocity dispersion, there is no relation between the star formation indicators and the temperature in the KS plane (Fig. A.7), which is consistent with the interpretation of the limited impact of feedback, even though the details, in particular on small scales, might be more complicated.
In terms of Mach number, the trends noted from the two underlying quantities (σ vel and T ) naturally combine to lead to a shift of the distributions of M towards low values with decreasing redshift, an effect which is significantly more pronounced for low-mass galaxies (Fig. 8).
A more detailed examination (Fig. A.5) reveals that the trends found in velocity dispersion and in temperature conspire to give rise to a clear evolution of M along the KS relation, with 'tighter correlation' with Σ gas and Σ SFR than the individual σ vel and T .This further demonstrates the paramount role of turbulence in the star formation activity, in particular in the KS plane.
Higher Mach numbers favor higher density contrasts in the ISM (i.e. a wider gas density PDF, see Federrath et al. 2008), and thus the formation of a larger faction of dense molecular gas.This explains that dwarf galaxies at low redshift, with a low Mach number, tend to have lower molecular gas fractions (Fig. A.2) than their massive counterparts, and thus appear below the canonical KS relation, in the so-called "break" (see Kraljic et al. 2014), when considering the total neutral gas (Fig. A.1), but are shifted toward the low gas densities when considering the molecular phase only (Fig. 2).Finally, the high Σ SFR of these galaxies implies that this shift toward low Σ H 2 places them above the canonical KS relation, which drives the flattening of the relation of these sub-populations (recall Fig. 4).
Furthermore, the histograms of Fig. 8 show that the distributions of velocity dispersion in low-mass galaxies become peaked toward the low-value end at low redshift, while the massive galaxies only exhibit a tail with only a few cases at such lowvelocity dispersion, and the bulk of their distribution remains centered around higher values (∼ 10 km s −1 ) with little evolution after z ≲ 2. In other words, the lower end of the distributions in σ vel gets more and more populated with low-mass galaxies with decreasing redshift, while the distribution of velocities dispersion of massive galaxies ceases to evolve (statistically).For the reasons discussed before, this transpires in the histograms of Mach number, and finally in the distribution of galaxies in the KS plane.Therefore, the mass-dependent evolution of the velocity dispersion explains the convergence of the slope of the KS relation at high mass after z ≈ 2, and the absence of the convergence in low-mass galaxies, noted in Fig. 4.

Scale and projection effects
So far, we have conducted our analysis using the observables Σ H 2 , or Σ H i+H 2 , and Σ SFR , and identified relations in the KS plane.However, by doing so, we effectively introduce an arbitrary choice for the spatial scale used in the measurement of both quantities, which necessarily impacts the values of the surface densities and possibly artificially distorts the distributions of galaxies in the KS plane.To establish whether our conclusions depend on our choice of examining projected quantities, and over the scale of the effective radius, we plot in Fig. 9 the distributions of galaxies of our sample in the plane of molecular gas mass vs. SFR, i.e. a deprojected version of the KS plot, colourcoded by mass (top) and size R eff (bottom).At all redshifts, more massive galaxies tend to have higher SFR and molecular gas mass, however, there is no obvious correlation between SFR, M H 2 , and effective radius of galaxies.Furthermore, all the trends (or the absence thereof) with molecular and total gas fractions, and M seen in the KS parameter space are retrieved with deprojected quantities, as shown in Fig. C.1.Therefore, the trends seen in the KS plane are not primarily driven by measuring the physical quantities in projection rather than in 3D, and the scale adopted does not introduce biases in the distributions.
The diversity of star formation activities seen in the wide distribution of Σ H 2 , Σ H i+H 2 , and Σ SFR , but also in the slopes and the scatters of the KS relation, results from the convolution of two other effects: (i) the diversity of galaxies in the sample, illustrated by the variations of the KS relations with stellar mass and redshift (Fig. 2), and (ii) the integration of the local, smallscale star formation law over entire galactic scales where not all the ISM is star-forming.Indeed, the shape (slope, offset, break, scatter) of the distributions of galaxies in the KS plane is driven by the underlying distribution of the physical properties of the star-forming regions within each galaxy, and it is the evolution of these distributions as functions of redshift, galactic mass, and other factors like the environment, that sets the evolution of the KS relations.

Sub-grid models
Our work points out the key role of turbulence (Mach number) in driving the KS relation at the galactic scale, already at high redshift.This confirms the analytical results of Renaud et al. (2012), and the numerical work of Kraljic et al. (2014) which conducted a similar study without cosmological context, and for galaxies in the nearby Universe only (i.e. at low gas fractions).We note that the star formation model in Kraljic et al. (2014) differs from those in NewHorizon, as it uses a fixed star formation efficiency per free-fall time and is applied at a higher resolution (∼ 1 pc).In other words, the right-hand side of Equation 1 reduces to a constant value in the former.The paramount role of turbulence in setting the KS relation found with both models (independently of an explicit inclusion of turbulence in the star formation model in the latter), further strengthens our conclusion.
Yet, it is important to keep in mind that the differences in the sub-grid models used do not necessarily correspond to different physics.Schemes that do not capture the formation of star-forming clouds (i.e. at resolutions ≳ 20 pc) ought to incorporate this aspect in their star formation prescription.This can be achieved by imposing a criterion on the instability of the gas, through e.g.converging flows and/or the virial parameter, as is done in NewHorizon.However, at cloud scale (≈ 1 − 10 pc), the fragmentation of the ISM into individual clouds is already captured by the simulations, such that an instability criterion is redundant: at high resolution, the dense gas has necessarily gone through the instability phase.As for any physical mechanism, it is crucial to identify which processes are not captured explicitly by the simulation, in order to construct and parametrize the "sub" aspect of the sub-grid models.
By conducting a spatially resolved study of the KS relation within the Fire (Hopkins et al. 2014) framework, Orr et al. (2018) pinpointed instead the central role of stellar feedback in regulating star formation at small scales (see also Dekel et al. 2019, for similar conclusion for the global KS relation).We note however that their sub-grid treatment de facto implies an important role for feedback, as it is required to regulate the star formation process set with an efficiency of 100 percent (as opposed to ∼ 0.1 − 10% in most of the literature from observations and simulations).It is, therefore, possible that a different conclusion, perhaps closer to ours, would be reached by adopting a lower star formation efficiency, and thus automatically decreasing the impact of feedback (see e.g.Brucy et al. 2020Brucy et al. , 2023;;Hu et al. 2022).For instance, in the analytical model of Renaud et al. (2012), stellar feedback can be introduced as a cap on the effective star formation efficiency.While this results in lowering the slope of the KS relation, a requirement to match observations, a KS-like power-law does exist without feedback, and originates solely from the log-normal shape of the distribution of gas density, itself known to be set by turbulence (e.g.Vazquez-Semadeni 1994).The fact that the KS relation can be modeled from such different underlying physics suggests that divergences should be sought in more fundamental quantities or behaviour, and possi- .Dash-dotted lines are fits in the logarithmic space at each redshift, with the slope a and the standard deviation of residuals σ.The coefficient R is the Pearson correlation coefficient.At all redshifts, more massive galaxies tend to have higher SFR and M H 2 compared to their lower mass counterparts.As galaxies grow in mass, they grow in size, however, no obvious correlation is seen between R eff , SFR and molecular gas mass.bly at smaller scales, i.e. before the differences between models get blended in integrated, projected, and galactic-scale averaged measurements.Considering other, more extreme environments where the relative contributions of the mechanisms involved vary could certainly provide interesting insights.

A variety of possible underlying physical mechanisms
Taking advantage of the broad diversity of galaxies in NewHorizon, we show here that the gas fraction does not strongly influence the star formation relation, at least when integrated over entire galaxies, and that the mild trends found between the gas fraction and its surface density (be it neutral or molecular) actually get reversed at z ≈ 2. This change entails an evolution of the morphology and size of the star-forming volume, likely connected to a more concentrated activity.Underlying physical reasons could be the intrinsic evolution of discs towards fueling more and more gas to the nuclei, and/or environmental effects due to gravitational torques exerted on disc material by more and more massive companion galaxies.Some of these aspects have been explored in the special case of the Milky Way-like galaxy, using the Vintergatan simulation (Agertz et al. 2021;Renaud et al. 2021b).These works highlight the necessity for the galactic disk to be in place (Park et al. 2021) for the galaxy to strongly react via large-scale wakes to interactions through a starburst activity (Segovia Otero et al. 2022).The redshift-dependence of tidal compression, both in terms of intensity and mass involved also appears as crucial in the cosmic evolution of starbursts (Renaud et al. 2022) as it controls energy input at the cascade injection scale.Exploring these points further and over an entire population of galaxies, like that in NewHorizon, requires a dedicated analysis of the diversity of individual evolution that builds the population statistics shown here, which we leave for a forthcoming paper (Kraljic et al. in preparation).
Galaxies with the highest global turbulence level are not only those which host the densest gas and form the most stars (as shown by their locations at high surface densities of gas and of SFR).They are also the galaxies with the shortest depletion time.This is particularly visible at high redshift (z ≳ 2) in Fig. 7, where the most turbulent galaxies lie around the starburst rela-tion from Daddi et al. (2010), about 1 dex above the canonical KS relation.This is in line with the conclusions of Renaud et al. (2014Renaud et al. ( , 2019aRenaud et al. ( , 2021a) ) which reported that the increase of the level of compressive turbulence in mergers can lead to a starburst activity.In this context, it is crucial to differentiate the production of many stars (high SFR) from the fast production of stars (short depletion time).While the two aspects are independent, our results show that high levels of turbulence make some high redshift galaxies reach both a high Σ SFR , and a short depletion time.This then changes at lower redshift, when the most turbulent galaxies of our sample do not necessarily yield the shortest depletion times.Such an evolution could be connected with the rarefaction of the major mergers at late epochs, and thus the statistical lowering of the tidal and turbulent compression (Renaud et al. 2022), and likely has implications on the evolution of the normalization of the main sequence of galaxy formation (e.g.Tacconi et al. 2018).

Conclusion
We have investigated the distributions of galaxies from the cosmological hydrodynamic simulation NewHorizon in the KS parameter space, as a function of their stellar mass and redshift, the emergence of the star formation scaling-laws at the galactic scale, and its physical drivers.Our main results are: -Both the stellar mass and redshift influence the overall location of the galaxy population in the KS plane.-A power-law relation of the form Σ SFR ∝ Σ a gas with a slope a ≈ 1.4 emerges at z ≈ 2 − 3 for the most massive half of the galaxy population (M ⋆ ≳ 10 8 M ⊙ at these redshifts) in agreement with observations up to these redshifts.However, the slope of the relation varies at earlier epochs, with an increased scatter.This indicates that the KS relation might not provide a robust calibration for star formation in galaxies at very high redshift.For the least massive galaxies, there is no sign of the convergence of the slope of their distribution in the KS plane, as it continues to get shallower at the last epochs.The slopes are systematically higher when considering the total neutral gas as opposed to the molecular gas.At all stellar masses, the dispersion around the best-fit relation decreases with the decreasing redshift.
-The gas fraction (neutral or molecular) does not correlate with the star formation activity as traced by Σ SFR and therefore does not play a primary role in establishing the KS relation.Similarly, neither the velocity dispersion of the gas nor its temperature alone can fully explain the star formation activity of galaxies as captured by the KS relation, pointing towards a limited impact of feedback.
-Conversely, the level of turbulence of the interstellar medium, as quantified by the Mach number, is found to drive the relation between gas and SFR densities at all redshifts, independently of stellar mass.More specifically, it is the ability of a galaxy to reach a supersonically turbulent regime that matters, with the Mach number (M > 1) being the driver of the KS relation independent of stellar mass.At high redshift, for a given gas density, the most turbulent galaxies yield short depletion times, characteristic of starburst galaxies.Their frequency decreases at low redshift.
The evolution reported here and a number of previous works on the star formation activity at galactic scales point toward an important role of inflow, interactions, mergers, and the proximity of the disk to marginal stability in driving the star formation relations and their scatters.The latter could act as a confounding factor for efficient turbulent cascade and star formation, explaining the emergence of tighter KS scaling relations, when secular dissipative processes take over.Exploring these aspects requires tracking individual galaxies along their merger histories, and seeking changes in the properties of the star-forming material during the starburst phases, both at cloud and galactic scales.We will cover these topics in the forthcoming papers of this series.

Fig. 1 .
Fig. 1.Projection of gas density, within the 10R eff thick slices, of representative galaxies at different redshifts (rows) and stellar masses (columns).Dashed circles show the effective radius of the stellar component (see Section 2.2 for the definition), and the white horizontal bars indicate a 1 kpc scale.

Fig. 2 .
Fig. 2. Distribution of galaxies in the KS plane at different redshifts (rows), and in four equally-populated stellar mass quartiles at each redshift (columns).The mass range of each quartile is shown in square brackets (in log M ⊙ ).The colours indicate the specific star formation rate of the galaxies, measured over a time scale of 10 Myr.Dash-dotted lines are fits at each redshift and mass bin, with the slope a and σ, the standard deviation of residuals of the best-fit relation, shown in the lower right corners.The coefficient R shown on the bottom right of each panel is the Pearson correlation coefficient.The solid black and dotted black lines show the sequence of discs and starbursts, respectively, from Daddi et al. (2010), for reference.Note, however, that the fitting method differs from the one adopted in this work.Regardless of stellar mass, the distributions of galaxies move within the KS plane towards lower values of Σ H 2 and Σ SFR with decreasing redshift.At each redshift and in each stellar mass bin, the sSFR of galaxies strongly correlates with Σ SFR .A version of this figure using the neutral gas is available in Fig. A.1.

Fig. 3 .
Fig.3.Same as Figs.2 and A.1, but without binning the stellar masses, and considering the molecular gas only (top), and the neutral gas (bottom).The correlation between sSFR and Σ SFR seen at in different stellar mass bins (Fig/2) is still apparent when stacking all galaxies.At all redshifts, the slope (a) and the dispersion (σ) around the best-fit relation (dash-dotted line) are larger for the neutral gas than for the molecular gas.At all redshifts, the correlation is stronger for molecular gas than for neutral gas.

Fig. 4 .
Fig. 4. Evolution of the slope of the fits of the KS relation from Figs. 2and A.1, i.e. using the Bayesian fitting method, for the four mass bins considered, from the low mass bin in a light color to the most massive one in dark color.The red points show the slope of the entire galaxy population at a given redshift (i.e.without accounting for their mass, as in Fig.3).The emergence of the KS relation is shown by the convergence of the slope of the massive galaxies (from the two most massive bins) near the observed relation at z ≈ 2 − 3. Low mass galaxies do not show signs of convergence toward a fixed slope: their KS relation gets continuously shallower after z ≈ 2 − 3.

Fig. 5 .
Fig. 5. Same as Fig. 4, but showing the dispersion around the best fit.Only the vertical dispersion, i.e. in log(Σ SFR ) is considered here.

Fig. 7 .PDFFig. 8 .
Fig. 7. Same as Fig. 3, but colour-coded by the Mach number M. At all redshifts, M increases monotonically with increasing both Σ H 2 and Σ SFR .At fixed Σ H 2 , M is correlated with Σ SFR .

Fig. 9 .
Fig.9.Distribution of galaxies in the M H 2 -SFR plane, i.e. deprojected version of the KS plane, at redshifts 4, 3, 2, 1 and 0.25, from left to right, respectively, as a function of galaxy mass (top) and effective radius (bottom).Dash-dotted lines are fits in the logarithmic space at each redshift, with the slope a and the standard deviation of residuals σ.The coefficient R is the Pearson correlation coefficient.At all redshifts, more massive galaxies tend to have higher SFR and M H 2 compared to their lower mass counterparts.As galaxies grow in mass, they grow in size, however, no obvious correlation is seen between R eff , SFR and molecular gas mass.