Dynamically cold disks in the early Universe: myth or reality?

Theoretical models struggle to reproduce dynamically cold disks with significant rotation-to-dispersion support($V_{\rm{rot}}/\sigma$) observed in star-forming galaxies in the early Universe, at redshift $z>4$. We aim to explore the possible emergence of dynamically cold disks in cosmological simulations and to understand if different kinematic tracers can help reconcile the tension between theory and observations. We use 3218 galaxies from the SERRA suite of zoom-in simulations, with $8<\log(M_*/M_{\odot})<10.3$ and SFR$<128\,M_{\odot}{yr}^{-1}$, within $4<z<9$ range. We generate hyper-spectral data cubes for 6436 synthetic observations of H$\alpha$ and [CII]. We find that the choice of kinematic tracer strongly influences gas velocity dispersion estimates. When using H$\alpha$ ([CII]) synthetic observations, we observe a strong (mild) correlation between $\sigma$ and $M_*$. Such a difference arises mostly for $M_*>10^9\,M_{\odot}$ galaxies, for which $\sigma_{H\alpha}>2\sigma_{CII}$ for a significant fraction of the sample. Regardless of the tracer, our predictions suggest the existence of massive ($M_*>10^{10}M_{\odot}$) galaxies with $V_{rot}/\sigma>10$ at $z>4$, maintaining cold disks for>10 orbital periods (200Myr). Furthermore, we do not find any significant redshift dependence for $V_{rot}/\sigma$ ratio in our sample. Our simulations predict the existence of dynamically cold disks in the early Universe. However, different tracers are sensitive to different kinematic properties. While [CII] effectively traces the thin, gaseous disk of galaxies, H$\alpha$ includes the contribution from ionized gas beyond the disk, characterized by prevalent vertical or radial motions that may be associated with outflows. The presence of H$\alpha$ halos could be a signature of such galactic outflows. This emphasizes the importance of combining ALMA and JWST/NIRspec studies of high-z galaxies.


Introduction
Disks are almost ubiquitous within the star-forming galaxies (SFGs) population in the local Universe.Somewhat surprisingly, recent observations (Ferreira et al. 2022;Kartaltepe et al. 2023;Robertson et al. 2023;Tohill et al. 2023) using the James Webb Space Telescope (JWST) have unveiled their presence also in the very early Universe, stretching as far back as  ∼ 9.However, these early disks might exhibit distinct dynamical characteristics when compared to their local counterparts.
Lastly, JWST/NIRspec multi-object spectroscopy of H and [O III] emission lines has started shedding light on the kinematics of ionized gas in EoR galaxies (de Graaff et al. 2023;Parlanti et al. 2023b), hinting that early galaxies can settle into dynamically cold disks (DCD).Despite this progress, the statistical relevance of DCDs in early SFGs is uncertain, given the limited number of deep spatially resolved observations.From a theoretical standpoint, most studies (Dekel & Burkert 2014;Zolotov et al. 2015;Hayward & Hopkins 2017;Pillepich et al. 2019) struggle to explain the existence of DCDs at  > 4. For instance, tng50 (Pillepich et al. 2019) shows an average  rot / < 3 at  = 4, similarly to most cosmological simulations.Nonetheless, a few studies predict that relatively massive galaxies ( ★ ∼ 10 10  ⊙ ) can temporarily sustain cold disks formed through intense accretion of co-planar, co-rotating gas at  > 3 (Kretschmer et al. 2022).Regarding EoR galaxies, Kohandel et al. 2020 show that moderate rotation support ( rot / ∼ 7) can be achieved in  ★ ∼ 10 10  ⊙ galaxies as far back as  ∼ 6, however, this conclusion is based on a relatively small sample of galaxies.
In this Letter, we exploit the serra simulations (Pallottini et al. 2022) to investigate the kinematic properties of normal SFGs at 4 <  < 9 through the analysis of mock observations of [C II] emission line at 158 m, a tracer of cold ( ∼ 100 K) molecular/neutral gas, and H emission line, a tracer of warm ( ∼ 10 4 K) ionized gas.

