© The Authors 2023

1. Introduction

The galaxy assembly and evolution in the first billion years of the Universe is still a hot topic of modern astrophysics. Primeval galaxies are thought to be subjected to galaxy-galaxy interactions, merging processes, gravitational instabilities, and feedback from bursts of star formation and active galactic nuclei (AGN; Dekel et al. 2009; Krumholz & Burkert 2010; Green et al. 2014; Somerville & Davé 2015; Dayal & Ferrara 2018; Krumholz et al. 2018; Ginzburg et al. 2022), which may affect the formation of the galactic disks and their properties. It is therefore fundamental to explore the kinematical properties of galaxies to understand the main processes that play a crucial role in the galaxy growth over cosmic time.

Several observing programs have focused on the gas kinematics of galaxies at 0 < z < 4 and provided thousands of spatially resolved observations of ionized, neutral, and molecular gas (Cresci et al. 2009; Förster Schreiber et al. 2009, 2018; Epinat et al. 2010; Gnerucci et al. 2011; Ianjamasimanana et al. 2012; Green et al. 2014; Wisnioski et al. 2015; Mogotsi et al. 2016; Di Teodoro et al. 2016; Harrison et al. 2017; Swinbank et al. 2017; Turner et al. 2017; Johnson et al. 2018; Übler et al. 2019; Girard et al. 2021). The kinematic analysis of these galaxies shows that the intrinsic velocity dispersion of gas (σ_gas) increases with redshift, suggesting that galaxies in the cosmic noon epoch are more turbulent than their local counterparts (Wisnioski et al. 2015, 2019; Johnson et al. 2018; Übler et al. 2019). These observations have differentiated between rotation-dominated and dispersion-dominated disks by comparing the ordered rotational velocities, V_max (the maximum rotational velocity of the galaxy), with the σ_gas. The distinction between the two kinematics is mainly based on the V_max/σ_gas, with a threshold value that spans the range between 1 and 3. Studies of the high-z Universe usually adopted the threshold based on the kinematics of the exponential disk¹. A ratio of indicates that the dynamical support from random motions of gas is stronger than the rotational support from ordered motions (Binney & Tremaine 2008; Wisnioski et al. 2019; Förster Schreiber & Wuyts 2020). In summary, high V_max/σ_gas ratios suggest a dynamically cold disk, and low values are usually associated with a turbulent gas disk that is mainly supported by random motions.

At 0 < z < 1, the bulk of the galaxy population with a stellar mass M_⋆ between 10⁹ M_⊙ and 10¹¹ M_⊙ shows a σ_gas ∼ 10 − 40 km s⁻¹, depending on the gas tracer, and a V_max/σ_gas ratio of ∼10 that indicates a rotationally supported gas disk similar to that of our Galaxy (Epinat et al. 2010; Ianjamasimanana et al. 2012; Di Teodoro et al. 2016; Mogotsi et al. 2016; Wisnioski et al. 2015; Harrison et al. 2017; Swinbank et al. 2017; Johnson et al. 2018; Übler et al. 2019). Observations in the distant Universe instead reveal a different scenario. At intermediate redshift (1 < z < 3), different studies found galaxies with σ_gas up to 70 km s⁻¹ and V_max/σ_gas values down to 1–2, suggesting that distant galaxies can be described as rotating turbulent disks (Cresci et al. 2009; Förster Schreiber et al. 2009, 2018; Wisnioski et al. 2015; Turner et al. 2017; Übler et al. 2019; Girard et al. 2021).

In the past decade, several models have been proposed to explain the observed evolution of the velocity dispersion with redshift (Krumholz & Burkert 2010; Forbes et al. 2014; Krumholz et al. 2018; Ginzburg et al. 2022; Ostriker & Kim 2022). A continuous energy injection into the system is necessary to maintain its high-velocity dispersion. This energy is thought to derive from three main mechanisms: feedback from star formation and AGN, gravitational instabilities within the disk, and gas accretion. Supernova explosions and stellar winds are expected to inject energy and momentum into the surrounding environment. These processes increase the temperature and the velocity dispersion of the gas in the galaxy (Hayward & Hopkins 2017; Orr et al. 2020; Ostriker & Kim 2022). Mergers and gravitational interactions between galaxies can also give rise to perturbations and generate gravitationally unstable regions in the galactic disk (Kohandel et al. 2020). In these regions, the axisymmetry of the disk is broken, promoting the formation of clumps (Zanella et al. 2015; Kohandel et al. 2019). The gravitational field will exert a torque on the clumps that will transport the mass within the disk from a larger to a smaller radius, which increases the turbulence and the star formation rate (SFR). The transport of mass would release gravitational energy and restore the stability of the disk (Krumholz & Burkert 2010; Forbes et al. 2014). Recently, even the accretion process of gas from the circumgalactic medium was taken into account as it increases the turbulence in galaxies. The gas accretes onto the galaxy via cold streams, and the streams fragment and turn into clumps. The energy of the accreting gas can turn into heat via shock or can remain kinetic energy that will turn into turbulence (Dekel et al. 2009; Forbes et al. 2023; Ginzburg et al. 2022).

The increase in velocity dispersion over cosmic time and the decrease in V_max/σ_gas also seem to be supported by cosmological simulations. These simulations report an increase in the velocity dispersion at fixed stellar masses (Dekel & Burkert 2014; Zolotov et al. 2015; Pillepich et al. 2019). The velocity dispersion of simulated high-redshift galaxies is higher by an order of magnitude than in their local counterparts. The simulations are consistent with the observations up to z = 4, supporting the scenario that the contribution from random motion increases across cosmic time.

Wisnioski et al. (2015) claimed that the evolution of the gas kinematical properties of galaxies only depends on the evolution of the gas fraction (f_gas), (Tacconi et al. 2020). Star-forming galaxies are thought to be in steady equilibrium between outflow, star formation, and gas inflow (Förster Schreiber et al. 2006; Mannucci et al. 2010; Genzel et al. 2011), and the balance between the heating and the cooling of gas maintains the disk in a quasi-stable state with the Toomre parameter Q ∼ 1 (Toomre 1964). This assumption leads to the relation

(1)

where 1 < a < 2, in particular, for a flat rotation curve, for a uniform disk, and a = 2 for a solid-body rotation. With this model based on marginally stable disks, the decrease in V_max/σ_gas and the increase in σ_gas at increasing redshift is due to the increment of the gas fraction at high redshift (Tacconi et al. 2020). This model would also explain the observed clumpy structure of high-redshift galaxies (Förster Schreiber et al. 2006, 2011; Genzel et al. 2008; Carniani et al. 2017, 2018) because Toomre’s stability criterion leads to the formation of star-forming clumps in gas-rich disk when Q ≤ 1.

With the advent of the Atacama Large Millimeter/sub-Millimeter Array (ALMA), we can exploit the brightest rest-frame far-infrared (FIR) lines, which fall in the millimeter bands (0.5–2 mm), to trace the dynamics of galaxies at z > 4. Among the rest-frame FIR lines, the [C II] (157.7 μm) and [O III] (88 μm) are the brightest lines (Carilli & Walter 2013). The [C II] transition is one of the principal coolants of the interstellar medium (ISM) and one of the most luminous FIR lines (Stacey et al. 1991; Malhotra et al. 1997; Luhman et al. 1998, 2003). The ionization potential of the carbon atom is ∼11.3 eV, which is lower than that of hydrogen (i.e., 13.6 eV), and can be used to trace different phases (neutral, molecular, and mildly ionized gas) of the ISM of a galaxy (Shibai et al. 1991; Heiles 1994; Stacey et al. 2010; Pineda et al. 2013; Vallini et al. 2015). On the other hand, the transition lines of [O III] are good tracers of H II regions because the oxygen atom has an ionization potential of 35.1 eV (Carilli & Walter 2013).

Through ALMA, a multitude of galaxies have been observed at z > 4 targeting FIR emission lines (Capak et al. 2015; Willott et al. 2015; Carniani et al. 2017, 2018; Smit et al. 2018; Harikane et al. 2020; Neeleman et al. 2020; Rizzo et al. 2020, 2021; Lelli et al. 2021; Fraternali et al. 2021; Jones et al. 2021; Herrera-Camus et al. 2021, 2022; Schouws et al. 2022; Bouwens et al. 2022). So far, only a few studies have focused on the gas kinematics of these high-redshift systems, and the results yield a contrasting picture of the kinematic properties of galaxies in the primeval Universe. On the one hand, Jones et al. (2021) and Herrera-Camus et al. (2021) found a range of the velocity dispersion from 20 km s⁻¹ to 116 km s⁻¹ for nine galaxies at 4 < z < 6, and V_max/σ_gas ∼ 3, which is consistent with what is observed at z = 1 − 3. On the other hand, other studies have found a lower velocity dispersion and a higher rotational velocity by studying the [C II] emission line for ten massive starburst dusty galaxies at z ∼ 4 − 5 (Sharda et al. 2019; Neeleman et al. 2020; Rizzo et al. 2020, 2021; Fraternali et al. 2021; Lelli et al. 2021). The measured velocity dispersion is as low as 15–20 km s⁻¹, with an average V_max/σ_gas of 10 and an extreme value of 20 for one of the selected galaxies (Fraternali et al. 2021). These estimates are comparable to the σ_gas and V_max/σ_gas obtained for local spiral galaxies (Epinat et al. 2010; Ianjamasimanana et al. 2012; Di Teodoro et al. 2016; Mogotsi et al. 2016). This implies that the internal ISM turbulence of these galaxies is similar to that of local galaxies, and despite their young age, the gas is already settled into a dynamically cold rotating disk. The existence of dynamically cold disks at z > 6 is predicted in cosmological simulations (Pallottini et al. 2019; Kohandel et al. 2020), where V_max/σ ∼ 7 was reported for a simulated Lyman-Break galaxy during its evolution when observed with the [CII] emission line.

