# The ALMA-ALPAKA survey II. Evolution of turbulence in galaxy disks across cosmic time: difference between cold and warm gas

F Rizzo, C Bacchini, M Kohandel, L Di Mascolo, F Fraternali, F

Roman-Oliveira, A Zanella, G Popping, F Valentino, G Magdis, et al.

► To cite this version:

F Rizzo, C Bacchini, M Kohandel, L Di Mascolo, F Fraternali, et al. The ALMA-ALPAKA survey II. Evolution of turbulence in galaxy disks across cosmic time: difference between cold and warm gas. Astronomy & Astrophysics - A&A, 2024, 689, pp.A273. 10.1051/0004-6361/202450455. hal-04661932

HAL Id: hal-04661932

https://hal.science/hal-04661932v1

Submitted on 16 Jan 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# The ALMA-ALPAKA survey

## II. Evolution of turbulence in galaxy disks across cosmic time: Difference between cold and warm gas

F. Rizzo<fnref n="1" />,<fnref n="2" />,<fnref n="3" />,*, C. Bacchini<fnref n="4" />, M. Kohandel<fnref n="5" />, L. Di Mascolo<fnref n="6" />, F. Fraternali<fnref n="3" />, F. Roman-Oliveira<fnref n="3" />, A. Zanella<fnref n="7" />, G. Popping<fnref n="8" />, F. Valentino<fnref n="8" />, G. Magdis<fnref n="1" />,<fnref n="9" />, and K. Whitaker<fnref n="10" />

Received 20 April 2024 / Accepted 4 July 2024

## ABSTRACT

The gas in the interstellar medium (ISM) of galaxies is supersonically turbulent. Measurements of turbulence typically rely on cold gas emission lines for low-z galaxies and warm ionized gas observations for $z > 0$ galaxies. Studies of warm gas kinematics at $z > 0$ conclude that the turbulence strongly evolves as a function of redshift, due to the increasing impact of gas accretion and mergers in the early Universe. However, recent findings suggest potential biases in turbulence measurements derived from ionized gas at high-z, impacting our understanding of turbulence origin, ISM physics and disk formation. We investigate the evolution of turbulence using velocity dispersion ($\sigma$) measurements from cold gas tracers (i.e., CO, [CI], [CII]). The initial dataset comprises 17 galaxy disks with high data quality from the ALPAKA sample, supplemented with galaxies from the literature, resulting in a sample of 57 galaxy disks spanning the redshift range $z = 0-5$. This extended sample consists of main-sequence and starburst galaxies with stellar masses $\ge 10^{10} M_{\odot}$. The comparison with current $H\alpha$ kinematic observations and existing models demonstrates that the velocity dispersion inferred from cold gas tracers differ by a factor of $\approx 3$ from those obtained using emission lines tracing the warm, ionized gas. We show that stellar feedback is the main driver of turbulence measured from cold gas tracers and the physics of turbulence driving does not appear to evolve with time. This is fundamentally different from the conclusions of studies based on warm gas, which had to consider additional turbulence drivers to explain the high values of $\sigma$. We present a model predicting the redshift evolution of turbulence in galaxy disks, attributing the increase of $\sigma$ with redshift to the higher energy injected by supernovae due to the elevated star-formation rate in high-z galaxies. This supernova-driven model suggests that turbulence is lower in galaxies with lower stellar mass compared to those with higher stellar mass. Additionally, it forecasts the evolution of $\sigma$ in Milky-Way like progenitors.

**Key words.** galaxies: evolution – galaxies: high-redshift – galaxies: ISM – galaxies: kinematics and dynamics – galaxies: star formation

## 1. Introduction

The gas within the interstellar medium (ISM) of both low- (e.g., Elmegreen & Scalo 2004; Utomo et al. 2019; Bacchini et al. 2020a) and high-z (e.g., see Förster Schreiber & Wuyts 2020, for a review) galaxies is known to be supersonically turbulent. The strongest observational evidence for the presence of turbulent motions within the ISM is that the gas velocity dispersion $\sigma$, which measures the broadening of emission lines, exceeds the expected broadening from thermal motions. This additional broadening allows for quantifying the random, chaotic motions of the gas in the ISM.

For $z \approx 0$ galaxies, turbulence is usually measured through spatially-resolved observations of emission lines tracing the cold molecular and atomic gas (e.g., CO, HI Ianjamasimanana et al.

2015; Mogotsi et al. 2016; Utomo et al. 2019; Bacchini et al. 2020a). For $z > 0$ galaxies, obtaining spatially-resolved observations of cold gas emission lines requires exceptionally long integration times (Rizzo et al. 2023), even with the Atacama Large Millimeter Array (ALMA). Consequently, studies of gas turbulence at $z > 0$ have mainly relied on Integral Field Unit (IFU) or long-slit observations of emission lines tracing warm, ionized gas (e.g., $H\alpha$, [OII], [OIII]; Epinat et al. 2009; Stott et al. 2016; Wisnioski et al. 2015; Di Teodoro et al. 2016; Simons et al. 2016; Turner et al. 2017), that were possible only up to $z \approx 3$ before the launch of JWST. Crucially, these studies based on ionized gas measurements have assumed that, for high-z galaxies, the ionized gas kinematics (i.e., $\sigma$ and rotation velocity $V$) is consistent with the ones derived from cold gas (e.g., Genzel et al. 2011; Wisnioski et al. 2015; Johnson et al. 2018).

One of the main results from warm gas kinematic studies is that the turbulence strongly evolves as a function of

<footnote n="1">¹ Cosmic Dawn Center (DAWN), Copenhagen, Denmark</footnote>

<footnote n="2">² Niels Bohr Institute, University of Copenhagen, Jagtvej 128, 2200 Copenhagen N, Denmark</footnote>

<footnote n="3">³ Kapteyn Astronomical Institute, University of Groningen, Landleven 12, 9747 AD Groningen, The Netherlands</footnote>

<footnote n="4">⁴ Dark, Niels Bohr Institute, University of Copenhagen, Jagtvej 128, DK-2200 Copenhagen, Denmark</footnote>

<footnote n="5">⁵ Scuola Normale Superiore, Piazza dei Cavalieri 7, I-56126 Pisa, Italy</footnote>

<footnote n="6">⁶ Laboratoire Lagrange, Université Côte d'Azur, Observatoire de la Côte d'Azur, CNRS, Blvd de l'Observatoire, CS 34229, 06304 Nice Cedex 4, France</footnote>

<footnote n="7">⁷ Istituto Nazionale di Astrofisica, Vicolo dell'Osservatorio 5, 35122 Padova, Italy</footnote>

<footnote n="8">⁸ European Southern Observatory, Karl-Schwarzschild-Str. 2, D-85748 Garching bei München, Germany</footnote>

<footnote n="9">⁹ DTU-Space, Technical University of Denmark, Elektrovej 327, DK-2800 Kgs. Lyngby, Denmark</footnote>

<footnote n="10">¹⁰ Department of Astronomy, University of Massachusetts, Amherst, MA 01003, USA</footnote>

<footnote>* Corresponding author; rizzo@astro.rug.nl</footnote>
redshift (e.g., Epinat et al. 2009; Förster Schreiber et al. 2009; Wisnioski et al. 2015; Simons et al. 2017; Turner et al. 2017; Johnson et al. 2018; Birkin et al. 2023; Puglisi et al. 2023). Galaxy disks at $z = 2–3$ are more turbulent and less rotationally supported (i.e., values of rotation-to-velocity dispersion ratio, $V/\sigma \le 5$) than galaxy disks at $z < 1$. These results have been recently challenged by an increasing number of kinematic studies obtained from [CII] observations of $z \ge 4$ galaxies (e.g., Neeleman et al. 2020; Rizzo et al. 2020, 2021; Lelli et al. 2021; Tsukui & Iguchi 2021; Pope et al. 2023; Posses et al. 2023; Roman-Oliveira et al. 2024) using ALMA. These works find that the gas turbulence within $z > 4$ galaxy disks is, on average, lower than those expected by extrapolating the evolution of the $\sigma-z$ relation obtained from warm gas studies at $z = 0–3$ (Übler et al. 2019) and that $z \sim 4$ disks are dynamically cold, with $V/\sigma$ values of $\approx 10$ (e.g., Neeleman et al. 2020; Rizzo et al. 2020, 2021; Fraternali et al. 2021; Lelli et al. 2021; Tsukui & Iguchi 2021; Pope et al. 2023; Roman-Oliveira et al. 2024).

The above discrepancy between the kinematics obtained from [CII] observations and those extrapolated using ionized gas could be due to intrinsic differences between the cold and warm gas components. At $z = 0$, the values of $\sigma$ from cold gas (e.g., CO, HI) are, in fact, systematically lower than the ones obtained from $H\alpha$ (e.g., Epinat et al. 2008; Varidel et al. 2020; Law et al. 2022). Moreover, Kohandel et al. (2024) use synthetic observations of zoom-in simulations of galaxies at $z > 4$ to highlight a potential additional bias in turbulence measurements derived from warm gas velocity dispersion. In their studies, unresolved non-circular motions induced by outflows contaminate the measurements of turbulence from warm gas. As a result, warm gas $\sigma$ values are 2–3 times higher than the cold gas ones. In other words, in addition to the thermal and turbulent motions, the shape of the emission lines tracing the warm gas can be affected by the presence of non-circular (radial or vertical) motions, induced by outflows. These motions can mimic the broadening due to turbulence when the spectral or angular resolutions of the data are suboptimal.

In addition to the different kinematic tracers, potential contamination from outflows, and the redshift range (i.e., 1–3 for warm gas and $>4$ for cold gas), there are two additional potential reasons that can explain the discrepancy between warm and cold gas kinematic studies at $z > 0$: intrinsic differences in the galaxy populations or biases due to observational limitations. First, the results obtained with warm gas rely on observations of main-sequence galaxies (Rodighiero et al. 2011; Whitaker et al. 2012) with a wide range of stellar masses ($M_\star = 10^9–10^{10} M_\odot$), while the [CII] studies at $z > 4$ have focused on massive ($M_\star \ge 10^{10} M_\odot$) upper main-sequence or starburst galaxies. Second, the typical angular resolution of ground-based telescopes used to observe $z = 1–3$ galaxies (e.g., $\ge 0.5''$ in seeing-limited mode Wisnioski et al. 2015; Harrison et al. 2017; Turner et al. 2017; Puglisi et al. 2023) poses challenges due to beam smearing, potentially leading to degeneracies between rotation velocity and velocity dispersion (Di Teodoro & Fraternali 2015; Rizzo et al. 2022). Additionally, even when the galaxies are well spatially resolved thanks to JWST, the spectral resolutions of current facilities (e.g., FWHM $\ge 100$ km/s for JWST/NIRSpec) complicate the kinematic analysis, hampering the possibility of both cleanly disentangling between turbulent and non-circular motions and measuring values of $\sigma \le \text{FWHM}/2.35 \approx 40$ km/s. On the other hand, cold gas emission lines usually fall in the sub-mm to radio range of the electromagnetic spectrum. In these wavelength ranges, the typical spectral resolutions of instru-

ments is of a few km/s for low-$z$ data and 30 km/s for high-$z$ data. Combined with the typical angular resolution of $\approx 0.2''$, this kind of observations allows one to identify or isolate the circular motions in the rare cases when cold gas outflows are present (Rizzo et al. 2023; Roman-Oliveira et al. 2024).

Understanding whether the origin of the discrepancy between studies based on warm and cold gas tracers is driven by different galaxy populations or biases due the kinematic tracers and observational limitations has significant implications on galaxy formation and evolution, even besides the theories of disk formation. For instance, the clumpy morphology of high-$z$ galaxies and the formation of galaxy bulges at high-$z$ have been explained within the framework of highly turbulent disks (Bournaud et al. 2010). In this scenario, mergers and vigorous gas accretion give rise to highly turbulent motions and instabilities within disks. These instabilities cause fragmentation of the gas disk, the formation of massive star-forming regions, and their migration into the inner regions (Dekel et al. 2013; Bournaud et al. 2014). Finally, inaccurate measurements of the velocity dispersion can significantly affect the derivation of the mass and distribution of dark matter in galaxies, as their rotation curves must be corrected for the asymmetric drift which depends on $\sigma$ (see Lelli 2022, for details).

Given the aforementioned potential biases on the $\sigma$ measured using the warm, ionized gas, key questions arise: Do the cold gas kinematics exhibit differences compared to warm counterparts when galaxies with similar physical properties are considered? Do the dynamical properties of galaxies depend on the gas tracers used to probe the gas kinematics? Does turbulence in disks evolve with cosmic time when measured from cold gas tracers?

