© The Authors 2024

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 σ, 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 ≈ 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 Millimiter 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α, [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 ≈ 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., σ 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 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/σ ≲ 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 ≳ 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 σ − z relation obtained from warm gas studies at z = 0 − 3 (Übler et al. 2019) and that z ∼ 4 disks are dynamically cold, with V/σ values of ≈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 σ from cold gas (e.g., CO, HI) are, in fact, systematically lower than the ones obtained from Hα (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 σ 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_⋆ = 10⁹ − 10¹⁰ M_⊙), while the [CII] studies at z > 4 have focused on massive (M_⋆ ≳ 10¹⁰ M_⊙) upper main-sequence or starburst galaxies. Second, the typical angular resolution of ground-based telescopes used to observe z = 1 − 3 galaxies (e.g., ≳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 ≳ 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 σ ≲ FHWM/2.35 ≈ 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 instruments 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 ≈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 σ (see Lelli 2022, for details).

Given the aforementioned potential biases on the σ 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 σ and V/σ 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 σ 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 capabilities, long integration time (i.e., ≳20 hours) are typically needed to spatially resolve the CO emission in normal main-sequence galaxies at z ∼ 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 (σ) profiles derived using ^3DBAROLO (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 ≲60%) estimates of their stellar masses M_⋆ 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_⋆ ≳ 10¹⁰ M_⊙. 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 Δ_MS = SFR/SFR_MS(M_⋆, z) using the prescriptions recently published by Popesso et al. (2023) and normalized to a Chabrier Initial Mass Function (IMF, Chabrier 2003):

$$\begin{array}{cc}\hfill \mathrm{MS}(M_{\star },z)\equiv log\left[\frac{{\mathrm{SFR}}_{\mathrm{MS}}(M_{\star },z)}{M_{\odot }{\mathrm{yr}}^{-1}}\right]=& (-0.034t(\mathrm{Gyr})+4.722)\times m\hfill \\ \hfill & -0.1925\times m^2\hfill \\ \hfill & +(-26.16+0.2t(\mathrm{Gyr})),\hfill \end{array}$$(1)

where t is the cosmic time and

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

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 defintion in Rizzo et al. 2023) and with stellar masses M_⋆ ≳ 10¹⁰ M_⊙ 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 Δ_MS versus redshift for all the galaxies in the extended sample. These galaxies span both the main-sequence and starburst regimes (i.e., Δ_MS > 4, Rodighiero et al. 2011). In Table 2, we list the 50th, 2nd and 98th percentiles of the Δ_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 Δ_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 σ 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/σ values are obtained as a ratio between the maximum value of the rotation velocity and the average velocity dispersion.

Fig. 1. Location of the extended sample in the plane of main-sequence offset ΔMS vs. redshift. The gray dotted line shows the location of main-sequence galaxies with ΔMS = 1 for reference (Popesso et al. 2023). The gray solid line shows the separation between main-sequence and starburst populations at ΔMS = 4.

2.1. Low-z subsample

In the low-z subsamples, we select 23 galaxies with kinematic data obtained from CO observations and Δ_MS ≳ 0.3, matching their properties with those of galaxies at z ≳ 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 ^3DBAROLO (Di Teodoro & Fraternali 2015). In our extended sample, we include the 4 galaxies from Bacchini et al. (2020a) with stellar masses M_⋆ ≳ 10¹⁰ M_⊙.

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 σ values used in this paper are derived after applying ^3DBAROLO (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_⋆ > 10¹⁰ M_⊙) disk galaxies at z ≈ 0.1 with properties akin to main-sequence galaxies at z = 1.5 − 2. DYNAMO galaxies are starbursts, exhibiting high SFR (≈10 M_⊙/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 ≳ 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 σ 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 ≈ 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 ≳ 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. However, 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_⋆ ≈ 10¹⁰ − 10¹¹ M_⊙) 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 σ and V/σ 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 σ (left panels) and V/σ (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 ≈2 − 3 in σ (V/σ) values in the redshift range z = 0 − 1. Interestingly, these trends become less pronounced at z > 1, where a pleatau in both σ and V/σ is visible.

Fig. 2. Redshift distribution of velocity dispersion (σ, left panel) and rotational support (V/σ, right panel) for our extended sample (colored markers) consisting of massive galaxies. The σ and V/σ 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-σ uncertainties (dark pink area) and the intrinsic scatter (pink area).

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

$$\begin{array}{cc}\hfill \sigma & =a_{\sigma }+b_{\sigma }exp(-z_{0,y}/z)+s_{\sigma },\hfill \end{array}$$(2)

$$\begin{array}{cc}\hfill V/\sigma & =a_{V/\sigma }+b_{V/\sigma }exp(-z/z_{0,y})+s_{V/\sigma }.\hfill \end{array}$$(3)

In both equations, the parameters s_σ and s_V/σ 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_σ and s_V/σ. The values of σ at z = 0, 2, and 4 derived from the best-fit relations are 9, 30, 37 km/s, while the corresponding values of V/σ are 18, 9, and 8. The high values of V/σ up to z ∼ 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 presciptions) 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/σ) 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α) for most galaxies in our extended sample. Therefore, we compare the best-fit relations derived in Sect. 3.1 with trends of σ − z and V/σ − z obtained from warm, ionized gas observations and models available in the literature.