To summarize, despite the large progress in studying the gas kinematics in the first billion years of the Universe (Sharda et al. 2019; Neeleman et al. 2020; Rizzo et al. 2020, 2021; Fraternali et al. 2021; Jones et al. 2021; Lelli et al. 2021), it is still unclear whether primeval galaxies are dominated by rotation or dispersion in the early phase of their evolution. Most of the studies in the literature have focused on a limited number of targets that were preselected to be starburst (SFR > 100 M_⊙), massive (M_⋆ > 10^10.5 M_⊙), and luminous in [C II]. As these galaxies might not represent the bulk of the galaxy population at high redshift, we need to analyze a larger sample of main-sequence star-forming galaxies to assess the dynamical state of early galaxies. Here we present a kinematic study of 22 galaxies observed with ALMA at z > 4, which are expected to represent the bulk of the galaxy population in the first billion years of the Universe. For the galaxies that were observed with both the [C II] and [O III] lines, we also explore the possibility of a bias in the kinematics traced with one FIR tracer alone.

This paper is structured as follows. In Sect. 2 we describe the sample selection and data reduction processes. In Sect. 3 we present the data analysis process and the algorithm we adopted to recover the kinematic properties of the sample. In Sects. 4 and 5 we present our results, compare them with the other literature findings, and discuss the evolution of the velocity dispersion and V_max/σ_gas as a function of redshift and other physical parameters of galaxies. In Sect. 6 we summarize the findings of this work and draw our conclusions. We adopt the cosmological parameters from Planck Collaboration XIII (2016): H₀ = 67.8 km s⁻¹ Mpc⁻¹, Ω_m = 0.308, and Ω_Λ = 0.685.

2. Sample selection and observations

We selected our galaxy sample from the public ALMA data archive². We queried the database (on 12 January 2022) to obtain only star-forming galaxies at z > 4 with observations of either [C II] or [O III] or both emission lines in bands 5, 6, 7, or 8 and angular resolutions lower than 1.5″ (Indeed, an angular size of 1.5″at z ∼ 5 corresponds to ∼7 kpc, which is ∼3 − 4 times larger than the typical size of a galaxy at this redshift). We aimed to explore the kinematic properties of high-redshift star-forming galaxies and therefore excluded known submillimeter galaxies (SMGs) and quasars, which do not represent the bulk of the galaxy population (e.g., Weiß et al. 2013). We thus selected galaxies whose previous measurements of stellar masses and SFR are compatible main-sequence of star-forming galaxies (see details in Sect. 5.2). We also excluded from the selected targets galaxies with clear merger features. Finally, we removed observations in which the far-infrared lines are detected with a signal-to-noise ratio (S/N) < 7, which is not sufficient to perform a kinematic analysis on the data. In this work we define the S/N as the median value across the map of the ratio of the integrated flux emission of the line (moment 0 map) and the pixel-by-pixel estimated error σ^flux defined in Eq. (2). Our final sample is composed of 22 galaxies at 4.2 < z < 7.6, some of which were observed in the [C II] and [O III] emission lines (see Tables 1 and A.1).

We calibrated the visibilities with the Common Astronomy Software Application (CASA; McMullin et al. 2007) by using the pipeline scripts delivered with the data that are available in the archive. We used the appropriate CASA version for each target as stated in the scripts. We then performed the ‘cleaning process on the calibrated visibilities by using the CASA task tclean. To obtain the best sensitivity from the ALMA data sets, we adopted a natural weighting scale to generate the final datacubes, yielding an angular resolution between 0.2″ and 1.5″ depending on the data sets. The spatial pixel size of each datacube was set up to be between ∼1/5 to ∼1/10 of the minor axis of the beam, which allowed us to both sample the beam profile and reduce the number of correlated pixels in the datacube (Remijan et al. 2019). The frequency channel spacing of the datacube was set up to exploit the maximum spectral resolution of each data set (10–20 km s⁻¹).

3. Data analysis

In this section, we perform the kinematic analysis of the galaxies of our sample with the goal of determining their kinematic properties. We present two different methods for the kinematic fitting procedure that we adopted to determine the rotational velocity curves and velocity dispersion from the [C II] and [O II]data on the basis of the angular resolution and sensitivity of the observations.

3.1. Moment maps

Initially, we performed a pixel-by-pixel Gaussian fitting on the datacube to generate the moment maps: flux (zeroth moment), velocity (first moment), and velocity dispersion (square root of the second moment). We also accounted for an additional constant to match any potential residual of continuum emission.

Before analyzing the moment and performing the kinematic fitting to the maps, we removed the parts of the maps that were dominated by noise fluctuations. We estimated the conservative pixel-by-pixel error of the flux map as

(2)

where σ_k is the ALMA noise level at the spectral channel k determined in a free-target region of the datacube, v_min and v_max are defined as the 16th and 84th percentiles in spectral channels of our best-fitting Gaussian model, respectively, and Δv is the frequency channel spacing of the datacube in velocity units. We then excluded all the pixels in the kinematic maps with [C II] integrated flux emission lower than 5σ^flux. For [O III] observations, we adopted a lower threshold of 3σ^flux because the observed flux density of the oxygen line is usually fainter than that of the [C II] line (Carilli & Walter 2013; Inoue et al. 2016; Carniani et al. 2017; Marrone et al. 2018; Walter et al. 2018; Hashimoto et al. 2019; Harikane et al. 2020; Witstok et al. 2022).

3.2. Kinematic models

To reproduce the kinematic maps of the FIR lines, we used a rotational velocity model, V_rot(r), defined by

(3)

where r is the distance from the galaxy center, and V_d(r) and V_b(r) are the velocity profile due to the disk and bulge components, respectively. We focus on high-redshift galaxies and used observations whose sensitivity is insufficient to probe the rotational velocity at the very large radii from the center. We therefore neglected the contribution of dark matter to the circular motion and focused on the disk and bulge components, which are dominant within two scale radii from the galaxy center (Sofue 2013). Our assumption is also supported by the results from Genzel et al. (2017), Fraternali et al. (2021), and Lelli et al. (2021), who have found that the rotation curves of 0.6 < z < 5 galaxies do not show the classical flat rotation profile due to dark matter, but the gas kinematics is well modeled with a baryon-dominated disk with a negligible mass of the dark matter in the inner regions.

Although most of the high-redshift studies (Neeleman et al. 2020; Fujimoto et al. 2020; Genzel et al. 2017; Übler et al. 2019; Gnerucci et al. 2011) show that the gas kinematics of galaxies at z = 1 − 4 can be modeled only with an exponential gas disk neglecting the stellar bulge component due to the increasing gas fraction (f_gas = M_gas/(M_gas + M_⋆)) with increasing redshift (Tacconi et al. 2020), recent [C II] observations (Rizzo et al. 2020; Lelli et al. 2021) have revealed that massive (M_⋆ > 10¹⁰ M_⊙) z > 4 galaxies might require a bulge component in addition to the exponential disk to reproduce their kinematics properties. In this work, we therefore fit the kinematics maps of massive galaxies with two models: one model that includes both components, and the other model with the exponential disk alone. Then we select the best-fitting model that better reproduces the moment maps by minimizing the χ².

In cases in which the sensitivity and resolution are not high enough to determine the presence of the bulge or M_⋆ < 10¹⁰ M_⊙, we adopt the simple exponential disk model.

For the exponential disk, we assumed that the light is emitted by a thin disk (Freeman 1970; Förster Schreiber et al. 2006; Binney & Tremaine 2008) whose surface brightness radial profile is

(4)

where r_d is the disk scale length. Assuming a constant mass-to-light ratio, we write the profile of the surface mass density as

(5)

The velocity contribution due to the exponential disk component can therefore be written as a function of the radius as

(6)

where , while ℐ₀, ℐ₁, 𝒦₀, and 𝒦₁ are modified Bessel functions. As this works focuses on galaxies at z > 4, which have a typical disk scale length of < 2 kpc (Shibuya et al. 2015), we defined the mass of the disk, M_d, as the mass computed within a radius R₀ = 5 kpc, which corresponds to 2–3 disk scale lengths: M_d = M_d(R₀) = M_d(r = 5 kpc).

On the other hand, the light profile of the bulge component is assumed to match a Sérsic profile (Sérsic 1963),

(7)

where I₀ is the surface brightness at the center of the bulge, r_e is the effective radius, and n is the Sérsic index. is a function of n (Prugniel & Simien 1997).

We can then describe the contribution of the bulge component to the circular velocity assuming a constant mass-to-light ratio as

(8)

(Terzić & Graham 2005), where M_b is the total stellar mass of the bulge, p is a function of the Sérsic index (p = 1 − 0.6097/n + 0.05563/n²; Lima Neto et al. 1999), and γ and Γ are the incomplete and complete gamma function, respectively. In our analysis, we opted for a de Vaucouleurs profile by assuming n = 4 for the bulge model (de Vaucouleurs 1958). Because of the sensitivity and angular resolution of our observations, different values of the Sérsic index do not change the results of the velocity dispersion. Future studies based on high-sensitivity and high angular resolution observations should leave n as a parameter in the kinematics fitting process.

A caveat of the present work is that we only considered rotating disk models. An alternative scenario would be that the observed velocity gradients are a consequence of two spatially unresolved sources at different systematic velocities that are located at a distance smaller than the half-beam size (i.e., < 2 kpc). Rizzo et al. (2022) showed that to perform a robust analysis and to distinguish disks from mergers in the early Universe, data with a resolution ∼0.2″ and a signal-to-noise > 10 are required. This means that the resolution of the ALMA observations is not sufficient to distinguish the two scenarios, and future high-angular resolutions will be fundamental to further this type of studies.

3.3. Beam-smearing effect

The limited angular resolution strongly affects the maps of flux, velocity, and velocity dispersion, especially in high-redshift observations, where the angular resolution is comparable to the galaxy size. From an analytical point of view, the beam-smearing effect introduces a convolution product in the equations of the moment maps. If the beam (or point spread function) of the telescope is greater than the characteristic size of the galaxies where its velocity profile changes steeply, the finite beam size smooths the velocity gradient in the velocity map and increases the velocity dispersion, especially in the central region, where the observed σ_gas is greater than the intrinsic velocity dispersion of the gas through the broadening of the emitting line. The broadening of the line can be mistaken as the result of random motions of the gas, yielding an overestimation of the gas velocity (Di Teodoro & Fraternali 2015; Kohandel et al. 2020; Concas et al. 2022). The net result of beam smearing might underestimate the rotation velocity and overestimate the velocity dispersion.