In this paper, we address these questions by analyzing dynamical properties obtained from cold gas tracers across a redshift range from 0 to 5. The paper is structured as follows: in Section 2, we describe the sample used in our analysis; in Section 3, we provide an empirical parametrization for the redshift evolution of $\sigma$ and $V/\sigma$ and we compare it with ionized gas studies from the literature. In Section 4, we briefly discuss the potential drivers of turbulence, while in Section 5, we present a physically-motivated model, calibrated on observations, predicting the turbulence evolution in star-forming galaxies. Caveats regarding the available values of $\sigma$ and star-formation rate (SFR) and a summary of this work are discussed in Sections 6 and 7, respectively.

## 2. ALPAKA and the extended sample

The primary objective of this paper is to study turbulence through the kinematics of cold gas across the widest possible redshift range. While previous works have explored turbulent motions and their driving mechanisms at $z > 0$, the focus has largely been on the range of $z = 4–5$ (e.g., Neeleman et al. 2020; Rizzo et al. 2021; Tsukui & Iguchi 2021; Roman-Oliveira et al. 2024). This particular redshift interval offers favorable conditions for conducting systematic, high-resolution observations crucial for accurate kinematic measurements, primarily due to the brightness of the [CII]–158 µm emission line (Carilli & Walter 2013). Observational campaigns targeting the intermediate-redshift range (i.e., $z = 0.5–4$) are significantly challenging. In fact, observations of [CII] are almost unfeasible because of the limited atmospheric transmission in the frequency range of the redshifted [CII] line and the restricted coverage of ALMA bands (Carilli & Walter 2013). CO or [CI] emission lines are promising alternatives for probing the kinematics of cold gas within this intermediate redshift range. Nonetheless, the faintness of CO and [CI] emission lines poses significant observational challenges. Even with ALMA's
<table><caption>**Table 1.** Datasets included in the extended sample.</caption><thead><tr><th>Name</th><th>Number of galaxies</th><th>Redshift</th><th>Kinematic tracer</th><th>Reference for the kinematic analysis</th></tr></thead><tbody><tr><td>HERACLES</td><td>4</td><td>0</td><td>CO(2–1)</td><td>Bacchini et al. (2020a)</td></tr><tr><td>PHANGS</td><td>10</td><td>0</td><td>CO(2–1)</td><td>This work</td></tr><tr><td>DYNAMO</td><td>9</td><td>0.07–0.13</td><td>CO(3–2), CO(4–3)</td><td>Girard et al. (2021)</td></tr><tr><td>ALPAKA</td><td>17</td><td>0.5–3.5</td><td>CO and [CI] lines</td><td>Rizzo et al. (2023)</td></tr><tr><td>Girard+19</td><td>2</td><td>1.03, 1.04</td><td>CO(4–3)</td><td>Girard et al. (2019)</td></tr><tr><td>Lelli+23</td><td>2</td><td>1.46, 2.23</td><td>CO(4–3)</td><td>Lelli et al. (2023)</td></tr><tr><td>Rizzo+20, 21</td><td>6</td><td>4–5</td><td>[CII]</td><td>Rizzo et al. (2021)</td></tr><tr><td>Fraternali+21</td><td>2</td><td>4.54, 4.57</td><td>[CII]</td><td>Fraternali et al. (2021)</td></tr><tr><td>Lelli+21</td><td>1</td><td>4.76</td><td>[CII]</td><td>Lelli et al. (2021)</td></tr><tr><td>Roman-Oliveira+23</td><td>4</td><td>4–5</td><td>[CII]</td><td>Roman-Oliveira et al. (2024)</td></tr></tbody></table>

capabilities, long integration time (i.e., ≥20 hours) are typically needed to spatially resolve the CO emission in normal main-sequence galaxies at $z \sim 2$. Consequently, despite a decade of ALMA operations, there have been no high-resolution systematic surveys of CO at $z = 0.5–4$, but only observations of single targets.

In response to this limitation, we initiated the “ALMA Archival Large Program to Advance Kinematic Analysis” (ALMA-ALPAKA) project (Rizzo et al. 2023). ALPAKA is based on the mining of the ALMA archive to gather high-quality observations of CO or [CI] emission lines from non-lensed, star-forming galaxies within the redshift range of $z = 0.5–4$. The ALPAKA project encompasses ≈0.25” observations for 28 massive star-forming galaxies, representing the most extensive $z = 0.5–4$ sample, to date, with spatially-resolved cold gas kinematics.

In Rizzo et al. (2023), we used the ALMA data to infer the dynamical state of the ALPAKA galaxies, and then divide the sample in three kinematic classes: rotating disks (19/28), interacting systems (2/28), uncertain (7/28). Rotation velocity ($V$) and velocity dispersion ($\sigma$) profiles derived using 3D BAROLO (Di Teodoro & Fraternali 2015) are presented and discussed in detail in Rizzo et al. (2023).

For the analysis presented in this paper, we selected the 17 ALPAKA disks with robust (i.e., uncertainties $\le 60\%$) estimates of their stellar masses $M_\star$ and SFR (Rizzo et al. 2023). Due to the aforementioned observational limitations, the 17 ALPAKA disks are representative of the star-forming galaxy population, with stellar masses $M_\star \gtrsim 10^{10} M_\odot$. Most of the galaxies that are included in the ALPAKA sample were discovered as bright sources in the infrared or submillimeter wavelength (see Rizzo et al. 2023, for details). The 17 ALPAKA galaxies have SFRs that are 1.3 to 28 times higher than the SFR of a main-sequence galaxy (SFR$_{MS}$). We quantify this offset $\Delta_{MS} = SFR/SFR_{MS}(M_\star, z)$ using the prescriptions recently published by Popesso et al. (2023) and normalized to a Chabrier Initial Mass Function (IMF, Chabrier 2003):

$$MS(M_{\star}, z) \equiv \log \left[ \frac{SFR_{MS}(M_{\star}, z)}{M_{\odot} \text{yr}^{-1}} \right] = (-0.034 t(\text{Gyr}) + 4.722) \times m \\ - 0.1925 \times m^2 \\ + (-26.16 + 0.2 t(\text{Gyr})), \quad (1)$$

where $t$ is the cosmic time and

$$m = \left[ \log \left( \frac{M_{\star}}{M_{\odot}} \right) + 0.025 \right].$$

To trace the evolution of the dynamical properties of galaxies across cosmic time, we supplement the ALPAKA dataset with kinematic measurements obtained from CO, [CI], or [CII] spatially resolved observations from the literature. Throughout the paper, we refer to this sample as the extended sample. This includes only galaxies that are classified as rotating disks (see the definition in Rizzo et al. 2023) and with stellar masses $M_\star \gtrsim 10^{10} M_\odot$ to ensure consistency with the stellar mass range covered by the ALPAKA sample. In Table 1, we summarize the properties of the extended sample, which is described in details in Sects. 2.1 and 2.2. Figure 1 displays the distribution of $\Delta_{MS}$ versus redshift for all the galaxies in the extended sample. These galaxies span both the main-sequence and starburst regimes (i.e., $\Delta_{MS} > 4$, Rodighiero et al. 2011). In Table 2, we list the 50th, 2nd and 98th percentiles of the $\Delta_{MS}$ distributions for the galaxies included in the extended sample (see Appendix A for further details). These summary statistics allow us to quantify the spread of the $\Delta_{MS}$ values in three redshift bins: 0–0.13, 0.5–3.5, 4–5. While the three bins contain both main-sequence and starburst galaxies, the redshift interval $z = 3–5$ is skewed towards the starburst region. This bias stems from a selection effect, as at $z > 3$ high-resolution observations are available mostly for bright dusty galaxies. The $\sigma$ used in this paper refer to the values obtained as a radial average of the velocity dispersion profiles for each galaxy in the extended sample (see Appendix A). The $V/\sigma$ values are obtained as a ratio between the maximum value of the rotation velocity and the average velocity dispersion.

## 2.1. Low-$z$ subsample

In the low-$z$ subsamples, we select 23 galaxies with kinematic data obtained from CO observations and $\Delta_{MS} \gtrsim 0.3$, matching their properties with those of galaxies at $z \gtrsim 0.5$. Throughout the paper, we will denote the galaxies in this subsample by the name of the surveys to which they belong, as outlined below.

**HERACLES:** Bacchini et al. (2020a) analyzed the CO(2–1) kinematics of 7 disk galaxies from the Heterodyne Receiver Array CO Line Extragalactic Survey (HERACLES, Leroy et al. 2009) using 3D BAROLO (Di Teodoro & Fraternali 2015). In our extended sample, we include the 4 galaxies from Bacchini et al. (2020a) with stellar masses $M_\star \gtrsim 10^{10} M_\odot$.

**PHANGS:** We include ten galaxies from the Physics at High Angular resolution in Nearby Galaxies (PHANGS)-ALMA survey (Leroy et al. 2021a). PHANGS-ALMA maps the CO(2–1) emission from a sample of nearby star-forming galaxies. The velocity $V$ and velocity dispersion $\sigma$ values used in this paper are derived after applying 3D BAROLO
<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Location of the extended sample in the plane of main-sequence offset $\Delta_{MS}$ vs. redshift. The gray dotted line shows the location of main-sequence galaxies with $\Delta_{MS} = 1$ for reference (Popesso et al. 2023). The gray solid line shows the separation between main-sequence and starburst populations at $\Delta_{MS} = 4$.</figcaption></figure>

(Di Teodoro & Fraternali 2015) to the publicly available data cubes (Leroy et al. 2021a,b, see Appendix A for details).

**DYNAMO:** The Dynamics of Newly Assembled Massive Objects (DYNAMO) survey (Green et al. 2014) includes massive ($M_{\star} > 10^{10} M_{\odot}$) disk galaxies at $z \approx 0.1$ with properties akin to main-sequence galaxies at $z = 1.5-2$. DYNAMO galaxies are starbursts, exhibiting high SFR ($\approx 10 M_{\odot}/\text{yr}$) and high gas fractions (20–40% Fisher et al. 2014, 2017). In our sample, we include the 9 DYNAMO galaxies with gas kinematics derived from CO(3–2) and CO(4–3) emission lines (Girard et al. 2021).

## 2.2. High-z subsamples

For galaxies at $z \ge 0.5$, in addition to the stellar mass-based selection criterion, we imposed a criterion based on the quality of observational data. We specifically selected galaxies from the literature covered by at least 3 independent resolution elements along the major axis of the gas disk. This ensures the reliability of $V$ and $\sigma$ measurements and allows us to minimize contamination from mergers and beam smearing (Rizzo et al. 2022). This quality criterion is not needed for the selection of the low-$z$ samples presented in Sect. 2.1 as their major axes are typically resolved with more than ten independent resolution elements.

In addition to the ALPAKA galaxies, the sample at $0.5 < z < 4$ comprises two gravitationally lensed galaxies at $z \approx 1$, whose kinematic analysis from CO(4–3) is detailed in Girard et al. (2019), and two main-sequence galaxies at $z = 1.46$ and 2.23, with CO(4–3) kinematics presented by Lelli et al. (2023).

At $z > 4$, we identified 6 lensed (Rizzo et al. 2020, 2021) and 7 non-lensed disk galaxies (Fraternali et al. 2021; Lelli et al. 2021; Roman-Oliveira et al. 2024) with spatially-resolved [CII] observations. To the best of our knowledge, there is currently no availability of resolved (i.e., at least 3 resolution elements along the major axis) [CII] observations for disk galaxies at $z \ge 5$.

## 2.3. Considerations on the extended sample

The studies of the dynamical properties of galaxies and the impact of different astrophysical mechanisms (e.g., environment, stellar mass, SFR) on the ISM properties requires the analysis of a sample spanning a large range of physical properties. How-

<table><caption>**Table 2.** Percentiles of the $\Delta_{MS}$ distributions for the galaxies in the extended sample, divided in three redshift bins.</caption><thead><tr><th>Subsample</th><th>Redshift range</th><th>50th</th><th>2nd</th><th>98th</th></tr></thead><tbody><tr><td>Low-z</td><td>0.5–3.5</td><td>3.6</td><td>0.4</td><td>15</td></tr><tr><td>Intermediate-z</td><td>0.5–3.5</td><td>2.8</td><td>0.5</td><td>20</td></tr><tr><td>High-z</td><td>4–5</td><td>9</td><td>1</td><td>46</td></tr></tbody></table>