3.2.1. Hα observations

The redshift evolution of turbulence and V/σ from warm gas tracers has been mainly studied using compilations of main-sequence galaxies with a wide range of stellar masses (log(M_⋆/M_⊙) = 9 − 11, see Appendices C.1 and C.2). However, to investigate the impact on the kinematic tracers on the derived σ and V/σ 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 σ and V/σ values of galaxies in our extended samples contained in the same redshift bins of the comparison Hα subsamples.

In our comparative analysis at z ≈ 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α or Hβ 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_⋆ ≳ 10¹⁰ M_⊙ and Δ_MS ≳ 0.3. Their Hα kinematics is presented by Varidel et al. (2020). Girard et al. (2021) have previously studied the comparison between CO and warm gas σ values for the 9 DYNAMO galaxies, revealing that, on average, warm gas σ values are a factor of 3 ± 1 greater than those derived from CO. Likewise, the median CO σ value within our extended sample at z ≈ 0 (i.e., PHANGS and HERACLES) is 8 ± 2 km/s, that is a factor of 3 lower compared to the median Hασ of 25 ± 6 km/s obtained for the SAMI galaxies. The V/σ values from CO are on average a factor of 3 higher with respect to the warm gas V/σ values for the 9 DYNAMO galaxies. A similar difference of ≈3 is obtained when we compare the median CO V/σ from our extended sample to the median HαV/σ for the SAMI galaxies.

In our comparative analysis at z > 0, we include the z ∼ 2 samples studied by Hogan et al. (2021) and Birkin et al. (2023). To our knowledge, there are no systematic Hα 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_⋆/M_⊙) > 10) and Δ_MS matching 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 σ of 160 ± 60 km/s and 87 ± 6 km/s for the ULIRG and DSFG samples, respectively. The median V/σ 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α with the median σ and V/σ 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 σ 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/σ from cold gas are a factor of 3 and 4 higher than the Hα values. This result clearly indicates that, when considering galaxies with similar stellar masses and Δ_MS, the values of σ and V/σ from Hα 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 σ aligns with observations from the few z > 1 galaxies for which both high-resolution CO (or [CII]) and Hα data are available (Girard et al. 2021; Parlanti et al. 2024).

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α (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 σ and V/σ derived from warm gas tracers (see Sect. 3.2).

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),

$$\begin{array}{c}\hfill Q=\frac{\kappa \sigma }{\pi G{\mathrm{\Sigma }}_{\mathrm{gas}}}·\end{array}$$(4)

In Eq. (4), Σ_gas 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_gas 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:

$$\begin{array}{cc}\hfill \sigma (z)& =\frac{Q}{\sqrt{2}}V{f}_{\mathrm{gas}}(z)\hfill \end{array}$$(5)

$$\begin{array}{cc}\hfill \frac{V}{\sigma }& =\frac{\sqrt{2}}{f_{\mathrm{gas}}(z)Q},\hfill \end{array}$$(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_gas 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^10.5 M_⊙ 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) (σ_exp, DI, (V/σ)_exp, DI) for the galaxies in our extended sample. The values of σ_exp, DI and (V/σ)_exp, DI represent, therefore, the values predicted by the disk-instability model for galaxies with V and f_gas 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 σ_exp, DI and (V/σ)_exp, DI.

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^1/2 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 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).

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 trasport+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.

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, Ė_SN, equal to the rate at which turbulent kinetic energy, E_kin is dissipated, that is, Ė_SN = E_kin/τ. The kinetic energy and supernova energy rate are defined as