The ratio V_max/σ_gas is then highly impacted by this effect and can be underestimated if no correction is applied. We thus take into account how the beam-smearing impacts the observations in the fitting process. In Appendix C we test in detail the beam-smearing effects on datacubes and our ability to recover the galaxy kinematics, even for low angular resolutions. We find that V_max/σ_gas is underestimated if the beam size is more than three times the galaxy size.

In addition to these effects, the beam smearing introduces a strong correlation between neighboring spatial pixels in the datacubes, and neglecting the pixel-pixel noise correlation might lead to an underestimation of the uncertainties of the kinematic properties. In the next sections and in Appendix B, we describe how we corrected for these effects.

3.4. 3D datacube model

We used the Python routine KINMS (Davis et al. 2013) to create the datacube models for the fitting process. KINMS initializes 5 × 10⁵ line-emitting particles given by the user. Each particle is assigned a position in the galaxy space (x, y, z) to reproduce the sky brightness distribution given as input by the user (e.g., the exponential profile in Eq. (4)). Each particle is also given a velocity (v_x, v_y, v_z) following the assigned velocity profile (e.g., rotational velocity in Eq. (3)) and a random velocity component σ that represents the intrinsic velocity dispersion of the gas.

The routine then projects the 6D space-velocity array of each line-emitting particle onto the sky plane based on a given inclination and rotated according to a position angle. After the rotation, every particle is assigned new coordinates onto the sky plane (x_s, y_s), and the velocity along the line of sight v_los. Finally, KINMS creates a datacube: Each particle with a given position and projected velocity along the line of sight is assigned to a pixel (i, j, k) in the datacube. The size of the spatial pixels is rebinned to match the size of the spatial pixels of the observation. Each spectral slice of the datacube is convoluted with a point spread function (PSF) as large as the ALMA beam of the data sets analyzed in this work. The result is a datacube with the same spatial and spectral resolution as the observations. The moment maps derived from the mock datacube have the same effect of the beam smearing as the real data. In this way, the fitting process already takes into account the effect of beam smearing without the need for further corrections.

3.5. Method I: 2D kinematic fitting

After we generated the moment maps for each target, we performed a kinematical fitting by using the mock datacube as the model to fit the data (see details in Appendix B). The parameter space was explored using a Bayesian approach based on the Monte Carlo Markov chains (MCMC) method, which allowed us to infer the posterior probability distribution of the free parameters. In particular, we used the package EMCEE (Foreman-Mackey et al. 2013).

The free parameters of our kinematic models were the scale radius of the exponential disk (r_d), the inclination (inc) and position angle of the disk (PA), the coordinates of the center (x₀, y₀), the rotational mass (M_rot) that is the mass of the disk (or disk + bulge) required to reproduce the velocity maps without taking into account the contributions from the turbulence support, the systemic velocity of the galaxy (v_sys), and the intrinsic velocity dispersion (σ_gas). For the galaxies that were fit with both bulge and disk components, we did not add additional free parameters, but we assumed that the bulge effective radius is equal to half of the exponential disk scale radius and that the bulge mass is as high as the stellar mass, which has been reported in the literature for our galaxies (see Table A.1).

The three main steps of the algorithm we adopted for the fitting procedure are shown in Fig. 1. We initially fit the flux map and determined the center and size of the galaxy. These best-fitting parameters were then kept fixed in the next steps. In the second step, we fit the velocity map and inferred the best-fitting results for the v_sys, inc, PA, and M_rot, which were used as priors for the last steps, in which all moment maps were fit simultaneously. As the neighbor spatial pixels of the moment maps are correlated due to the beam-smearing effect (see Sect. 3.3), performing the standard kinematical fitting on the maps might underestimate the parameter uncertainties because in the computation of the likelihood, all pixel are considered as independent values. Therefore, we opted to estimate the likelihood by using only independent pixels in the moment maps that were randomly extracted from the maps at a distance greater than the semi-major axis of the beam. This procedure was performed 100 times for each target to ensure that the convergence of the free parameters did not depend on the pixel selection. We then combined the walkers for every trial (after a 50% burn-in phase) to derive the best parameters (median values) and error estimates (16th and 84th percentile) of the free parameter distribution for each target. We discuss this method in more detail in Appendix B. An example of the best-fitting results from method I is shown in Fig. 2.

Fig. 1. Summary of the method I fitting procedure to derive the model parameters.

Fig. 2. Best-fitting results of HZ4 obtained with method I. The left panels show the moment [C II] maps: the normalized flux map (top), the velocity map (middle), and the velocity dispersion (bottom). From left to right, the panels illustrate the data, models, and residuals. The color bars of the residual range between −3σ and +3σ, and the black lines indicate ±1σ. The ALMA beam is shown as the gray ellipse in the flux map. The multidimensional parameter space explored by step 3 of the method I algorithm is shown in the corner plot on the right.

3.6. Method II: 1D kinematic fitting

Some of ALMA observations of our sample do not show a clear velocity gradient in the kinematic maps. The reasons why we could not see a clear velocity gradient could be: a low angular resolution of the observations, the low sensitivity, or the disk not being dominated by rotation, but being dominated by random motions. In these cases, the fitting process does not converge and thus returns a flat posterior distribution for the free parameters. To overcome this problem and gain some information regarding the kinematics of these galaxies, we determined the kinematic properties from the spatially integrated spectrum. The spatially integrated spectrum is mainly used in radio observations of the HI emission line because it can show the kinematical properties of spatially unresolved galaxies (Roberts 1978; Yu et al. 2020; Stewart et al. 2014). The basic idea is to determine the velocity dispersion by comparing the integrated spectrum of the observations with the integrated spectrum of the model datacube. For this purpose, we need to limit the number of free parameters. In particular, we need to estimate the radius and mass of the galaxy a priori because these parameters have a strong impact on the shape of the line profile (Roberts 1978; de Blok & Walter 2014) and are very degenerate with the other disk inclination and intrinsic velocity dispersion parameters (Kohandel et al. 2019).

We initially fit the flux map (step 1 of method I) and inferred both the scale length radius and a first constraint on the disk inclination. We then derive the mass as

(9)

where M_* is the stellar mass, based on observations reported in the literature (Capak et al. 2015; Willott et al. 2015; Matthee et al. 2017; Smit et al. 2018; Harikane et al. 2020; Herrera-Camus et al. 2021; Fudamoto et al. 2021; Bouwens et al. 2022; Dayal et al. 2022; Schouws et al. 2022), and M_gas is the gas mass estimated by using the empirical relation by Zanella et al. (2018),

(10)

with α_[CII] = 30 M_⊙/L_⊙. The uncertainty on the conversion factor between [C II] luminosity and gas mass is about 0.2 dex.

After the scale radius and mass of the galaxy were estimated, the only free parameters of our datacube model are the disk inclination inc and the intrinsic velocity dispersion σ_gas.

In the fitting procedure, the model spectrum, F^datacube, was extracted in the model datacube from the same region defined in the observations. In addition to the disk inclination and σ_gas, we also added two other free parameters to match the normalization of the spectrum (a) and any offset due to the continuum emission (b). In conclusion, the spectrum of each galaxy was fit with a spectrum S^model defined as

(11)

We fit the data spectrum by using an MCMC procedure and adopting a flat prior distribution for the velocity dispersion in the range 0 km s⁻¹ and 300 km s⁻¹. Since the inclination plays a major role in the global profile shape and the subsequent determination of velocity dispersion of the galaxy, we left as parameter the inclination assuming flat priors between ±3σ of the best-fitting inclination estimated from the flux map. Figure 3 illustrates an example of the best-fitting results from method II.

Fig. 3. Best-fitting results from method II for the target HZ9. The top left panels illustrate the observed flux map, model, and residuals. The ALMA beam is shown as the gray ellipse in the bottom left corner of the observed flux map. The bottom right panel shows the observed integrated spectrum and the best-fitting model. The multidimensional parameter space explored by the method II algorithm is shown in the corner plot on the right.

We note that by applying method II to the galaxies that were fit using method I, the best-fitting velocity dispersion is systematically higher by 1.6 times than that obtained from method I (Appendix D). Hence, we applied a correction factor of 1.6 to the results inferred from method II.

4. Results

In this section, we report the results from the kinematic fitting by adopting the rotating models described in Sect. 3.2. Table 1 summarizes the results provided by method I and II. In addition to Figs. 2 and 3, the best-fitting result and posterior distribution of free parameter for each target are presented in Appendices E and F. For the galaxies that were fit by an exponential disk alone and by the model composed of bulge and exponential components, we present the model parameters for the model with the smaller χ². Due to the poor angular resolution of some observations, we were not able to constrain the disk inclination, dynamical mass, and velocity dispersion of all galaxies. In these cases, we inferred upper limits to the kinematics parameters.

The values of the velocity dispersion obtained for galaxies at z > 5 that were studied with method I range between 37 km s⁻¹ and 132 km s⁻¹, with a median velocity dispersion of km s⁻¹. For the rotational support, we find a median value of V_max/σ_gas = 2 that is compatible with turbulent but still rotationally supported disk galaxies.

Galaxies analyzed with method II show a velocity dispersion, corrected by a factor 1.6 (see Sect. 3.6), between 37 and 126 km s⁻¹ and a median value of km s⁻¹, which is compatible with the median value obtained with the other method. The median value of V_max/σ_gas is 3.2, which is slightly higher than the result found with method I. This can be due to the uncertainties on the mass inferred from the luminosity of the [C II] emission and the rest-frame UV observations, to the uncertainty of the scale radius of the exponential disk obtained with the fitting of the flux map, or to the uncertainties on the inclination estimates, which can broaden the line.

Interestingly, we note that for the galaxies that are observed in the two FIR lines, the velocity dispersion of [C II] and [O III] are similar, but they trace different gas phases of the interstellar medium. In only one out of four galaxies is the [O III] velocity twice higher than the velocity dispersion obtained by fitting the [C II] map. The [O III]and [C II] observations in this sample predict similar values for the velocity dispersion of the galaxies. Simulations and observations showed that different gas tracers give rise to different velocity dispersions (Kretschmer et al. 2022; Ejdetjärn et al. 2022).