ever, the current limitations in sample size and statistics hinder such investigations at present. Only with the next generation of facilities (e.g., ALMA upgrade, ngVLA), these studies can be conducted with greater precision (Carpenter et al. 2020; Kadler et al. 2023), allowing for a more comprehensive understanding of the evolution of galaxy dynamics across cosmic time. Nevertheless, leveraging the sample compiled in this work offers a promising avenue to initiate these studies. Despite comprising galaxies selected through various methods (e.g., sub-mm, infrared, or optical observations), the extended sample covers a narrow range in stellar mass ($M_{\star} \approx 10^{10}-10^{11} M_{\odot}$) but 4 orders of magnitude in SFR. These characteristics allow us to investigate the evolution of the dynamics of the massive galaxy population. Furthermore, the extended sample allows us to minimize observational biases (e.g., beam-smearing effect, inclusion of non-resolved interacting systems) as it includes only observations with high data quality.

# 3. Redshift evolution of galaxy dynamics using cold gas tracers

In this section, we use the extended sample to study how $\sigma$ and $V/\sigma$ evolve as a function of redshift (Sect. 3.1). We then compare the dynamical properties of our extended sample with previous results from the literature (Sect. 3.2).

## 3.1. Empirical relations using the extended sample

In Fig. 2, we present the distribution of $\sigma$ (left panels) and $V/\sigma$ (right panel) with respect to redshift for our extended sample, represented by colored markers. Despite the scatter, discernible trends emerge in both panels. Specifically, we observe an increase (decrease) by a factor of $\approx 2-3$ in $\sigma$ ($V/\sigma$) values in the redshift range $z = 0-1$. Interestingly, these trends become less pronounced at $z > 1$, where a pleatau in both $\sigma$ and $V/\sigma$ is visible.

To quantify these observed trends, we employ the following functional forms:

$$ \sigma = a_{\sigma} + b_{\sigma} \exp(-z_{0,y}/z) + s_{\sigma}, \quad (2) $$

$$ V/\sigma = a_{V/\sigma} + b_{V/\sigma} \exp(-z/z_{0,y}) + s_{V/\sigma}. \quad (3) $$

In both equations, the parameters $s_{\sigma}$ and $s_{V/\sigma}$ represent the intrinsic scatter around the median scaling relation. The fits are performed adopting a hierarchical Bayesian method and sampling the parameter space by means of an Hamiltonian Monte Carlo algorithm (see Appendix B for details). Table 3 reports the best-fit parameters, while the pink solid lines in Fig. 2 show the best fit relation with the two dark and shaded regions representing the scatters due to the uncertainties on the parameters and the intrinsic scatters parametrized by $s_{\sigma}$ and $s_{V/\sigma}$. The values of $\sigma$ at $z = 0, 2$, and 4 derived from the best-fit relations are 9, 30, 37 km/s, while the corresponding values of $V/\sigma$ are 18, 9, and 8.
<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Redshift distribution of velocity dispersion ($\sigma$, left panel) and rotational support ($V/\sigma$, right panel) for our extended sample (colored markers) consisting of massive galaxies. The $\sigma$ and $V/\sigma$ values are derived using emission lines tracing cold gas (i.e., CO, [CI], [CII]). The pink solid lines show the best-fit empirical relations, with the $1-\sigma$ uncertainties (dark pink area) and the intrinsic scatter (pink area).</figcaption></figure>

<table><caption>**Table 3.** Best-fit parameters defining the relations between $\sigma-z$, Eq. (2) and $V/\sigma-z$, Eq. (3).</caption><thead><tr><th>Parameter</th><th>y = σ</th><th>y = V/σ</th></tr></thead><tbody><tr><td>a<sub>y</sub></td><td>9 ± 1 km/s</td><td>7 ± 2</td></tr><tr><td>b<sub>y</sub></td><td>39 ± 5 km/s</td><td>13 ± 2</td></tr><tr><td>z<sub>0,y</sub></td><td>1.1 ± 0.3</td><td>1.5 ± 0.9</td></tr><tr><td>s<sub>y</sub></td><td>3.3 ± 0.9 km/s</td><td>4.5 ± 0.6</td></tr></tbody></table>

The high values of $V/\sigma$ up to $z \sim 5$ indicate that massive galaxy disks are dynamically cold across cosmic time, despite having turbulent motions higher than their local counterparts.

Our results corroborate and strengthen the conclusions of previous works based on samples covering smaller redshift ranges (e.g., Neeleman et al. 2020; Rizzo et al. 2020; Fraternali et al. 2021; Lelli et al. 2021, 2023; Tsukui & Iguchi 2021; Roman-Oliveira et al. 2024): star-forming massive disk at $z = 0-5$ are well settled and rotation-dominated, rather than being highly perturbed by mergers or intense gas accretion. It is noteworthy that despite these results being obtained from a sample comprising both main-sequence and starburst galaxies, this conclusion remains conservative, as turbulent motions would be even lower if only main-sequence galaxies were considered (see Sect. 5 for details).

The discovery of dynamically cold disks at ~4–5 (Neeleman et al. 2020; Rizzo et al. 2020, 2021; Lelli et al. 2021) prompted the exploitation of zoom-in simulations (Pallottini et al. 2022) to identify galaxies with similar dynamical properties (Kohandel et al. 2020; Kretschmer et al. 2022; Kohandel et al. 2024). For instance, Kretschmer et al. (2022) observed that these dynamically cold disks can form as transient phenomena in their simulations following the accretion of co-planar, co-rotating gas via cold cosmic-web streams. Similarly, Kohandel et al. (2024) identified dynamically cold disks among their galaxies at $z > 4$. However, owing to the intrinsic characteristics of zoom-in simulations (e.g., limited volumes, resolutions, implementation of feedback prescriptions) and the absence of galaxies with

stellar masses and SFR comparable to those in the extended sample, it remains challenging to quantify whether state-of-the-art cosmological simulations can produce dynamically cold massive disks with properties (e.g., SFR, $V/\sigma$) similar to those in our extended sample.

## 3.2. Comparison with warm gas kinematics

Ideally, one would like to compare the cold and warm gas kinematics in the same sample of galaxies. However, this comparison is hampered by the lack of warm gas observations (e.g., $H\alpha$) for most galaxies in our extended sample. Therefore, we compare the best-fit relations derived in Sect. 3.1 with trends of $\sigma-z$ and $V/\sigma-z$ obtained from warm, ionized gas observations and models available in the literature.

### 3.2.1. $H\alpha$ observations

The redshift evolution of turbulence and $V/\sigma$ from warm gas tracers has been mainly studied using compilations of main-sequence galaxies with a wide range of stellar masses ($\log(M_\star/M_\odot) = 9-11$, see Appendices C.1 and C.2). However, to investigate the impact on the kinematic tracers on the derived $\sigma$ and $V/\sigma$ and minimize potential biases due to differences in physical properties of galaxies, we make comparisons using samples that match the stellar masses and SFRs of our extended sample. In the central columns of Table 4, we summarize the main properties of the different subsamples. For comparison, we report in the right columns the number and the median $\sigma$ and $V/\sigma$ values of galaxies in our extended samples contained in the same redshift bins of the comparison $H\alpha$ subsamples.

In our comparative analysis at $z \approx 0$, we incorporate two distinct subsets. The first subset consists of the 9 DYNAMO galaxies, which are part of our extended sample and for which $H\alpha$ or $H\beta$ kinematics are presented in Girard et al. (2021). In the second subset, we included galaxies from the Sydney-AAO Multiobject Integral field Spectrograph (SAMI) survey (Croom et al. 2012; Scott et al. 2018). In particular, we selected the 67 SAMI galaxies with $M_\star \gtrsim 10^{10} M_\odot$ and $\Delta_{MS} \gtrsim 0.3$. Their $H\alpha$
<table><caption>**Table 4.** Main properties of the comparison and extended sample in four redshift bins.</caption><thead><tr><th rowspan="3">Redshift</th><th rowspan="3">Name</th><th colspan="5">Comparison samples (warm gas)</th><th colspan="3">Extended sample (cold gas)</th></tr><tr><th rowspan="2">Number</th><th colspan="3">Δ<sub>MS</sub></th><th rowspan="2">σ (km/s)</th><th rowspan="2">V/σ</th><th rowspan="2">Number</th><th rowspan="2">σ (km/s)</th><th rowspan="2">V/σ</th></tr><tr><th>50th</th><th>2nd</th><th>98th</th></tr></thead><tbody><tr><td>≈0</td><td>SAMI <sup>(a)</sup></td><td>67</td><td>1.7</td><td>0.3</td><td>6</td><td>25 ± 6</td><td>7<sup>+2</sup><sub>-1</sub></td><td>15</td><td>8 ± 2</td><td>21<sup>+8</sup><sub>-6</sub></td></tr><tr><td>0.07–0.13</td><td>DYNAMO <sup>(b)</sup></td><td>9</td><td>9</td><td>4</td><td>17</td><td>36<sup>+10</sup><sub>-9</sub></td><td>5<sup>+4</sup><sub>-2</sub></td><td>9</td><td>12<sup>+3</sup><sub>-4</sub></td><td>17<sup>+9</sup><sub>-2</sub></td></tr><tr><td>2.1–2.4</td><td>Hogan+21 <sup>(c)</sup></td><td>6</td><td>2.6</td><td>0.8</td><td>6.2</td><td>160 ± 60</td><td>2.6 ± 0.9</td><td>4</td><td>36<sup>+34</sup><sub>-11</sub></td><td>9<sup>+5</sup><sub>-3</sub></td></tr><tr><td>1.3–2.6</td><td>KAOSS <sup>(d)</sup></td><td>31</td><td>2.9</td><td>0.9</td><td>8.2</td><td>87 ± 6</td><td>2.4 ± 0.5</td><td>13</td><td>40<sup>+9</sup><sub>-11</sub></td><td>10 ± 4</td></tr></tbody></table>

**Notes.** Kinematics from: <sup>(a)</sup>Varidel et al. (2020); <sup>(b)</sup>Girard et al. (2021); <sup>(c)</sup>Hogan et al. (2021); <sup>(d)</sup>Birkin et al. (2023).

<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Comparison between cold and warm gas dynamical evolution. The gray points and pink lines show our extended sample and the best-fit relations presented in Fig. 2. The two green boxes indicate the redshift coverages and the 16th and 84th percentiles for samples with galaxy dynamics derived from $H\alpha$ (Varidel et al. 2020; Girard et al. 2021; Hogan et al. 2021; Birkin et al. 2023). The dashed and dash-dotted black lines show the expected relations for samples of main-sequence massive galaxies and for the galaxies in our extended sample, respectively. These relations are derived using the disk-instability model proposed by Wisnioski et al. (2015), which is used to reproduce the redshift evolution of $\sigma$ and $V/\sigma$ derived from warm gas tracers (see Sect. 3.2).</figcaption></figure>

kinematics is presented by Varidel et al. (2020). Girard et al. (2021) have previously studied the comparison between CO and warm gas $\sigma$ values for the 9 DYNAMO galaxies, revealing that, on average, warm gas $\sigma$ values are a factor of 3 ± 1 greater than those derived from CO. Likewise, the median CO $\sigma$ value within our extended sample at $z \approx 0$ (i.e., PHANGS and HERACLES) is 8 ± 2 km/s, that is a factor of 3 lower compared to the median $H\alpha \sigma$ of 25 ± 6 km/s obtained for the SAMI galaxies. The $V/\sigma$ values from CO are on average a factor of 3 higher with respect to the warm gas $V/\sigma$ values for the 9 DYNAMO galaxies. A similar difference of ≈3 is obtained when we compare the median CO $V/\sigma$ from our extended sample to the median $H\alpha V/\sigma$ for the SAMI galaxies.

In our comparative analysis at $z > 0$, we include the $z \sim 2$ samples studied by Hogan et al. (2021) and Birkin et al. (2023). To our knowledge, there are no systematic $H\alpha$ studies at $z > 4$. The sample analyzed by Hogan et al. (2021) consists of 6 ultra luminous infrared galaxies (ULIRG) at $z = 2.1–2.4$, while the one analyzed by Birkin et al. (2023) consists of 31 dusty star-forming galaxies (DSFGs) at $z = 1.3–2.6$. Galaxies in both studies have high stellar masses ($\log(M_{\star}/M_{\odot}) > 10$) and $\Delta_{MS}$ match-

