ALMA hints at the presence of turbulent disk galaxies at z > 5

Context. High-redshift galaxies are expected to be more turbulent than local galaxies because of their smaller size and higher star formation and thus stronger feedback from star formation, frequent mergers events, and gravitational instabilities. However, this scenario has recently been questioned by the observational evidence of a few galaxies at z ∼ 4 − 5 with a gas velocity dispersion similar to what is observed in the local population. Aims. Our goal is to determine whether galaxies in the ﬁrst billion years of the Universe have already formed a dynamically cold rotating disk similar to the local counterparts. Methods. We studied the gas kinematic of 22 main-sequence star-forming galaxies at z > 5 and determined their dynamical state by estimating the ratio of the rotational velocity and of the gas velocity dispersion. We mined the ALMA public archive and exploited the [C II ] and [O III ] observations to perform a kinematic analysis of the cold and warm gas of z > 5 main-sequence galaxies. We compared our results with what was found in the local and distant Universe and investigated the evolution of the gas velocity dispersion with redshift. We also compared the observations with theoretical expectations to assess the main driver of the gas turbulence at z > 5. Results. The gas kinematics of the high-z galaxy population observed with ALMA is consistent within the errors with rotating but turbulent disks. We indeed infer a velocity dispersion that is systematically higher by 4-5 times than the local galaxy population and the z ∼ 5 dust-obscured galaxies reported in the literature. The di ﬀ erence between our results and those reported at similar redshift can be ascribed to the systematic di ﬀ erence in the galaxy properties in the two samples: the disks of massive dusty galaxies are dynamically colder than the disks of dust-poor galaxies. The comparison with the theoretical predictions suggests that the main driver of the velocity dispersion in high-redshift galaxies is the gravitational energy that is released by the transport of mass within the disk. Finally, we stress that future deeper ALMA high-angular resolution observations are crucial to constrain the kinematic properties of high-z galaxies and to distinguish rotating disks from kiloparsec-scale mergers.


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;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;Förster Schreiber et al. 2018;Johnson et al. 2018;Übler et al. 2019;Girard et al. 2021).The kinematic anal-ysis 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;Johnson et al. 2018;Wisnioski et al. 2019;Ü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 V max /σ gas = √ 3.36 based on the kinematics of the exponential disk (For an exponential disk, V c = (V 2 rot + 3.36σ 2 gas ) 0.5 Eq. 1 Förster Schreiber et al. (2018)).A ratio of V max /σ gas < √ 3.36 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 9 M and 10 11 M shows a σ gas ∼ 10 − 40 km s −1 , 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 −1 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;Wisnioski et al. 2015;Turner et al. 2017;Förster Schreiber et al. 2018;Ü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. 2022;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 where 1 < a < 2, in particular, a = √ 2 for a flat rotation curve, a = √ 3 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;Genzel et al. 2008;Förster Schreiber et al. 2011;Carniani et al. 2017Carniani et al. , 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. 1998Luhman et al. , 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 HII 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. 2017Carniani et al. , 2018;;Smit et al. 2018;Harikane et al. 2020;Neeleman et al. 2020;Rizzo et al. 2020Rizzo et al. , 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 −1 to 116 km s −1 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. 2020Rizzo et al. , 2021;;Fraternali et al. 2021;Lelli et al. 2021).The measured velocity dispersion is as low as 15-20 km s −1 , 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. 2020Rizzo et al. , 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 (S FR > 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 starforming 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 Section 2 we describe the sample selection and data reduction processes.In Section 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 section 6 we summarize the findings of this work and draw our conclusions.We adopt the cosmological parameters from Planck Collaboration et al. ( 2016): H 0 = 67.8km s −1 Mpc −1 , Ω m = 0.308, and Ω Λ = 0.685.

Sample selection and observations
We selected our galaxy sample from the public ALMA data archive 1 .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.52 .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 Sec.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 ratio3 (S/N)< 7, which is not sufficient to perform a kinematic analysis on the data.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 Table 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 −1 ).

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 III] data on the basis of the angular resolution and sensitivity of the observations.

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