The rotation-to-random motion ratios estimated in our sample are lower on average than the V_max/σ_gas ratios observed in the local galaxies and in the dusty massive galaxies at z > 4 reported in the literature. To exclude any possible bias due to the different fitting algorithms, we also analyzed four high-redshift galaxies that were studied by Lelli et al. (2021), Neeleman et al. (2020) Jones et al. (2021), Herrera-Camus et al. (2022), and Posses et al. (2023), and we compare their results with ours.

Lelli et al. (2021) studied the kinematic properties of the massive dusty galaxy ALESS 073.1 using the software ^3DBAROLO (Di Teodoro et al. 2016), which directly takes the beam-smearing effect into account on the datacube and uses a 3D fitting algorithm with tilted-ring models to infer the kinematical properties of galaxies. Our best-fitting results based on the moment maps agrees with the results found by Lelli et al. (2021), who reported low values of σ_gas ranging from ∼60 km s⁻¹ in the internal regions to ∼10 km s⁻¹ in the outer region. Our results in Fig. E.13 show a bimodal distribution of the velocity dispersion that converges at ∼30 km s⁻¹ when most of the selected pixels were in the outer regions, and to ∼56 km s⁻¹ when the central pixels were selected. Despite the two different methods of fitting, our results agree with the kinematical cold disk with values of V_max/σ_gas comparable to those obtained for local galaxies. De Breuck et al. (2014) analyzed the same galaxy using a 2D fitting algorithm of data with a poorer resolution and found values of V_max/σ_gas ∼ 3.1, but this is mostly due to the underestimated rotational velocity (∼120 km s⁻¹) and not to the overestimated velocity dispersion (σ_gas = 40 ± 10 km s⁻¹).

Another galaxy that was observed with ALMA and was recently analyzed by Neeleman et al. (2020) and Jones et al. (2021) is DLA0817g. Neeleman et al. (2020) exploited the higher-resolution observations (∼0.2″ as used in this work) and used an approach similar to ours, in which the flux was fit assuming a rotating thin disk with an exponential brightness profile. They compared it with the results from ^3DBAROLO. Jones et al. (2021) studied the kinematics by exploiting lower-resolution data (∼1″) using ^3DBAROLO. Their results agree with those by Neeleman et al. (2020), who they found km s⁻¹ and . We note that we obtain a higher value for the maximum rotational velocity because we also considered the bulge component in our fitting because this model better reproduced the kinematic maps. When we use the exponential disk alone, we obtain a , while the velocity dispersion that we obtain is the same using the bulge+exponential disk or the exponential disk model. Our results disagree with the results found by Jones et al. (2021), who reported a velocity dispersion that is lower by a factor of two.

Herrera-Camus et al. (2022) presented the kinematical study of the target HZ4 by modeling the galaxy as a turbulent disk plus a dark matter halo. They found values of the velocity dispersion of km s⁻¹ and a value of V_max/σ_gas = 2.2. These results are comparable within the uncertainties to our best-fitting values.

Finally, the target COS-29 was analyzed by Posses et al. (2023) using the algorithm ^3DBAROLO. They found that the average velocity dispersion has an upper limit of 30 km s⁻¹ and a ratio of the velocity and the velocity dispersion greater than 1.4, which is comparable to the values found in this work.

Overall, we find that the parameters inferred from our algorithms agree with the results found by other high-redshift studies, and we found no evidence of a bias in the recovery of the velocity dispersion. In the next sections, we therefore compare our results to those of other studies from z = 0 to z = 8 to understand the evolution of turbulence with redshift and galaxy properties better.

5. Discussion

5.1. Velocity dispersion and V_max/σ_gas evolution with redshift

The left panel of Fig. 4 shows the evolution of the gas velocity dispersion as a function of the redshift from z = 0 to z = 8. In the figure we show previous results from the literature, and the best-fitting estimates from our work. Differently from the conclusion reported by Lelli et al. (2021), Rizzo et al. (2020, 2021), and Fraternali et al. (2021), who reported that high-redshift (z > 4) galaxies have low values of the velocity dispersion (i.e., σ_gas ∼ 30 − 60 km s ⁻¹; blue region at 4 < z < 4.5 in Fig. 4) with respect to theoretical models, our sample at z ∼ 5 − 7 seems to be consistent with the prediction by Wisnioski et al. (2015) of the velocity dispersion evolution with redshift,

(12)

Fig. 4. Evolution of the velocity dispersion (left) and the Vmax/σgas ratio (right) with redshift. Blue and purple marks are the values obtained from our sample by analyzing the [C II] and [O III] observations, respectively. Diamonds are obtained by fitting the three moment maps (method I), and circles represent the results found for the marginally resolved galaxies (method II). Previous σ and Vrot/σ estimates reported in the literature are represented as color-coded rectangles: the observations with the [C II] emission line at 4 < z < 5 by Fraternali et al. (2021), Rizzo et al. (2020, 2021), Lelli et al. (2021), Sharda et al. (2019), and Neeleman et al. (2020) at 4 < z < 6 from Jones et al. (2021) are plotted in blue. Red rectangles show the results of the studies of ionized gas (Hα, [O III] λ 5007) in z < 4 galaxies (Epinat et al. 2010; Green et al. 2014; Turner et al. 2017; Übler et al. 2019). Green rectangles represent the measurements derived from molecular gas (Mogotsi et al. 2016; Girard et al. 2021), (CO). The pink symbols were derived with traced of neutral gas (H I) at z = 0 (Mogotsi et al. 2016). The boxes extend from the 25th to the 75th percentile of the data in the samples, and solid black lines show the median of the value. The light blue area represents the predictions by Wisnioski et al. (2015). The upper and lower boundaries of the shaded light blue region represent the predicted Vmax/σ (see Eq. (13)) for Q = 0.67 and Q = 2 and the solid light blue line represents the Vmax/σgas values predicted for Q = 1. The dashed black line represents the evolution of the best-fitting velocity dispersion trend by Übler et al. (2019). The gray and orange areas show the results obtained by the TNG50 (Pillepich et al. 2019) and SERRA simulations (Kohandel et al. in prep.), respectively.

The trend of Eq. (12) with Q = 1 and 100 km s⁻¹ < V_max < 250 km s⁻¹ is shown with the light blue shaded curve in the left panel of Fig. 4. As the gas fraction increases with redshift (Tacconi et al. 2020), the velocity dispersion also increases, reaching values of ∼40 − 100 km s⁻¹ at z ∼ 6 − 8, independently of the tracer that is used to map the gas kinematics, which seems to be confirmed by our results based on two different FIR lines, [C II] and [O III].

The evolution of the V_max/σ_gas ratio as a function of the redshift from z = 0 to z = 8 is shown in the left panel of Fig. 4. We also superpose the expected decrease in V_max/σ_gas with redshift based on the Toomre instability criterion (Wisnioski et al. 2015; Genzel et al. 2011),

(13)

In the figure we also report the results obtained from the TNG50 simulations (Pillepich et al. 2019), and from the SERRA simulations (Pallottini et al. 2022, Kohandel et al. in prep.) for the kinematic studied with the [C II] emission line. While our V_max/σ_gas ratios agreed with the estimates found for galaxies in the ALPINE survey (Jones et al. 2021), our analysis differs from the results determined in the sample of galaxies studied by Sharda et al. (2019), Rizzo et al. (2020, 2021), Fraternali et al. (2021) and Lelli et al. (2021), hereafter the dusty star-forming galaxy (DSFG) sample (blue box at 4 < z < 5). The DSFG sample indeed shows a level of rotational support similar to what is observed in the local Universe (i.e., V_max/σ ∼ 10). In contrast, the sample targeted in this work and in the ALPINE survey has values of V_max/σ_gas between 1 and 4 that are consistent with the values expected by Wisnioski et al. (2015) and the trend in the TNG50 simulations (Pillepich et al. 2019) and is slightly lower than what was observed in the SERRA simulations (V_max/σ_gas ∼ 8 − 4, Kohandel et al. in prep.)

The discrepancy in the velocity dispersion and in V_max/σ_gas between our sample and the DSFGs suggests that the gas kinematic properties do not only depend on the redshift alone, but also on the properties of the galaxies themselves. In the next sections, we therefore analyze the properties of the galaxies and verify whether the gas kinematics depend on these properties.

5.2. Sample properties

To understand whether the physical properties of the galaxies can impact the galaxy kinematics, we analyzed the difference between the three samples (i.e, ALPINE, DSFG, and our sample) in the SFR-stellar mass diagram and dust content. The values of stellar mass, SFR, and dust mass were obtained from the literature (see Table A.1 for details).

We analyzed the SFR-stellar mass diagram (Fig. 5), in which our galaxies, the ALPINE, and the DSFG sample are indicated with diamonds, triangles, and squares, respectively. We also plot the relation for main-sequence galaxies as parameterized by Schreiber et al. (2015),

(14)

Fig. 5. SFR as a function of the stellar mass for our sample of galaxies and the results found in the literature. The colors represent the redshift on the galaxy (see the color bar). The diamonds are the galaxies targeted in this work. The squares show the results found by Sharda et al. (2019), Rizzo et al. (2020, 2021), Lelli et al. (2021), and Fraternali et al. (2021). The triangles show the ALPINE sample (Jones et al. 2021). In green we plot the starburst region at z ∼ 4 − 5 as predicted by Caputi et al. (2017). In light blue we plot the main-sequence region extrapolated between z ∼ 4 − 8 by Schreiber et al. (2015).

extrapolated for 4 < z < 8³, and the empirical relation for starburst galaxies at z ∼ 4 − 5 by Caputi et al. (2017),

(15)

The ALPINE sample and our sample lie above the z = 5 main-sequence relation, and the DSFGs are about 0.5 dex more massive and more consistent with the high-redshift starburst populations. Although at low redshift, starburst galaxies tend to higher velocity dispersion values than main-sequence galaxies (Übler et al. 2019; Perna et al. 2022), our analysis reveals the opposite scenario in the distant Universe. The DSFG sample seems to be less turbulent than the high-z main-sequence galaxies, suggesting that the gas is settled in a dynamically cold rotating disk despite the high star formation activity. The reason for the discrepancy in the velocity dispersion and V_max/σ between the two samples may therefore be ascribed to the different galaxy population properties of the selected galaxies.