Galaxy formation and evolution
The serra suite of simulations focuses on studying the formation and evolution of galaxies during the EoR (Pallottini et al. 2022).Gas and dark matter are evolved using a customized version of the adaptive mesh refinement code ramses (Teyssier 2002).krome (Grassi et al. 2014) is employed to model the nonequilibrium chemical network, that includes H, H + , H − , He, He + , He ++ , H2, H2 + , electrons, and metals, encompassing ∼ 40 reactions (Bovino et al. 2016;Pallottini et al. 2017a).The tracking of metallicity () involves summing heavy elements, assuming solar abundance ratios for different metal species (Asplund et al. 2009).Dust is approximated to scale with a fixed dust-to-metal ratio, denoted as D = D ⊙ (/Z ⊙ ), where D ⊙ /Z ⊙ ≃ 0.3 for the Milky Way (Hirashita & Ferrara 2002).A Milky Way-like grain size distribution is adopted (Weingartner & Draine 2001).An initial metallicity floor of  floor = 10 −3 Z ⊙ is adopted, as expected from pre-enrichment of the intergalactic medium around density peaks (Madau et al. 2001;Pallottini et al. 2014a,b).
The conversion of molecular hydrogen into stars follows a Schmidt (1959)-Kennicutt (1998)-like relation (Pallottini et al. 2017a).These stars, in turn, act as sources of metals, mechanical energy, and radiation (Pallottini et al. 2017b).The interstellar radiation field (ISRF) is dynamically evolved on the fly using the moment-based solver from ramses-rt (Rosdahl et al. 2013), which is linked to the chemical evolution (Pallottini et al. 2019;Decataldo et al. 2019).To efficiently model radiation propagation, the speed of light is reduced by a factor of 10 3 in serra, leading to negligible deviations compared to a 102 reduction (Pallottini et al. 2019;Lupi et al. 2020).Simulations track radiation in 5 energy bins, with one bin partially covering the Habing band (6.0 < ℎ < 11.2), one dedicated to the Lyman-Werner band (11.2 < ℎ < 13.6) to address H 2 photoevaporation, and the remaining three bins covering ionization processes from H to the first ionization level of He (13.6 < ℎ < 24.59).
Each run in the serra suite is initialized at  = 100 from cosmological conditions generated with music (Hahn & Abel 2011), and zooms in on target DM halos selected at around  ≃ 6.The cosmic volume considered is (20, Mpc/h) 3 , evolved with a base grid of 8 levels.The zoom-in region has a volume of about (2.1 Mpc/h) 3 (approximately 5 times the virial radius of the target DM halo), and has 3 additional levels, resulting in a gas mass resolution of   = 1.2 × 10 4 M ⊙ .Additionally, 6 extra levels of refinement are enabled in the zoom-in region based on a Lagrangian-like criterion, allowing the simulation to reach scales of  res ≃ 30 pc at  = 6 in the densest regions, similar to Galactic molecular clouds (Federrath & Klessen 2013).
Due to the coarse nature of the chemical network used in hydrodynamical simulations, precise emission computation requires post-processing of the data to extract kinematic information.The line luminosity ( em−line ) for each gas cell are obtained using the spectral synthesis code cloudy (Ferland et al. 2017).This process takes into account the interstellar radiation field, the turbulent and clumpy structure of the interstellar medium (ISM), which is parameterized as a function of the local gas Mach number (Vallini et al. 2018;Pallottini et al. 2022).
With information on  em−line , position (x), velocity (v), and thermal+turbulent line broadening (( 2 th +  2 nt )1/2 ) for each gas cell within a specified field of view (FOV)1 and along a line of sight direction, we construct 3-dimensional HDCs.These cubes comprise two spatial dimensions and one spectral dimension, effectively mapping the 6-dimensional data to coordinates (, ,   ).In HDCs, the surface brightness of the emission line is recorded for each voxel, providing valuable insights into the spatial and spectral distribution of the emission.The contribution of all gas cells within the FOV can be directly summed for optically thin emission lines.However, for optically thick lines, radiative transfer through dust needs to be considered when comparing with pre-dereddened observations (Behrens et al. 2018).
In this study, we model the [C II] 158m emission line as a tracer of cold neutral/molecular gas and the nebular H emission line as a tracer of warm ionized gas2.Different gas phase tracers can yield different values of  rot and  (Kohandel et al. 2020;Ejdetjärn et al. 2022).In this work, we estimate the velocity dispersion of a galaxy separately using two tracers: [C II] observations ( [CII] ) and H ( H ).Following a similar approach to Kohandel et al. ( 2020),  em−line represents the luminosity-weighted average velocity dispersion, calculated using the moment-2 and moment-0 maps of the corresponding emission line.For the rotational velocity, we estimate3  rot with the circular velocities of the galaxy, denoted as  c = (  dyn /  ) 1/2 , where  dyn =  g +  ★ is the dynamical mass within the desired FOV, and   is the disk effective radius, i.e.where 50 % of gas mass is contained.In other words, we keep  rot constant regardless of the tracer we use.Con-Table 1: Relevant properties of the simulated sample.Shown are the number of galaxies in each mass bin, their average  rot / [CII] and  rot / H ratios.In this study, we focus on serra galaxies that exhibit a stellar mass of 10 8  ⊙ and higher and covers the redshift range 4 ≤  ≤ 9. Our sample includes 3218 galaxies with SFRs ranging from 0.04 to 128 M ⊙ yr −1 and stellar masses in 8 ≤ log  ★ / ⊙ ≤ 10.3.We refer to Tab. 1 for the overview of our sample.