Kinematic models
To reproduce the kinematic maps of the FIR lines, we used a rotational velocity model, V rot (r), defined by 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 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 χ 2 .
In cases in which the sensitivity and resolution are not high enough to determine the presence of the bulge or M < 10 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 where r d is the disk scale length.Assuming a constant mass-tolight ratio, we write the profile of the surface mass density as The velocity contribution due to the exponential disk component can therefore be written as a function of the radius as where y ≡ r 2r d , while I 0 , I 1 , K 0 , and K 1 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 0 = 5 kpc, which corresponds to 2-3 disk scale lengths: On the other hand, the light profile of the bulge component is assumed to match a Sérsic profile (Sérsic 1963), where I 0 is the surface brightness at the center of the bulge, r e is the effective radius, and n is the Sérsic index.b is a function of n b = 2n -1/3 + 0.009876/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 (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 2 ; 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.

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.

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 5 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.

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 0 , y 0 ), the rotational mass 4 (M rot ), 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 Figure 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 bestfitting 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 Sec. 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 4 M rot is the mass of the disk (or disk + bulge) that is required to reproduce the velocity maps without taking into account the contributions from the turbulence support as well.
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.

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: low angular resolution of the observations; low sensitivity; the disk is not dominated by rotation, but it is instead dominated by random motion, hence we cannot see a velocity gradient.
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 where M * is the stellar mass, based on observations reported in the literature (Capak et al. 2015;Willott et al. 2015; Matthee with α [CII] = 30M /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 S model (a, b, σ gas , inc) = b + aF datacube (σ gas , inc). (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 −1 and 300 km s −1 .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.
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.
Fig. 3 illustrates an example of the best-fitting results from method II.

Results
In this section, we report the results from the kinematic fitting by adopting the rotating models described in Sec.3.2.Table 1 summarizes the results provided by method I and II.In addition to Fig. 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 χ 2 .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 −1 and 132 km s −1 , with a median velocity dispersion of 67 +15 −13 km s −1 .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 Sec. 3.6), between 37 and 126 km s −1 and a median value of 60 +19 −11 km s −1 , 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)  ), 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 −1 in the internal regions to ∼ 10 km s −1 in the outer region.Our results in Figure E.13 show a bimodal distribution of the velocity dispersion that converges at ∼ 30 km s −1 when most of the selected pixels were in the outer regions, and to ∼ 56 km s −1 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 −1 ) and not to the overestimated velocity dispersion (σ gas = 40 ± 10 km s −1 ).
Another galaxy that was observed with ALMA and was recently analyzed by Neeleman et al. (2020)  −0.3 .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 V max /σ gas = 4.01 +0.21  −0.73 , 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 σ gas = 65.8 +2.9 −3.3 km s −1 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. (2022) using the algorithm 3D Barolo.They found that the average velocity dispersion has an upper limit of 30 km s −1 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.

Velocity dispersion and V max /σ gas evolution with redshift
The left panel of Figure 4 shows the evolution of the gas velocity dispersion as a function of the redshift from z = 0 to z = 8.The best-fitting velocity dispersion estimates of our sample are shown with blue ([C II]) and purple ([O III]) marks (diamonds for the results found with method I, and circles for the results found with method II).The median values inferred for the sample of galaxies analyzed with methods I and II are reported as black marks.In the figure, we also report previous results from the literature that are shown as boxes that extend from the 25th to 75th percentile of the data in their samples and are color-coded as a function of the tracer adopted to study the gas kinematics.Differently from the conclusion reported by Lelli et al. (2021), Rizzo et al. (2020Rizzo et al. ( , 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 −1 ; blue region at 4 < z < 4.5 in figure 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, The trend of eq. 12 with Q = 1 and 100 km s −1 < V max < 250 km s −1 is shown with the light blue shaded curve in the left panel of Figure 4.As the gas fraction increases with redshift (Tacconi et al. 2020), the velocity dispersion also increases, reaching values of ∼ 40−100 km s −1 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 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 Figure 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), The upper and lower boundaries of the shaded light blue region represent the predicted V max /σ for Q = 0.67 and Q = 2, respectively, and the solid light blue line represents the values predicted for Q = 1 (Wisnioski et al. 2015).The shaded gray region illustrates the results from the TNG50 simulations (Pillepich et al. 2019), and the shaded orange region reports the results 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 agreet with the estimates found for galaxies in the ALPINE survey (Jones et al. 2021) 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.

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 (Figure 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),  (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 Figure 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.While for the dust mass of the sample targeted in this work see Table A.1.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.
The figure shows that while our sample and the ALPINE samples have a star-to-dust mass ratio of 10 −3 on average, the DSFGs are characterized by M dust /M > 10 −2 , 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.

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 10 have V max /σ gas > 3, while M < 10 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.(Wisnioski et al. 2019).Circles show the V max /σ 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.
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 (equation 1; see also Wisnioski et al. 2015).In Figure 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.
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-Fig.8. V max /σ 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. 2020Rizzo et al. , 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.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.

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 Figure 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.
-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, starforming 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 highredshift observations are not able to explore the region at low SFR, so that we are not able to distinguish between the feed-back+transport and the transport-only models.
In Figure 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 where f gas,p is the gas fraction in the midplane of the galaxy ( f gas,p = 1.5 f 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 Figure 11 S FR/ f gas,p as a function of σ gas v 2 .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.
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 feed-back+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).

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 −1 , 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 3D Barolo.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 starforming 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 explo-sions) 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 (V max /σ gas √ 3.36) 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(Rizzo et al. , 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 Sec. 5.1).However, we note that these massive (M * = 10 11 M ) galaxies are highly starforming galaxies (S FR > 100M yr −1 ), as shwon in figure 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 Sec. 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 section 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 Section 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 highredshift 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.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-tonoise 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.
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 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 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 −1 ).For the models with low S/N and 1 beam, especially with the dispersiondominated 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.   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. 1 .
Fig. 1.Summary of the method I fitting procedure to derive the model parameters.