In Fig. 6 we report the distributions for the three samples of log(M_dust/M_*) and log(M_dust/M_tot), where M_tot is the total baryonic mass computed as M_gas + M_*, and M_dust is the dust mass of the galaxy. The values of M_dust were obtained from Pozzi et al. (2021) and Béthermin et al. (2020) for the ALPINE sample, and from Aravena et al. (2016) for the DSFG. see Table A.1 for the dust mass of the sample targeted in this work. For some galaxies, the dust temperature was assumed to be a fixed value. We therefore note that the correlation between the dust temperature and the dust mass may have overestimated some of the dust mass estimates.

Fig. 6. Dust mass content in the three galaxy samples targeted in this work. In the upper panel, we plot the distribution of the ratio of the dust mass and stellar mass for the sample targeted in this work (in blue) and the sample of DSFG targeted by Rizzo et al. (2020, 2021), the starburst dusty galaxies (Sharda et al. 2019; Fraternali et al. 2021; in red), and the ALPINE sample (Jones et al. 2021; in green). In the lower panel we present the distribution of the ratio of the dust mass and total baryonic mass (Mgas + M*) for the sample targeted in this work, the sample of DSFG targeted by other high-z studies, and the ALPINE sample.

The figure shows that while our sample and the ALPINE samples have a star-to-dust mass ratio of 10⁻³ on average, the DSFGs are characterized by M_dust/M_⋆ > 10⁻², similar to what was observed in other high-z DSFGs and in submillimeter galaxies (Casey et al. 2014).

The low amount of dust and the high velocity dispersion in our sample with respect to DSFGs suggests that the mechanisms that cause the ISM to become turbulent might also have driven shocks that affect the dust distribution.

5.3. V_max/σ_gas evolution with galaxy properties

By combining the three samples studied at z > 4, we can also assess whether the V_max/σ_gas ratio of our galaxies is correlated with stellar mass, dust mass, or gas fraction.

Lower-redshift studies have shown a dependence on the V_max/σ_gas and the stellar mass (Kassin et al. 2012; Wisnioski et al. 2019; Förster Schreiber & Wuyts 2020), where more massive galaxies have higher values of V_max/σ_gas, while evidence of rotating but turbulent disks are common in less massive galaxies. The cosmological simulation TNG50 (Pillepich et al. 2019) shows a slight but similar dependence on the stellar mass of the galaxy with greater V_max/σ with increasing stellar mass at all redshifts, while SERRA simulations predict a stronger dependence on the stellar mass, where V_max/σ at a fixed stellar mass remains almost constant for 6 < z < 8.

Figure 7 shows the evolution of V_max/σ_gas with redshift, color-coded with stellar mass bins. The dependence of V_max/σ_gas on stellar mass is still valid at z > 5, where galaxies with M_⋆ > 10¹⁰ have V_max/σ_gas > 3, while M_⋆ < 10¹⁰ are characterized by lower values (V_max/σ_gas ∼ 2). As DSFGs and our sample cover two different ranges of stellar masses, the discrepancy in V_max/σ_gas values between the two populations might be due to a bias in the sample selection.

Fig. 7. Vmax/σgas as a function of redshift and stellar mass. Squares represent the results from KMOS3D (Wisnioski et al. 2019). Circles show the Vmax/σgas for the whole z > 4 sample including the results of this work. Solid lines show the results obtained by TNG50 simulations (Pillepich et al. 2019). Dashed lines show the results obtained by SERRA simulations (Kohandel et. al. in prep.). Different stellar mass bins are represented with different colors.

Analyzing the parameter space V_max/σ_gas function of the ratio M_dust/M_⋆, we note a clear bimodality distribution between dusty galaxies and galaxies with low dust content. Figure 8 shows that the kinematics of gas in high-z galaxies is slightly related to the presence of dust in the interstellar medium. The ISM of high-z dusty galaxies seems to be less affected by star formation feedback and gravitational torques. This might imply that the processes that increase the gas velocity dispersion in the galaxies might also create shocks in the gas that destroy the dust.

Fig. 8. Vmax/σgas as a function of redshift, and the ratio of the dust mass and stellar mass. In red we show the DSFG sample (Sharda et al. 2019; Rizzo et al. 2020, 2021; Lelli et al. 2021; Fraternali et al. 2021), in green we plot the ALPINE sample (Jones et al. 2021), and in blue we show the sample targeted in this work.

Finally, we explored whether the V_max/σ ratio is driven by the gas fraction, as expected for a disk that is in equilibrium between gas heating and cooling, and assuming a Toomre parameter Q = 1 (Eq. (1); see also Wisnioski et al. 2015). In Fig. 9 we report the V_max/σ_gas as a function of gas fraction for z > 4 galaxies. Our data (blue marks) show a flat distribution between f_gas = 0.3 and f_gas = 1.0, which is marginally consistent with the relation by Wisnioski et al. (2015). However, taking the distribution of V_max/σ_gas values of all high-z galaxies into account, we note that the dependence of V_max/σ_gas on f_gas drops. Several galaxies show similar f_gas but quite different V_max/σ_gas. This suggests that the state (i.e., V_max/σ_gas) of the rotating disk does not depend on the gas fraction alone.

Fig. 9. Vmax/σgas as a function of gas fraction. The black line represents the model by Wisnioski et al. (2015) with the assumption of Q = 1 and . In dark blue, we show the data obtained from the sample of galaxies targeted in this work.

Another possible explanation for the differences between the DSFG sample, this work, and the ALPINE sample is suggested by the theoretical work by Kretschmer et al. (2022). The zoom-in simulations developed by the authors show that when a galaxy is in a phase of constructive gas accretion (i.e., the cold gas streams are coplanar and corotate with the disk), the galaxy disk is then primarily supported by rotation with V_max/σ_gas ∼ 5. This phase only lasts for about five rotational periods (∼400 Myr), then the disk can be disrupted by merging episodes. Therefore, we might see a difference because of the different merging rates for high- and low-mass galaxies. Massive galaxies undergo fewer major mergers than less massive galaxies (Dekel et al. 2020), which means that the ordered disk can survive longer. When this phase of low-velocity dispersion has a limited duration, however, this effect cannot explain why all of the galaxies in our reference sample are in this state while the majority of our galaxies are not.

5.4. Drivers of the velocity dispersion

We analyzed the dependence of the velocity dispersion as a function of SFR. The velocity dispersion is expected to increase with increasing SFR due to two main mechanisms: feedback, and gravitational instabilities.

In Fig. 10 we report the velocity dispersion as a function of the SFR of the data as well as the empirical models derived by Krumholz et al. (2018), assuming the high-z fiducial values for the model (see Table 3 from Krumholz et al. 2018). These analytical models are based on the vertical hydrostatic equilibrium of the disk, on the marginal gravitational stability (Q ∼ 1), and on the balance between the energy injected into the ISM by the star formation feedback and the gas transport from the outer to the inner region of the disk, and the dissipation of the energy from the velocity dispersion.

Fig. 10. Dependence of the velocity dispersion on the SFR and the rotational velocity. In the left panel, we show the model by Krumholz et al. (2018) with the feedback as the only driver of the velocity dispersion. In the central panel, we show the model with the release of gravitational energy due to the transport of the gas across the disk as the only driver of the velocity dispersion. In the right panel, we show the model with the combination of feedback and transport as drivers of the velocity dispersion.

The left panel shows the prediction from the model called “Feedback” here, which is based on the model that Krumholz et al. (2018) called “No-transport fixed Q”. In this scenario, the authors assumed that galaxy-scale star-formation activity helps to maintain the Toomre parameter Q close to the value of 1. When the Toomre parameter drops below 1, star-forming clumps begin to form, which increases the SFR and therefore the feedback effects. Since the Toomre parameter is proportional to the velocity dispersion, Q increases due to the turbulence induced by supernova explosions (i.e., feedback). When Q > 1, the star formation is instead not triggered, and when no other drivers of the velocity dispersion exist, the energy is dissipated, and the Toomre parameter decreases.

The central panel depicts the model “No-Feedback”. In this model, the gravitational energy released by the mass transport within the disk is the only driver of the velocity dispersion. The assumption for this model is that the Toomre parameter Q is always equal to 1, and all the gravitational energy released is dissipated by turbulence.

The right panel describes the dependence of the velocity dispersion as a function of the SFR when the velocity dispersion is driven by both the star formation feedback and the mass transport within the disk (i.e., galaxy-scale gravitational instabilities). At high SFR, the trend is similar to that of the model “No-Feedback”, while at low SFRs, there is a flat plateau in the velocity-dispersion-SFR trend that is based on the assumption that the star formation efficiency is constant for the entire disk, while the Toomre parameter is left free to vary. These curves match the trend observed in local star-forming galaxies with a low SFR, which show a plateau in the velocity dispersion values independently of their SFR (Epinat et al. 2010; Ianjamasimanana et al. 2012; Mogotsi et al. 2016).

The gas velocity dispersion of our sample of high-redshift galaxies cannot be totally explained with the model in which only the feedback from star formation injects energy into the interstellar medium (left panel). To obtain high values of the velocity dispersion by using only the feedback effects, our galaxies should have an SFR higher by one order of magnitude than observed. Finally, models that include the gravitational processes for galaxy-scale turbulence driving (central and right panels) agree better with our data. On the other hand, the DSFG sample shows a velocity dispersion that agrees with the feedback-only driven turbulence (Rizzo et al. 2021). Unfortunately, our high-redshift observations are not able to explore the region at low SFR, so that we are not able to distinguish between the feedback+transport and the transport-only models.

In Fig. 10 we study the dependence of σ_gas on SFR and rotational velocity. However, the Krumholz et al. (2018) model expects that the gas velocity dispersion driven by the mass transport also depends on the gas fraction. In particular, for the “No Feedback” model, the dependence is

(16)