Rotation support in early galaxies
In Fig. 1 and Fig. 2, we present the kinematic characteristics of our sample within the  ★ - and  rot /-4 planes, respectively.In both cases,  values are derived from [C II] 158 m and H synthetic observations.In Fig. 2, along with serra galaxies, we plot predictions from tng50 simulations (Pillepich et al. 2019) Regarding the  ★ - relation, we find two distinct behaviours for the two tracers.Specifically, we find that  H exhibits a steeper increase with  ★ compared to  [CII] .Interestingly, this trend is particularly pronounced in ∼ 40% of our high-mass bin galaxies and 30% of those in the intermediate-mass bin, where  H > 2 [CII] .Such behaviour is not entirely unexpected, as even surveys of local galaxies (Levy et al. 2018;Girard et al. 2021) have highlighted systematic differences between hot ionized and neutral/molecular cold gas velocity dispersions.This disparity leads to a consistent difference between the results obtained using different tracers in serra, indicating that cold gas exhibits a higher degree of rotational support than warm ionized gas (see Tab. 1 and Fig. 2).We will explore this point in more detail in Sec. 4.
As for the  rot /-, regardless of the kinematic tracer employed, we do not find any notable correlation even when different mass bins are considered.Nevertheless, a clear trend appears with the stellar mass of galaxies, indicating that most massive galaxies exhibit greater rotation support than the less massive ones.A milder version of such a correlation was also predicted in tng50 for  < 4 galaxies.
The most noteworthy finding of our analysis is that when we classify our galaxies based on their  rot / ratio (as outlined in Tab. 2), not only we identify Super Cold disks ( rot / > 10 similar to SFGs observed at  ∼ 4 (Rizzo et al. 2021;Fraternali et al. 2021) and at  ∼ 7 at Rowland et al., in prep. 2023) within our massive sub-sample ( ★ > 10 10  ⊙ ), but we also ascertain that [C II] emitting gas in ∼ 60% of the whole sample is dynamically cold, having 4 <  rot / < 10.This finding suggests that galactic disks can form as early as the EoR and if deep ALMA observations targeting SFGs at  > 4 become available, more dynamically cold disks will likely be uncovered.
For the majority of galaxies with [C II] observations,  rot / values at  > 5 fall below the predicted mean values from serra simulations.This discrepancy can primarily be attributed to the marginal resolution of these observations (∼ 0.1 − 1.5 arcsecs), compared to our synthetic datacubes featuring a higher spatial resolution of 0.005 arcsecs.In fact, in Kohandel et al. 2020, we showed how the beam-smearing effect in low-resolution observations can lead to a substantial overestimate of velocity dispersion, reaching up to ∼ 100% (see also Ejdetjärn et al. 2022;Rizzo et al. 2022).Another potential cause of such discrepancy is the challenge of accurately estimating disk inclination in high- kinematic observations.This becomes crucial when determining kinematic properties, such as velocity dispersion, from integrated spectra, where the shape heavily correlates with the disk inclination (Kohandel et al. 2019).

Why [CII] and H𝛼 kinematics are different?
In Fig. 3, we show the synthetic [C II] and H maps and kinematics observables for face-on and edge-on views of Hibiscus.The [C II] spectrum is characterized by a narrow and prominently Gaussian-shaped profile (FWHM∼ 185 km s −1 ), while the H spectrum appears more complex, broader (FWHM∼ 437 km s −1 ), and exhibits high-velocity wings.Such broad wings in the spectrum could be indicative of outflowing gas.Indeed, comparing the moment maps, we see that the [C II] emission line effectively traces Hibiscus thin gaseous disk, while H traces ionized gas that lies beyond the disk plane, including gas that might be in an in/outflowing state.Such a difference between various phases of the ISM could conceivably arise due to distinct effects of stellar feedback influencing them, as suggested by isolated disk galaxy simulations (Ejdetjärn et al. 2022).This illustrates that the observed velocity dispersion in a given galaxy using H data may not solely arise from turbulence within the galactic disks.Instead, a substantial contribution from outflows may be in effect, introducing an additional layer of complexity to the data interpretation.
We note that the spatial extent of [C II] and H emission in Hibiscus differs significantly.The [C II] is 4× more extended than the stellar effective radius, similar to observed high- galaxies (Fujimoto et al. 2019;Carniani et al. 2020;Fudamoto et al. 2022).The H distribution is even more far-flung.This happens since H originates from the  ∼ 10 4 K photoionized regions outside the disk that are part of an expanding, cooling outflow through which LyC photons percolate.As carbon in these regions is ionized to higher states (e.g.CIII), [C II] emission is limited to denser, more confined regions where recombination rates are higher.Interestingly, this shows that H halos are intimately linked to the presence and morphology of these outflows, offering intriguing prospects for their detection with JWST.
Lastly, it is important to highlight that [C II] kinematics can also be challenging and merits deeper exploration.In particular, observations of the so-called [C II] halos (Gallerani et al. 2018;Fujimoto et al. 2019Fujimoto et al. , 2021;;Ginolfi et al. 2020) have been inter- preted in the framework of outflow models (Pizzati et al. 2020(Pizzati et al. , 2023)), but they have not yet been reproduced by cosmological simulations (Fujimoto et al. 2019;Arata et al. 2020).