ing the properties of the intermediate-$z$ subsample (see the percentile values in Table 4). Hogan et al. (2021) and Birkin et al. (2023) report median values of $\sigma$ of 160±60 km/s and 87±6 km/s for the ULIRG and DSFG samples, respectively. The median $V/\sigma$ values are 2.6 ± 0.9 for the ULIRGs and 2.4 ± 0.5 for the DSFGs (green shaded area in Fig. 3). We now compare these values from $H\alpha$ with the median $\sigma$ and $V/\sigma$ values of our extended sample in the two redshift bins $z = 2.1–2.4$ and $z = 1.3–2.6$ (see Table 4). The values of $\sigma$ from cold gas are a factor of 2 and 4 lower than those found by Hogan et al. (2021) and Birkin et al. (2023), while the corresponding values of $V/\sigma$ from cold gas are a factor of 3 and 4 higher than the $H\alpha$ values. This result clearly indicates that, when considering galaxies with similar stellar masses and $\Delta_{MS}$, the values of $\sigma$ and $V/\sigma$ from $H\alpha$ are, on average, systematically higher and lower, respectively, than those obtained using cold gas tracers. Providing further support to the findings in this section, we highlight that the systematic disparity of approximately a factor of 3 between cold and warm gas $\sigma$ aligns with observations from the few $z > 1$ galaxies for which both high-resolution CO (or [CII]) and $H\alpha$ data are available (Girard et al. 2021; Parlanti et al. 2024).
Due to the data quality of current Hα observations (see details in Sect. 6), we cannot determine what is causing the discrepancy between warm and cold gas kinematics. As a result, we cannot discriminate whether this difference is mainly due to: (i) differences in the disk thickness of the two gas phases; (ii) outflow contaminating the measurements of turbulence from warm gas; (iii) insufficient angular and/or spectral resolution which are biasing the measurements of warm gas σ.

### 3.2.2. Models calibrated on Hα observations

In addition to the comparison with the Hα samples, we consider the expected evolution of σ and V/σ according to the disk-instability model proposed by Wisnioski et al. (2015). This model is calibrated to reproduce the redshift trends of σ and V/σ from warm gas observations (e.g., Wisnioski et al. 2015; Turner et al. 2017; Johnson et al. 2018) and it is based on two assumptions: cold gas disks are marginally unstable to perturbations; and the cold gas σ values are similar to the warm gas ones. The presence of instabilities is commonly quantified through measurements of the Toomre parameter Q (Toomre 1964),

$$Q = \frac{\kappa\sigma}{\pi G \Sigma_{\text{gas}}}. \qquad (4)$$

In Eq. (4), Σ<sub>gas</sub> is the surface density of the gas and κ is the epicyclic frequency which is a function of the angular velocity Ω = V/R. A disk is unstable if Q ≲ 1, that is, when the self-gravity of the gas overcomes the repelling forces due to gas pressure (i.e., σ) and rotation support (Ω). In Wisnioski et al. (2015) model, the evolution of the dynamical properties of the star-forming galaxy population is explained in terms of the redshift increase of gas fraction f<sub>gas</sub> within marginally unstable disks. In particular, the evolution of σ and V/σ, defined by approximating Eq. (4) (Genzel et al. 2011), are given by the following equations:

$$\sigma(z) = \frac{Q}{\sqrt{2}} V f_{\text{gas}}(z) \qquad (5)$$

$$\frac{V}{\sigma} = \frac{\sqrt{2}}{f_{\text{gas}}(z)Q}, \qquad (6)$$

Following Wisnioski et al. (2015), we use Eqs. (5) and (6) with Q = 1, which is the typical value for marginal instability. We compute the expected evolution of σ and V/σ in the two cases:

1. massive main-sequence galaxies. We compute f<sub>gas</sub> using the empirical relation in Eq. (3) in Wisnioski et al. (2015), which is a function of redshift and stellar mass. We fix the stellar mass at 10<sup>10.5</sup> M<sub>☉</sub> and use V = 200 km/s, the average rotation velocity typically measured for massive main-sequence galaxies (Harrison et al. 2017; Förster Schreiber et al. 2018). The resulting relations are shown in Fig. 3 (bottom panels) by a black dashed line.

2. galaxies from our extended sample. We compute the expected values of σ and V/σ using Eqs. (5) and (6) (σ<sub>exp,DI</sub>, (V/σ)<sub>exp,DI</sub>) for the galaxies in our extended sample. The values of σ<sub>exp,DI</sub> and (V/σ)<sub>exp,DI</sub> represent, therefore, the values predicted by the disk-instability model for galaxies with V and f<sub>gas</sub> equal to the measured values of our extended sample. We then fitted the resulting distributions using Eqs. (2) and (3) (see Appendix C.3 for details). The dash-dotted lines in Fig. 3 are the best-fit relations to the resulting redshift distributions of σ<sub>exp,DI</sub> and (V/σ)<sub>exp,DI</sub>.

Both the dashed and dash-dotted curves in Fig. 3 systematically overestimate the values of σ from cold gas by factors of 2–4 at z ≳ 1, while simultaneously underestimating V/σ by the same factors. This discrepancy implies that the disk-instability model fails to accurately reproduce the evolution of the dynamical properties inferred from cold gas tracers. It is important to note that, although the instability-model curves are consistent with the observational results from Hα kinematics (green areas), this agreement does not imply that the warm gas disks are marginally unstable. First, the instability-model proposed by Wisnioski et al. (2015) is valid for the cold gas, under the assumption that cold and warm gas have similar kinematics. As shown in Sect. 3.2, this assumption is not valid as, on average, the σ values from cold gas appear to be a factor of 3 lower than the warm values. Second, Bacchini et al. (2024) show that the approximation of the Toomre parameter used by Wisnioski et al. (2015) to derive Eqs. (5) and (6) significantly underestimates the values of Q.

## 4. Drivers of turbulence

Theoretical models and simulations indicate that turbulence decays rapidly relative to the dynamical timescale of the galaxy (e.g., Mac Low et al. 1998; Kim & Ostriker 2015). Therefore, continuous driving mechanisms are necessary to explain the observed level of turbulence. Stellar feedback is expected to be the primary source of turbulence energy in star-forming galaxies (Mac Low & Klessen 2004; Ostriker & Shetty 2011; Hayward & Hopkins 2017; Orr et al. 2020). However, a number of additional processes (e.g., radial gas flows due to gravitational instabilities Bournaud et al. 2010; Krumholz & Burkhart 2016; Krumholz et al. 2018; Ginzburg et al. 2022), gas accretion from the cosmic web (Klessen & Hennebelle 2010; Forbes et al. 2023; Jiménez et al. 2023) are considered when the impact of stellar feedback alone is insufficient to account for the observed turbulence.

Krumholz & Burkhart (2016) and Krumholz et al. (2018, K18, hereafter) developed analytical predictions for the dependence of the cold gas velocity dispersion on SFR considering two driving mechanisms of turbulence, stellar feedback and radial gas flows (also called transport model). The transport model predicts values of σ which are more than 4 times higher than those predicted by the feedback model at fixed SFR, due to the different scaling between σ and SFR in the transport and feedback models. In the transport model, σ ∝ SFR, while σ ∝ SFR<sup>1/2</sup> for the feedback model with a star formation efficiency free to vary as a function of σ; no relationship between SFR and σ is predicted for the feedback model with fixed star formation efficiency (Appendix D). Comparisons of K18 models with warm gas σ of z = 1–3 galaxies suggest that energy injected by stellar feedback is not sufficient to explain the high values of the observed velocity dispersion. In contrast, the transport model is in good agreement with high-z observations derived from warm gas (see also Johnson et al. 2018; Übler et al. 2019). The feedback model alone reproduces the observations at low-z. The fiducial model proposed by K18 is, therefore, a combination of the feedback and transport models and it can reproduce both low and high-z galaxies.

In this context, we compare the values of cold gas σ with the fiducial and feedback models proposed by K18, in order to understand what is the potential origin of the cold gas turbulence in our extended sample. In Fig. 4, we show the distribution of the extended sample in the σ–SFR plane and the feedback and transport+feedback models. In both models, we assume the fiducial
<figure><img src="image_5.png" /><figcaption>Fig. 4. Distribution of σ vs. SFR for our extended sample divided in three redshift bins, as indicated in the legend, compared with the feedback model (left panel) and the transport+feedback model (right panel) from K18. The gray dashed and colored solid lines are obtained using the prescriptions for low and high-z galaxies, as reported in K18. For the high-z models, we vary the rotation velocity and orbital times (parenthesis in the legend) in the range covered by the extended sample (shaded area). In the left panel, the solid dark blue and dashed light blue lines overlap.</figcaption></figure>

parameters (i.e., rotation velocities and orbital times) for high-z galaxies reported by K18 and we use the range of values of our extended sample (shaded area).

The σ-SFR distribution of our sample aligns with the predictions of the stellar feedback model proposed by K18. Conversely, the transport+feedback model predicts that high-z galaxies with SFR in the range of the extended sample have σ values which are, on average, 3–10 times higher than the values of σ of our galaxies. The findings presented in this section corroborate previous studies showing that the turbulence measured from cold gas kinematics can be explained by stellar feedback (Girard et al. 2021; Rizzo et al. 2021; Roman-Oliveira et al. 2024). This underscores that the warm gas kinematics is not the ideal observable to investigate the mechanisms driving turbulence.

# 5. A predictive model for the evolution of turbulence

## 5.1. Feedback-driven model

In Sect. 4, we demonstrate that the σ-SFR distribution for galaxies in our extended sample aligns with the K18 feedback-driven models. These models, however, are based on assumptions about star-formation efficiency and the Toomre parameter (see Appendix D for details). Additionally, the K18 models rely on a series of free parameters that are adjusted for different redshift ranges and specific galaxy populations. Given these complexities, in this section, we adopt the simpler feedback-driven model proposed by Fraternali et al. (2021, hereafter F21). The F21 model has fewer free parameters, enabling us to better control and understand their impact on the results. Despite its simplicity, the F21 model effectively describes galaxies at both low (Bacchini et al. 2020a) and high redshift (Rizzo et al. 2021; Roman-Oliveira et al. 2024). This approach allows us to maintain accuracy while reducing the uncertainties inherent in more complex models. F21 assumes that supernovae must provide energy at a rate, $\dot{E}_{SN}$, equal to the rate at which turbulent kinetic energy, $E_{kin}$, is dissipated, that is, $\dot{E}_{SN} = E_{kin}/\tau$. The kinetic

energy and supernova energy rate are defined as

$$E_{\text{kin}} = \frac{3}{2} M_{\text{gas}} \sigma^2 \qquad (7)$$

$$\dot{E}_{\text{SN}} = \epsilon_{\text{SN}} \text{SFR} \eta_{\text{SN}} E_{\text{SN}} \qquad (8)$$

where $\epsilon_{SN}$ is the efficiency at which the kinetic energy is transferred to the ISM, $\eta_{SN}$ is the supernova rate, $E_{SN}$ is the supernova energy equal to $10^{51}$ erg. The timescale of turbulence dissipation $\tau$ can be defined as $\tau = 2h/\sigma$ (Mac Low et al. 1998) where $h$ is the scale height of the gas disk. For a gas disk in hydrostatic equilibrium, the scale height can be written, to first-order approximation, as:

$$h(R) = \sigma(R) \left[ \frac{\partial \phi(R, 0)}{\partial Z} \right]^{-1/2} \qquad (9)$$

where $\phi(R, 0)$ is the galactic gravitational potential in the plane of the disk ($Z = 0$). Using $\dot{E}_{SN} = E_{kin}/\tau$, we derive the relation between $\sigma$ and SFR:

$$\sigma = \text{SFR}^{1/3} \left( \frac{4 \epsilon_{\text{SN}} \eta_{\text{SN}} E_{\text{SN}} h}{3 M_{\text{gas}}} \right)^{1/3} \qquad (10)$$