Fig. 2 .
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.

Fig. 3 .
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.
and Jones et al. (2021) is DLA0817g.Neeleman et al. (2020) exploited the higherresolution 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 3D Barolo.Jones et al. (2021) studied the kinematics by exploiting lower-resolution data (∼ 1 ) using 3D Barolo.Their results agree with those byNeeleman et al. (2020), who they found σ gas = 80 +13 −11 km s −1 and V max /σ gas =3.4+1.1

Fig. 4 .
Fig. 4. Evolution of the velocity dispersion (left) and the V max /σ 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 V rot/σ 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), Rizzo et al. (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 byWisnioski et al. (2015), and the solid light blue line represents the V max /σ 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.

Fig. 5 .
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), Rizzo et al. (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 byCaputi et al. (2017).In light blue we plot the main-sequence region extrapolated between z ∼ 4 − 8 bySchreiber et al. (2015)

Fig. 6 .
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 (M gas + M * ) for the sample targeted in this work, the sample of DSFG targeted by other high-z studies, and the ALPINE sample.

Fig. 7 .
Fig. 7. V max /σ gas as a function of redshift and stellar mass.Squares represent the results from KMOS 3D(Wisnioski et al. 2019).Circles show the V max /σ 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.

Fig. 9 .
Fig. 9. V max /σ 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 a= √2.In dark blue, we show the data obtained from the sample of galaxies targeted in this work.

Fig. 10 .
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.

Fig. 11 .
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 Notes.(1): target name; (2) and (3): stellar masses computed from the UV luminosity and references; (4) and (5): SFR computed as SFR= SFR UV + SFR IR (if not indicated otherwise with a symbol) and references; (6) and (7): dust masses computed by SED fitting of the rest-frame FIR emission temperatures and references.[CII] Estimated from the [C II] luminosity SED Estimated from SED fitting IR Estimated from the IR luminosity UV Estimated from the UV luminosity the galaxy size is well recovered within the uncertainties in all observations selected for this work.

Fig
Fig. C.1.Recovered galaxy size (2r d ) 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). E.5).

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

Fig
Fig. C.3.Best-fit values of σ gas (left) and V max /σ gas (right) for the mock datacubes with different S/N and different beams for the dispersiondominated 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.

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

Fig. D. 2 .
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.

Fig. E. 1 .
Fig. E.1.Kinematical fitting using method I. On the left, we show moment maps 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
Fig. E.2.Target COS-29 observed with the [C II] emission line.See the caption of Fig. E.1.

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

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

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

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

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

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

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

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

Fig
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 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
Fig. F.8. Target J0235 observed with the [C II] emission line, see the caption of F.1.

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

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

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

Table 1 .
Results of the kinematical fitting for our sample of galaxies.
Notes.(1) target name; (2) observed line; (3) beam FWHM; (4) best-fitting method; (5) redshift of the FIR line; (6) scale radius of the exponential disk; (7) disk inclination in degrees; (8) position angle of the galaxy; (9) logarithm of the dynamical mass in solar masses; (10) maximum rotational velocity as computed for an exponential disk with the dynamical mass and the scale radius obtained from the fitting; (11) best-fit velocity dispersion; (12) ratio of the maximum rotational velocity and the velocity dispersion.b best-fit results with the model exponential disk + stellar bulge.
Lelli et al. (2021)2021)rom the results determined in the sample of galaxies studied bySharda et al. (2019),Rizzo et al. (2020),Rizzo et al. (2021),Fraternali et al. (2021)andLelli 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