where f_gas, p is the gas fraction in the midplane of the galaxy (f_gas, p = 1.5f_gas; Übler et al. 2019), Q is the Toomre parameter, and v and σ are the rotational velocity and the velocity dispersion, respectively. With this dependence, galaxies at fixed σ_gas and SFR can have a low rotational velocity but higher gas fraction, or conversely, a low gas fraction and high rotational velocity. To take the different behavior at different gas fractions into account, we plot in Fig. 11SFR/f_gas, p as a function of σ_gasv². Despite the large uncertainties and low statistics, our sample scatters around Q = 1 compatible with the lower redshift study (Übler et al. 2019), while the majority of the DSFG sample does not follow this trend because gravitational instabilities are not the dominant driver of turbulence in these galaxies. As even at lower redshift, the galaxies lie on the same relation, we conclude that main-sequence galaxies self-regulate over the cosmic time (from z = 8 to z = 0) and evolve along roughly the lines of constant Q.

Fig. 11. Relation of the kinematical properties and SF-related properties for our sample of galaxies following Eq. (16). The dashed line corresponds to the model with Q = 1, the dotted line shows the model with Q = 0.1, and the dash-dotted line shows the model assuming Q = 10

Figures 10 and 11 show that our sample of galaxies is consistent with the models by Krumholz et al. (2018), who assumed a velocity dispersion that is mainly driven by gravitational instabilities. The lower-redshift observations by Übler et al. (2019) and Girard et al. (2021) also agree with the model of feedback+transport. Rizzo et al. (2021) instead found that the velocity dispersion of their galaxy can be explained only by the energy injected by the stellar feedback due to their higher SFR. This scenario is also consistent with the fact that the majority of the galaxies analyzed by Rizzo et al. (2021) are in a starburst phase (see Fig. 5).

6. Conclusions

Based on the kinematics analysis of the selected 22 galaxies observed with ALMA and studied through the [C II] and [O III] emission lines, we found:

• The median velocity dispersion of the gas in the sample is ∼65 km s⁻¹, even though we have large uncertainties due to the low signal-to-noise ratio and the low resolution of the observations. This value is about four to five times higher than what is observed in local galaxies (z ∼ 0) and two to three times higher than the velocity dispersion measured in galaxies at z ∼ 1 − 2. We therefore conclude that the velocity dispersion increases with redshift, which is consistent with the predictions of Wisnioski et al. (2015) and the cosmological simulations of Pillepich et al. (2019) and Kohandel et al. in prep.

• The 2D kinematic analysis does not show a large bias in the recovery of the velocity dispersion and the value of V_max/σ_gas due to beam smearing in comparison to the 3D fitting algorithm ^3DBAROLO.

• The velocity dispersion computed with both the gas tracer [C II], which traces the multiphase medium, and [O III], which only traces the ionized gas, does not show a significant difference for three out of four galaxies, while a galaxy shows a velocity dispersion traced with ionized gas greater than a factor of 2, as it does for lower-redshift galaxies (Girard et al. 2021).

• The high-velocity dispersion of our sample of normal star-forming galaxies when compared with the Krumholz et al. (2018) models can be explained by the release of gravitational energy from the transport of mass across the disk. This high-velocity dispersion cannot be sustained by similar models in which star formation feedback (i.e., supernova explosions) is the only source of the velocity dispersion. This is supported by the Kohandel et al. (2020) simulation, in which they showed that 90% of the observed [C II] velocity dispersion is powered by bulk motions due to gravitational processes, such as merging and accretion events.

• The ratio V_max/σ_gas does not show a clear dependence on the gas fraction, as was assumed by Wisnioski et al. (2015). However, we also found a dependence on the stellar mass and the dust mass. We therefore conclude that at fixed redshift, galaxies with higher stellar mass have a more regularly rotating disk () than those with lower stellar mass.

Current studies focusing on z > 5 galaxies show a difference between rotation-dominated disks and dispersion-dominated disks. Our analysis of a sample of normal star-forming galaxies at z > 4 suggests that high-redshift galaxies are more chaotic and less rotationally supported than their local counterparts. Our results contradict some previous studies (Rizzo et al. 2020, 2021; Lelli et al. 2021; Fraternali et al. 2021; Sharda et al. 2019) that analyzed z ∼ 4 massive galaxies with low values of the velocity dispersion and higher values of the ratio V_max/σ_gas comparable to local star-forming galaxies (see Sect. 5.1). However, we note that these massive (M_* = 10¹¹ M_⊙) galaxies are highly star-forming galaxies (SFR > 100 M_⊙ yr⁻¹), as shwon in Fig. 5, and they appear to be more evolved than the galaxies in our sample, with the evidence of an already formed stellar bulge just 1.2 Gyr after the Big Bang (Lelli et al. 2021) and a greater dust content.

Because of the difference in the physical properties (SFR, stellar mass, and dust content) of the samples targeted (see Sect. 5.2), we can explain the contrast between their kinematical behavior, with these samples being part of two different galaxy populations that both exist at high redshift. We showed in Sect. 5.3 that the gas fraction seems not the most important parameter for assessing the rotational support of a galaxy, as was shown by Wisnioski et al. (2015), and other physical quantities such as stellar mass and dust content may also affect the value of V_max/σ_gas.

The discrepancy between our results and those found in massive galaxies at similar redshift may also be explained by comparing the measured dispersion with the theoretical prediction by Krumholz et al. (2018) to investigate the drivers of velocity dispersion. As discussed in Sect. 5.4, while the velocity dispersion of massive galaxies in the literature can be explained with star formation feedback alone, the high-velocity dispersion of the galaxies targeted in this work requires the contribution of gravitational instabilities as well, according to the models considered.

The challenge of inferring the kinematical properties of high-redshift data is particularly influenced by the angular resolution of the observations and the sample selection of high-redshift galaxies. The selection of high-redshift galaxies is biased toward DSFGs or UV-bright galaxies. To understand and explain the evolution of the velocity dispersion and the rotational support of galaxies with redshift and explain the presence of dynamically cold but extremely star-forming galaxies and dynamically hot main-sequence galaxies, we need more high-resolution observations for high-redshift galaxies and a higher S/N.

Acknowledgments

This paper makes use of the following ALMA data: ADS/JAO.ALMA #2017.1.01052.S, #2017.1.01471.S, #2018.1.01359.S, #2017.1.00508.S, #2017.1.01451.S, #2018.1.01359.S, #2017.1.00604.S, #2018.1.00429.S, #2018.1.00085.S, #2019.1.01634.L. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), MOST 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. EP is grateful to F. Rizzo for useful comments and discussions. SC is supported by European Union’s HE ERC Starting Grant No. 101040227 – WINGS. MK and AP acknowledge support from the ERC Advanced Grant INTERSTELLAR H2020/740120.

References

Appendix A: Galaxy properties

Table A.1 lists the properties collected from the literature of the galaxies analyzed in this work.

Appendix B: Method I uncertainties and pixel correlation

Our fitting process, similar to others used in the literature (Neeleman et al. 2020; Pensabene et al. 2020; Neeleman et al. 2021), determined the best-fitting parameters by maximizing the likelihood function ℒ = exp(−0.5χ²), and χ² is defined as

(B.1)

where the D_ij and M_ij are the data and model moment maps, respectively, and is the error associated with the data maps. We note that neighboring pixels within a beam are strongly correlated and not independent because of the angular resolution of the data. This means that the assumption that individual pixels in the moment maps are independent can lead to an underestimation of the free parameter uncertainties.

To better understand the impact of the pixel-pixel correlation on the kinematic analysis, we simulated an ALMA mock observation of a rotating gas disk with the CASA task tool simobserve and generated four different datacubes with the task tclean with the same angular resolution, but different pixel sizes: 1/20, 1/10, 1/5, and 1/3 of the beam. Finally, we fit the 2D maps with the same algorithm as explained in Section 3.3 and left only the inclination as a free parameter, which is one of the parameters that is more affected by the beam smearing (Pensabene et al. 2020). Figure B.1 shows the uncertainties on the disk inclination parameter obtained by assuming all the pixels in the map as independent pixels as the red line. Although the datacubes were generated from the same visibilities dataset, the uncertainties on the free parameter depend on the pixel size that is adopted to create the final image.

Fig. B.1. Error on the inclination parameter for different datacubes with different pixel sizes. The errors derived from our method in which we estimated the likelihood by using only independent pixels is shown in blue, and the uncertainties obtained by assuming all pixel of a moment map as independent pixels are shown in red.

To better constrain the uncertainties on our free parameter and avoid pixel correlation, we estimated the likelihood by using only uncorrelated pixels. In each map, we randomly selected pixels that were separated from each other by a distance greater than the semi-major axis of the beam (see an example of the randomly selected pixel for the target COS-30 in figure B.2). We repeated this procedure 100 times to ensure that the likelihood estimate does not depend on the pixel selection. The blue line in Fig. B.1 shows the uncertainties obtained with this method for the four datacubes, and we note that the error does not depend on the pixel size. This confirms that the selected pixels are uncorrelated. This method therefore returns more conservative and physical uncertainties on the free parameters than the standard fitting procedure.

Fig. B.2. Example of a randomly selected pixel for each fitting iteration.

Appendix C: Impact of beam smearing on the disk size estimate

In this section, we test the reliability of method I when the kinematic properties of the galaxy sample are estimated even in cases of poor angular resolutions and a low signal-to-noise ratio. As described in section 3.3, the beam size can strongly affect the recovery of the kinematical parameters, leading to an underestimation of V_max/σ_gas. This impact is more predominant when the radius of the galaxy is much smaller than the beam size.

We initially assessed whether the angular resolution of our observations might lead to a bias in the estimate of the galaxy disk sizes. Figure C.1 shows the disk size (2r_d) from the fit of flux maps as a function of the full width at half maximum of the ALMA beam, FWHM_beam. Although the ratio of the beam and the disk size spans a wide range from 0.3 to 3.5, we do not find any clear relation between the two parameters, suggesting that the galaxy size is well recovered within the uncertainties in all observations selected for this work.

Fig. C.1. Recovered galaxy size (2rd) as a function of the beam size. Circles are galaxies fit with method I, and stars are galaxies fit with method II. The colors represent the redshift of the galaxy (see the color bar).