Establishing the relationship between the local values of the SFR and $\sigma$ using Eq. (10) is non-trivial due to the dependence of $h(R)$ not only on $\sigma(R)$ but also on the SFR via $M_{gas}$, which affects the gravitational potential (Bacchini et al. 2020b; Ostriker & Kim 2022). As a consequence the ratio $h/M_{gas}$ may depend on both $\sigma$ and SFR. To investigate this, we used the scale height profiles $h(R)$ for 17 galaxies (4 at low-z and 13 at high-z) from our extended sample. These profiles are presented and described in Bacchini et al. (2019, 2024), and they are derived using a method which assumes vertical hydrostatic equilibrium and takes into account the full galactic potential, including contributions from stars, dark matter, and gas (Bacchini et al. 2019). We averaged these profiles radially and analyzed them in relation to $\sigma$, SFR, and $\sigma \times$ SFR. Additionally, we examined the
distribution of $\epsilon_{SN}h/M_{\text{gas}}$ as a function of $\sigma$, SFR, and $\sigma \times \text{SFR}$. This analysis revealed that the dependencies of $h$ and $\epsilon_{\text{SN}}h/M_{\text{gas}}$ on $\sigma$ and SFR cancel out when considering the average values of $h(R)$ and $\sigma(R)$ along with the global SFRs and $M_{\text{gas}}$. Moreover, the factor $CC = [\epsilon_{\text{SN}}(h/\text{pc})/(M_{\text{gas}}/M_{\odot})]^{1/3}$ remains approximately constant over cosmic time. The individual parameters in CC evolve in such a way that they cancel out any redshift dependence. For low-$z$ galaxies, $CC \sim 10^{-3}$ considering: $h \sim 200\,\text{pc}$ (Bacchini et al. 2020b), $\epsilon_{\text{SN}} \sim 0.01$ when global values of $h$ and $\sigma$ are considered (see discussion in Bacchini et al. 2020a), $M_{\text{gas}} \sim 10^9 M_{\odot}$ (e.g., Leroy et al. 2021b). For high-$z$ galaxies, $CC \sim 10^{-3}$, considering: $h \sim 200\,\text{pc}$ (Bacchini et al. 2024), $\epsilon_{\text{SN}} \sim 0.01-0.1$ (Rizzo et al. 2020, 2021; Roman-Oliveira et al. 2024), $M_{\text{gas}} \sim 10^{10} M_{\odot}$. Throughout the rest of the paper, we adopt CC as a constant value, leading to $\sigma \propto \text{SFR}^{1/3}$, which differs from the K18 feedback models, which predict that $\sigma \propto \text{SFR}^{1/2}$ or $\sigma$ remains constant.

## 5.2. Calibrating the model

In this section, we calibrate the physically-motivated model described in Sect. 5.1 using our extended sample. The final aim is to provide a relationship that predicts the average turbulence, derived from cold gas tracers, for galaxies with known global SFRs. To quantify the relationship between $\sigma$ and SFR, we fit a linear model between $\log(\sigma)$ and $\log(\text{SFR})$ (see Appendix B for details, Fig. 5), and we fixed its slope at 0.33, the value expected in the case of F21 feedback-driven turbulence. The analysis yields the following median relation:

$$ \log \sigma = 0.33 \times \log \text{SFR} + (0.77 \pm 0.03), \quad (11) $$

with an intrinsic scatter of 0.13 dex. In Eq. (11), $\sigma$ is in units of km s$^{-1}$ and SFR is in units of $M_{\odot}$ yr$^{-1}$. We note that, if the slope is left free to vary, we find a value of $0.27 \pm 0.03$, consistent with the fixed-slope case within the 2-$\sigma$ uncertainties and $\ge 8-\sigma$ larger than the slope of $\ge 0.5$ predicted by K18 models (see Appendix D for details). Furthermore, the leave-one-out cross-validation and the widely applicable information criterion methods (Vehtari et al. 2017) indicate that the model with the fixed slope is favored over the one with the free value. This is another confirmation of the supernova-driven nature of the turbulence. Based on Eq. (10), the relative small scatter in this relation confirms that the factor $CC = [\epsilon_{\text{SN}}(h/\text{pc})/(M_{\text{gas}}/M_{\odot})]^{1/3}$ has small variations that do not correlate with the galaxy SFR and the galaxy mass. In Appendix E, we also perform a diagnostic check to ensure that the residuals of the linear relation (Eq. 11) do not correlate with $\log \sigma$ and $\log \text{SFR}$. Additionally, we estimate the partial correlation between these two quantities to isolate their relationship against residuals.

## 5.3. Turbulence across cosmic time

In this section, we present a semi-empirical model predicting the redshift evolution of turbulence when considering populations of main-sequence and starburst galaxies, at fixed $M_{\star}$ and $\Delta_{\text{MS}}$. The main assumption of this model is that the evolution of $\sigma$ is mainly driven by the injection of energy due to the supernova explosions (see Sect. 4) and that, at fixed stellar mass, high-$z$ galaxies have higher SFRs than the local ones. By assuming the prescription of the redshift evolution of the SFR proposed by Popesso et al. (2023, Eq. (1)), and combining it with Eq. (11), we derive the following evolution of $\sigma$ for galaxies at $z = 0-7$:

$$ \log \sigma(M_{\star}, z, \Delta_{\text{MS}}) = 0.33 \times [\text{MS}(M_{\star}, z) + \log \Delta_{\text{MS}}] + 0.77. \quad (12) $$

<figure><img src="image_6.png" /><figcaption>Fig. 5. Distribution of $\sigma$ vs. SFR for our extended sample divided in three redshift bins, as indicated in the legend. The black line shows the best-fit relation with the 1-$\sigma$ uncertainties (dark gray area) and the intrinsic scatter (gray area).</figcaption></figure>

It is important to emphasize that this evolution of $\sigma$ as a function of redshift is entirely due to the evolution of the main-sequence relation. The physics of turbulence driving does not appear to evolve with time (see Fig. 5).

Fig. 6 shows the redshift evolution for populations of main-sequence galaxies with $M_{\star}$ of $10^{11} M_{\odot}$ and $10^{10} M_{\odot}$ in the left panel, while the right panel showcases the evolution for populations of starburst galaxies with SFR ten times higher than that of a main-sequence galaxy. We highlight the $\sigma-z$ relations in Fig. 6 are not evolutionary tracks for individual galaxies, but profiles showing the evolution of $\sigma$ when galaxy populations with similar properties (i.e., $M_{\star}$, $\Delta_{\text{MS}}$) are investigated. Interestingly, the supernova-driven model predicts that, at fixed redshift, main-sequence galaxies with lower stellar masses have lower velocity dispersions compared to galaxies with higher stellar masses, due to their lower SFRs. We also show the evolution of $\sigma$ for galaxies with $M_{\star} = 10^9 M_{\odot}$ with the dotted line. This model should be taken with caution as it relies on the extrapolation of the $\sigma$-SFR relation below the mass range probed by the extended sample. The values of $\sigma$ for the main-sequence galaxies at $z \ge 4$ are consistent with the [CII]-derived values of $\sigma$ reported for galaxies from the zoom-in SERRA simulation (Kohandel et al. 2024).

The models shown in Fig. 6 suggest only a mild redshift evolution of turbulence motion within the range $z = 0-2$ and no significant evolution at $z > 2$. For instance, a main-sequence galaxy with $M_{\star} = 10^{10} M_{\odot}$ has $\sigma \approx 5$ km/s at $z = 0$ to $\sigma \approx 16, 20$, and $21$ km/s at $z = 2, 4$, and 6. On the other hand, for a starburst galaxy with $\Delta_{\text{MS}} = 10$, the model predicts $\sigma$ values of $\approx 10, 35, 46$ km/s at $z = 0, 2$ and 6. For galaxies with smooth build-up history that grow in mass following the main-sequence relation (e.g., the Milky Way; van Dokkum et al. 2013), our model predicts that there is only a moderate evolution of turbulence over a large fraction of cosmic time. For instance, the Milky Way is expected to grow from $\approx 10^{10} M_{\odot}$ at $z = 2$ to a stellar mass of
<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Predicted evolution of $\sigma$ as a function of redshift for populations of main-sequence ($\Delta_{\text{MS}}$ = 1, left panel) and starburst ($\Delta_{\text{MS}}$ = 10, right panel) galaxies. Different lines show the evolution at different stellar masses, as indicated in the legend.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Fig. 7.** Evolution of the $V/\sigma$ ratio for populations of main-sequence galaxies. Note that at $z > 0$, these profiles should be considered as lower limits.</figcaption></figure>

$\approx 5 \times 10^{10}$ at $z = 0$ (Bland-Hawthorn & Gerhard 2016). During this time range, the average $\sigma$ values are expected to evolve from 16 to 6 km/s if the galaxy evolves following the main-sequence relation. Remarkably, the value of 6 km/s predicted by our model for a Milky-Way like galaxy at $z = 0$ is consistent with the average values of the velocity dispersion measured from molecular ($\approx 4$ km/s) and atomic gas ($\approx 8$ km/s) (Kramer & Randall 2016; Marasco et al. 2017).

Predicting the evolution of the $V/\sigma$ across cosmic time is not straightforward, as it requires determining the evolution of the rotation velocity. The latter depends on the dynamical mass and the size of the galaxy. However, we can set lower limits for the evolution of $V/\sigma$ at fixed stellar mass by employing the local stellar Tully-Fisher relation, a scaling relation between the velocity and the stellar mass of a galaxy (Tully & Fisher 1977; Reyes et al. 2011). The evolution of the stellar Tully-Fisher relation (TFR) has not been established within the redshift range of our extended sample; nonetheless, we can use the local stellar TFR to determine the lower limits of the velocity of galaxies at

fixed stellar masses. This is because, within the $\Lambda$CMD cosmology framework, the stellar TFR is expected to evolve such that at fixed stellar mass, high-$z$ galaxies have higher velocity than $z = 0$ galaxies (Dutton et al. 2011). In Fig. 7, we show the profiles obtained by dividing the velocity from the local stellar TFR from Reyes et al. (2011) with respect to $\sigma$ from the supernova-driven model. This combination of the TFR and the supernova-driven model predicts that at $z > 1$, galaxies are rotationally supported with $V/\sigma$ of $\approx 10$, even for the low-mass ($M_{\star} = 10^{9} M_{\odot}$) galaxy population.

As a consistency check, in Fig. 8, we compare the distribution of the main-sequence (upper panels) and starburst (bottom panels) galaxies in the extended sample with the predictions from our model, by dividing the galaxies in three stellar mass bins. The model is in good agreement with observations, but the typical uncertainties on the $\sigma$ values and the low number of galaxies in each bin of stellar mass and $\Delta_{\text{MS}}$ inhibit a detailed comparison. In the future, high-data quality observations (e.g., high S/N, angular resolutions) and larger sample sizes will facilitate the validation and refinement of this supernova-driven evolution model. For reference, in each panel of Fig. 8, we show the empirical relation of $\sigma-z$ presented in Sect. 3.1 with a gray dotted line. Its comparison with the supernova-driven models highlights how this relation reflects the heterogeneity of the extended sample, which consists of galaxies spanning a wide range of $\Delta_{\text{MS}}$.

# 6. Caveats

The results presented in this work rely on three observables: cold gas $\sigma$, warm gas $\sigma$, and SFR. In this section, we address potential caveats that may influence their measurements.

## 6.1. Uncertainties on the velocity dispersion

As discussed in Sect. 2, deriving the kinematics of high-$z$ galaxies poses significant challenges. A primary limitation arises from angular resolution, often implying that the emission of high-$z$ galaxies is resolved by only few independent resolution elements. Consequently, the kinematics of high-$z$ galaxies can be influenced by the beam-smearing effect, leading to a degen-
<figure><img src="image_9.png" /><figcaption>**Fig. 8.** Comparison between the distribution of the extended sample in the σ–z plane and the model predictions presented in Section 5. The upper and bottom panels show the main-sequence and starbursts galaxies, respectively. From left to right, the panels show observations and model curves in three stellar mass bins, as indicated in the legend. The data and the curves are color-coded according to the ΔMS. For comparison, we show the best-fit empirical relation (see left panel in Fig. 2) with a gray dotted line.</figcaption></figure>

eracy between *V* and σ and artificially inflated values of σ (Di Teodoro & Fraternali 2015; Rizzo et al. 2022). For example, for about half of the extended sample at *z* > 0.5, the major axis is sampled by 3–4 independent resolution elements. While the kinematic analysis for all galaxies of our extended sample was performed using state-of-the-art forward-modeling techniques to mitigate the beam-smearing effect, overcoming this issue in the case of low-resolution data remains difficult. Consequently, velocity dispersion values may be slightly overestimated. Moreover, the velocity dispersion at *z* ≳ 1 are ≈30 km/s, a value similar to the typical spectral resolution of the ALMA data. Thus, we cannot exclude that certain σ values within our high-z extended subsample may be slightly overestimated. This makes our conclusion even more robust. In fact, if the σ values from cold gas were smaller, the discrepancy with the σ values of warm gas would be even stronger.

