11institutetext: Observatoire Astronomique de Strasbourg, Université de
Strasbourg, CNRS, UMR 7550, F-67000 Strasbourg, France
11email<EMAIL_ADDRESS>22institutetext: Lund Observatory,
Division of Astrophysics, Department of Physics, Lund University, Box 43,
SE-221 00 Lund, Sweden 33institutetext: University of Strasbourg Institute for
Advanced Study, 5 allée du Général Rouvillois, F-67083 Strasbourg, France
44institutetext: Institut d’Astrophysique de Paris, CNRS and Sorbonne
Université, UMR 7095, 98 bis Boulevard Arago, F-75014 Paris, France
55institutetext: IPhT, DRF-INP, UMR 3680, CEA, L’Orme des Merisiers, Bât 774,
91191 Gif-sur-Yvette, France 66institutetext: Korea Institute for Advanced
Study, 85 Hoegi-ro, Dongdaemun-gu, Seoul 02455, Republic of Korea
77institutetext: Department of Astrophysics, American Museum of Natural
History, 200 Central Park West, New York, NY 10024, USA 88institutetext:
Department of Physics, University of Oxford, Keble Road, Oxford OX1 3RH,
United Kingdom 99institutetext: Centre for Astrophysics Research, Department
of Physics, Astronomy and Mathematics, University of Hertfordshire, Hatfield
AL10 9AB, UK 1010institutetext: Department of Astronomy and Yonsei University
Observatory, Yonsei University, Seoul 03722, Republic of Korea
1111institutetext: Korea Astronomy and Space Science Institute, 776
Daedeokdae-ro, Yuseong-gu, Daejeon 34055, Republic of Korea 1212institutetext:
Steward Observatory, University of Arizona, 933 N. Cherry Ave, Tucson, AZ
85719, USA 1313institutetext: Observatoire de la Côte d’Azur, CNRS,
Laboratoire Lagrange, Bd de l’Observatoire, CS 34229, F-06304 Nice Cedex 4,
France

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

Katarina Kraljic 11 0000-0001-6180-0245 Florent Renaud 2233
0000-0001-5073-2267 Yohan Dubois 44 0000-0003-0225-6387 Christophe Pichon
445566 Oscar Agertz 22 0000-0002-4287-1088 Eric Andersson 77
0000-0003-3479-4606 Julien Devriendt 88 Jonathan Freundlich 11
0000-0002-5245-7796 Sugata Kaviraj 99 0000-0002-5601-575X Taysun Kimm 1010
Garreth Martin 11111212 Sébastien Peirani 1313 Álvaro Segovia Otero 22
0000-0002-0415-3077 Marta Volonteri 44 Sukyoung K. Yi 1010 0000-0002-4556-2619