To verify the impact of the beam smearing on the inferred kinematic parameters, we simulated a set of ALMA observations, performed a kinematic analysis based on method I, and evaluated the discrepancy between the best-fits results and input parameters. In particular, we simulated ALMA observations with different angular resolutions of a disk galaxy at z = 6 with a scale radius of 1 kpc. For each selected angular resolution, we produced a set of ten ALMA mock datacubes with a signal-to-noise ratio of ∼10 so that we can determine the uncertainties and biases associated with our inferred kinematic parameters. In Figure C.2 we show the results of the kinematical fitting of four important kinematical parameters (i.e, scale radius, velocity dispersion, inclination, and V_max/σ_gas) as a function of the FWHM_beam/2r_d ratio. The uncertainties on the kinematic parameters increase with decreasing angular resolution, but we do not find any clear bias or trend. However, we note that the V_max/σ_gas ratio is underestimated in all ten ALMA simulations with FWHM_beam/2r_d = 3.5, but this does not change the result of our work on the real data because we have only one dataset with a poor angular resolution like this.

Fig. C.2. Best-fitting results of the scale radius (upper left), inclination (bottom left), σgas (top right), and Vmax/σgas (bottom right) for the mock datacubes with different FWHMbeam/2rd. The dashed red line represents the input parameters of the model we used to create the mock ALMA datasets.

We furthermore tested the reliability of this method in determining the velocity dispersion and the ratio V_max/σ_gas even for data with low S/N. We created two different mock datacubes with infinite S/N and angular resolution, but different values of the kinematical parameters, such that one is dispersion dominated, and the other is rotation dominated. In particular, the dispersion-dominated model had a scale radius of 1.8 kpc, corresponding to 0.31″at redshift 6, an inclination of 40°, a mass of 10¹⁰, and a velocity dispersion of 80 km/s, resulting in a V_max/σ_gas = 1.4. The rotation-dominated model had the same values for the radius and the inclination, but a mass of 3 × 10¹⁰ M_⊙ and a velocity dispersion of 40 km/s, resulting in a V_max/σ_gas = 5. We then simulated ALMA observations by changing both the angular resolution and the S/N. For the angular resolution, we tested two values of 0.8 and 1.5 times the mock galaxy radius that represents the bulk of our data, while for the S/N, we test four different levels of S/N that comprise the range of S/N of the data. Using method I, we then fit the mock cubes and recovered the kinematical parameters.

In Fig. C.3 we report the results of the velocity dispersion and the ratio V_max/σ_gas for the dispersion-dominated model and the rotation-dominated model, respectively. We do not find any significant bias in the recovery of these parameters. The median value of the velocity dispersion is always within the limit of the velocity resolution of the datacube (10 km s⁻¹). For the models with low S/N and 1″beam, especially with the dispersion-dominated datacube, we find that the degeneration between the rotational mass and the inclination of the disk does not allow us to constrain the value of the velocity, but we only find the upper limits, similar to what happens for our data (see Fig. E.5).

Fig. C.3. Best-fit values of σgas (left) and Vmax/σgas (right) for the mock datacubes with different S/N and different beams for the dispersion-dominated model (upper panel) and rotation-dominated model (lower panel). In green we plot the results obtained from the model with 0.5″beam, and in blue we plot the results obtained from the model smeared with a beam of 1″. The dashed red line represents the values of the mock datacube with infinite angular resolution.

Appendix D: Method II tests and correction factor

In this section, we test the results of method II, first by using the results of method I to understand the correctness in the estimation of σ_gas, and then we try to understand the limits of this method, which is principally based on an a priori determination of some kinematical parameters. To understand the correctness of the value of velocity dispersion extrapolated with method II, we exploited the galaxies analyzed with method I, and we also analyzed them with method II. The results of the comparison are shown in Figure D.1. This analysis shows an overestimation of σ obtained from method II in comparison with the value of σ obtained with method I for the same galaxy. Analyzing the difference between the two different values, we found that the velocity dispersion estimated with method II is overestimated in the median by a factor of 1.6. Therefore, we applied this correction factor to the velocity dispersion obtained from galaxies analyzed with the second method.

Fig. D.1. Comparison between the velocity dispersion of the same galaxies obtained with the two different methods.

As already mentioned in Sec. 3.6, for this method, we had to determine the radius and the galaxy mass a priori while we constrained the inclination between 3σ to the results of the flux map fitting. To understand the effect of the radius and inclination in the determination of the velocity dispersion, we performed a kinematic fitting on simulated ALMA observations with method II and also left the radius and the inclination as free parameters The mock datacube simulated the ALMA observations of a disk galaxy with a radius of 1 kpc, an inclination of 45°, a mass of 3 × 10¹⁰ M_⊙, and a σ_gas = 80 km s⁻¹. The S/N of the datacube was about 10. In the kinematic fitting, all free parameters had flat priors, and in particular, the prior within [5 degrees, 85 degrees], [0.3 kpc, 5 kpc], and [0 km s⁻¹, 300 km s⁻¹ ] was set for disk inclination, scale radius, and velocity dispersion, respectively. In Figure D.2 we show the corner plot of the posterior distribution from the kinematic fitting on the mock datacube. We note that the posterior distributions of the disk inclination and scale radius are almost flat, while the posterior of the velocity dispersion parameters clearly peaks around 100 km⁻¹ and has a width of 80 km/s. On the other hand, the posterior distribution that we obtain with a fixed scale radius peaks at 87 km⁻¹, which is more consistent with the input value, and has a width of 70km⁻¹. The combination of the two results indicates that the shape of the posterior distribution of the velocity dispersion parameter is only slightly modulated by the disk inclination and scale radius parameter.

Fig. D.2. Corner plot distribution of the free parameters for the mock datacube analyzed with method II, in which the radius was left free to vary. In red we show the flat priors.

Appendix E: 2D kinematical fitting results

In this section, we show the moment maps derived from the pixel-by-pixel Gaussian fitting, the best model, and the 1D distribution of the free parameter recovered with the 2D fitting procedure explained in Sec. 3.5. To recover the best values of the free parameters, we combined the results of 100 trials, in which each time, we randomly selected pixels at a distance greater than half of the beam size to take the correlation between spatial pixels into account. To create the final distributions of the free parameters, we combined the results of each distribution for every fitting.

Fig. E.1. Kinematic fitting using method I for COS-30 observed with the [C II] emission line. In the left column, we present the flux, velocity, and velocity dispersion map derived with the Gaussian fitting. The beam is shown as the gray ellipse in the flux map. In the central panels, we present from top to bottom the best-result maps of the zeroth, first, and second moment. In the right panels, we show the residuals of the fitting, with color bars ranging between ±3σ, the black line on the color bar indicates 1σ. The flux maps are normalized for a pixel with a maximum flux = 1. The x and y values of the kinematical maps are the measure in arcseconds from the central pixel of the image. In the center, we show the corner plot distribution of the free parameters. In the right panel, we show the position-velocity diagram along the major and minor axes in the upper and lower panel, respectively. In black we plot the data, in red we plot the model, and the contours are computed at 3 and 6 σ, respectively.

Fig. E.2. Target COS-29 observed with the [C II] emission line. See the caption of Fig. E.1.

Fig. E.3. Target J1211 observed with the [C II] emission line. See the caption of Fig. E.1.

Fig. E.4. Target J1211 observed with the [O III] emission line. See the caption of Fig. E.1.

Fig. E.5. Target J0217 observed with the [O III] emission line. See the caption of Fig. E.1.

Fig. E.6. Target HZ7 observed with the [C II] emission line. See the caption of Fig. E.1.

Fig. E.7. Target UVISTA-Z-349 observed with the [C II] emission line. See the caption of Fig. E.1.

Fig. E.8. Target UVISTA-Y-001 observed with the [C II] emission line. See the caption of Fig. E.1.

Fig. E.9. Target UVISTA-Y-004 observed with the [C II] emission line. See the caption of Fig. E.1.

Fig. E.10. Target UVISTA-Z-001 observed with the [C II] emission line. See the caption of Fig. E.1.

Fig. E.11. Target UVISTA-Y-879 observed with the [C II] emission line. See the caption of Fig. E.1.

Fig. E.12. Target DLA0817g observed with the [C II] emission line. See the caption of Fig. E.1.

Fig. E.13. Target ALESS 073.1 observed with the [C II] emission line. See the caption of Fig. E.1.

Appendix F: 1D kinematical fitting results

In this section, we present the results of the 1D fitting procedure for galaxies with unresolved kinematics. We recovered the radius and the galaxy center by fitting the flux map that we present for each target in the following figures. We recovered the total baryonic mass by summing the stellar masses (estimated in UV) and the gas mass recovered from the [C II] luminosity. We then performed the fitting on the spatially integrated spectrum. We present the spectrum, the model spectrum, and the corner plot to show the distribution of the free parameters for each target in the following figures.

Fig. F.1. 1D fitting results for the target VR7. In the left panel, we present the flux map at the top. From left to right, we show the observed flux, the model, and the residuals. The beam size is shown as the gray ellipse. At the bottom, we show the best-fit and the observed integrated spectrum. On the right, we show the corner plot of the free parameter of the fitting of the integrated spectrum.

Fig. F.2. Target CLM1, see the caption of F.1.

Fig. F.3. Target UVISTA-Z-004, see the caption of F.1.

Fig. F.4. Target SUPER-8, see the caption of F.1.

Fig. F.5. Target UVISTA-Z-019, see the caption of F.1.

Fig. F.6. Target UVISTA-Y-003, see the caption of F.1

Fig. F.7. Target UVISTA-Z-049, see the caption of F.1.

Fig. F.8. Target J0235 observed with the [C II] emission line, see the caption of F.1.

Fig. F.9. Target J0235 observed with the [O III] emission line, see the caption of F.1.

Fig. F.10. Target J0217 observed with the [C II] emission line, see the caption of F.1.

Fig. F.11. Target COS-30 observed with the [O III] emission line, see the caption of F.1.

All Tables

All Figures

	Fig. 1. Summary of the method I fitting procedure to derive the model parameters.
In the text

Fig. 2. Best-fitting results of HZ4 obtained with method I. The left panels show the moment [C II] maps: the normalized flux map (top), the velocity map (middle), and the velocity dispersion (bottom). From left to right, the panels illustrate the data, models, and residuals. The color bars of the residual range between −3σ and +3σ, and the black lines indicate ±1σ. The ALMA beam is shown as the gray ellipse in the flux map. The multidimensional parameter space explored by step 3 of the method I algorithm is shown in the corner plot on the right.