The Hα measurements used in Sect. 3.2 are potentially affected by more severe issues compared to the CO and [CII] data. Both in the case of low- and high-z galaxies, the spectral resolutions of the instruments used for Hα are worse than the ones used for cold gas observations (see also Sect. 1). For instance, the spectral resolution (FWHM) of the SAMI data is ≈70 km/s, while the one for the KMOS data used by Birkin et al. (2023) is ≈150 km/s. In addition, the IFU data at z ∼ 2 are taken in seeing-limited mode (resolution of ≈0.6", a factor of 3 worse than the ALMA data), resulting in Hα emission which are only marginally resolved in some cases (Hogan et al. 2021; Birkin et al. 2023). As a result of these observational limitations, the warm gas σ may be overestimated by some factor. This latter is not straightforward to quantify, as it would require the modelling of several mock datacubes at different spatial and spectral resolutions and for galaxies with different properties (e.g., size, inclination), a task beyond the scope of this work.

## 6.2. Different cold gas tracers

In this study, we investigate the kinematics of cold gas ($T < 10^4$ K) using a range of emission lines, including multiple CO transitions, [CI], and [CII]. While we expect minimal differences in the kinematics among various CO transitions and [CI] lines (as discussed in Rizzo et al. 2023), the [CII] emission line may exhibit variations compared to the kinematics traced by CO and [CI] lines. Having a lower ionization potential than HI (11.3 eV instead of 13.6 eV), [CII] can be generated in multiple phases of the ISM, including molecular, atomic, and ionized phases (e.g., Rigopoulou et al. 2014; Díaz-Santos et al. 2017; Wolfire et al. 2022), although several studies suggest that the majority (60–80%) of [CII] emission arises from the molecular gas phase (e.g., Pineda et al. 2013; Vallini et al. 2015; De Breuck et al. 2019). As a rough estimate, if [CII] originates from warmer gas ($T \approx 10^2-10^3$ K) than that traced by CO ($T \lesssim 10^2$ K), we estimate a thermal broadening speed of ≈2 km/s. Consequently, similar to CO at $z < 4$, the dominant contribution to the [CII] broadening (≈30–40 km/s) at $z = 4–5$ stems from gas turbulence. Conversely, if [CII] originates from the diffuse atomic and ionized phases, we expect [CII] velocity dispersion values that are, on average, larger than CO values. The absence of kinematic studies on both CO and [CII] for the same galaxies, both locally and at high-z, precludes us from providing quantitative estimates at present. However, future observations can potentially complement ALMA observations of [CII] at $z > 3$ with CO or [CI] observations using the high-frequency bands of the Square Kilometre Array (SKA2, Braun et al. 2019) or the next generation Very Large Array (ngVLA, Kadler et al. 2023). Considering that the average [CII] σ at $z > 4$ is only about 1.2 times larger than those measured at $z \sim 2$, any disparities between CO and [CII] kinematics would strengthen our
ADS/JAO.ALMA#2017.1.01471.S; ADS/JAO.ALMA#2015.1.00456.S;
ADS/JAO.ALMA#2011.0.00124.S; ADS/JAO.ALMA#2012.1.00978.S;
ADS/JAO.ALMA#2015.1.01564.S. ALMA is a partnership of ESO (repre-
senting its member states), NSF (USA) and NINS (Japan), together with NRC
(Canada), NSC and ASIAA (Taiwan), and KASI (Republic of Korea), in coop-
eration with the Republic of Chile. The Joint ALMA Observatory is operated by
ESO, AUI/NRAO and NAOJ. We acknowledge usage of the Python program-
ming language (Van Rossum & Drake 2009), Astropy (Astropy Collaboration
2013), Matplotlib (Hunter 2007), NumPy (van der Walt et al. 2011), and SciPy
(Virtanen et al. 2020).

References

Astropy Collaboration (Robitaille, T. P., et al.) 2013, A&A, 558, A33
Bacchini, C., Fraternali, F., Iorio, G., & Pezzulli, G. 2019, A&A, 622, A64
Bacchini, C., Fraternali, F., Iorio, G., et al. 2020a, A&A, 641, A70
Bacchini, C., Fraternali, F., Pezzulli, G., & Marasco, A. 2020b, A&A, 644, A125
Bacchini, C., Nipoti, C., Iorio, G., et al. 2024, A&A, 687, A115
Birkin, J. E., Smail, I., Swinbank, A. M., et al. 2023, ArXiv e-prints
[arXiv:2301.05720]
Bland-Hawthorn, J., & Gerhard, O. 2016, ARA&A, 54, 529
Bournaud, F., Elmegreen, B. G., Teyssier, R., Block, D. L., & Puerari, I. 2010,
MNRAS, 409, 1088
Bournaud, F., Perret, V., Renaud, F., et al. 2014, ApJ, 780, 57
Braun, R., Bonaldi, A., Bourke, T., Keane, E., & Wagg, J. 2019, ArXiv e-prints
[arXiv:1912.12699]
Carilli, C. L., & Walter, F. 2013, ARA&A, 51, 105
Carpenter, J., Iono, D., Kemper, F., & Wootten, A. 2020, ArXiv e-prints
[arXiv:2001.11076]
Chabrier, G. 2003, PASP, 115, 763
Cioffi, D. F., McKee, C. F., & Bertschinger, E. 1988, ApJ, 334, 252
Croom, S. M., Lawrence, J. S., Bland-Hawthorn, J., et al. 2012, MNRAS, 421,
872
Curran, P. A. 2014, ArXiv e-prints [arXiv:1411.3816]
De Breuck, C., Weiß, A., Béthermin, M., et al. 2019, A&A, 631, A167
Dekel, A., Zolotov, A., Tweed, D., et al. 2013, MNRAS, 435, 999
Díaz-Santos, T., Armus, L., Charmandaris, V., et al. 2017, ApJ, 846, 32
Di Mascolo, L. 2024, https://doi.org/10.5281/zenodo.10528327
Di Teodoro, E. M., & Fraternali, F. 2015, MNRAS, 451, 3021
Di Teodoro, E. M., Fraternali, F., & Miller, S. H. 2016, A&A, 594, A77
Dutton, A. A., van den Bosch, F. C., Faber, S. M., et al. 2011, MNRAS, 410,
1660
Elmegreen, B. G., & Scalo, J. 2004, ARA&A, 42, 211
Epinat, B., Amram, P., & Marcelin, M. 2008, MNRAS, 390, 466
Epinat, B., Contini, T., Le Fèvre, O., et al. 2009, A&A, 504, 789
Fisher, D. B., Glazebrook, K., Bolatto, A., et al. 2014, ApJ, 790, L30
Fisher, D. B., Glazebrook, K., Damjanov, I., et al. 2017, MNRAS, 464, 491
Forbes, J. C., Emami, R., Somerville, R. S., et al. 2023, ApJ, 948, 107
Förster Schreiber, N. M., & Wuyts, S. 2020, ARA&A, 58, 661
Förster Schreiber, N. M., Genzel, R., Bouché, N., et al. 2009, ApJ, 706, 1364
Förster Schreiber, N. M., Renzini, A., Mancini, C., et al. 2018, ApJS, 238, 21
Fraternali, F., Karim, A., Magnelli, B., et al. 2021, A&A, 647, A194
Genzel, R., Newman, S., Jones, T., et al. 2011, ApJ, 733, 101
Ginzburg, O., Dekel, A., Mandelker, N., & Krumholz, M. R. 2022, MNRAS,
513, 6177
Girard, M., Dessauges-Zavadsky, M., Combes, F., et al. 2019, A&A, 631, A91
Girard, M., Fisher, D. B., Bolatto, A. D., et al. 2021, ApJ, 909, 12
Green, A. W., Glazebrook, K., McGregor, P. J., et al. 2014, MNRAS, 437, 1070
Harrison, C. M., Johnson, H. L., Swinbank, A. M., et al. 2017, MNRAS, 467,
1965
Hayward, C. C., & Hopkins, P. F. 2017, MNRAS, 465, 1682
Hoffman, M. D., & Gelman, A. 2011, ArXiv e-prints [arXiv:1111.4246]
Hogan, L., Rigopoulou, D., Magdis, G. E., et al. 2021, MNRAS, 503, 5329
Hunter, J. D. 2007, Comput. Sci. Eng., 9, 90
Ianjamasimanana, R., de Blok, W. J. G., Walter, F., et al. 2015, AJ, 150, 47
Jiménez, E., Lagos, C. D. P., Ludlow, A. D., & Wisnioski, E. 2023, MNRAS,
524, 4346
Johnson, H. L., Harrison, C. M., Swinbank, A. M., et al. 2018, MNRAS, 474,
5076
Kadler, M., Riechers, D. A., Baczko, A. K., et al. 2023, ArXiv e-prints
[arXiv:2311.10056]
Kassin, S. A., Weiner, B. J., Faber, S. M., et al. 2012, ApJ, 758, 106
Kelly, B. C. 2007, ApJ, 665, 1489
Kendall, M., & Stuart, A. 1973, The Advanced Theory of Statistics. Vol. 2:
Inference and: Relationship (Griffin)
Kim, C.-G., & Ostriker, E. C. 2015, ApJ, 815, 67
Kim, C.-G., & Ostriker, E. C. 2017, ApJ, 846, 133
Klessen, R. S., & Hennebelle, P. 2010, A&A, 520, A17
Kohandel, M., Pallottini, A., Ferrara, A., et al. 2020, MNRAS, 499, 1250
Kohandel, M., Pallottini, A., Ferrara, A., et al. 2024, A&A, 685, A72
Kramer, E. D., & Randall, L. 2016, ApJ, 829, 126
Kretschmer, M., Dekel, A., & Teyssier, R. 2022, MNRAS, 510, 3266
Krumholz, M. R., & Burkhart, B. 2016, MNRAS, 458, 1671
Krumholz, M. R., & Burkhart, B., Forbes, J. C., & Crocker, R. M. 2018, MNRAS,
477, 2716
Law, D. R., Belfiore, F., Bershady, M. A., et al. 2022, ApJ, 928, 58
Lelli, F. 2022, Nat. Astron., 6, 35
Lelli, F., Di Teodoro, E. M., Fraternali, F., et al. 2021, Science, 371, 713
Lelli, F., Zhang, Z.-Y., Bisbas, T. G., et al. 2023, A&A, 672, A106
Leroy, A. K., Walter, F., Bigiel, F., et al. 2009, AJ, 137, 4670
Leroy, A. K., Schinnerer, E., Hughes, A., et al. 2021a, ApJS, 257, 43
Leroy, A. K., Hughes, A., Liu, D., et al. 2021b, ApJS, 255, 19
Mac Low, M.-M., & Klessen, R. S. 2004, Rev. Mod. Phys., 76, 125
Mac Low, M.-M., Klessen, R. S., Burkert, A., & Smith, M. D. 1998, Phys. Rev.
Lett., 80, 2754
Marasco, A., Fraternali, F., van der Hulst, J. M., & Oosterloo, T. 2017, A&A,
607, A106
Martizzi, D., Faucher-Giguère, C.-A., & Quataert, E. 2015, MNRAS, 450, 504
Mogotsi, K. M., de Blok, W. J. G., Caldú-Primo, A., et al. 2016, AJ, 151, 15
Neal, R. 2011, Handbook of Markov Chain Monte Carlo (Chapman &
Hall/CRC), 113
Neeleman, M., Prochaska, J. X., Kanekar, N., & Rafelski, M. 2020, Nature, 581,
269
Orr, M. E., Hayward, C. C., Medling, A. M., et al. 2020, MNRAS, 496, 1620
Ostriker, E. C., & Kim, C.-G. 2022, ApJ, 936, 137
Ostriker, E. C., & Shetty, R. 2011, ApJ, 731, 41
Pacifi, C., Iyer, K. G., Mobasher, B., et al. 2023, ApJ, 944, 141
Pallottini, A., Ferrara, A., Gallerani, S., et al. 2022, MNRAS, 513, 5621
Parlanti, E., Carniani, S., Übler, H., et al. 2024, A&A, 684, A24
Pillepich, A., Nelson, D., Springel, V., et al. 2019, MNRAS, 490, 3196
Pineda, J. L., Langer, W. D., Velusamy, T., & Goldsmith, P. F. 2013, A&A, 554,
A103
Pope, A., McKinney, J., Kamieneski, P., et al. 2023, ApJ, 951, L46
Pomposo, P., Concas, A., Cresci, G., et al. 2023, MNRAS, 519, 1526
Posses, A. C., Aravena, M., González-López, J., et al. 2023, A&A, 669, A46
Puglisi, A., Dudzevičiūtė, U., Swinbank, M., et al. 2023, MNRAS, 524, 2814
Reyes, R., Mandelbaum, R., Gunn, J. E., Pizagno, J., & Lackner, C. N. 2011,
MNRAS, 417, 2347
Rigopoulou, D., Hopwood, R., Magdis, G. E., et al. 2014, ApJ, 781, L15
Rizzo, F., Vegetti, S., Powell, D., et al. 2020, Nature, 584, 201
Rizzo, F., Vegetti, S., Fraternali, F., Stacey, H. R., & Powell, D. 2021, MNRAS,
507, 3952
Rizzo, F., Kohandel, M., Pallottini, A., et al. 2022, A&A, 667, A5
Rizzo, F., Roman-Oliveira, F., Fraternali, F., et al. 2023, A&A, 679, A129
Rodighiero, G., Daddi, E., Baronchelli, I., et al. 2011, ApJ, 739, L40
Roman-Oliveira, F., Fraternali, F., & Rizzo, F. 2024, MNRAS, 521, 1045
Sawicki, M. 2012, PASP, 124, 1208
Scott, N., van de Sande, J., Croom, S. M., et al. 2018, MNRAS, 481, 2299
Sereno, M. 2016, MNRAS, 455, 2149
Simons, R. C., Kassin, S. A., Trump, J. R., et al. 2016, ApJ, 830, 14
Simons, R. C., Kassin, S. A., Weiner, B. J., et al. 2017, ApJ, 843, 46
Stott, J. P., Swinbank, A. M., Johnson, H. L., et al. 2016, MNRAS, 457, 1888
Tacconi, L. J., Neri, R., Genzel, R., et al. 2013, ApJ, 768, 74
Thornton, K., Gaudlitz, M., Janka, H. T., & Steinmetz, M. 1998, ApJ, 500, 95
Toomre, A. 1964, ApJ, 139, 1217
Tsukui, T., & Iguchi, S. 2021, Science, 372, 1201
Tully, R. B., & Fisher, J. R. 1977, A&A, 54, 661
Turner, O. J., Cirasuolo, M., Harrison, C. M., et al. 2017, MNRAS, 471, 1280
Übler, H., Genzel, R., Tacconi, L. J., et al. 2018, ApJ, 854, L24
Übler, H., Genzel, R., Wisnioski, E., et al. 2019, ApJ, 880, 48
Utomo, D., Blitz, L., & Falgarone, E. 2019, ApJ, 871, 17
Vallini, L., Gallerani, S., Ferrara, A., Pallottini, A., & Yue, B. 2015, ApJ, 813,
36
van der Walt, S., Colbert, S. C., & Varoquaux, G. 2011, Comput. Sci. Eng., 13,
22
van Dokkum, P. G., Leja, J., Nelson, E. J., et al. 2013, ApJ, 771, L35
Van Rossum, G., & Drake, F. L. 2009, Python 3 Reference Manual (Scotts
Valley: CreateSpace)
Varidel, M. R., Croom, S. M., Lewis, G. F., et al. 2020, MNRAS, 495, 2265
Vehtari, A., Melton, A., & Gabry, J. 2017, Stat. Comput., 27, 1413
Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nat. Meth., 17, 261
Whitaker, K. E., van Dokkum, P. G., Brammer, G., & Franx, M. 2012, ApJ, 754,
L29
Wisnioski, E., Förster Schreiber, N. M., Wuyts, S., et al. 2015, ApJ, 799, 209
Wolfire, M. G., Vallini, L., & Chevance, M. 2022, ARA&A, 60, 247
## Appendix A: Catalogue and table for the extended sample