Are high-𝑧 dynamically cold disks a transient feature?
To explore the stability of cold disks, we take a closer look at the evolutionary path undertaken by individual galaxies.In Fig. 4 Regarding the  rot / H , there is an interesting difference between the two systems.As witnessed in Fig. 3, the H emitting gas is found to be dynamically warm in the last ∼ 200 Myr for Hibiscus, potentially attributable to the presence of outflows.However, in the case of Narcissus, the gas traced by both emission lines follows a similar evolutionary path.Despite a slightly lower  [CII] than  H , this galaxy remains dynamically cold according to both tracers.Considering the comparable stellar masses of these galaxies, differences in their star formation histories, feedback effects, or other global properties may be responsible for the different behaviour of the tracers.

Conclusions
In this Letter, we investigate the possible existence of dynamically cold disks (with significant rotation support) in the early Universe using a sample of 3218 Lyman Break Galaxies (LBG) from the serra zoom-in cosmological simulations.We analyze the kinematic of both [C II] and H in the redshift range of 4 ≤  ≤ 9 for LBGs with 8 ≤ log ( ★ / ⊙ ) ≤ 10.3 and 0 < SFR ≤ 128.Our main conclusions are the following.
• A strong (mild) correlation is present between stellar mass and gas velocity dispersion when H ([C II]) synthetic observations are considered.The difference arises mostly for  ★ > 10 9  ⊙ galaxies where  H > 2 [CII] .• Irrespective of galaxy mass and the chosen kinematic tracer, our analysis reveals no significant redshift dependence in the ratio  rot /.• Massive ( ★ ≥ 10 10  ⊙ ) galaxies in serra settle into dynamically Super Cold disks with  rot / > 10 at  > 4. Such cold disks are not transient features as they last for more than 10 galaxy orbital times (∼ 200 Myr).
We have shown that in serra galaxies, [C II] effectively traces the thin gaseous disks within galaxies, while H emission can also trace the ionized gas outside the disk.The differences in the kinematics of [C II] and H may be attributed to galactic outflows, although further exploration is necessary to substantiate and statistically quantify this point.We show that the identification of H halos could be a signature of such galactic outflows.We foresee that more high- dynamically cold disks will be found with the increasing availability of deep ALMA observations targeting [C II] 158 m in galaxies with stellar masses exceeding 10 9  ⊙ .
In view of the essential role of multiple tracers in gaining a comprehensive understanding of early galaxy kinematics, we emphasize that the ALMA-JWST/NIRspec synergy will be essential.

Fig. 1 :
Fig. 1: The relationship between velocity dispersion () and stellar mass ( ★ ) in the serra galaxy sample: Blue (red) contours show the 1, 2, and 3-sigma probability density function levels for the  ★ −  relationship derived from synthetic [C II] (H) observations.Individual data points are represented by crosses.

Fig. 3 :
Fig. 3: Multi-wavelength kinematics of "Hibiscus" at  4.5.In the top (last) two rows we display FIR [C II] 158m (nebular H) syntactic observation for line spectra, 0, 1, and 2 moment maps in different columns, for both face-on and edge-on views.

Fig. 4 :
Fig. 4: Example of dynamical evolution of two serra galaxies traced by [CII] and H; at  = 4.5, Hibiscus (top) has a [CII] dynamically cold disk and a turbulent H emitting gas possibly featuring galactic outflows (see also Fig. 3); at  = 6.8,Narcisuss (bottom) appears as a super cold galaxy both in [C II] and H.

Table 2 :
5, as well as observed data of  > 4 galaxies through FIR [C II] emission line by ALMA and nebular H and [O III] lines by JWST/NIRspec.Dynamical categorization of serra galaxies.Shown are the number of simulated galaxies in each dynamical stage depending on the adopted tracer.
, we show the evolution of  rot / as a function of lookback time for Hibiscus and Narcissus.We see that, considering only [C II] emitting gas, both galaxies exhibit a consistent  rot / around 2 − 4 up until roughly 200 Myr.During this interval, the rotation support rises due to the effective accretion of gas and the efficient transfer of angular momentum into the disk.Indeed, if we estimate the disk orbital time  orb = 2  / rot , for Hibiscus it is ∼ 16 Myr, and ∼ 21 Myr for Narcissus.Therefore, these high- DCDs survive for more than 10 orbital times.