In the text

Fig. 3. Best-fitting results from method II for the target HZ9. The top left panels illustrate the observed flux map, model, and residuals. The ALMA beam is shown as the gray ellipse in the bottom left corner of the observed flux map. The bottom right panel shows the observed integrated spectrum and the best-fitting model. The multidimensional parameter space explored by the method II algorithm is shown in the corner plot on the right.

In the text

Fig. 4. Evolution of the velocity dispersion (left) and the Vmax/σgas ratio (right) with redshift. Blue and purple marks are the values obtained from our sample by analyzing the [C II] and [O III] observations, respectively. Diamonds are obtained by fitting the three moment maps (method I), and circles represent the results found for the marginally resolved galaxies (method II). Previous σ and Vrot/σ estimates reported in the literature are represented as color-coded rectangles: the observations with the [C II] emission line at 4 < z < 5 by Fraternali et al. (2021), Rizzo et al. (2020, 2021), Lelli et al. (2021), Sharda et al. (2019), and Neeleman et al. (2020) at 4 < z < 6 from Jones et al. (2021) are plotted in blue. Red rectangles show the results of the studies of ionized gas (Hα, [O III] λ 5007) in z < 4 galaxies (Epinat et al. 2010; Green et al. 2014; Turner et al. 2017; Übler et al. 2019). Green rectangles represent the measurements derived from molecular gas (Mogotsi et al. 2016; Girard et al. 2021), (CO). The pink symbols were derived with traced of neutral gas (H I) at z = 0 (Mogotsi et al. 2016). The boxes extend from the 25th to the 75th percentile of the data in the samples, and solid black lines show the median of the value. The light blue area represents the predictions by Wisnioski et al. (2015). The upper and lower boundaries of the shaded light blue region represent the predicted Vmax/σ (see Eq. (13)) for Q = 0.67 and Q = 2 and the solid light blue line represents the Vmax/σgas values predicted for Q = 1. The dashed black line represents the evolution of the best-fitting velocity dispersion trend by Übler et al. (2019). The gray and orange areas show the results obtained by the TNG50 (Pillepich et al. 2019) and SERRA simulations (Kohandel et al. in prep.), respectively.

In the text

Fig. 5. SFR as a function of the stellar mass for our sample of galaxies and the results found in the literature. The colors represent the redshift on the galaxy (see the color bar). The diamonds are the galaxies targeted in this work. The squares show the results found by Sharda et al. (2019), Rizzo et al. (2020, 2021), Lelli et al. (2021), and Fraternali et al. (2021). The triangles show the ALPINE sample (Jones et al. 2021). In green we plot the starburst region at z ∼ 4 − 5 as predicted by Caputi et al. (2017). In light blue we plot the main-sequence region extrapolated between z ∼ 4 − 8 by Schreiber et al. (2015).

In the text

Fig. 6. Dust mass content in the three galaxy samples targeted in this work. In the upper panel, we plot the distribution of the ratio of the dust mass and stellar mass for the sample targeted in this work (in blue) and the sample of DSFG targeted by Rizzo et al. (2020, 2021), the starburst dusty galaxies (Sharda et al. 2019; Fraternali et al. 2021; in red), and the ALPINE sample (Jones et al. 2021; in green). In the lower panel we present the distribution of the ratio of the dust mass and total baryonic mass (Mgas + M*) for the sample targeted in this work, the sample of DSFG targeted by other high-z studies, and the ALPINE sample.

In the text

Fig. 7. Vmax/σgas as a function of redshift and stellar mass. Squares represent the results from KMOS3D (Wisnioski et al. 2019). Circles show the Vmax/σgas for the whole z > 4 sample including the results of this work. Solid lines show the results obtained by TNG50 simulations (Pillepich et al. 2019). Dashed lines show the results obtained by SERRA simulations (Kohandel et. al. in prep.). Different stellar mass bins are represented with different colors.

In the text

	Fig. 8. Vmax/σgas as a function of redshift, and the ratio of the dust mass and stellar mass. In red we show the DSFG sample (Sharda et al. 2019; Rizzo et al. 2020, 2021; Lelli et al. 2021; Fraternali et al. 2021), in green we plot the ALPINE sample (Jones et al. 2021), and in blue we show the sample targeted in this work.
In the text

	Fig. 9. Vmax/σgas as a function of gas fraction. The black line represents the model by Wisnioski et al. (2015) with the assumption of Q = 1 and . In dark blue, we show the data obtained from the sample of galaxies targeted in this work.
In the text

Fig. 10. Dependence of the velocity dispersion on the SFR and the rotational velocity. In the left panel, we show the model by Krumholz et al. (2018) with the feedback as the only driver of the velocity dispersion. In the central panel, we show the model with the release of gravitational energy due to the transport of the gas across the disk as the only driver of the velocity dispersion. In the right panel, we show the model with the combination of feedback and transport as drivers of the velocity dispersion.

In the text

	Fig. 11. Relation of the kinematical properties and SF-related properties for our sample of galaxies following Eq. (16). The dashed line corresponds to the model with Q = 1, the dotted line shows the model with Q = 0.1, and the dash-dotted line shows the model assuming Q = 10
In the text

	Fig. B.1. Error on the inclination parameter for different datacubes with different pixel sizes. The errors derived from our method in which we estimated the likelihood by using only independent pixels is shown in blue, and the uncertainties obtained by assuming all pixel of a moment map as independent pixels are shown in red.
In the text

	Fig. B.2. Example of a randomly selected pixel for each fitting iteration.
In the text

	Fig. C.1. Recovered galaxy size (2rd) as a function of the beam size. Circles are galaxies fit with method I, and stars are galaxies fit with method II. The colors represent the redshift of the galaxy (see the color bar).
In the text

	Fig. C.2. Best-fitting results of the scale radius (upper left), inclination (bottom left), σgas (top right), and Vmax/σgas (bottom right) for the mock datacubes with different FWHMbeam/2rd. The dashed red line represents the input parameters of the model we used to create the mock ALMA datasets.
In the text

Fig. C.3. Best-fit values of σgas (left) and Vmax/σgas (right) for the mock datacubes with different S/N and different beams for the dispersion-dominated model (upper panel) and rotation-dominated model (lower panel). In green we plot the results obtained from the model with 0.5″beam, and in blue we plot the results obtained from the model smeared with a beam of 1″. The dashed red line represents the values of the mock datacube with infinite angular resolution.

In the text

	Fig. D.1. Comparison between the velocity dispersion of the same galaxies obtained with the two different methods.
In the text

	Fig. D.2. Corner plot distribution of the free parameters for the mock datacube analyzed with method II, in which the radius was left free to vary. In red we show the flat priors.
In the text

Fig. E.1. Kinematic fitting using method I for COS-30 observed with the [C II] emission line. In the left column, we present the flux, velocity, and velocity dispersion map derived with the Gaussian fitting. The beam is shown as the gray ellipse in the flux map. In the central panels, we present from top to bottom the best-result maps of the zeroth, first, and second moment. In the right panels, we show the residuals of the fitting, with color bars ranging between ±3σ, the black line on the color bar indicates 1σ. The flux maps are normalized for a pixel with a maximum flux = 1. The x and y values of the kinematical maps are the measure in arcseconds from the central pixel of the image. In the center, we show the corner plot distribution of the free parameters. In the right panel, we show the position-velocity diagram along the major and minor axes in the upper and lower panel, respectively. In black we plot the data, in red we plot the model, and the contours are computed at 3 and 6 σ, respectively.

In the text

	Fig. E.2. Target COS-29 observed with the [C II] emission line. See the caption of Fig. E.1.
In the text

	Fig. E.3. Target J1211 observed with the [C II] emission line. See the caption of Fig. E.1.
In the text

	Fig. E.4. Target J1211 observed with the [O III] emission line. See the caption of Fig. E.1.
In the text

	Fig. E.5. Target J0217 observed with the [O III] emission line. See the caption of Fig. E.1.
In the text

	Fig. E.6. Target HZ7 observed with the [C II] emission line. See the caption of Fig. E.1.
In the text

	Fig. E.7. Target UVISTA-Z-349 observed with the [C II] emission line. See the caption of Fig. E.1.
In the text

	Fig. E.8. Target UVISTA-Y-001 observed with the [C II] emission line. See the caption of Fig. E.1.
In the text

	Fig. E.9. Target UVISTA-Y-004 observed with the [C II] emission line. See the caption of Fig. E.1.
In the text

	Fig. E.10. Target UVISTA-Z-001 observed with the [C II] emission line. See the caption of Fig. E.1.
In the text

	Fig. E.11. Target UVISTA-Y-879 observed with the [C II] emission line. See the caption of Fig. E.1.
In the text

	Fig. E.12. Target DLA0817g observed with the [C II] emission line. See the caption of Fig. E.1.
In the text

	Fig. E.13. Target ALESS 073.1 observed with the [C II] emission line. See the caption of Fig. E.1.
In the text

Fig. F.1. 1D fitting results for the target VR7. In the left panel, we present the flux map at the top. From left to right, we show the observed flux, the model, and the residuals. The beam size is shown as the gray ellipse. At the bottom, we show the best-fit and the observed integrated spectrum. On the right, we show the corner plot of the free parameter of the fitting of the integrated spectrum.

In the text

	Fig. F.2. Target CLM1, see the caption of F.1.
In the text

	Fig. F.3. Target UVISTA-Z-004, see the caption of F.1.
In the text

	Fig. F.4. Target SUPER-8, see the caption of F.1.
In the text

	Fig. F.5. Target UVISTA-Z-019, see the caption of F.1.
In the text

	Fig. F.6. Target UVISTA-Y-003, see the caption of F.1
In the text

	Fig. F.7. Target UVISTA-Z-049, see the caption of F.1.
In the text

	Fig. F.8. Target J0235 observed with the [C II] emission line, see the caption of F.1.
In the text

	Fig. F.9. Target J0235 observed with the [O III] emission line, see the caption of F.1.
In the text

	Fig. F.10. Target J0217 observed with the [C II] emission line, see the caption of F.1.
In the text

	Fig. F.11. Target COS-30 observed with the [O III] emission line, see the caption of F.1.
In the text