With this paper we release<fnref n="1" /> a catalogue containing the physical properties for the 57 galaxies in the extended sample. The kinematic analysis of these galaxies is presented and discussed in previous papers (see Table 1), with the exception of the PHANGS galaxies. For the latter, we used the combined 12-m and 7-m publicly available data cubes<fnref n="2" /> (Leroy et al. 2021b,a) to fit their kinematics using ³DBAROLO. The parameter files used for the fitting and the corresponding outputs are available online<fnref n="3" />.

## Appendix B: Fitting empirical relations

To fit Eqs. (2), (3) and (11), we use a hierarchical Bayesian approach<fnref n="4" /> for the estimation of the likelihood function of the measured data (Kelly 2007; Sereno 2016). This hierarchical Bayesian method accounts for measurement uncertainties, intrinsic scatter, and upper limits. The sampling is performed by means of a Hamiltonian Monte Carlo (HMC; Neal 2011) inference approach, with No U-Turn Sampler (Hoffman & Gelman 2011) adaptation of the HMC step length.

### B.1. Fitting equations (2) and (3)

<table><caption>**Table B.1.** Ranges of the uniform priors used to perform the fitting of the $\sigma$-z and $V/\sigma$-z distributions using Eqs. (2) and (3).</caption><thead><tr><th>Parameter</th><th>σ</th><th>V/σ</th></tr></thead><tbody><tr><td>a</td><td>[0, 100]</td><td>[0, 100]</td></tr><tr><td>b</td><td>[0, 100]</td><td>[0, 100]</td></tr><tr><td>z<sub>0</sub></td><td>[0, 10]</td><td>[0, 10]</td></tr><tr><td>s<sub>y</sub></td><td>[0, 100]</td><td>[0, 100]</td></tr></tbody></table>

We aim at constraining the set of free parameters $\theta = \{a, b, z_0\}$ of a given functional form $f(z; \theta)$ (i.e., Eq. (2) or (3)) described by a set of independent variables $z$ and observed quantities $y$ (i.e., $\sigma$, $V/\sigma$). Here, we assume $z$ to represent a set of measurements not affected by intrinsic scatter and characterized by negligible uncertainties. This is related to the corresponding observed quantity $y$ through the relation

$$y = \eta_y(z) + \epsilon_y, \qquad (B.1)$$

where $\epsilon_y$ is the corresponding measurement uncertainty that follows a normal probability distribution with known variance $u_y^2$. The parameter $\eta_y$ instead represents the true value of the observed quantity $y$. To account for any potential scatter intrinsic to the measurement process or the measured quantities themselves, we assume $\eta$ to be drawn from the conditional distribution $p(\eta|z, \theta)$, assumed to be a normal density distribution $N\{f(z; \theta), s_y^2\}$ with mean equal to $f(z; \theta)$ and variance $s_y^2$. The

<figure><img src="image_10.png" /><figcaption>**Fig. B.1.** Posterior distributions for the parameter defining the empirical relation in Eq. (2).</figcaption></figure>

latter term specifically quantifies the intrinsic scatter of the analysed data around the mean relation $f(z; \theta)$, and we include it as free parameter in the fit. The likelihood function of each measured data point is thus given by

$$p(y|z, \theta) = \int p(y|z, \eta)p(\eta|z, \theta)d\eta. \qquad (B.2)$$

Since the data points are statistically independent, the full likelihood is the product of the likelihood function for the individual data points. For the prior probability of each model parameter, we assume wide uniform priors (see Table B.1).

### B.2. Fitting equation (11)

In this case, both the independent ($\sigma$) and dependent variable (SFR) sets comprise observed quantities. As such, we should account for any associated measurement uncertainties and potential intrinsic scatter in both cases. To do so, we extend the hierarchical approach presented in the previous section to both the dependent and independent variables $x \equiv \sigma$ and $y \equiv \text{SFR}$:

$$ \begin{aligned} x &= \eta_x(\xi) + \epsilon_x, & \eta_x &= N\{f_x(\xi; \theta_x), s_x^2\} \\ y &= \eta_y(\xi) + \epsilon_y, & \eta_y &= N\{f_y(\xi; \theta_y), s_y^2\}. \end{aligned} \qquad (B.3) $$

In accordance with the notation in Sec. B.1, the variables $\eta_i$, $\epsilon_i$, denote the true values and measurement uncertainties of each variable $i = \{x, y\}$, respectively, and the corresponding functional model $f_i$, parameters $\theta_i$ and intrinsic scatter $s_i^2$. Following Kelly (2007), Sereno (2016), we further define $\xi$ as a reference independent variable drawn from a mixture of Gaussian probability distributions and related to the true values $\eta_i$ through the normal probability distributions $N(f_i; s_i^2)$. In our specific case,

$$ \begin{aligned} f_x &= f_x(\xi) = 10^{\xi}, \\ f_y &= f_y(\xi; b) = 10^{0.33\xi + b}. \end{aligned} \qquad (B.4) $$

<footnote n="1">¹ The information on the ALPAKA sources and the extended sample is available in electronic format at https://alpaka-survey.github.io, where later releases and updates will also be placed.</footnote>

<footnote n="2">² https://www.canfar.net/storage/list/phangs</footnote>

<footnote n="3">³ https://doi.org/10.5281/zenodo.10968307</footnote>

<footnote n="4">⁴ The analysis code — scattr (Di Mascolo 2024) — is available at: https://github.com/lucadimascolo/scattr.</footnote>
<figure><img src="image_11.png" /><figcaption>Fig. B.2. Posterior distributions for the parameter defining the empirical relation in Eq. (3).</figcaption></figure>

For the prior probability on the offset parameter b in Eq. (B.4), we assume a wide uniform prior in the range [-100, 100]. The joint likelihood function can instead be rewritten as

$$p(x, y|\theta) = \iint p(x, y|\eta_x, \eta_y)p(\eta_y|\xi, \theta)p(\eta_x|\xi, \theta)p(\xi)d\eta_x d\eta_y, \quad (B.5)$$

where $p(\xi)$ is the underlying Gaussian mixture model. The term $p(x, y|\eta_x, \eta_y)$ represents the joint likelihood distribution for the {x, y} variables at given {$\eta_x, \eta_y$} and, such, could be interpreted as the probability function accounting for the uncertainties on each observed quantity. For most of the data in our sample, this can be described as a normal distribution with mean equal to the observed quantity i and variance $u_i^2$ ($i$ = {x, y}). However, for a subset of the analysed sample, we have access only to marginal constrain on the corresponding observed y values. For these data points, we thus encode the available upper limits in the probability function $p(x, y|\eta_x, \eta_y)$ following the prescription by Sawicki (2012) and considering the cumulative distribution function of a normal probability distribution integrated up to the given upper limit. For a given set of n independent data points {x, y} with m available upper limits on y, we can express the joint likelihood distribution as follows:

$$p(x, y|\eta_x, \eta_y) = \prod_{k=1}^{n} \frac{1}{\sqrt{2\pi} u_{x,k}^2} \exp \left\{ -\frac{1}{2} \left( \frac{x_k - \eta_{x,k}}{u_{x,k}} \right)^2 \right\} \\ \times \prod_{k=1}^{n-m} \frac{1}{\sqrt{2\pi} u_{y,k}^2} \exp \left\{ -\frac{1}{2} \left( \frac{y_k - \eta_{y,k}}{u_{y,k}} \right)^2 \right\} \\ \times \prod_{k=m}^{n} \frac{1}{2} \left[ 1 + \text{erf} \left( \frac{y_k - \eta_{y,k}}{\sqrt{2} u_{y,k}} \right) \right]. \quad (B.6)$$

# Appendix C: Further comparison with previous works

## C.1. Warm gas kinematics from main-sequence galaxies

In Sect. 3.2, we compare the redshift distribution of $\sigma$ from cold gas with the warm gas observations and models. Such comparison was conducted using observational samples that match the physical properties (i.e., stellar masses, SFR) of our extended sample. In this Section, we compare the cold gas $\sigma$ with the warm gas $\sigma$ derived using samples of main-sequence galaxies with a wider range of stellar masses than our extended sample.

<figure><img src="image_12.png" /><figcaption>Fig. C.1. Comparison between cold gas dynamics and previous results from the literature. The gray points and pink lines show the extended sample and the best-fit relations presented in Fig. 2. The black and gray dashed lines show the best-fit relations found by Übler et al. (2019) and Simons et al. (2017) using samples of main-sequence galaxies with stellar masses within a range larger ($\log(M_*/M_\odot) = 9-11.5$) than the ones probed by the extended sample. The relations derived by Übler et al. (2019) and Simons et al. (2017) are obtained by fitting data points in the redshift range [0, 3.8] and [0.1-2.5], respectively. The dashed curves show the extrapolations of these relations up to $z=5$. Similarly, the purple dotted line shows the best-fit relation obtained by fitting molecular gas data for 9 main-sequence galaxies at $z=0.7-2.4$ (Übler et al. 2019) and its extrapolation up to $z=5$.</figcaption></figure>