$$\begin{array}{cc}\hfill E_{\mathrm{kin}}& =\frac{3}{2}{M}_{\mathrm{gas}}{\sigma }^2\hfill \end{array}$$(7)

$$\begin{array}{cc}\hfill {\dot{E}}_{\mathrm{SN}}& ={ϵ}_{\mathrm{SN}}\mathrm{SFR}{\eta }_{\mathrm{SN}}{E}_{\mathrm{SN}}\hfill \end{array}$$(8)

where ϵ_SN is the efficiency at which the kinetic energy is transferred to the ISM, η_SN is the supernova rate, E_SN is the supernova energy equal to 10⁵¹ erg. The timescale of turbulence dissipation τ can be defined as τ = 2h/σ (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:

$$\begin{array}{c}\hfill h(R)=\sigma (R){\left[\frac{\partial \varphi (R,0)}{\partial Z}\right]}^{-1/2}\end{array}$$(9)

where ϕ(R, 0) is the galactic gravitational potential in the plane of the disk (Z = 0). Using Ė_SN = E_kin/τ, we derive the relation between σ and SFR:

$$\begin{array}{c}\hfill \sigma ={\mathrm{SFR}}^{1/3}{\left(\frac{4{ϵ}_{\mathrm{SN}}{\eta }_{\mathrm{SN}}{E}_{\mathrm{SN}}h}{3M_{\mathrm{gas}}}\right)}^{1/3}.\end{array}$$(10)

Establishing the relationship between the local values of the SFR and σ using Eq. (10) is non-trivial due to the dependence of h(R) not only on σ(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 σ 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 σ, SFR, and σ × SFR. Additionally, we examined the distribution of ϵ_SNh/M_gas as a function of σ, SFR, and σ × SFR. This analysis revealed that the dependencies of h and ϵ_SNh/M_gas on σ and SFR cancel out when considering the average values of h(R) and σ(R) along with the global SFRs and M_gas. Moreover, the factor CC = [ϵ_SN(h/pc)/(M_gas/M_⊙]^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 ∼ 10⁻³ considering: h ∼ 200 pc (Bacchini et al. 2020b), ϵ_SN ∼ 0.01 when global values of h and σ are considered (see discussion in Bacchini et al. 2020a), M_gas ∼ 10⁹ M_⊙ (e.g., Leroy et al. 2021b). For high-z galaxies, CC ∼ 10⁻³, considering: h ∼ 200 pc (Bacchini et al. 2024), ϵ_SN ∼ 0.01 − 0.1 (Rizzo et al. 2020, 2021; Roman-Oliveira et al. 2024), M_gas ∼ 10¹⁰ M_⊙. Throughout the rest of the paper, we adopt CC as a constant value, leading to σ ∝ SFR^1/3, which differs from the K18 feedback models, which predict that σ ∝ SFR^1/2 or σ 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 σ and SFR, we fit a linear model between log(σ) and log(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:

$$\begin{array}{c}\hfill log\sigma =0.33\times log\mathrm{SFR}+(0.77\pm 0.03),\end{array}$$(11)

Fig. 5. Distribution of σ 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-σ uncertainties (dark gray area) and the intrinsic scatter (gray area).

with an intrinsic scatter of 0.13 dex. In Eq. (11), σ is in units of km s⁻¹ and SFR is in units of M_⊙ yr⁻¹. We note that, if the slope is left free to vary, we find a value of 0.27 ± 0.03, consistent with the fixed-slope case within the 2-σ uncertainties and ≳8-σ larger than the slope of ≳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 = [ϵ_SN(h/pc)/(M_gas/M_⊙]^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 σ and logSFR. 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_⋆ and Δ_MS. The main assumption of this model is that the evolution of σ 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 σ for galaxies at z = 0 − 7:

$$\begin{array}{c}\hfill log\sigma (M_{\star },z,{\mathrm{\Delta }}_{\mathrm{MS}})=0.33\times [\mathrm{MS}(M_{\star },z)+log{\mathrm{\Delta }}_{\mathrm{MS}}]+0.77.\end{array}$$(12)

It is important to emphasize that this evolution of σ 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_⋆ of 10¹¹ M_⊙ and 10¹⁰ M_⊙ 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 σ − z relations in Fig. 6 are not evolutionary tracks for individual galaxies, but profiles showing the evolution of σ when galaxy populations with similar properties (i.e., M_⋆, Δ_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 σ for galaxies with M_⋆ = 10⁹ M_⊙ with the dotted line. This model should be taken with caution as it relies on the extrapolation of the σ–SFR relation below the mass range probed by the extended sample. The values of σ for the main-sequence galaxies at z ≳ 4 are consistent with the [CII]-derived values of σ reported for galaxies from the zoom-in SERRA simulation (Kohandel et al. 2024).

Fig. 6. Predicted evolution of σ as a function of redshift for populations of main-sequence (ΔMS = 1, left panel) and starburst (ΔMS = 10, right panel) galaxies. Different lines show the evolution at different stellar masses, as indicated in the legend.

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_⋆ = 10¹⁰ M_⊙ has σ ≈ 5 km/s at z = 0 to σ ≈ 16, 20, and 21 km/s at z = 2, 4, and 6. On the other hand, for a starburst galaxy with Δ_MS = 10, the model predicts σ values of ≈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 ≈10¹⁰ M_⊙ at z = 2 to a stellar mass of ≈5 × 10¹⁰ at z = 0 (Bland-Hawthorn & Gerhard 2016). During this time range, the average σ 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 (≈4 km/s) and atomic gas (≈8 km/s) (Kramer & Randall 2016; Marasco et al. 2017).

Predicting the evolution of the V/σ 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/σ 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 Λ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 σ 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/σ of ≈10, even for the low-mass (M_⋆ = 10⁹ M_⊙) galaxy population.

Fig. 7. Evolution of the V/σ ratio for populations of main-sequence galaxies. Note that at z > 0, these profiles should be considered as lower limits.

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 σ values and the low number of galaxies in each bin of stellar mass and Δ_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 σ − 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 Δ_MS.

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.

6. Caveats

The results presented in this work rely on three observables: cold gas σ, warm gas σ, 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 degeneracy 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⁴ 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 ≈ 10² − 10³ K) than that traced by CO (T ≲ 10² 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 ∼ 2, any disparities between CO and [CII] kinematics would strengthen our findings, further attenuating the redshift evolution of turbulent motion.

6.3. Considerations on SFR

The model presented in Sect. 5 relies on the σ–SFR relation. However, the SFR values of the extended sample are derived from diverse methodologies and assumptions as they are taken from the literature. We normalized the SFR values to a common IMF, but variations in assumptions and methodologies may introduce scatter (Pacifici et al. 2023).

Although radial profiles of σ are available for the majority of galaxies within our extended sample, the absence of SFR surface density maps impedes an investigation of the efficiency of the turbulence drivers as a function of the galactocentric distance, as has been done for local galaxy samples (Utomo et al. 2019; Bacchini et al. 2020a). Analyzing global σ and SFR values enables an assessment of the primary mechanisms driving turbulence. Future high-quality observations mapping the Hα or UV emission will allow us to obtain spatially-resolved maps of the SFR distribution and finally investigate the relative role of different driving mechanisms in the different regions of galaxies.

7. Summary and conclusions

In this study, we analyze the redshift evolution of galaxy dynamics (i.e., σ, V/σ) using kinematics from cold gas tracers (CO, [CI], [CII]). Our initial sample consists of 17 galaxy disks from the ALPAKA sample with high data quality. We complement ALPAKA with σ and V/σ measurements from the literature for subsamples of galaxies with spatially-resolved observations of emission lines tracing cold gas. This extended sample consists of 57 main-sequence and starburst disks at z = 0 − 5 and it covers a narrow range in stellar masses (M_⋆ ≈ 10¹⁰ − 10¹¹ M_⊙), but a wide range of 4 order of magnitude in SFRs. Our main results are the following.

• There is a trend of σ increasing by around 2 to 3 times across the redshift range z = 0 − 1, while V/σ decreases by a similar amount. However, these trends become less apparent at z > 1, where both σ and V/σ exhibit a plateau. We quantify the observed evolution using empirical relations, even though the derived trends reflect the heterogeneity of our extended sample, consisting of galaxies with a large range of SFRs.

• We compare the kinematics of cold gas from our extended sample with the ones of warm gas. For the latter, we selected galaxies from the literature with physical properties (i.e., M_⋆, SFR) matching the ones of our extended sample. The values of σ (V/σ) from Hα are a factor of ≈3 higher (lower) than the ones from cold gas both for low and high-z galaxies.

• We investigate the potential drivers of turbulence and find that stellar feedback can reproduce the observed turbulence, traced by velocity dispersion measurements obtained from cold gas tracers. In contrast, gravitational instabilities tend to overpredict turbulence measured from cold gas. This finding diverges from previous studies based on warm gas kinematics, which suggested that turbulence is mainly driven by gravitational instabilities in z > 1 galaxies.

• We present a physically-motivated model predicting that turbulence is fed by the energy injected into the ISM by supernova explosions, and that the physics of turbulence driving does not to evolve with time. We calibrated this model on observations and developed a semi-empirical framework to predict the redshift evolution of turbulence in star-forming galaxies based on the evolution of the main-sequence relation. This model suggests a mild increase in turbulence within the redshift range 0−2 and no significant evolution at z > 2 for populations of galaxies at fixed stellar mass. Within this framework, for a Milky-Way-like progenitor, the velocity dispersion from cold gas evolve from ≈16 km/s to 6 km/s in the last 10 Gyr (z ≲ 2) of cosmic history. Interestingly, the supernova-driven model indicates that galaxies with lower stellar masses have lower velocity dispersion compared to more massive counterparts. By combining this model with the local stellar Tully–Fisher relation, we set lower limits for the redshift evolution of V/σ at fixed stellar mass. We predict that high-z, low mass (M_⋆ ≳ 10⁹ M_⊙) galaxy disks should be rotationally supported, with V/σ ∼ 10.

The results presented in this paper indicate that the dynamical properties of galaxies significantly differ when derived from warm or cold gas. This work highlights the importance of studying the evolution of galaxy dynamics, the turbulence within the ISM and its origin using cold gas tracers. The difference between cold and warm gas kinematics evidenced in our analysis may be due to (i) diffuse nature of warm, ionized gas compared to the cold one; (ii) bias of the velocity dispersion measurements to higher values due to contamination from non-circular motions due to outflows; (iii) insufficient spectral and angular resolution. In the future, studies of the ionized gas kinematics using high angular and spectral resolution observations (FWHM ≲ 40 km/s, e.g., VLT/ERIS; ELT/HARMONI) will be able to discriminate between these three hypotheses.

The supernova-driven model presented in this paper makes testable predictions on the evolution of turbulence in galaxy disks across cosmic time. As more cold gas observations of star-forming galaxies with a wide range of stellar masses will be collected over the next years, we should be able to gain a more complete picture on the formation of the progenitors of disks in Milky-Way-like and local spiral galaxies.

Acknowledgments

The author thank the anonymous referee for constructive comments and suggestions. FR is grateful to Pavel Mancera Piña for useful comments and discussions. FR acknowledges support from the European Union’s Horizon 2020 research and innovation program under the Marie Sklodowska-Curie grant agreement No. 847523 ‘INTERACTIONS’. This work has been supported by the Cosmic Dawn Center (DAWN), funded by the Danish National Research Foundation under grant No. 140. CB acknowledges support from the Carlsberg Foundation Fellowship Programme by Carlsbergfondet. LDM acknowledges support from the French government, through the UCA^J.E.D.I. Investments in the Future project managed by the National Research Agency (ANR) with the reference number ANR-15-IDEX-01. FRO acknowledges support from the Dutch Research Council (NWO) through the Klein-1 Grant code OCEN2.KLEIN.088. This paper makes use of the following ALMA data: ADS/JAO.ALMA#2017.1.01659.S; ADS/JAO.ALMA#2017.1.00471.S; ADS/JAO.ALMA#2017.1.01674.S; ADS/JAO.ALMA#2015.1.00862.S; ADS/JAO.ALMA#2017.1.01228.S; ADS/JAO.ALMA#2018.1.00974.S; ADS/JAO.ALMA#2017.1.00413.S; ADS/JAO.ALMA#2019.1.01362.S; ADS/JAO.ALMA#2017.1.01045.S; ADS/JAO.ALMA#2013.1.00059.S; ADS/JAO.ALMA#2018.1.00543.S; ADS/JAO.ALMA#2015.1.00723.S; ADS/JAO.ALMA#2018.1.01146.S; ADS/JAO.ALMA#2016.1.01155.S; ADS/JAO.ALMA#2017.1.01677.S; ADS/JAO.ALMA#2018.1.01306.S; ADS/JAO.ALMA#2017.1.00908.S; ADS/JAO.ALMA#2012.1.00650.S; ADS/JAO.ALMA#2013.1.00803.S; ADS/JAO.ALMA#2013.1.01161.S; ADS/JAO.ALMA#2015.1.00121.S; ADS/JAO.ALMA#2015.1.00782.S; ADS/JAO.ALMA#2015.1.00925.S; ADS/JAO.ALMA#2015.1.00956.S; ADS/JAO.ALMA#2016.1.00386.S; ADS/JAO.ALMA#2017.1.00392.S; ADS/JAO.ALMA#2017.1.00766.S; ADS/JAO.ALMA#2017.1.00886.L; ADS/JAO.ALMA#2018.1.00484.S; ADS/JAO.ALMA#2018.1.01321.S; ADS/JAO.ALMA#2018.1.01651.S; ADS/JAO.ALMA#2018.A.00062.S; ADS/JAO.ALMA#2019.1.01235.S; ADS/JAO.ALMA#2019.2.00129.S; ADS/JAO.ALMA#2016.1.01499.S; ADS/JAO.ALMA#2015.1.00330.S; ADS/JAO.ALMA#2017.1.00127.S; ADS/JAO.ALMA#2017.1.00394.S; ADS/JAO.ALMA#2017.1.01052.S; ADS/JAO.ALMA#2018.1.00001.S; ADS/JAO.ALMA#2017.1.01020.S; ADS/JAO.ALMA#2019.1.00862.S; 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 (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), NSC and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ. We acknowledge usage of the Python programming 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

Appendix A: Catalogue and table for the extended sample

With this paper we release¹ 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² (Leroy et al. 2021b, 2021a) to fit their kinematics using ^3DBAROLO. The parameter files used for the fitting and the corresponding outputs are available online³.

Appendix B: Fitting empirical relations

To fit Eqs. (2), (3) and (11), we use a hierarchical Bayesian approach⁴ 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)

We aim at constraining the set of free parameters θ = {a, b, z₀} of a given functional form f(z; θ) (i.e., Eq. (2) or (3)) described by a set of independent variables z and observed quantities y (i.e., σ, V/σ). 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

$$\begin{array}{c}\hfill y={\eta }_y(z)+{ϵ}_y,\end{array}$$(B.1)

where ϵ_y is the corresponding measurement uncertainty that follows a normal probability distribution with known variance $u_y^2$. The parameter η_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 η to be drawn from the conditional distribution p(η|z, θ), assumed to be a normal density distribution $\mathcal{N}\{f(z;\mathit{\theta }),s_y^2\}$ with mean equal to f(z; θ) and variance $s_y^2$. The latter term specifically quantifies the intrinsic scatter of the analysed data around the mean relation f(z; θ), and we include it as free parameter in the fit. The likelihood function of each measured data point is thus given by

$$\begin{array}{c}\hfill p(y|z,\mathit{\theta })={\displaystyle \int p}(y|z,\eta )p(\eta |z,\mathit{\theta })d\eta .\end{array}$$(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).

Fig. B.1. Posterior distributions for the parameter defining the empirical relation in Eq. (2).

Fig. B.2. Posterior distributions for the parameter defining the empirical relation in Eq. (3).

B.2. Fitting equation (11)

In this case, both the independent (σ) 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 ≡ σ and y≡ SFR:

$$\begin{array}{cc}\hfill x& ={\eta }_x(\xi )+{ϵ}_x,{\eta }_x=\mathcal{N}\{f_x(\xi ;{\mathit{\theta }}_x),s_x^2\}\hfill \\ \hfill y& ={\eta }_y(\xi )+{ϵ}_y,{\eta }_y=\mathcal{N}\{f_y(\xi ;{\mathit{\theta }}_y),s_y^2\}.\hfill \end{array}$$(B.3)

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

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

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

$$\begin{array}{c}\hfill p(x,y|\mathit{\theta })={\displaystyle \int \int p}(x,y|{\eta }_x,{\eta }_y)p({\eta }_y|\xi ,\mathit{\theta })p({\eta }_x|\xi ,\mathit{\theta })p(\xi )d{\eta }_{x}d{\eta }_y,\end{array}$$(B.5)

where p(ξ) is the underlying Gaussian mixture model. The term p(x, y|η_x, η_y) represents the joint likelihood distribution for the {x, y} variables at given {η_x, η_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|η_x, η_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:

$$\begin{array}{cc}\hfill p(x,y|{\eta }_x,{\eta }_y)& ={\displaystyle \prod _{k=1}^n}\frac{1}{\sqrt{2\pi u_{x,k}^2}}exp\{-\frac{1}{2}{\left(\frac{x_k-{\eta }_{x,k}}{u_{x,k}}\right)}^2\}\hfill \\ \hfill & \times {\displaystyle \prod _{k=1}^{n-m}}\frac{1}{\sqrt{2\pi u_{y,k}^2}}exp\{-\frac{1}{2}{\left(\frac{y_k-{\eta }_{y,k}}{u_{y,k}}\right)}^2\}\hfill \\ \hfill & \times {\displaystyle \prod _{k=m}^n}\frac{1}{2}[1+\mathrm{erf}\left(\frac{y_k-{\eta }_{y,k}}{\sqrt{2}{u}_{y,k}}\right)].\hfill \end{array}$$(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 σ 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 σ with the warm gas σ derived using samples of main-sequence galaxies with a wider range of stellar masses than our extended sample.

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⊙) = 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.

Fig. C.2. Redshift distribution of σexp, DI and (V/σ)exp, DI for the low-z starbursts and the z ≳ 0.5 galaxies in the extended sample. The black dash-dotted lines show the best-fit relations.

In Fig. C.1, the black dashed line represents the σ-redshift relation derived by Übler et al. (2019). This relation was established by fitting warm gas σ 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_⊙) = 9 − 11.5. The dashed gray line in Fig. C.1 represents the σ-z relation derived by Simons et al. (2017) using Hα 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_⊙) = 9 − 11.5 at z = 0.1 − 2.5. Interestingly, Simons et al. (2017) find that the σ-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 σ due to the beam-smearing effect (Sect. 6). This effect might be amplified in 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] σ 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 (≳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 σ 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 σ values, most of which are systematically lower (∼2×) 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 σ-z at z > 2.4. As noted already by Rizzo et al. (2021), this extrapolation tends to overestimate, on average, the values of σ measured at z > 2.

C.3. Disk instability model and Hα 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 σ_exp, DI and (V/σ)_exp, DI, we use Eq. (5) and Eq. (6) adopting our f_gas measurements. Given that z = 0 galaxies are not unstable (i.e., Q ≫ 1), we compute σ_exp, DI and (V/σ)_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 σ_exp, DI and (V/σ)_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_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

$$\begin{array}{c}\hfill {\dot{E}}_{\mathrm{SN},\mathrm{K}18}={\mathrm{\Sigma }}_{\mathrm{SFR}}\left(\frac{p_{\star }}{m_{\star }}\right)\sigma ,\end{array}$$(D.1)

where Σ_SFR is the SFR per unit area. Similarly to F21, the kinetic energy dissipated per unit area is defined as

$$\begin{array}{c}\hfill {\dot{E}}_{\mathrm{kin},\mathrm{K}18}=\frac{3}{2}\frac{{\mathrm{\Sigma }}_{\mathrm{gas}}{\sigma }^3}{h}.\end{array}$$(D.2)

In K18, Eq. (D.2) is also expressed in terms of other model parameters, such as Q, Ω 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:

$$\begin{array}{c}\hfill {\mathrm{\Sigma }}_{\mathrm{SFR}}={ϵ}_{\mathrm{ff}}\frac{{\mathrm{\Sigma }}_{\mathrm{gas}}}{t_{\mathrm{ff}}},\end{array}$$(D.3)

where ϵ_ff is the star formation efficiency and t_ff is the free-fall time and Σ_gas is expressed as a function of σ and Q, Eq. (4).

In the K18 model, equations for ϵ_ff and Q are obtained by equating Eq. (D.1) and (D.2). These two equations are then used to derive the relation between σ and Σ_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 ϵ_ff ∝ σ, which, when substituted into Eq. (D.3), results in $\sigma \propto {\mathrm{\Sigma }}_{\mathrm{SFR}}^{1/2}\propto {\mathrm{SFR}}^{1/2}$. In the second scenario, the disk is assumed to be Toomre stable, with a fixed ϵ_ff and variable Q. By equating Eqs. (D.1) and (D.2), it is determined that σ can only take one specific value, resulting in no direct relationship between σ and SFR ∝ Σ_SFR.

D.2. Energy vs momentum

In this section, we derive the relation between σ 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

$$\begin{array}{c}\hfill {\dot{E}}_{\mathrm{SN},\mathrm{p}}=\mathrm{SFR}\left(\frac{p_{\star }}{m_{\star }}\right)\sigma .\end{array}$$(D.4)

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

$$\begin{array}{c}\hfill {\sigma }^2=\mathrm{SFR}\left(\frac{p_{\star }}{m_{\star }}\right)\frac{4h}{3M_{\mathrm{gas}}}\end{array}$$(D.5)

The primary difference from the F21 model lies in the slope. Here, σ ∝ SFR^1/2 when assuming a constant value for $\left(\frac{p_{\star }}{m_{\star }}\right)$ as in K18, or σ ∝ SFR when assuming p_⋆ ∝ σ. In contrast, F21 finds that σ ∝ SFR^1/3.

For a consistency check, we fit the relation log σ = m × logSFR + 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.

Fig. D.1. Same as Fig. 5 but with the slope of the log σ − logSFR fixed at 0.5 (dashed, dark green line) and 1 (dotted green line).

Appendix E: Statistical analysis of the residuals

In Sect. 5, we derived the expected relation between σ and SFR (Eq. (10) and discuss the potential dependence of h and M_gas on σ 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 σ and logSFR (Fig. E.1). Any correlations between the residuals and log σ and logSFR 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 σ_i, logSFR_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 ρ (Kendall & Stuart 1973) between SFR and residuals and between σ 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_Pearson and the p-values are plotted in Fig. E.2. We find that the correlation coefficients are $r_{\mathrm{Pearson}}=-0.{23}_{-0.31}^{+0.15}$ for the logSFR vs residual and $r_{\mathrm{Pearson}}=0.{13}_{-0.03}^{+0.23}$ for the log σ 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 logSFR, and between residuals and log σ, respectively.

As an additional consistency check, we calculate the partial correlation coefficients (Kendall & Stuart 1973) between log σ and logSFR 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 σ and logSFR 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_gas and, therefore, on σ and SFR. In Fig. E.3, we present the distributions of the partial correlation coefficients (r_partial) and the corresponding p-values. The analysis yields $r_{\mathrm{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 σ and logSFR 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 σ and logSFR is not merely an artifact of the residual effects of the assumed model.

Fig. E.1. SFR and σ 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.

Fig. E.2. Distribution of the Pearson rank correlation coefficient (upper panels) and the p-value (bottom panels) between logSFR and residuals (left panels) and log σ and residuals (right panels).

Fig. E.3. Distribution of the partial rank correlation coefficient (left panel) and the p-value (right panel) between logSFR and log σ, while excluding the residuals.

All Tables

All Figures

	Fig. 1. Location of the extended sample in the plane of main-sequence offset ΔMS vs. redshift. The gray dotted line shows the location of main-sequence galaxies with ΔMS = 1 for reference (Popesso et al. 2023). The gray solid line shows the separation between main-sequence and starburst populations at ΔMS = 4.
In the text

Fig. 2. Redshift distribution of velocity dispersion (σ, left panel) and rotational support (V/σ, right panel) for our extended sample (colored markers) consisting of massive galaxies. The σ and V/σ 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-σ uncertainties (dark pink area) and the intrinsic scatter (pink area).

In the text

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α (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 σ and V/σ derived from warm gas tracers (see Sect. 3.2).

In the text

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 trasport+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.

In the text

	Fig. 5. Distribution of σ 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-σ uncertainties (dark gray area) and the intrinsic scatter (gray area).
In the text

	Fig. 6. Predicted evolution of σ as a function of redshift for populations of main-sequence (ΔMS = 1, left panel) and starburst (ΔMS = 10, right panel) galaxies. Different lines show the evolution at different stellar masses, as indicated in the legend.
In the text

	Fig. 7. Evolution of the V/σ ratio for populations of main-sequence galaxies. Note that at z > 0, these profiles should be considered as lower limits.
In the text

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.

In the text

	Fig. B.1. Posterior distributions for the parameter defining the empirical relation in Eq. (2).
In the text

	Fig. B.2. Posterior distributions for the parameter defining the empirical relation in Eq. (3).
In the text

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⊙) = 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.

In the text

	Fig. C.2. Redshift distribution of σexp, DI and (V/σ)exp, DI for the low-z starbursts and the z ≳ 0.5 galaxies in the extended sample. The black dash-dotted lines show the best-fit relations.
In the text

	Fig. D.1. Same as Fig. 5 but with the slope of the log σ − logSFR fixed at 0.5 (dashed, dark green line) and 1 (dotted green line).
In the text

	Fig. E.1. SFR and σ 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.
In the text

	Fig. E.2. Distribution of the Pearson rank correlation coefficient (upper panels) and the p-value (bottom panels) between logSFR and residuals (left panels) and log σ and residuals (right panels).
In the text

	Fig. E.3. Distribution of the partial rank correlation coefficient (left panel) and the p-value (right panel) between logSFR and log σ, while excluding the residuals.
In the text