The gas in the interstellar medium ISM of galaxies is supersonically turbulent. Measurements of turbulence typically rely on cold gas emission lines for lowz galaxies and warm ionized gas observations for z gt 0 galaxies. Studies of warm gas kinematics at z gt 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 highz impacting our understanding of turbulence origin ISM physics and disk formation. We investigate the evolution of turbulence using velocity dispersion  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 mainsequence and starburst galaxies with stellar masses 10SUP10SUP MSUBSUB. The comparison with current H kinematic observations and existing models demonstrates that the velocity dispersion inferred from cold gas tracers differ by a factor of 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 . We present a model predicting the redshift evolution of turbulence in galaxy disks attributing the increase of  with redshift to the higher energy injected by supernovae due to the elevated starformation rate in highz galaxies. This supernovadriven 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  in MilkyWay like progenitors.
['galaxies: evolution', 'galaxies: high-redshift', 'galaxies: ISM', 'galaxies: kinematics and dynamics', 'galaxies: star formation', 'Astrophysics - Astrophysics of Galaxies', 'Astrophysics - Cosmology and Nongalactic Astrophysics']
The ALMAALPAKA survey II. Evolution of turbulence in galaxy disks across cosmic time Difference between cold and warm gas
{'The ALMA-ALPAKA survey II. Evolution of turbulence in galaxy disks across cosmic time: Difference between cold and warm gas': "F. Rizzo 1 , 2 , 3 , C. Bacchini 4 , M. Kohandel 5 , L. Di Mascolo 6 , F. Fraternali 3 , F. Roman-Oliveira 3 A. Zanella 7 , G. Popping 8 , F. Valentino 8 , G. Magdis 1 , 9 , and K. Whitaker 10  \n- 1 Cosmic Dawn Center (DAWN)\n- 2 Niels Bohr Institute, University of Copenhagen, Jagtvej 128, 2200 Copenhagen N, Denmark\n- 3 Kapteyn Astronomical Institute, University of Groningen, Landleven 12, 9747 AD, Groningen, The Netherlands e-mail: rizzo@astro.rug.nl\n- 4 Dark, Niels Bohr Institute, University of Copenhagen, Jagtvej 128, DK-2200 Copenhagen, Denmark\n- 5 Scuola Normale Superiore, Piazza dei Cavalieri 7, I-56126 Pisa, Italy\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\n- 7 Istituto Nazionale di Astrofisica, Vicolo dell'Osservatorio 5, 35122, Padova, Italy\n- 8 European Southern Observatory, Karl-Schwarzschild-Str. 2, D-85748 Garching bei München, Germany\n- 9 DTU-Space, Technical University of Denmark, Elektrovej 327, DK-2800, Kgs. Lyngby, Denmark\n- 10 Department of Astronomy, University of Massachusetts, Amherst, MA 01003, USA  \nJuly 16, 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 lowz 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 highz , impacting our understanding of turbulence origin, ISM physics and disk formation. We investigate the evolution of turbulence using velocity dispersion ( σ ) 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 ≳ 10 10 M ⊙ . The comparison with current H α kinematic observations and existing models demonstrates that the velocity dispersion inferred from cold gas tracers di ff er by a factor of ≈ 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 di ff erent from the conclusions of studies based on warm gas, which had to consider additional turbulence drivers to explain the high values of σ . We present a model predicting the redshift evolution of turbulence in galaxy disks, attributing the increase of σ with redshift to the higher energy injected by supernovae due to the elevated star-formation rate in highz 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 σ in Milky-Way like progenitors.  \nKey words. Galaxies: evolution - Galaxies: high-redshift - Galaxies: ISM - Galaxies: kinematics and dynamics', '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 highz (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.  \nFor 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 observa-  \nons 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 highz 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).  \nOne 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.  \n2018; 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. 2023) 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; Fraternali et al. 2021; Rizzo et al. 2021; Lelli et al. 2021; Tsukui & Iguchi 2021; Pope et al. 2023; Roman-Oliveira et al. 2023).  \nThe above discrepancy between the kinematics obtained from [CII] observations and those extrapolated using ionized gas could be due to intrinsic di ff erences 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. (2023) 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 a ff ected 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.  \nIn addition to the di ff erent 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 di ff erences 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 9 -10 10 M ⊙ ), while the [CII] studies at z > 4 have focused on massive ( M ⋆ ≳ 10 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 lowz data and 30 km / s for highz data. Combined with the typical angular resolution of ≈ 0.2", this kind of observations allows one to identify or isolate the circular motions in  \nthe rare cases when cold gas outflows are present (Rizzo et al. 2023; Roman-Oliveira et al. 2023).  \nUnderstanding whether the origin of the discrepancy between studies based on warm and cold gas tracers is driven by di ff erent 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 highz galaxies and the formation of galaxy bulges at highz 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 a ff ect 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).  \nGiven the aforementioned potential biases on the σ measured using the warm, ionized gas, key questions arise: Do the cold gas kinematics exhibit di ff erences 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?  \nIn 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. 2023). This particular redshift interval o ff ers 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 intermediateredshift 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  \nin normal main-sequence galaxies at z ∼ 2. Consequently, despite a decade of ALMA operations, there have been no highresolution systematic surveys of CO at z = 0 . 5 -4, but only observations of single targets.  \nIn 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, starforming 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.  \nIn 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 3D B arolo (Di Teodoro & Fraternali 2015) are presented and discussed in detail in Rizzo et al. (2023).  \nFor 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 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 (SFRMS). We quantify this o ff set ∆ MS = SFR / SFRMS( M ⋆, z ) using the prescriptions recently published by Popesso et al. (2023) and normalized to a Chabrier Initial Mass Function (IMF, Chabrier 2003):  \nMS ( M ⋆, z ) ≡ log " SFRMS( M ⋆, z ) M ⊙ yr -1 # = ( -0 . 034 t (Gyr) + 4 . 722) × -0 . 1925 × m 2 + + ( -26 . 16 + 0 . 2 t (Gyr)) , (1)  \nwhere t is the cosmic time and  \nm = " log M ⋆ M ⊙ ! + 0 . 025 # .  \nTo 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 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). This summary statistics allow us to quantify the spread of the ∆ MS values in three redshift bins: 0 -0 . 13, 0 . 5 -3 . 5,  \n4 -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 e ff ect, 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.  \nFig. 1. Location of the extended sample in the plane of main-sequence o ff set ∆ MS vs redshift. The gray dotted line shows the location of mainsequence 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.  \n<!-- image -->', '2.1. Low-z subsample': 'In the lowz 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.  \nHERACLES: 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 B arolo (Di Teodoro & Fraternali 2015). In our extended sample, we include the 4 galaxies from Bacchini et al. (2020a) with stellar masses M ⋆ ≳ 10 10 M ⊙ .  \nPHANGS: We include ten galaxies from the Physics at High Angular resolution in Nearby Galaxies (PHANGS)-ALMA survey (Leroy et al. 2021b). 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 3D B arolo (Di Teodoro & Fraternali 2015) to the publicly available data cubes (Leroy et al. 2021b,a, see Appendix A for details).  \nDYNAMO: The DYnamics of Newly Assembled Massive Objects (DYNAMO) survey (Green et al. 2014) includes massive ( M ⋆ > 10 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).  \nm  \nTable 1. Datasets included in the extended sample.Table 2. Percentiles of the ∆ MS distributions for the galaxies in the extended sample, divided in three redshift bins.', '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 lowz samples presented in Sect. 2.1 as their major axes are typically resolved with more than ten independent resolution elements.  \nIn 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).  \nAt 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. 2023) 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 di ff erent 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,  \nleveraging the sample compiled in this work o ff ers 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 -10 11 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 e ff ect, 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.  \nTo quantify these observed trends, we employ the following functional forms:  \nσ = a σ + b σ exp( -z 0 , y / z ) + s σ, (2)  \nV /σ = aV /σ + bV /σ exp( -z / z 0 , y ) + sV /σ. (3)  \nIn both equations, the parameters s σ and sV /σ 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  \nTable 3. Best-fit parameters defining the relations between σ -z , Eq (2) and V /σ -z , Eq (3).  \nbest 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 sV /σ . 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.  \nOur 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; Tsukui & Iguchi 2021; RomanOliveira et al. 2023; Lelli et al. 2023): 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).  \nThe 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. 2023). 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. (2023) 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  \n(log( M ⋆/ M ⊙ ) = 9 -11, see Appendix 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 di ff erences 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 di ff erent 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.  \nIn 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 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 di ff erence of ≈ 3 is obtained when we compare the median CO V /σ from our extended sample to the median H α V /σ for the SAMI galaxies.  \nIn 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 starforming 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 intermediatez 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 lower than their 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. 2023).  \nDue 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 di ff erence is mainly due to: (i) di ff erences in the disk thickness of the two gas phases; (ii) out-  \nFig. 2. Redshift distribution of velocity dispersion ( σ , left panels) 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).  \n<!-- image -->  \nFig. 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).  \n<!-- image -->  \nlow contaminating the measurements of turbulence from warm gas; (iii) insu ffi cient 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  \nTable 4. Main properties of the comparison and extended sample in four redshift bins.  \nNotes. Kinematics from: ( a ) Varidel et al. (2020); ( b ) Girard et al. (2021); ( c ) Hogan et al. (2021); ( d ) Birkin et al. (2023) .  \net 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),  \nQ = κσ π G Σ gas . (4)  \nIn 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 selfgravity 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 starforming 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:  \nσ ( z ) = Q √ 2 V f gas( z ) (5)  \nV σ = √ 2 f gas( z ) Q , (6)  \nFollowing 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:  \n- 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.\n- 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.  \nBoth 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 insu ffi cient to account for the observed turbulence.  \nKrumholz & 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 e ffi ciency free to vary as a function of σ ; no relationship between SFR and σ is predicted for the feedback model with fixed star formation e ffi -ciency (Appendix D). Comparisons of K18 models with warm gas σ of z = 1 -3 galaxies suggest that energy injected by stellar feedback is not su ffi cient to explain the high values of the observed velocity dispersion. In contrast, the transport model is in  \ngood agreement with highz observations derived from warm gas (see also Johnson et al. 2018; Übler et al. 2019). The feedback model alone reproduces the observations at lowz . The fiducial model proposed by K18 is, therefore, a combination of the feedback and transport models and it can reproduce both low and highz galaxies.  \nIn 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 highz galaxies reported by K18 and we use the range of values of our extended sample (shaded area).  \nThe σ -SFR distribution of our sample aligns with the predictions of the stellar feedback model proposed by K18. Conversely, the transport + feedback model predicts, that highz 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.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 e ffi ciency and the Toomre parameter (see Appendix D for details). Additionally, the K18 models rely on a series of free parameters that are adjusted for di ff erent 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 e ff ectively 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, ˙ E SN, equal to the rate at which turbulent kinetic energy, E kin is dissipated, that is, ˙ E SN = E kin /τ . The kinetic energy and supernova energy rate are defined as  \nE kin = 3 2 M gas σ 2 (7)  \n˙ E SN = ϵ SN SFR η SN E SN (8)  \nwhere ϵ SN is the e ffi ciency at which the kinetic energy is transferred to the ISM, η SN is the supernova rate, E SN is the supernova energy equal to 10 51 erg. The timescale of turbulence dissipation τ can be defined as τ = 2 h /σ (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:  \nh ( R ) = σ ( R ) " ∂ϕ ( R , 0) ∂ Z # -1 / 2 (9)  \nArticle number, page 8 of 22  \nwhere ϕ ( R , 0) is the galactic gravitational potential in the plane of the disk ( Z = 0). Using ˙ E SN = E kin /τ , we derive the relation between σ and SFR:  \nσ = SFR 1 / 3 4 ϵ SN η SN E SN h 3 M gas ! 1 / 3 . (10)  \nEstablishing 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 lowz and 13 at highz ) from our extended sample. These profiles are presented and described in Bacchini et al. (2019) and Bacchini et al. (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 ϵ SN h / M gas as a function of σ , SFR, and σ × SFR. This analysis revealed that the dependencies of h and ϵ SN h / 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 lowz galaxies, CC ∼ 10 -3 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 9 M ⊙ (e.g., Leroy et al. 2021a). For highz galaxies, CC ∼ 10 -3 , considering: h ∼ 200 pc (Bacchini et al. 2024), ϵ SN ∼ 0 . 01 -0 . 1 (Rizzo et al. 2020, 2021; RomanOliveira et al. 2024), M gas ∼ 10 10 M ⊙ . Throughout the rest of the paper, we adopt CC as a constant value, leading to σ ∝ SFR 1 / 3 , which di ff ers 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:  \nlog σ = 0 . 33 × log SFR + (0 . 77 ± 0 . 03) , (11)  \nwith an intrinsic scatter of 0.13 dex. In Eq. (11), σ is in units of km s -1 and SFR is in units of M ⊙ yr -1 . 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. 2015) indicate that the model with the fixed slope is favored over the one with the free value. This is  \nFig. 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 highz galaxies, as reported in K18. For the highz 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.  \n<!-- image -->  \nanother 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 log 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 ⋆ 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, highz 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:  \nlog σ ( M ⋆, z , ∆ MS) = 0 . 33 × [MS( M ⋆, z ) + log ∆ MS] + 0 . 77 . (12)  \nIt is important to emphasize that this evolution of σ as a function of redshift is entirely due to the evolution of the mainsequence relation. The physics of turbulence driving does not appear to evolve with time (see Fig. 5).  \nFig. 6 shows the redshift evolution for populations of mainsequence galaxies with M ⋆ of 10 11 M ⊙ and 10 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,  \nFig. 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).  \n<!-- image -->  \nthe supernova-driven model predicts that, at fixed redshift, mainsequence 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 9 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.  \nFig. 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. Di ff erent lines show the evolution at di ff erent stellar masses, as indicated in the legend.  \n<!-- image -->  \nThe 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. 2023).  \nThe 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 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 10 M ⊙ at z = 2 to a stellar mass of ≈ 5 × 10 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).  \nPredicting 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 Λ CMDcosmology framework, the stellar TFR is expected to evolve such that at fixed stellar mass, highz 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  \nReyes 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 9 M ⊙ ) galaxy population.  \nFig. 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.  \n<!-- image -->  \nAs 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 evo-  \nion 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.', '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 highz galaxies poses significant challenges. A primary limitation arises from angular resolution, often implying that the emission of highz galaxies is resolved by only few independent resolution elements. Consequently, the kinematics of highz galaxies can be influenced by the beam-smearing e ff ect, 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 e ff ect, overcoming this issue in the case of low-resolution data remains di ffi cult. 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 highz 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.  \nThe 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 highz 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 di ff erent spatial and spectral resolutions and for galaxies with di ff erent 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 di ff erences 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  \ninstead 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 2 -10 3 K) than that traced by CO (T ≲ 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 di ff use 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 highz , 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).  \nAlthough 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 e ffi ciency 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 di ff erent driving mechanisms in the di ff erent 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 -10 11 M ⊙ ), but a wide range of 4 order of magnitude in SFRs. Our main results are the following.  \n- -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  \nFig. 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.  \n<!-- image -->  \n- 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.\n- -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 highz galaxies.\n- -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.\n- -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  \ngalaxies 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 highz , low mass ( M ⋆ ≳ 10 9 M ⊙ ) galaxy disks should be rotationally supported, with V /σ ∼ 10.  \nThe results presented in this paper indicate that the dynamical properties of galaxies significantly di ff er 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 di ff erence between cold and warm gas kinematics evidenced in our analysis may be due to (i) di ff use 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) insu ffi cient spectral and angular resolution. In the future, studies of the ionized gas kinematics using high an-  \ngular and spectral resolution observations (FWHM ≲ 40 km / s, e.g., VLT / ERIS; ELT / HARMONI) will be able to discriminate between these three hypotheses.  \nThe 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 starforming 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.  \nAcknowledgements. 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 SklodowskaCurie 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  \ndata: ADS / JAO.ALMA#2017.1.01659.S; ADS / JAO.ALMA#2017.1.00471.S;  \nADS / JAO.ALMA#2017.1.01674.S;  \nADS / JAO.ALMA#2015.1.00862.S;  \nADS / JAO.ALMA#2017.1.01228.S;  \nADS  \n/  \nJAO.ALMA#2018.1.00974.S;  \nADS / JAO.ALMA#2017.1.00413.S;  \nADS  \n/  \nJAO.ALMA#2019.1.01362.S;  \nADS / JAO.ALMA#2017.1.01045.S;  \nADS / JAO.ALMA#2013.1.00059.S;  \nADS / JAO.ALMA#2018.1.00543.S;  \nADS / JAO.ALMA#2015.1.00723.S;  \nADS  \n/  \nJAO.ALMA#2018.1.01146.S;  \nADS / JAO.ALMA#2016.1.01155.S;  \nADS / JAO.ALMA#2017.1.01677.S;  \nADS / JAO.ALMA#2018.1.01306.S;  \nADS / JAO.ALMA#2017.1.00908.S;  \nADS / JAO.ALMA#2012.1.00650.S;  \nADS / JAO.ALMA#2013.1.00803.S;  \nADS / JAO.ALMA#2013.1.01161.S;  \nADS / JAO.ALMA#2015.1.00121.S;  \nADS / JAO.ALMA#2015.1.00782.S;  \nADS / JAO.ALMA#2015.1.00925.S;  \nADS / JAO.ALMA#2015.1.00956.S;  \nADS  \n/  \nJAO.ALMA#2016.1.00386.S;  \nADS  \n/  \nJAO.ALMA#2017.1.00392.S;  \nADS  \n/  \nJAO.ALMA#2017.1.00766.S;  \nADS  \n/  \nJAO.ALMA#2017.1.00886.L;  \nADS / JAO.ALMA#2018.1.00484.S;  \nADS  \n/  \nJAO.ALMA#2018.1.01321.S;  \nADS / JAO.ALMA#2018.1.01  \n651.S;  \nADS / JAO.ALMA#2018.A.00062.S;  \nADS / JAO.ALMA#2019.1.01235.S;  \nADS  \n/  \nJAO.ALMA#2019.2.00129.S;  \nADS / JAO.ALMA#2016.1.01499.S;  \nADS  \n/  \nJAO.ALMA#2015.1.00330.S;  \nADS / JAO.ALMA#2017.1.00127.S;  \nADS  \n/  \nJAO.ALMA#2017.1.00394.S;  \nADS  \n/  \nJAO.ALMA#2017.1.01052.S;  \nADS  \n/  \nJAO.ALMA#2018.1.00001.S;  \nADS / JAO.ALMA#2017.1.01020.S;  \nADS  \n/  \nJAO.ALMA#2019.1.00862.S;  \nADS / JAO.ALMA#2017.1.01471.S;  \nADS / JAO.ALMA#2015.1.00456.S;  \nADS / JAO.ALMA#2011.0.00124.S;  \nADS  \n/  \nJAO.ALMA#2012.1.00978.S;  \nADS / 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.  \nWe acknowledge usage of the Python programming language (Van Rossum & Drake 2009), Astropy (Astropy Collaboration et al. 2013), Matplotlib (Hunter 2007), NumPy (van der Walt et al. 2011), and SciPy (Virtanen et al. 2020).", 'References': 'Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33  \n- Bacchini, C., Fraternali, F., Iorio, G., & Pezzulli, G. 2019, A&A, 622, A64\n- Bacchini, C., Fraternali, F., Iorio, G., et al. 2020a, A&A, 641, A70\n- Bacchini, C., Fraternali, F., Pezzulli, G., & Marasco, A. 2020b, A&A, 644, A125\n- Bacchini, C., Nipoti, C., Iorio, G., et al. 2024, arXiv e-prints, arXiv:2405.00103\n- Birkin, J. E., Smail, I., Swinbank, A. M., et al. 2023, arXiv e-prints, arXiv:2301.05720\n- Bland-Hawthorn, J. & Gerhard, O. 2016, ARA&A, 54, 529\n- Bournaud, F., Elmegreen, B. G., Teyssier, R., Block, D. L., & Puerari, I. 2010, MNRAS, 409, 1088\n- Bournaud, F., Perret, V., Renaud, F., et al. 2014, ApJ, 780, 57\n- Braun, R., Bonaldi, A., Bourke, T., Keane, E., & Wagg, J. 2019, arXiv e-prints, arXiv:1912.12699\n- Carilli, C. L. & Walter, F. 2013, ARA&A, 51, 105\n- Carpenter, J., Iono, D., Kemper, F., & Wootten, A. 2020, arXiv e-prints, arXiv:2001.11076\n- Chabrier, G. 2003, PASP, 115, 763 Cio ffi , 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\n- Curran, P. A. 2014, arXiv e-prints, arXiv:1411.3816\n- De Breuck, C., Weiß, A., Béthermin, M., et al. 2019, A&A, 631, A167\n- Dekel, A., Zolotov, A., Tweed, D., et al. 2013, MNRAS, 435, 999\n- Di Mascolo, L. 2024, lucadimascolo / scattr\n- Di Teodoro, E. M. & Fraternali, F. 2015, MNRAS, 451, 3021\n- Di Teodoro, E. M., Fraternali, F., & Miller, S. H. 2016, A&A, 594, A77\n- Díaz-Santos, T., Armus, L., Charmandaris, V., et al. 2017, ApJ, 846, 32\n- Dutton, A. A., van den Bosch, F. C., Faber, S. M., et al. 2011, MNRAS, 410, 1660\n- Elmegreen, B. G. & Scalo, J. 2004, ARA&A, 42, 211\n- Epinat, B., Amram, P., & Marcelin, M. 2008, MNRAS, 390, 466\n- Epinat, B., Contini, T., Le Fèvre, O., et al. 2009, A&A, 504, 789\n- Fisher, D. B., Glazebrook, K., Bolatto, A., et al. 2014, ApJ, 790, L30\n- Fisher, D. B., Glazebrook, K., Damjanov, I., et al. 2017, MNRAS, 464, 491\n- Forbes, J. C., Emami, R., Somerville, R. S., et al. 2023, ApJ, 948, 107\n- Förster Schreiber, N. M., Genzel, R., Bouché, N., et al. 2009, ApJ, 706, 1364\n- Förster Schreiber, N. M., Renzini, A., Mancini, C., et al. 2018, ApJS, 238, 21\n- Förster Schreiber, N. M. & Wuyts, S. 2020, ARA&A, 58, 661\n- Fraternali, F., Karim, A., Magnelli, B., et al. 2021, A&A, 647, A194\n- Genzel, R., Newman, S., Jones, T., et al. 2011, ApJ, 733, 101\n- Ginzburg, O., Dekel, A., Mandelker, N., & Krumholz, M. R. 2022, MNRAS, 513, 6177\n- 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\n- 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\n- Hayward, C. C. & Hopkins, P. F. 2017, MNRAS, 465, 1682\n- Ho ff man, M. D. & Gelman, A. 2011, arXiv e-prints, arXiv:1111.4246\n- Hogan, L., Rigopoulou, D., Magdis, G. E., et al. 2021, MNRAS, 503, 5329  \nHunter, J. D. 2007, Computing in Science Engineering, 9, 90  \n- Ianjamasimanana, R., de Blok, W. J. G., Walter, F., et al. 2015, AJ, 150, 47\n- Jiménez, E., Lagos, C. d. P., Ludlow, A. D., & Wisnioski, E. 2023, MNRAS, 524, 4346\n- Johnson, H. L., Harrison, C. M., Swinbank, A. M., et al. 2018, MNRAS, 474, 5076\n- Kadler, M., Riechers, D. A., Baczko, A. K., et al. 2023, arXiv e-prints, arXiv:2311.10056\n- Kassin, S. A., Weiner, B. J., Faber, S. M., et al. 2012, ApJ, 758, 106  \nKelly, B. C. 2007, ApJ, 665, 1489  \n- Kendall, M. & Stuart, A. 1973, The Advanced Theory of Statistics. Vol. 2: Inference and: Relationsship (Gri ffi n)\n- Kim, C.-G. & Ostriker, E. C. 2015, ApJ, 815, 67\n- Kim, C.-G. & Ostriker, E. C. 2017, ApJ, 846, 133\n- Klessen, R. S. & Hennebelle, P. 2010, A&A, 520, A17\n- Kohandel, M., Pallottini, A., Ferrara, A., et al. 2020, MNRAS, 499, 1250\n- Kohandel, M., Pallottini, A., Ferrara, A., et al. 2023, arXiv e-prints, arXiv:2311.05832\n- Kramer, E. D. & Randall, L. 2016, ApJ, 829, 126\n- Kretschmer, M., Dekel, A., & Teyssier, R. 2022, MNRAS, 510, 3266\n- Krumholz, M. R. & Burkhart, B. 2016, MNRAS, 458, 1671\n- Krumholz, M. R., Burkhart, B., Forbes, J. C., & Crocker, R. M. 2018, MNRAS, 477, 2716  \nLaw, D. R., Belfiore, F., Bershady, M. A., et al. 2022, ApJ, 928, 58  \n- Lelli, F. 2022, Nature Astronomy, 6, 35\n- Lelli, F., Di Teodoro, E. M., Fraternali, F., et al. 2021, Science, 371, 713\n- Lelli, F., Zhang, Z.-Y., Bisbas, T. G., et al. 2023, A&A, 672, A106\n- Leroy, A. K., Hughes, A., Liu, D., et al. 2021a, ApJS, 255, 19\n- Leroy, A. K., Schinnerer, E., Hughes, A., et al. 2021b, ApJS, 257, 43\n- Leroy, A. K., Walter, F., Bigiel, F., et al. 2009, AJ, 137, 4670\n- Mac Low, M.-M. & Klessen, R. S. 2004, Reviews of Modern Physics, 76, 125\n- Mac Low, M.-M., Klessen, R. S., Burkert, A., & Smith, M. D. 1998, Phys. Rev. Lett., 80, 2754\n- Marasco, A., Fraternali, F., van der Hulst, J. M., & Oosterloo, T. 2017, A&A, 607, A106\n- 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\n- Neal, R. 2011, in Handbook of Markov Chain Monte Carlo, 113-162\n- Neeleman, M., Prochaska, J. X., Kanekar, N., & Rafelski, M. 2020, Nature, 581, 269\n- Orr, M. E., Hayward, C. C., Medling, A. M., et al. 2020, MNRAS, 496, 1620\n- Ostriker, E. C. & Kim, C.-G. 2022, ApJ, 936, 137\n- Ostriker, E. C. & Shetty, R. 2011, ApJ, 731, 41\n- Pacifici, C., Iyer, K. G., Mobasher, B., et al. 2023, ApJ, 944, 141\n- Pallottini, A., Ferrara, A., Gallerani, S., et al. 2022, MNRAS, 513, 5621\n- Parlanti, E., Carniani, S., Übler, H., et al. 2023, arXiv e-prints, arXiv:2309.05713  \nPillepich, 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  \nPope, A., McKinney, J., Kamieneski, P., et al. 2023, ApJ, 951, L46 Popesso, 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ˇci¯ut˙e, 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  \nRigopoulou, D., Hopwood, R., Magdis, G. E., et al. 2014, ApJ, 781, L15  \nRizzo, F., Kohandel, M., Pallottini, A., et al. 2022, A&A, 667, A5  \n- Rizzo, F., Roman-Oliveira, F., Fraternali, F., et al. 2023, A&A, 679, A129\n- Rizzo, F., Vegetti, S., Fraternali, F., Stacey, H. R., & Powell, D. 2021, MNRAS, 507, 3952  \nRizzo, F., Vegetti, S., Powell, D., et al. 2020, Nature, 584, 201  \nRodighiero, G., Daddi, E., Baronchelli, I., et al. 2011, ApJ, 739, L40  \nRoman-Oliveira, F., Fraternali, F., & Rizzo, F. 2023, MNRAS, 521, 1045  \nRoman-Oliveira, F., Fraternali, F., & Rizzo, F. 2024, MNRAS, 521, 1045  \nSawicki, M. 2012, PASP, 124, 1208  \nScott, N., van de Sande, J., Croom, S. M., et al. 2018, MNRAS, 481, 2299 Sereno, M. 2016, MNRAS, 455, 2149  \nSimons, R. C., Kassin, S. A., Trump, J. R., et al. 2016, ApJ, 830, 14  \nSimons, R. C., Kassin, S. A., Weiner, B. J., et al. 2017, ApJ, 843, 46  \nStott, J. P., Swinbank, A. M., Johnson, H. L., et al. 2016, MNRAS, 457, 1888  \nTacconi, L. J., Neri, R., Genzel, R., et al. 2013, ApJ, 768, 74  \nThornton, K., Gaudlitz, M., Janka, H. T., & Steinmetz, M. 1998, ApJ, 500, 95  \nToomre, A. 1964, ApJ, 139, 1217  \nTsukui, T. & Iguchi, S. 2021, Science, 372, 1201  \nTully, R. B. & Fisher, J. R. 1977, A&A, 54, 661  \nTurner, O. J., Cirasuolo, M., Harrison, C. M., et al. 2017, MNRAS, 471, 1280  \nÜbler, H., Genzel, R., Tacconi, L. J., et al. 2018, ApJ, 854, L24  \nÜbler, H., Genzel, R., Wisnioski, E., et al. 2019, ApJ, 880, 48  \nUtomo, D., Blitz, L., & Falgarone, E. 2019, ApJ, 871, 17  \nVallini, L., Gallerani, S., Ferrara, A., Pallottini, A., & Yue, B. 2015, ApJ, 813, 36 van der Walt, S., Colbert, S. C., & Varoquaux, G. 2011, Computing in Science Engineering, 13, 22  \nvan Dokkum, P. G., Leja, J., Nelson, E. J., et al. 2013, ApJ, 771, L35  \nVan Rossum, G. & Drake, F. L. 2009, Python 3 Reference Manual (Scotts Valley,  \nCA: CreateSpace)  \nVaridel, M. R., Croom, S. M., Lewis, G. F., et al. 2020, MNRAS, 495, 2265 Vehtari, A., Gelman, A., & Gabry, J. 2015, arXiv e-prints, arXiv:1507.04544 Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261 Whitaker, K. E., van Dokkum, P. G., Brammer, G., & Franx, M. 2012, ApJ, 754, L29  \nWisnioski, 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 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 2 (Leroy et al. 2021a,b) to fit their kinematics using 3D B arolo . The parameter files used for the fitting and the corresponding outputs are available online 3 .', 'Appendix B: Fitting empirical relations': 'To fit Eqs. (2), (3) and (11), we use a hierarchical Bayesian approach 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 (Ho ff man & Gelman 2011) adaptation of the HMC step length.', 'Appendix B.1: Fitting equations (2) and (3)': 'We aim at constraining the set of free parameters θ = { a , b , z 0 } 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 a ff ected by intrinsic scatter and characterized by negligible uncertainties. This is related to the corresponding observed quantity y through the relation  \ny = η y ( z ) + ϵ y , (B.1)  \nwhere ϵ y is the corresponding measurement uncertainty that follows a normal probability distribution with known variance u 2 y . 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 N{ f ( z ; θ ) , s 2 y } with mean equal to f ( z ; θ ) and variance s 2 y . 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  \np ( y | z , θ ) = Z p ( y | z , η ) p ( η | z , θ ) d η. (B.2)  \nSince 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).  \nTable B.1. Ranges of the uniform priors used to perform the fitting of the σ -z and V /σ -z distributions using Eqs. (2) and (3).  \nFig. B.1. Posterior distributions for the parameter defining the empirical relation in Eq. (2).  \n<!-- image -->', 'Appendix 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:  \nx = η x ( ξ ) + ϵ x , η x = N{ fx ( ξ ; θ x ) , s 2 x } y = η y ( ξ ) + ϵ y , η y = N{ fy ( ξ ; θ y ) , s 2 y } . (B.3)  \nIn 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 fi , parameters θ i and intrinsic scatter s 2 i . 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 N ( fi ; s 2 i ). In our specific case,  \nfx = fx ( ξ ) = 10 ξ , fy = fy ( ξ ; b ) = 10 0 . 33 ξ + b . (B.4)  \nFor the prior probability on the o ff set parameter b in Eq. (B.4), we assume a wide uniform prior in the range  \nFig. B.2. Posterior distributions for the parameter defining the empirical relation in Eq. (3)  \n<!-- image -->  \n.  \n[ -100 , 100]. The joint likelihood function can instead be rewritten as  \np ( x , y | θ ) = " p ( x , y | η x , η y ) p ( η y | ξ, θ ) p ( η x | ξ, θ ) p ( ξ ) d η xd η y , (B.5)  \nwhere 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 2 i ( 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:  \np ( x , y | η x , η y ) = n Y k = 1 1 q 2 π u 2 x , k exp        -1 2 xk -η x , k ux , k ! 2        × n -m Y k = 1 1 q 2 π u 2 y , k exp        -1 2 yk -η y , k uy , k ! 2        × n Y k = m 1 2        1 + erf        yk -η y , k √ 2 uy , k               . (B.6)  \nArticle number, page 16 of 22', 'Appendix C: Further comparison with previous works': 'Appendix C.1: Warm gas kinematics from main-sequence galaxies  \nIn 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.  \nFig. 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.  \n<!-- image -->  \nIn 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. 2023). However, this discrepancy between observations and simulations might be ex-  \nstematic bias towards high σ due to the beamsmearing e ff ect (Sect. 6). This e ff ect might be amplified in lowmass 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.', 'Appendix 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 spatiallyresolved 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.', 'Appendix 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: lowz 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.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 (Cio ffi 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  \nis, therefore, defined, as  \n˙ E SN , K18 = Σ SFR p ⋆ m ⋆ ! σ, (D.1)  \nwhere Σ SFR is the SFR per unit area. Similarly to F21, the kinetic energy dissipated per unit area is defined as  \n˙ E kin , K18 = 3 2 Σ gas σ 3 h . (D.2)  \nIn 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 di ff erence between F21 and K18 is that, in K18, the rate of star formation per unit area is defined by the star formation law:  \nΣ SFR = ϵ ff Σ gas t ff , (D.3)  \nwhere ϵ ff is the star formation e ffi ciency and t ff is the free-fall time and Σ gas is expressed as a function of σ and Q , Eq. (4).  \nIn 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 e ffi ciency 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 σ ∝ Σ 1 / 2 SFR ∝ 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.', 'Appendix 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  \n˙ E SN , p = SFR p ⋆ m ⋆ ! σ. (D.4)  \nBy equating this expression with Eq. (7), we obtain the following relationship:  \nσ 2 = SFR p ⋆ m ⋆ ! 4 h 3 M gas (D.5)  \nThe primary di ff erence from the F21 model lies in the slope. Here, σ ∝ SFR 1 / 2 when assuming a constant value for GLYPH<16> p ⋆ m ⋆ GLYPH<17> as in K18, or σ ∝ SFR when assuming p ⋆ ∝ σ . In contrast, F21 finds that σ ∝ SFR 1 / 3 .  \nFor a consistency check, we fit the relation log σ = m × log 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.  \nFig. C.2. Redshift distribution of σ exp , DI and ( V /σ )exp , DI for the lowz starbursts and the z ≳ 0 . 5 galaxies in the extended sample. The black dashdotted lines show the best-fit relations.  \n<!-- image -->  \nFig. 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).  \n<!-- image -->', '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 log SFR (Fig. E.1). Any correlations between the residuals and log σ and log 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,  \nwe compute the orthogonal distance between each data point (log σ i , log SFRi) 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 coe ffi cient ρ (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 coe ffi cients are r Pearson = -0 . 23 + 0 . 15 -0 . 31 for the log SFR vs residual and r Pearson = 0 . 13 + 0 . 23 -0 . 03 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 log SFR, and between residuals and log σ , respectively.  \nAs an additional consistency check, we calculate the partial correlation coe ffi cients (Kendall & Stuart 1973) between log σ and log 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 σ and log 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 gas and, therefore, on σ and SFR. In Fig. E.3, we present the distributions of the partial correlation coe ffi cients ( r partial) and the corresponding p-values. The analysis yields r partial = 0 . 86 + 0 . 04 -0 . 06 . Based on the p-value distribution, where 100% of pvalues are less than 0.05, we conclude that the correlation between log σ and log 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  \nthat the relationship between log σ and log SFR is not merely an artifact of the residual e ff ects of the assumed model.  \n4  \nFig. 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.  \n<!-- image -->  \nF. Rizzo et al.: Dynamics from warm and cold gas tracersFig. E.2. Distribution of the Pearson rank correlation coe ffi cient (upper panels) and the p-value (bottom panels) between log SFR and residuals (left panels) and log σ and residuals (right panels).  \n<!-- image -->', 'A & A proofs: manuscript no. main': 'Fig. E.3. Distribution of the partial rank correlation coe ffi cient (left panel) and the p-value (right panel) between log SFR and log σ , while excluding the residuals.  \n<!-- image -->'}