In Fig. C.1, the black dashed line represents the $\sigma$-redshift relation derived by Übler et al. (2019). This relation was established by fitting warm gas $\sigma$ measurements from various IFU studies spanning $z=0-3.8$ and comprising around 500 galaxies, while the line extension beyond $z=3.8$ is based on an extrapolation. The compilation by Übler et al. (2019) includes main-sequence galaxies with stellar masses in the range $\log(M_*/M_\odot) = 9-11.5$. The dashed gray line in Fig. C.1 represents the $\sigma$-z relation derived by Simons et al. (2017) using H$\alpha$ slit observations from the DEEP2 (Kassin et al. 2012) and SIGMA (Simons et al. 2016) surveys, which include around 500 main-sequence galaxies with $\log(M_*/M_\odot) = 9-11.5$ at $z=0.1-2.5$. Interestingly, Simons et al. (2017) find that the $\sigma$-redshift relation does not exhibit strong dependence on stellar mass, although simulations suggest that, on average, low-mass galaxies have lower velocity dispersions than high-mass galaxies (Pillepich et al. 2019; Kohandel et al. 2024). However, this discrepancy between observations and simulations might be explained by a systematic bias towards high $\sigma$ due to the beam-smearing effect (Sect. 6). This effect might be amplified in
<figure><img src="image_13.png" /><figcaption>**Fig. C.2.** Redshift distribution of $\sigma_{\text{exp,DI}}$ and $(V/\sigma)_{\text{exp,DI}}$ for the low-$z$ starbursts and the $z \ge 0.5$ galaxies in the extended sample. The black dash-dotted lines show the best-fit relations.</figcaption></figure>

low-mass galaxies, which are typically only marginally resolved (see discussion in Übler et al. 2019). Similarly to the results described in Sect. 3.2, both the relations by Übler et al. (2019) and Simons et al. (2017), derived from warm gas, systematically overestimate the values of the majority of CO/[CII] $\sigma$ of our extended sample.

### C.2. Cold gas kinematics

To date, comprehensive investigations of the evolution of galaxy dynamics, leveraging kinematic data derived from cold gas tracers, remain largely unexplored due to the dearth of spatially-resolved observations for large samples of galaxies at redshifts beyond $z > 0$. An initial endeavor in this direction was undertaken by Übler et al. (2019), who assembled a dataset comprising 9 main-sequence galaxies spanning the redshift range $z = 0.7 - 2.4$, with low-angular resolution ($\ge 0.6''$) CO observations (Tacconi et al. 2013; Übler et al. 2018). The dotted purple line in Fig. C.1 is from Übler et al. (2019) and represents the first attempt to establish an empirical relation describing the evolution of $\sigma$ with redshift using cold gas kinematics. The best-fit relation found by Übler et al. (2019) is in overall agreement with the one obtained in Sect. 3.1 in the redshift range $z = 0.7 - 2$. However, Übler et al. (2019) do not provide an explanation for the data re-processing procedure used to infer the adopted CO $\sigma$ values, most of which are systematically lower ($\sim 2\times$) than the original values reported in the first paper describing their kinematic analysis (Tacconi et al. 2013). In Fig. C.1, we show the extrapolation of the $\sigma-z$ at $z > 2.4$. As noted already by Rizzo et al. (2021), this extrapolation tends to overestimate, on average, the values of $\sigma$ measured at $z > 2$.

### C.3. Disk instability model and H$\alpha$ observations

In Sect. 3.2, we compare the evolution of galaxy dynamics from cold gas with the predictions of the disk-instability model (Wisnioski et al. 2015). To estimate the expected values of $\sigma_{\text{exp,DI}}$ and $(V/\sigma)_{\text{exp,DI}}$, we use Eq. (5) and Eq. (6) adopting our $f_{\text{gas}}$ measurements. Given that $z = 0$ galaxies are not unsta-

ble (i.e., $Q \gg 1$), we compute $\sigma_{\text{exp,DI}}$ and $(V/\sigma)_{\text{exp,DI}}$ only for the galaxies that are expected to be unstable according to the Toomre instability criterion: low-$z$ starbursts (e.g., see Girard et al. 2021, for the analysis of $Q$ for the DYNAMO galaxies) and galaxies at $z > 0$ (Wisnioski et al. 2015; Turner et al. 2017; Johnson et al. 2018). In Fig. C.2, we show the best-fit relations, obtained by fitting Eqs. (2) and (3) to the redshift distributions of $\sigma_{\text{exp,DI}}$ and $(V/\sigma)_{\text{exp,DI}}$.

## Appendix D: Assumptions of feedback models

### D.1. Krumholz et al. 2018

The feedback models in K18 are derived under the assumption that the energy injected into the ISM by supernova explosions drives the turbulence. Unlike the approach taken by F21 (see Sect. 5.1), K18 define the rate at which star formation adds energy to the gas not in terms of the energy injected by each supernova $E_{\text{SN}}$ but in terms of the mean momentum injected per unit mass of stars formed ($p_*/m_*$), which is assumed to be 3000 km/s, based on numerical simulations (Cioffi et al. 1988; Thornton et al. 1998; Martizzi et al. 2015; Kim & Ostriker 2017). The rate of energy injected by supernovae into the ISM per unit area is, therefore, defined, as

$$ \dot{E}_{\text{SN,K18}} = \Sigma_{\text{SFR}} \left( \frac{p_{\star}}{m_{\star}} \right) \sigma, \qquad (\text{D.1}) $$

where $\Sigma_{\text{SFR}}$ is the SFR per unit area. Similarly to F21, the kinetic energy dissipated per unit area is defined as

$$ \dot{E}_{\text{kin,K18}} = \frac{3}{2} \frac{\Sigma_{\text{gas}} \sigma^3}{h}. \qquad (\text{D.2}) $$

In K18, Eq. (D.2) is also expressed in terms of other model parameters, such as $Q$, $\Omega$ and the fractional contribution of gas to $Q$. A key difference between F21 and K18 is that, in K18, the rate of star formation per unit area is defined by the star formation law:

$$ \Sigma_{\text{SFR}} = \epsilon_{\text{ff}} \frac{\Sigma_{\text{gas}}}{t_{\text{ff}}}, \qquad (\text{D.3}) $$
<figure><img src="image_14.png" /><figcaption>**Fig. D.1.** Same as Fig. 5 but with the slope of the log σ – log SFR fixed at 0.5 (dashed, dark green line) and 1 (dotted green line).</figcaption></figure>

where $\epsilon_{\text{ff}}$ is the star formation efficiency and $t_{\text{ff}}$ is the free-fall time and $\Sigma_{\text{gas}}$ is expressed as a function of $\sigma$ and $Q$, Eq. (4).

In the K18 model, equations for $\epsilon_{\text{ff}}$ and $Q$ are obtained by equating Eq. (D.1) and (D.2). These two equations are then used to derive the relation between $\sigma$ and $\Sigma_{\text{SFR}}$ and its integral SFR in two scenarios. In the first scenario, the disk is marginally unstable (i.e., value of $Q$ is fixed) and the star-formation efficiency is allowed to vary. By equating equations (D.1) and (D.2), it is found that $\epsilon_{\text{ff}} \propto \sigma$, which, when substituted into Eq. (D.3), results in $\sigma \propto \Sigma_{\text{SFR}}^{1/2} \propto \text{SFR}^{1/2}$. In the second scenario, the disk is assumed to be Toomre stable, with a fixed $\epsilon_{\text{ff}}$ and variable $Q$. By equating Eqs. (D.1) and (D.2), it is determined that $\sigma$ can only take one specific value, resulting in no direct relationship between $\sigma$ and SFR $\propto \Sigma_{\text{SFR}}$.

## D.2. Energy vs momentum

In this section, we derive the relation between $\sigma$ and SFR by adopting the approach of F21, but using the momentum injected by supernovae rather than their energy. The energy injection rate is defined as

$$ \dot{E}_{\text{SN,p}} = \text{SFR} \left( \frac{p_{\star}}{m_{\star}} \right) \sigma. \qquad (\text{D.4}) $$

By equating this expression with Eq. (7), we obtain the following relationship:

$$ \sigma^2 = \text{SFR} \left( \frac{p_{\star}}{m_{\star}} \right) \frac{4h}{3 M_{\text{gas}}} \qquad (\text{D.5}) $$

The primary difference from the F21 model lies in the slope. Here, $\sigma \propto \text{SFR}^{1/2}$ when assuming a constant value for $(p_{\star}/m_{\star})$ as in

K18, or $\sigma \propto \text{SFR}$ when assuming $p_{\star} \propto \sigma$. In contrast, F21 finds that $\sigma \propto \text{SFR}^{1/3}$.

For a consistency check, we fit the relation $\log \sigma = m \times \log \text{SFR} + c$ with $m$ fixed at 0.5 (dark green dashed curve in Fig. D.1) and 1 (light green dotted curve in Fig. D.1). For comparison, we also include the fit with $m = 0.33$ (see Sect. 5 and Fig. 5. As discussed in Sect. 5, the relation with a slope of $m = 0.33$ is the only one that accurately reproduces the data.

## Appendix E: Statistical analysis of the residuals

In Sect. 5, we derived the expected relation between $\sigma$ and SFR (Eq. (10)) and discuss the potential dependence of $h$ and $M_{\text{gas}}$ on $\sigma$ and SFR. As a diagnostic check, in this section we study the correlation between the residuals of the best-fit linear relations Eq. (11) as a function of $\log \sigma$ and $\log \text{SFR}$ (Fig. E.1). Any correlations between the residuals and $\log \sigma$ and $\log \text{SFR}$ would suggest that the linear model we are assuming with a slope fixed at 0.33 might not be appropriate. To perform this test, we compute the orthogonal distance between each data point ($\log \sigma_i, \log \text{SFR}_i$) and the best-fit linear relation. To accurately propagate the uncertainties of the data points, we simulate multiple (N = 10000) datasets through bootstrapping. We then compute the distributions of the Pearson rank correlation coefficient $\rho$ (Kendall & Stuart 1973) between SFR and residuals and between $\sigma$ and the residuals, and their corresponding p-value distributions using a Monte Carlo approach (Curran 2014). This methodology allows for robustly evaluating the correlation between quantities that have associated uncertainties. The resulting distributions for $r_{\text{Pearson}}$ and the p-values are plotted in Fig. E.2. We find that the correlation coefficients are $r_{\text{Pearson}} = -0.23_{-0.31}^{+0.15}$ for the log SFR vs residual and $r_{\text{Pearson}} = 0.13_{-0.03}^{+0.23}$ for the log $\sigma$ and residual distributions, respectively. The fractions of the two distributions with p-values < 0.05 are 36% and 11%, respectively. In other words, there is only a weak and no statistically significant correlation between residuals and log SFR, and between residuals and log $\sigma$, respectively.

As an additional consistency check, we calculate the partial correlation coefficients (Kendall & Stuart 1973) between $\log \sigma$ and $\log \text{SFR}$ while controlling for spurious or extra correlations beyond the linear relation by using the residuals as control variables. This analysis allows us to determine whether the correlation between $\log \sigma$ and $\log \text{SFR}$ is statistically significant or if it can be attributed to any correlations not accounted for in the adopted model. This is easily measured as a non-null correlation with the residuals, whose scatter and distribution could depend on $h/M_{\text{gas}}$ and, therefore, on $\sigma$ and SFR. In Fig. E.3, we present the distributions of the partial correlation coefficients ($r_{\text{partial}}$) and the corresponding p-values. The analysis yields $r_{\text{partial}} = 0.86_{-0.06}^{+0.04}$. Based on the p-value distribution, where 100% of p-values are less than 0.05, we conclude that the correlation between $\log \sigma$ and $\log \text{SFR}$ is statistically significant. This finding supports the robustness of the observed relationship, indicating that it persists even when accounting for the potential influence of residuals. Consequently, our results provide strong evidence that the relationship between $\log \sigma$ and $\log \text{SFR}$ is not merely an artifact of the residual effects of the assumed model.
<figure><img src="image_15.png" /><figcaption>**Fig. E.1.** SFR and $\sigma$ as a function of the residuals from the best-fit relation defined by Eq. (11). The gray shaded regions show the intrinsic scatter of the relation.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Fig. E.2.** Distribution of the Pearson rank correlation coefficient (upper panels) and the p-value (bottom panels) between log SFR and residuals (left panels) and log $\sigma$ and residuals (right panels).</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>**Fig. E.3.** Distribution of the partial rank correlation coefficient (left panel) and the p-value (right panel) between log SFR and log σ, while excluding the residuals.</figcaption></figure>