(Received July 1, 2023; accepted )

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 $\mbox{$\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.

###### Key Words.:

galaxies: evolution – galaxies: star formation – galaxies: ISM – turbulence –
methods: numerical

##  1 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 $\mbox{$\Sigma_{\rm
SFR}$}\propto\Sigma_{\rm gas}^{a}$ with an index $a\approx 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), sub-millimeter 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 $\sim 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 $\sim 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 ($\approx 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 starbursts 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., 2008, 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 galaxy-averaged (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-H2 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-H2 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 non-starbursting 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 $\sim$
30%.

Nowadays, this variety of results seems to be acknowledged as being mostly due
to systematics related to the above-mentioned 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.

##  2 Methods

Figure 1: Projection of gas density, within the 10$R_{\rm{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.

###  2.1 The NewHorizon simulation

This work makes use of the NewHorizon 111http://new.horizon-simulation.org
simulation (Dubois et al., 2021), a large-scale zoom-in simulation of a sub-
volume 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 $\Lambda$CDM cosmology compatible with the WMAP-7 data
(Komatsu et al., 2011). The mass resolution is $1.2\times 10^{6}\,\rm
M_{\odot}$ for the dark matter and $1.3\times 10^{4}\,\rm M_{\odot}$ 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 $\approx 10^{4}\,\rm K$ is allowed through collisional
ionization, excitation, recombination, Bremsstrahlung, and Compton cooling.
Further cooling of metal-enriched gas down to $0.1\,\rm K$ follows tabulated
rates from Dalgarno & McCray (1972) and Sutherland & Dopita (1993).

Gas above the density threshold of 10 $\rm{H}\,\rm{cm}^{-3}$ is converted into
stars following the Schmidt relation
$\dot{\rho}_{\star}=\epsilon_{\star}{\rho/t_{\rm ff}}$, where
$\dot{\rho}_{\star}$ is the star formation rate density, $\rho$ the gas mass
density, $t_{\rm ff}=\sqrt{3\pi/(32G\rho)}$ the local free-fall time of the
gas, $G$ the gravitational constant, and $\epsilon_{\star}$ is a varying star
formation efficiency (see Krumholz & McKee, 2005; Padoan & Nordlund, 2011;
Hennebelle & Chabrier, 2011; Federrath & Klessen, 2012; Kimm et al., 2017;
Trebitsch et al., 2017, 2021). $\epsilon_{\star}$ is a function of the local
turbulence Mach number ${\cal M}$, and the virial parameter $\alpha_{\rm
vir}=2E_{\rm kin}/E_{\rm grav}$ ($E_{\rm kin}$ and $E_{\rm grav}$ are
respectively the turbulent and gravitational energies):

$\epsilon_{\star}=\epsilon_{\star}({\cal M},\alpha_{\rm
vir})=\frac{\epsilon}{2\phi_{\rm
t}}\exp\left(\frac{3}{8}\sigma_{s}^{2}\right)\left[1+\mathrm{erf}\left(\frac{\sigma_{s}^{2}-s_{\rm
crit}}{\sqrt{2\sigma_{s}^{2}}}\right)\right],$ (1)

where $s_{\rm crit}(\alpha_{\rm vir},\mathcal{M})$ is the critical logarithmic
density contrast of the gas density probability distribution function with
variance $\sigma_{\rm s}^{2}(\mathcal{M})$ (see Dubois et al., 2021, for
details). The parameter $\phi_{\rm t}$ is set to the best-fit value between
the theory and the numerical experiments (Federrath & Klessen, 2012) and
$\epsilon$, 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.

NewHorizon includes feedback from type II supernovae (Thornton et al., 1998)
following the mechanical supernova feedback scheme of Kimm & Cen (2014, ) 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).

###  2.2 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.222Given that NewHorizon is a zoom
simulation embedded in a larger cosmological volume filled with lower DM
resolution particles, some halos of the zoom regions can be polluted with low-
resolution DM particles.

Following the convention adopted in Dubois et al. (2021), we identify the
neutral gas component (atomic and molecular), noted H i \+ H2, as denser than
0.1 $\rm{H}\,\rm{cm}^{-3}$ and colder than $2\leavevmode\nobreak\
\times\leavevmode\nobreak\ 10^{4}\,\rm\leavevmode\nobreak\ K$, and the H2
molecular component denser than 10 $\rm{H}\,\rm{cm}^{-3}$ and colder than
$2\times 10^{4}\,\rm 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
($\Sigma_{\rm\textsc{H\,i}+\rm{H}_{2}}$) and molecular gas ($\Sigma_{\rm
H_{2}}$), and of the star formation rate ($\Sigma_{\rm SFR}$) are computed
within the (three-dimensional) effective radius $R_{\rm 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$\alpha$-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 $\mbox{$\Sigma_{\rm gas}$}-\mbox{$\Sigma_{\rm SFR}$}$
relation (e.g. Feldmann et al., 2012; Andersson et al., 2021), in particular
at higher redshifts.

The turbulence Mach number $\mathcal{M}$ of each gas cell is computed as
$\mbox{$\mathcal{M}$}=\sigma_{\rm g}/(\sqrt{3}c_{\rm s})$, where $\sigma_{\rm
g}$ and $c_{\rm 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 $\mathcal{M}$ of every neutral gas (H
i $+$ H2) cell333However, considering the molecular gas alone yields
qualitatively similar results (not shown). within the effective radius $R_{\rm
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_{\star}$
above $10^{7}\,\mbox{$\rm{M}_{\odot}$}$, and having hosted star formation in
the last 10 Myr. At $z=4$, this corresponds to $\sim$ 90% of the entire sample
of galaxies with $M_{\star}\geq 10^{7}\,\mbox{$\rm{M}_{\odot}$}$, while with
decreasing redshift this fraction decreases to $\sim$ 40% at $z=0.25$. This is
essentially due to the lack of star formation activity during the last 10 Myr
and is limited to galaxies with $M_{\star}\lesssim
10^{8.5}\,\mbox{$\rm{M}_{\odot}$}$ at $z\geq 2$, while below $z\sim 2$, more
and more massive galaxies are concerned. Only $\lesssim 2\%$ of galaxies at
$z=4-1$ and $\sim 7\%$ at $z=0.25$ do not host any neutral gas within their
effective radius and these are limited to low mass range ($\lesssim 10^{8}$
$\rm{M}_{\odot}$) 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.

Table 1: Number of galaxies and median of their stellar mass (in log $\rm{M}_{\odot}$) used in our analysis, at each redshift. Note that our sample is limited to galaxies of stellar mass $\geq 10^{7}$ $\rm{M}_{\odot}$, with SFR $\geq 10^{-3}\,\mbox{$\rm{M}_{\odot}$}\,\rm yr^{-1}$ and containing neutral gas. $z$ | Number of galaxies | Median $M_{\star}$
---|---|---
4 | 535 | 7.55
3 | 558 | 7.74
2 | 582 | 7.99
1 | 303 | 8.63
0.25 | 153 | 9.07

###  2.3 Fitting method

To quantify the correlation between the surface densities of gas and SFR, we
fit the distributions with the relation $\log\Sigma_{\rm
SFR}\leavevmode\nobreak\ =\leavevmode\nobreak\ a(\log\Sigma_{\rm 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.

##  3 Results

###  3.1 Distributions of galaxies in the KS plane

Figure 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
$\rm{M}_{\odot}$). 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 $\sigma$, 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 $\Sigma_{\rm H_{2}}$
and $\Sigma_{\rm SFR}$ with decreasing redshift. At each redshift and in each
stellar mass bin, the sSFR of galaxies strongly correlates with $\Sigma_{\rm
SFR}$. A version of this figure using the neutral gas is available in Fig. 10.

Figure 3: Same as Figs. 2 and 10, but without binning the stellar masses, and
considering the molecular gas only (top), and the neutral gas (bottom). The
correlation between sSFR and $\Sigma_{\rm 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 ($\sigma$) 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.

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. Figure 2 shows the distribution of galaxies of the NewHorizon
simulation in the $\Sigma_{\rm H_{2}}$-$\Sigma_{\rm SFR}$ plane, in different
stellar mass bins (columns) at different redshifts (rows), colour-coded by
their specific star formation rate (sSFR $=\leavevmode\nobreak\
\mbox{SFR}/\mbox{$M_{\star}$}$), computed for the entire galaxy with the same
timescale as $\Sigma_{\rm 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
$\Sigma_{\rm H_{2}}$ and $\Sigma_{\rm 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 parameter space towards lower values of
$\Sigma_{\rm H_{2}}$ and $\Sigma_{\rm SFR}$ with decreasing redshift: fewer
and fewer galaxies are found far above the canonical KS relation444With the
term canonical, we refer to the sequence of normal, star-forming discs, as
defined in e.g. Daddi et al. (2010), i.e. $\Sigma_{\rm SFR}$ $\propto$
$\Sigma_{\rm{gas}}^{1.4}$. (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\leq 2$. It also decreases with increasing stellar mass at
each redshift. The sSFR of galaxies strongly correlates with $\Sigma_{\rm
SFR}$ at each redshift and in each stellar mass bin, essentially because the
sSFR is computed using stars with the same age as $\Sigma_{\rm SFR}$ ($<10$
Myr). This strong correlation vanishes when considering older stars (e.g.
$<100$ Myr). As a consequence, at fixed $\Sigma_{\rm 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 mass555As such, lines of constant depletion time have a slope
of unity in the KS plane. A slight difference exists with the observed
sequence of starbursts which yields a slope of 1.4 in the KS plane (Daddi et
al., 2010). This nuance is yet to be understood., 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
($\Sigma_{\rm\textsc{H\,i}+\rm{H}_{2}}$) is considered instead of the
molecular gas alone ($\Sigma_{\rm H_{2}}$, see Fig. 10), 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 $\Sigma_{\rm H_{2}}$– $\Sigma_{\rm
SFR}$ and $\Sigma_{\rm\textsc{H\,i}+\rm{H}_{2}}$– $\Sigma_{\rm SFR}$ planes,
but now for the entire population of galaxies at different redshifts, by
stacking all galaxies from Figs. 2 and 10, 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 shallower slope between the SFR density
and the molecular gas, than with the neutral gas.

The shallower slope of the correlation with H2 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\sim 1.5$; Daddi et al., 2010),
while the inclusion of dwarf galaxies tends to produce a shallower slope of
$\sim 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.

###  3.2 Emergence of the KS relation

Figure 4: Evolution of the slope of the fits of the KS relation from Figs. 2
and 10, 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\approx 2\\--3$. Low mass galaxies do not show signs of convergence
toward a fixed slope: their KS relation gets continuously shallower after
$z\approx 2\\--3$.

Figure 5: Same as Fig. 4, but showing the dispersion around the best fit. Only
the vertical dispersion, i.e. in $\log(\mbox{$\Sigma_{\rm SFR}$})$ is
considered here.

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 ($\Sigma_{\rm SFR}\propto\Sigma_{\rm gas}^{a}$)
fitted with the Bayesian linear regression method, 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 $\Sigma_{\rm
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\sim 2$), this regime disappears
at all masses but reappears at $z\lesssim 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 \+ H2), 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
star-forming 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\approx 2$ for the overall
population (red thick line on Fig. 4) and the most massive galaxies
($M_{\star}\gtrsim 10^{8}{\rm M}_{\odot}$ 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\approx 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 $\gtrsim 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\approx 2-3$ followed by a
relatively flat plateau. Hence, the relation for the low mass galaxies becomes
simultaneously tighter and shallower below $z\approx 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 low-mass 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)666As mentioned in
the previous section, a quantitative agreement cannot be reached due to the
diversity of fitting methods employed in the literature.. 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 $\Sigma_{\rm SFR}$ and the neutral gas
surface density $\Sigma_{\rm\textsc{H\,i}+\rm{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 $\Sigma_{\rm\textsc{H\,i}+\rm{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 $\Sigma_{\rm 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 $\Sigma_{\rm 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 metallicities777The metallicity is
computed within the $R_{\rm eff}$. For the gas metallicity, only the neutral
phase is considered. 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\gtrsim
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.

###  3.3 Physical drivers of the KS relation

####  3.3.1 Molecular and total gas content

Figure 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 $\Sigma_{\rm
H_{2}}$ at all redshifts. Baryonic molecular fraction (middle panel)
correlates with $\Sigma_{\rm H_{2}}$ at high redshifts, below z$\sim$2 this
correlation weakens and is essentially carried by low-mass galaxies. The
correlation between the neutral gas fraction and $\Sigma_{\rm H_{2}}$ is
apparent only at $z=4$. At $z\lesssim$2 the trend reverses such that this
fraction decreases with increasing $\Sigma_{\rm H_{2}}$.

We now investigate whether the relation between $\Sigma_{\rm H_{2}}$ and
$\Sigma_{\rm SFR}$ is driven by the gas content of galaxies.

Figure 6 (top panel) shows the molecular gas fraction, defined as the fraction
of H2 mass over the neutral gas, i.e.
M${}_{\rm{H}_{2}}$/(M${}_{\textsc{H\,i}}+\rm{M}_{\rm{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 $\Sigma_{\rm 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. 11, 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${}_{\rm{H}_{2}}$/($\rm{M}_{\textsc{H\,i}}+\rm{M}_{\rm{H}_{2}}+\rm{M}_{\star}$).
The correlation between the baryonic molecular gas fraction and $\Sigma_{\rm
H_{2}}$ is maintained at all redshifts (Fig. 6, middle panel), although it
weakens at $z\lesssim 2$ when it is only carried by the low-mass galaxies –
the massive galaxies having very low baryonic molecular gas fraction,
independently of $\Sigma_{\rm H_{2}}$ (see Fig. 12). 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.
(M${}_{\rm\textsc{H\,i}}+\rm{M}_{\rm{H}_{2}})/(\rm{M}_{\textsc{H\,i}}+\rm{M}_{\rm{H}_{2}}+\rm{M}_{\star})$.
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.
13). At high redshift ($z=4$), the neutral gas fraction increases with
$\Sigma_{\rm H_{2}}$ for all stellar masses, but this correlation disappears
at later epochs, where this fraction varies weakly with $\Sigma_{\rm H_{2}}$
at fixed stellar mass but varies strongly with stellar mass. The combination
of these relations between the neutral gas fraction and $\Sigma_{\rm H_{2}}$
at high $z$, and $M_{\star}$ 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 $\Sigma_{\rm 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 ($\Sigma_{\rm H_{2}}$), but have close to no
influence on the other ($\Sigma_{\rm 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.

####  3.3.2 Turbulence

Figure 7: Same as Fig. 3, but colour-coded by the Mach number $\mathcal{M}$.
At all redshifts, $\mathcal{M}$ increases monotonically with increasing both
$\Sigma_{\rm H_{2}}$ and $\Sigma_{\rm SFR}$. At fixed $\Sigma_{\rm H_{2}}$,
$\mathcal{M}$ is correlated with $\Sigma_{\rm SFR}$.

Figure 8: Normalised distributions of gas velocity dispersion (left),
temperature (middle), and Mach number (right) in different stellar mass bins
(colored lines, as indicated in the legend, reporting the logarithm of
$M_{\star}$ in units of $\rm{M}_{\odot}$) and redshifts (indicated in the left
upper part of each panel). Vertical lines represent medians of distributions.

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 $\mathcal{M}$. At all redshifts,
$\mathcal{M}$ evolves monotonically along the best-fit relation in the KS
parameter space, by increasing with $\Sigma_{\rm H_{2}}$ and $\Sigma_{\rm
SFR}$. Furthermore, at fixed $\Sigma_{\rm H_{2}}$, $\Sigma_{\rm SFR}$ is
positively correlated with $\mathcal{M}$.

As shown in Fig. 14, 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 $\mbox{$\mathcal{M}$}\propto\sigma_{\rm
vel}/\sqrt{T}$, and the underlying quantities which are the gas velocity
dispersion ($\sigma_{\rm vel}$) and the temperature ($T$)888As in the case of
velocity dispersion and Mach number, the temperature is also computed as a
mass-weighted average of gas cells within the $R_{\rm eff}$. Therefore, the
temperature values are not directly comparable to typical values within
individual molecular clouds..

While the velocity dispersion decreases with redshift for all mass bins, only
massive galaxies maintain high values at low $z$ ($\sim 10\ \rm{km\ s}^{-1}$),
leading to significantly higher median values. This is confirmed by
observations at low redshifts of lower dispersion in dwarf galaxies ($\sim 1\
\rm{km\ s}^{-1}$) than in massive galaxies ($\sim 10\ \rm{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 $\sigma_{\rm vel}$ in the KS plane
(Fig. 15) reveal complex correlations with the indicators of star formation:
while galaxies with short depletion times tend to have high $\sigma_{\rm
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\lesssim 1$). Contrary to the velocity
dispersion, there is no relation between the star formation indicators and the
temperature in the KS plane (Fig. 16), 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
($\sigma_{\rm vel}$ and $T$) naturally combine to lead to a shift of the
distributions of $\mathcal{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. 14) reveals that the trends found in
velocity dispersion and in temperature conspire to give rise to a clear
evolution of $\mathcal{M}$ along the KS relation, with ‘tighter correlation’
with $\Sigma_{\rm gas}$ and $\Sigma_{\rm SFR}$ than the individual
$\sigma_{\rm 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. 11) 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. 10), but are shifted toward the
low gas densities when considering the molecular phase only (Fig. 2). Finally,
the high $\Sigma_{\rm SFR}$ of these galaxies implies that this shift toward
low $\Sigma_{\rm 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 low-velocity dispersion, and the bulk of their distribution remains
centered around higher values ($\sim 10\ \rm{km\ s}^{-1}$) with little
evolution after $z\lesssim 2$. In other words, the lower end of the
distributions in $\sigma_{\rm 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\approx 2$, and the
absence of the convergence in low-mass galaxies, noted in Fig. 4.

##  4 Discussion

###  4.1 Scale and projection effects

Figure 9: Distribution of galaxies in the $M_{\rm 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 $\sigma$.
The coefficient $R$ is the Pearson correlation coefficient. At all redshifts,
more massive galaxies tend to have higher SFR and $M_{\rm 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_{\rm eff}$, SFR and
molecular gas mass.

So far, we have conducted our analysis using the observables $\Sigma_{\rm
H_{2}}$, or $\Sigma_{\rm\textsc{H\,i}+\rm{H}_{2}}$, and $\Sigma_{\rm 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, colour-coded by
mass (top) and size $R_{\rm 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_{\rm H_{2}}$, and effective radius of
galaxies. Furthermore, all the trends (or the absence thereof) with molecular
and total gas fractions, and $\mathcal{M}$ seen in the KS parameter space are
retrieved with deprojected quantities, as shown in Fig. 22. 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
$\Sigma_{\rm H_{2}}$, $\Sigma_{\rm\textsc{H\,i}+\rm{H}_{2}}$, and $\Sigma_{\rm
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, small-scale
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.

###  4.2 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 ($\sim$ 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
$\gtrsim 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 ($\approx 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 $\sim 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., 2020,
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 possibly 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.

###  4.3 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\approx 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\gtrsim 2$) in Fig. 7, where the most turbulent galaxies lie around the
starburst relation from Daddi et al. (2010), about 1 dex above the canonical
KS relation. This is in line with the conclusions of Renaud et al. (2014,
2019a, 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 $\Sigma_{\rm 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).

##  5 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 $\mbox{$\Sigma_{\rm SFR}$}\propto\Sigma_{\rm
gas}^{a}$ with a slope $a\approx 1.4$ emerges at $z\approx 2\\--3$ for the
most massive half of the galaxy population ($\mbox{$M_{\star}$}\gtrsim
10^{8}\,\mbox{$\rm{M}_{\odot}$}$ 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 $\Sigma_{\rm 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 ($\mathcal{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.

###### Acknowledgements.

This work was granted access to the high-performance computing resources of
CINES under the allocations c2016047637 and A0020407637 from GENCI, and KISTI
(KSC-2017-G2-0003). Large data transfer was supported by KREONET, which is
managed and operated by KISTI. This work relied on the HPC resources of the
Horizon Cluster hosted by Institut d’Astrophysique de Paris. We warmly thank
S. Rouberol for running the cluster on which the simulation was post-
processed. FR, OA, EA and ASO acknowledge support from the Knut and Alice
Wallenberg Foundation and the Swedish Research Council (grant 2019-04659). FR
acknowledges the support provided by the University of Strasbourg Institute
for Advanced Study (USIAS), within the French national programme Investment
for the Future (Excellence Initiative) IdEx-Unistra. EA acknowledges support
from US NSF grant AST18-1546. SK acknowledges support from the STFC [grant
numbers ST/S00615X/1 and ST/X001318/1] and a Senior Research Fellowship from
Worcester College Oxford. TK was supported by the National Research Foundation
of Korea (NRF) grant funded by the Korea government (No. 2020R1C1C1007079).
S.K.Y. acknowledges support from the Korean National Research Foundation
(2020R1A2C3003769). This work was supported in part by the Korean National
Research Foundation (2022R1A6A1A03053472). This work is partially supported by
the grant Segal ANR-19-CE31-0017 of the French Agence Nationale de la
Recherche and by the National Science Foundation under Grant No. NSF
PHY-1748958.

## References

  * Agertz & Kravtsov (2015) Agertz, O. & Kravtsov, A. V. 2015, ApJ, 804, 18
  * Agertz et al. (2021) Agertz, O., Renaud, F., Feltzing, S., et al. 2021, MNRAS, 503, 5826
  * Agertz et al. (2011) Agertz, O., Teyssier, R., & Moore, B. 2011, MNRAS, 410, 1391
  * Alves et al. (2007) Alves, J., Lombardi, M., & Lada, C. J. 2007, A&A, 462, L17
  * Andersson et al. (2021) Andersson, E. P., Renaud, F., & Agertz, O. 2021, MNRAS, 502, L29
  * André et al. (2010) André, P., Men’shchikov, A., Bontemps, S., et al. 2010, A&A, 518, L102
  * Appleby et al. (2020) Appleby, S., Davé, R., Kraljic, K., Anglés-Alcázar, D., & Narayanan, D. 2020, MNRAS, 494, 6053
  * Aubert et al. (2004) Aubert, D., Pichon, C., & Colombi, S. 2004, MNRAS, 352, 376
  * Bigiel et al. (2008) Bigiel, F., Leroy, A., Walter, F., et al. 2008, AJ, 136, 2846
  * Bigiel et al. (2011) Bigiel, F., Leroy, A. K., Walter, F., et al. 2011, ApJ, 730, L13
  * Bolatto et al. (2013) Bolatto, A. D., Wolfire, M., & Leroy, A. K. 2013, ARA&A, 51, 207
  * Bouché et al. (2007) Bouché, N., Cresci, G., Davies, R., et al. 2007, ApJ, 671, 303
  * Bournaud et al. (2010) Bournaud, F., Elmegreen, B. G., Teyssier, R., Block, D. L., & Puerari, I. 2010, MNRAS, 409, 1088
  * Brucy et al. (2020) Brucy, N., Hennebelle, P., Bournaud, F., & Colling, C. 2020, ApJ, 896, L34
  * Brucy et al. (2023) Brucy, N., Hennebelle, P., Colman, T., & Iteanu, S. 2023, arXiv e-prints, arXiv:2305.18012
  * Ciesla et al. (2023) Ciesla, L., Gómez-Guijarro, C., Buat, V., et al. 2023, A&A, 672, A191
  * Cormier et al. (2014) Cormier, D., Madden, S. C., Lebouteiller, V., et al. 2014, A&A, 564, A121
  * Daddi et al. (2010) Daddi, E., Elbaz, D., Walter, F., et al. 2010, ApJ, 714, L118
  * Dalgarno & McCray (1972) Dalgarno, A. & McCray, R. A. 1972, ARA&A, 10, 375
  * de los Reyes & Kennicutt (2019) de los Reyes, M. A. C. & Kennicutt, Robert C., J. 2019, ApJ, 872, 16
  * Dekel et al. (2019) Dekel, A., Sarkar, K. C., Jiang, F., et al. 2019, MNRAS, 488, 4753
  * Dubois et al. (2021) Dubois, Y., Beckmann, R., Bournaud, F., et al. 2021, A&A, 651, A109
  * Dubois et al. (2010) Dubois, Y., Devriendt, J., Slyz, A., & Teyssier, R. 2010, MNRAS, 409, 985
  * Dubois et al. (2012) Dubois, Y., Devriendt, J., Slyz, A., & Teyssier, R. 2012, MNRAS, 420, 2662
  * Dubois et al. (2014) Dubois, Y., Pichon, C., Welker, C., et al. 2014, MNRAS, 444, 1453
  * Federrath & Klessen (2012) Federrath, C. & Klessen, R. S. 2012, ApJ, 761, 156
  * Federrath et al. (2008) Federrath, C., Klessen, R. S., & Schmidt, W. 2008, ApJ, 688, L79
  * Feldmann et al. (2012) Feldmann, R., Gnedin, N. Y., & Kravtsov, A. V. 2012, ApJ, 758, 127
  * Filho et al. (2016) Filho, M. E., Sánchez Almeida, J., Amorín, R., et al. 2016, ApJ, 820, 109
  * Freundlich et al. (2013) Freundlich, J., Combes, F., Tacconi, L. J., et al. 2013, A&A, 553, A130
  * Freundlich et al. (2019) Freundlich, J., Combes, F., Tacconi, L. J., et al. 2019, A&A, 622, A105
  * Genzel et al. (2010) Genzel, R., Tacconi, L. J., Gracia-Carpio, J., et al. 2010, MNRAS, 407, 2091
  * Genzel et al. (2013) Genzel, R., Tacconi, L. J., Kurk, J., et al. 2013, ApJ, 773, 68
  * Gómez-Guijarro et al. (2022) Gómez-Guijarro, C., Elbaz, D., Xiao, M., et al. 2022, A&A, 659, A196
  * Grisdale et al. (2022) Grisdale, K., Hogan, L., Rigopoulou, D., et al. 2022, MNRAS, 513, 3906
  * Haardt & Madau (1996) Haardt, F. & Madau, P. 1996, ApJ, 461, 20
  * Hennebelle & Chabrier (2011) Hennebelle, P. & Chabrier, G. 2011, ApJ, 743, L29
  * Hogg et al. (2010) Hogg, D. W., Bovy, J., & Lang, D. 2010, arXiv e-prints, arXiv:1008.4686
  * Hopkins et al. (2014) Hopkins, P. F., Kereš, D., Oñorbe, J., et al. 2014, MNRAS, 445, 581
  * Hu et al. (2022) Hu, C.-Y., Smith, M. C., Teyssier, R., et al. 2022, arXiv e-prints, arXiv:2208.10528
  * Hunter et al. (2021) Hunter, D. A., Elmegreen, B. G., Archer, H., Simpson, C. E., & Cigan, P. 2021, AJ, 161, 175
  * Isobe et al. (1990) Isobe, T., Feigelson, E. D., Akritas, M. G., & Babu, G. J. 1990, ApJ, 364, 104
  * Jackson et al. (2021a) Jackson, R. A., Kaviraj, S., Martin, G., et al. 2021a, MNRAS, 502, 1785
  * Jackson et al. (2021b) Jackson, R. A., Martin, G., Kaviraj, S., et al. 2021b, MNRAS, 502, 4262
  * Kaviraj et al. (2017) Kaviraj, S., Laigle, C., Kimm, T., et al. 2017, MNRAS, 467, 4739
  * Kelly (2007) Kelly, B. C. 2007, ApJ, 665, 1489
  * Kennicutt (1989) Kennicutt, Robert C., J. 1989, ApJ, 344, 685
  * Kennicutt (1998) Kennicutt, Robert C., J. 1998, ApJ, 498, 541
  * Kennicutt et al. (2007) Kennicutt, Robert C., J., Calzetti, D., Walter, F., et al. 2007, ApJ, 671, 333
  * Kennicutt & De Los Reyes (2021) Kennicutt, Robert C., J. & De Los Reyes, M. A. C. 2021, ApJ, 908, 61
  * Kennicutt & Evans (2012) Kennicutt, R. C. & Evans, N. J. 2012, ARA&A, 50, 531
  * Kimm & Cen (2014) Kimm, T. & Cen, R. 2014, ApJ, 788, 121
  * Kimm et al. (2015) Kimm, T., Cen, R., Devriendt, J., Dubois, Y., & Slyz, A. 2015, MNRAS, 451, 2900
  * Kimm et al. (2017) Kimm, T., Katz, H., Haehnelt, M., et al. 2017, MNRAS, 466, 4826
  * Komatsu et al. (2011) Komatsu, E., Smith, K. M., Dunkley, J., et al. 2011, ApJS, 192, 18
  * Kraljic et al. (2014) Kraljic, K., Renaud, F., Bournaud, F., et al. 2014, ApJ, 784, 112
  * Kravtsov (2003) Kravtsov, A. V. 2003, ApJ, 590, L1
  * Krumholz & McKee (2005) Krumholz, M. R. & McKee, C. F. 2005, ApJ, 630, 250
  * Leroy et al. (2008) Leroy, A. K., Walter, F., Brinks, E., et al. 2008, AJ, 136, 2782
  * Leroy et al. (2013) Leroy, A. K., Walter, F., Sandstrom, K., et al. 2013, AJ, 146, 19
  * Liu et al. (2015) Liu, L., Gao, Y., & Greve, T. R. 2015, ApJ, 805, 31
  * Martin et al. (2021) Martin, G., Jackson, R. A., Kaviraj, S., et al. 2021, MNRAS, 500, 4937
  * Matzner & McKee (2000) Matzner, C. D. & McKee, C. F. 2000, ApJ, 545, 364
  * Narayanan et al. (2011) Narayanan, D., Krumholz, M., Ostriker, E. C., & Hernquist, L. 2011, MNRAS, 418, 664
  * Onodera et al. (2010) Onodera, S., Kuno, N., Tosaki, T., et al. 2010, ApJ, 722, L127
  * Orr et al. (2018) Orr, M. E., Hayward, C. C., Hopkins, P. F., et al. 2018, MNRAS, 478, 3653
  * Padoan & Nordlund (2011) Padoan, P. & Nordlund, Å. 2011, ApJ, 730, 40
  * Park et al. (2021) Park, M. J., Yi, S. K., Peirani, S., et al. 2021, ApJS, 254, 2
  * Rawle et al. (2014) Rawle, T. D., Egami, E., Bussmann, R. S., et al. 2014, ApJ, 783, 59
  * Renaud et al. (2021a) Renaud, F., Agertz, O., Andersson, E. P., et al. 2021a, MNRAS, 503, 5868
  * Renaud et al. (2021b) Renaud, F., Agertz, O., Read, J. I., et al. 2021b, MNRAS, 503, 5846
  * Renaud et al. (2019a) Renaud, F., Bournaud, F., Agertz, O., et al. 2019a, A&A, 625, A65
  * Renaud et al. (2019b) Renaud, F., Bournaud, F., Daddi, E., & Weiß, A. 2019b, A&A, 621, A104
  * Renaud et al. (2014) Renaud, F., Bournaud, F., Kraljic, K., & Duc, P. A. 2014, MNRAS, 442, L33
  * Renaud et al. (2012) Renaud, F., Kraljic, K., & Bournaud, F. 2012, ApJ, 760, L16
  * Renaud et al. (2022) Renaud, F., Segovia Otero, Á., & Agertz, O. 2022, MNRAS, 516, 4922
  * Rosdahl & Blaizot (2012) Rosdahl, J. & Blaizot, J. 2012, MNRAS, 423, 344
  * Roychowdhury et al. (2017) Roychowdhury, S., Chengalur, J. N., & Shi, Y. 2017, A&A, 608, A24
  * Salmi et al. (2012) Salmi, F., Daddi, E., Elbaz, D., et al. 2012, ApJ, 754, L14
  * Schmidt (1959) Schmidt, M. 1959, ApJ, 129, 243
  * Schruba et al. (2011) Schruba, A., Leroy, A. K., Walter, F., et al. 2011, AJ, 142, 37
  * Segovia Otero et al. (2022) Segovia Otero, Á., Renaud, F., & Agertz, O. 2022, MNRAS, 516, 2272
  * Semenov et al. (2019) Semenov, V. A., Kravtsov, A. V., & Gnedin, N. Y. 2019, ApJ, 870, 79
  * Sharon et al. (2013) Sharon, C. E., Baker, A. J., Harris, A. I., & Thomson, A. P. 2013, ApJ, 765, 6
  * Shetty et al. (2013) Shetty, R., Kelly, B. C., & Bigiel, F. 2013, MNRAS, 430, 288
  * Sun et al. (2023) Sun, J., Leroy, A. K., Ostriker, E. C., et al. 2023, ApJ, 945, L19
  * Sutherland & Dopita (1993) Sutherland, R. S. & Dopita, M. A. 1993, ApJS, 88, 253
  * Tacconi et al. (2010) Tacconi, L. J., Genzel, R., Neri, R., et al. 2010, Nature, 463, 781
  * Tacconi et al. (2018) Tacconi, L. J., Genzel, R., Saintonge, A., et al. 2018, ApJ, 853, 179
  * Tacconi et al. (2013) Tacconi, L. J., Neri, R., Genzel, R., et al. 2013, ApJ, 768, 74
  * Teng et al. (2023) Teng, Y.-H., Sandstrom, K. M., Sun, J., et al. 2023, arXiv e-prints, arXiv:2304.04732
  * Teyssier (2002) Teyssier, R. 2002, A&A, 385, 337
  * Teyssier et al. (2011) Teyssier, R., Moore, B., Martizzi, D., Dubois, Y., & Mayer, L. 2011, MNRAS, 414, 195
  * Thornton et al. (1998) Thornton, K., Gaudlitz, M., Janka, H. T., & Steinmetz, M. 1998, ApJ, 500, 95
  * Trebitsch et al. (2017) Trebitsch, M., Blaizot, J., Rosdahl, J., Devriendt, J., & Slyz, A. 2017, MNRAS, 470, 224
  * Trebitsch et al. (2021) Trebitsch, M., Dubois, Y., Volonteri, M., et al. 2021, A&A, 653, A154
  * Vazquez-Semadeni (1994) Vazquez-Semadeni, E. 1994, ApJ, 423, 681
  * Volonteri et al. (2020) Volonteri, M., Pfister, H., Beckmann, R. S., et al. 2020, MNRAS, 498, 2219

##  Appendix A Galaxy distribution in KS plane

###  A.1 Dependence on the total gas

Figure 10 shows the distribution of galaxies in the
$\Sigma_{\rm\textsc{H\,i}+\rm{H}_{2}}$-$\Sigma_{\rm SFR}$ parameter space at
different redshifts (rows) and in four equally-populated stellar mass
quartiles at each redshift (columns). Trends with the sSFR (given by the
colour coding) seen when considering the molecular gas (see Fig. 2) persist.
At each redshift and each stellar mass bin, the sSFR of galaxies strongly
correlated with $\Sigma_{\rm SFR}$, and regardless of the stellar mass, the
distributions of galaxies move within the KS parameter space towards lower
values of $\Sigma_{\rm\textsc{H\,i}+\rm{H}_{2}}$ and $\Sigma_{\rm SFR}$ with
decreasing redshift.

Figure 10: Same as Fig. 2, but for the neutral gas surface density
$\Sigma_{\rm\textsc{H\,i}+\rm{H}_{2}}$.

###  A.2 Molecular and total gas fractions

Figure 11 shows the molecular gas fraction, defined as the fraction of H2 mass
over the neutral gas, i.e.
M${}_{\rm{H}_{2}}$/(M${}_{\textsc{H\,i}}+\rm{M}_{\rm{H}_{2}}$), and its
evolution within the KS plane with the cosmic time for different stellar mass
bins.

Figure 12 shows the molecular gas content of galaxies defined in terms of
baryonic fraction, i.e.
M${}_{\rm{H}_{2}}/(\rm{M}_{\textsc{H\,i}}+\rm{M}_{\rm{H}_{2}}+\rm{M}_{\star})$,
as a function of cosmic time and stellar mass.

Figure 13 shows the neutral baryonic gas fraction, i.e.
(M${}_{\rm\textsc{H\,i}}+\rm{M}_{\rm{H}_{2}})/(\rm{M}_{\textsc{H\,i}}+\rm{M}_{\rm{H}_{2}}+\rm{M}_{\star})$,
at different redshifts and in different stellar mass bins.

Figure 11: Same as Fig. 2, but colour-coded by the molecular over cold gas
fraction.  Figure 12: Same as Fig. 2, but colour-coded by the baryonic
molecular gas fraction.  Figure 13: Same as Fig. 2, but colour-coded by the
neutral gas fraction.

###  A.3 Turbulence

Figure 14 shows the distributions of galaxies in the $\Sigma_{\rm
H_{2}}$-$\Sigma_{\rm SFR}$ parameter space as a function of Mach number
$\mathcal{M}$. At all redshifts and in all stellar mass bins $\mathcal{M}$
monotonically increases with increasing $\Sigma_{\rm H_{2}}$ and $\Sigma_{\rm
SFR}$. In addition, at fixed $\Sigma_{\rm H_{2}}$, $\Sigma_{\rm SFR}$
correlates with $\mathcal{M}$ at all redshifts, driving the offset of the most
turbulent galaxies above the observed KS relation onto the starburst regime
(dotted line) and beyond, as visible at high redshift.

At $z<2$, this regime is very sparsely populated. 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.

Figures 15 and 16 correspond to the distribution of galaxies in the
$\Sigma_{\rm H_{2}}$-$\Sigma_{\rm SFR}$ parameter space color-coded by their
gas velocity dispersion and temperature, respectively.

Figure 14: Same as Fig. 2, but colour-coded by the Mach number.  Figure 15:
Same as Fig. 2, but colour-coded by the gas velocity dispersion.  Figure 16:
Same as Fig. 2, but colour-coded by the gas temperature.

##  Appendix B Fitting methods

We fit the distribution of galaxies within the $\log\Sigma_{\rm
gas}-\log\Sigma_{\rm SFR}$ plane with the linear relation $\log\Sigma_{\rm
SFR}=a(\log\Sigma_{\rm gas})+b$, with the best-fit values for the slope $a$
and intercept $b$ given in each panel. We compared three different fitting
methods, largely used and advocated in the literature. These are the ordinary
least square (OLS) technique, OLS bisector (OLS-bis) technique (Isobe et al.
1990), and Bayesian linear regression (BLR) with the Student’s
t-distribution999We use the probabilistic programming package for Python
PyMC3, available at https://docs.pymc.io/en/v3/index.html.  for the likelihood
to minimise the impact of outliers. Normal likelihood gives the data used in
this work identical results as OLS. The same is true for the hierarchical
Bayesian model linmix (Kelly 2007) that we have tested as well101010We use the
Python package implementing the linmix algorithm available on github at
https://github.com/jmeyers314/linmix.. Overall, comparing the three above-
mentioned methods, for the data used in this work, we conclude that:

  * •

For the fits involving the entire population of galaxies (e.g. Fig. 3), the
derived slopes show the following hierarchy $a_{\rm{OLS}}\leq
a_{\rm{BLR}}<a_{\rm{OLS-bis}}$ at all redshifts (with $a_{\rm{OLS}}$ and
$a_{\rm{BLR}}$ being often consistent within the error-bars) and follows the
same trend with redshift. The dispersion around the best fit (i.e. the rms in
$\log(\mbox{$\Sigma_{\rm SFR}$})$) follows the same trend. This is true for
both $\Sigma_{\rm H_{2}}$\- and
$\Sigma_{\rm\textsc{H\,i}+\rm{H}_{2}}$-$\Sigma_{\rm SFR}$ fits.

  * •

For the fits at different mass bins and redshifts (e.g. Fig. 2) the same
hierarchy is followed by slopes with $a_{\rm{OLS}}$ and $a_{\rm{BLR}}$ being
most of the time within uncertainties of each other, for both $\Sigma_{\rm
H_{2}}$ and $\Sigma_{\rm\textsc{H\,i}+\rm{H}_{2}}$. The same applies to the
dispersion around the best fit.

All three explored fitting methods provide qualitatively similar trends for
the slopes and dispersions around the best fit as a function of redshift and
stellar mass. There are however notable quantitative differences, with slopes
being systematically higher for the bisector OLS method compared to the two
other methods. Similarly, the dispersion around the best fit is always larger
for the bisector OLS method compared to the standard OLS and Bayesian methods.

###  B.1 Distribution of galaxies in the KS plane

Figure 17 shows the distribution of galaxies in the $\Sigma_{\rm
SFR}$-$\Sigma_{\rm H_{2}}$ plane, as in Fig. 2, but reporting the results of
the OLS bisector fits for comparison (red solid lines) to the Bayesian fit
(black dash-dotted lines). The corresponding best-fit values for the OLS fit
slope $a$ and the standard deviation of residuals $\sigma$ are shown in the
bottom right corners of each panel.

Figure 18 is the same as Fig. 17, but showing instead the distribution of
galaxies in the $\Sigma_{\rm SFR}$-$\Sigma_{\rm\textsc{H\,i}+\rm{H}_{2}}$
parameter space.

Figure 17: Same as Fig. 2, but showing in addition, OLS bisector fits for
comparison (red solid lines). The corresponding best-fit values for slope $a$
and the standard deviation of residuals $\sigma$ are shown in the bottom right
corners of each panel.  Figure 18: Same as Fig. 17, but for the neutral gas
surface density $\Sigma_{\rm\textsc{H\,i}+\rm{H}_{2}}$.

###  B.2 Integrated KS

Figure 19 is equivalent to Fig. 3, but showing the results of the OLS bisector
fit (solid red lines), together with its best-fit values that can be compared
to the Bayesian fit (black dash-dotted lines).

Figure 19: Same as Fig. 3, but showing in addition, OLS bisector fits for
comparison (red solid lines). The corresponding best-fit values for slope $a$
and the standard deviation of residuals $\sigma$ are shown in the bottom right
corners of each panel.

###  B.3 Emergence of the KS relation

Figures 20 and 21 are the same as Figs. 4 and 5, but for the slopes and
dispersions around the best-fit relation using the OLS technique.

Figure 20: Same as Fig. 4, but using the OLS bisector fitting method (values
from Figs. 17 and 18).

Figure 21: Same as Fig. 5, but using the OLS bisector fitting method (values
from Fig. 17 and 18).

To conclude, our analysis reveals that the choice of the fitting method
impacts the quantitative conclusions for all the measurements and diagnostics
presented, but that the overall qualitative trends hold.

##  Appendix C Deprojected KS relation

Figure 22 displays the dependence of molecular gas fraction(top row), the
neutral gas fraction (middle row), and Mach number (bottom row) of galaxies at
different redshifts (columns) on the $M_{\rm H_{2}}$-SFR plane. The trends
seen for the projected quantities $\Sigma_{\rm H_{2}}$ and $\Sigma_{\rm SFR}$
are reproduced, i.e., the only physical parameter that evolves with
$\Sigma_{\rm H_{2}}$ and $\Sigma_{\rm SFR}$ is the Mach number $\mathcal{M}$.

We find that the results identified in the projected, i.e. observable,
versions of the KS plane are also present in their deprojected counterparts,
which means that our conclusions are not affected by projection artefacts.

Figure 22: Same as Fig. 9, but colour-coded by the molecular gas mass over the
total neutral mass (top row), i.e.
M${}_{\rm{H}_{2}}$/(M${}_{\textsc{H\,i}}+\rm{M}_{\rm{H}_{2}}$), the neutral
baryonic gas fraction (middle row), i.e.
(M${}_{\rm\textsc{H\,i}}+\rm{M}_{\rm{H}_{2}})/(\rm{M}_{\textsc{H\,i}}+\rm{M}_{\rm{H}_{2}}+\rm{M}_{\star})$,
and the Mach number (bottom row). For comparison, dotted-dash lines correspond
to the OLS bisector fits. The corresponding best-fit values for slope $a$ and
the standard deviation of residuals $\sigma$ are shown in the top left corners
of each panel.