During the past decades, there has been outstanding progress in studying the internal dynamics of high-z galaxies. Nearinfrared (NIR) spectroscopy with integral field units (IFUs) allowed for the kinematics of warm (T ' 10 4 K) ionized gas to be traced using the H↵ emission line at z ' 1.0 2.5 (e.g., Förster Schreiber et al. 2009;Wisnioski et al. 2015;Stott et al. 2016) and the [O III] 5007 Å line up to z ' 3.5 (e.g., Gnerucci et al. 2011;Turner et al. 2017). Radio and submillimeter observations with the Jansky Very Large Array (JVLA) and the NOrthern Extended Millimeter Array (NOEMA) allowed for the kinematics of cold neutral gas (T ' 10 100 K) to be traced using CO transitions at z ' 1 4 (e.g., Hodge et al. 2012;Übler et al. 2018). Moreover, the Atacama Large Millimeter Array (ALMA) made it possible to study gas dynamics using [C I] lines at z ' 1 3 (Lelli et al. 2018;Dye et al. 2022;Gururajan et al. 2022), the [C II] line at z ' 4 7 (De Breuck et al. 2014;Jones et al. 2017;Smit et al. 2018), and high-J CO lines (Tadaki et al. 2017;Talia et al. 2018).
The first IFU surveys of massive galaxies (M ? & 10 10 M ) at z ' 1 3 suggested that about one-third of star-forming galaxies were rotation-dominated disks, one-third were dispersiondominated objects, and another third were merging systems (Förster Schreiber et al. 2009;Gnerucci et al. 2011). Subsequent observations, however, led to a drastically di↵erent picture: it is now clear that at least 80% of star-forming galaxies at z ' 1 3 have rotating gas disks, whereas dispersion-dominated and merging systems constitute a minority of the star-forming population (Wisnioski et al. 2015(Wisnioski et al. , 2019;;Stott et al. 2016). The reasons for such a changing view are complex, but an important role has been played by the limited spatial resolution and the resulting beam-smearing e↵ects (Warner et al. 1973;Bosma 1978;Begeman 1989). When a rotating disk is spatially resolved with only a few resolution elements, di↵erent line-of-sight velocity projections are flux-averaged within the resolution element creating two main observational e↵ects: (1) the line profiles are artificially broadened, so a rotation-supported disk may appear as a dispersion-dominated object when observed at low spatial resolution, and (2) the line profiles become asymmetric, typically with long tails of emission toward the systemic velocity, so the gas kinematics may appear more complex than they really are. These two observational e↵ects can then lead to an artificially high fraction of dispersion-dominated and merging systems with respect to rotation-dominated ones. The same issues may occur in low-resolution [C II] surveys of galaxies at z ' 4 7 (Le Fèvre et al. 2020;Jones et al. 2021;Neeleman et al. 2021).
While there is now overall agreement about the existence of rotating disks at high z, their kinematic properties remain debated. Kinematic studies of warm ionized gas led to the common view that high-z disks are more turbulent than their local analogs (Förster Schreiber et al. 2009;Lehnert et al. 2009;Gnerucci et al. 2011). The gas velocity dispersion ( V ) is thought to increase systematically with z, while the degree of rotation support (V rot / V ) decreases (Wisnioski et al. 2015;Stott et al. 2016;Übler et al. 2019). Low values of V rot / V would imply that the gas disk is not fully supported by rotation: if one aims to trace the circular velocity of a test particle in the equilibrium gravitational potential, the pressure support would need to be taken into account (the so-called asymmetric-drift correction). For example, Genzel et al. (2017) and Lang et al. (2017) found that H↵ rotation curves decline in the outer region steeper than the Newtonian expectation for a thin disk, and argued that such a super-Keplerian decline is due to pressure gradients from strong, turbulent gas motions.
Kinematic studies of cold neutral gas from ALMA are painting a di↵erent picture. High-resolution observations of CO, [C I] and [C II] lines revealed the existence of gas disks with low turbulence and high V rot / values at z ' 2 7 (Lelli et al. 2018(Lelli et al. , 2021;;Kaasinen et al. 2020;Fraternali et al. 2021;Xiao et al. 2022;Posses et al. 2023), so the observed rotation curves do not require corrections for pressure support. Dynamically cold gas disks have also been found in gravitationally lensed galaxies, in which the e↵ective spatial resolution is exceptionally high (⇠100-500 pc) due to the magnification e↵ect (Di Teodoro et al. 2018;Rizzo et al. 2020Rizzo et al. , 2021;;Dye et al. 2022).
The cause of the disagreement between optical and submillimeter studies may be twofold. On the one hand, beam-smearing e↵ects and their modeling may again play a role: even when a rotating disk is clearly recognized, the finite spatial resolution has nontrivial e↵ects in assessing its kinematic properties. The artificial broadening of the emission line can lead to severe overestimates of the gas velocity dispersion, while the artificial skewness can lead to severe underestimates of the rotation velocity (e.g., Swaters et al. 2009;Di Teodoro & Fraternali 2015;Di Teodoro et al. 2016). For example, one may derive a slowly rising rotation curve even though the intrinsic rotation curve rises steeply (e.g., Lelli et al. 2010). In addition to getting observations with the highest possible spatial resolution, the best approach to account for beam-smearing e↵ects is modeling the three-dimensional (3D) emission-line datacube in a nonparametric fashion (e.g., Corbelli & Schneider 1997;Sicking 1997;Swaters 1999;Fraternali et al. 2001Fraternali et al. , 2002;;Gentile et al. 2003Gentile et al. , 2004Gentile et al. , 2007;;Józsa et al. 2007;Swaters et al. 2009;Lelli et al. 2010Lelli et al. , 2012aLelli et al. ,b, 2014b;;Di Teodoro & Fraternali 2015).
On the other hand, the di↵erent gas phases probed by di↵erent emission lines may also play a role. At z ' 0, ionized gas disks tend to display higher velocity dispersions and more complex noncircular motions than atomic and molecular gas disks (e.g., Lelli 2022). A similar behavior may occur at high z. The velocity dispersions from CO lines, indeed, appear systematically smaller than those from the H↵ line (Übler et al. 2019). In turn, H↵ velocity dispersions seem to be smaller than those from the [O II] line and stellar absorption lines (Übler et al. 2022).
To shed new light on galaxy dynamics at cosmic noon, we present an in-depth study of two galaxies with kinematic and photometric data of the highest quality. We used the Atacama Large Millimeter/submillimeter Array (ALMA) to obtain high-resolution data of the CO(3 2) and CO(4 3) lines for a main-sequence galaxy at z ' 2.24: zC-400569 (Mancini et al. 2011) also known as COSMOS 0488950 (Capak et al. 2007;Mobasher et al. 2007) or COSMOS2015 0300906 (Laigle et al. 2016). This galaxy is one of the best observed object from the SINS IFU survey (Förster Schreiber et al. 2009, 2018;Genzel et al. 2017). Luckily, next to zC-400569, there is another main-sequence galaxy at z = 1.47 (COSMOS 0488879 or COSMOS2015 0301356) so we set up our ALMA observations to simultaneously target the CO(2 1) and CO(3 2) lines of this second object (see Fig. 1). For brevity, we refer to this second galaxy as zC-488879. Essentially, we got two birds with one stone. Liu et al. (2019a) provide the stellar mass (M ? ) and starformation rate (SFRs) of both galaxies by fitting their spectral energy distributions (SEDs), combining optical and NIR photometry from COSMOS with the ALMA submillimeter continuum. According to their work, zC-400569 has M ? = 2⇥10 11 M and SFR = 81 M yr 1 while zC-488879 has M ? = 5 ⇥ 10 10 M and SFR = 115 M yr 1 . Thus, both galaxies lie on the starformation main sequence within the observed scatter (cf. with Schreiber et al. 2015).
Throughout this paper, we assume a flat ⇤ cold dark matter (⇤CDM) cosmology with H 0 = 67.4 km s 1 Mpc 1 , ⌦ m = 0.315, and ⌦ ⇤ = 0.685 (Planck Collaboration VI 2020). In this cosmology, 1 arcsec corresponds to 8.45 kpc at z = 2.24 and 8.68 kpc at z = 1.47. The age of the Universe and the light travel time are, respectively, 2.9 Gyr and 10.9 Gyr at z = 2.24, and 4.3 Gyr and 9.5 Gyr at z = 1.47.
ALMA observations were obtained through two di↵erent projects (2017.1.01020.S and 2019.1.00862.S; PI: T. Bisbas): the first one provides low-resolution observations in band 4, while the second one is a high-resolution follow-up in both bands 4 and 3. The two galaxies (zC-400569 and zC-488879) were simultaneously observed in the same field of view, pointing toward the center of zC-400569. In both bands, we used a mixed spectral setup with four spectral windows (SPWs) with a bandwidth of 1.875 GHz each. For band-3 observations, one SPW was used to cover the CO(2 1) line of zC-488879, another one for the CO(3 2) line of zC-400569, and the remaining two SPWs for the 3-mm continuum. For band-4 observations, the four SPWs cover the CO(3 2) line of zC-488879, the CO(4 3) line of zC-400569, the [C I](1 0) line of zC-400569, and the 2-mm continuum. All line SPWs were covered with 1920 channels giving a channel width of 4 5 km s 1 , while the continuum SPWs were covered with 128 channels giving a channel width of 60 100 km s 1 .
In project 2017.1.01020.S, band-4 observations used the nominal C43-3 array configuration with 41-45 antennas, providing minimum and maximum baselines of 14 and 740 m, respectively. Three execution blocks were obtained in May 2018 giving a total on-source integration time of about 2.5 h. In project 2019.1.00862.S, bands-3 and -4 observations used the nominal C43-6 configuration with 40-45 antennas, providing minimum and maximum baselines of 14 and 3638 m, respectively. Five A106, page 2 of 21 Fig. 1. Galaxies zC-488879 at z ' 1.47 (left) and zC-400569 at z ' 2.24 (right). The HST/WFC3 image in the NIR F160W filter (gray scale) is overlaid with CO(2 1) emission (green contours), CO(3 2) emission (blue contours), and CO(4 3) emission (cyan contours). Contours are at (3,6,12,24) map with the values of 3 map given in Table 1. The band 4 observations have lower quality than the band 3 ones, so the di↵erent spatial extent of CO(3 2) and CO(4 3) emission in zC-400569 may be due to sensitivity rather than being physical. Notes. All cubes have channel width of 15 km s 1 and velocity resolution of 30 km s 1 after Hanning smoothing.
executions blocks were obtained in band 4 in May 2021, while six execution blocks were obtained in band 3 between May and July 2021. Unfortunately, the execution blocks taken on 11 May 2021 (band 4) and 26 May 2021 (band 3) did not reach the expected quality and are therefore not used in this study. The final on-source integration time is ⇠4 h in band 3 and ⇠3 h in band 4.
The data reduction was performed using the Common Astronomy Software Applications (Casa) package (McMullin et al. 2007;The CASA Team 2022). The Fourierplane data were flagged and calibrated using the appropriate pipeline version provided by the ALMA team, which varies for di↵erent observing runs. The imaging was performed with the tclean task in Casa (version 6.1.2.7) using a Briggs' robust parameter of 1.5 and interactive cleaning. Emission-line cubes were derived using a channel width of ⇠15 km s 1 to have adequate signal-to-noise ratio (S/N). In the individual SPWs, the continuum emission is undetected, so no continuum subtraction was performed. For our kinematic analysis (Sect. 4), we use cubes without primary-beam correction, so the noise structure is uniform and well behaved. A primary-beam map is used to correct moment-one maps and calculate total fluxes (as we describe below) but the flux correction is entirely negligible for zC-400569 and only 1-2% for zC-488879 because both galaxies are well within the half-power beam-width (HPBW) of 41 00 in band 3 and 62 00 in band 4.
After imaging, the emission-line cubes were analyzed using the 3D Barolo software (version 1.6;Di Teodoro & Fraternali 2015). To further boost the S/N, the cubes were Hanning smoothed over three channels using the Smoothspec task, giving a velocity resolution of 30 km s 1 and the rms noise ( cube ) listed in Table 1. Subsequently, moment maps were obtained with the Makemask task, considering the signal inside a Boolean mask that was created with the Smooth & Search task. When using a mask, the noise in the moment-zero map ( map ) varies from pixel to pixel. Following Verheijen & Sancisi (2001) and Lelli et al. (2014c), we build a S/N map considering channel dependencies, then we computed a pseudo 3 map value taking the median intensity of pixels with S/N between 2.9 and 3.1. Total line fluxes are measured summing pixels with S /N > 3 in the moment-zero maps; uncertainties are estimated repeating the sum for pixels with S /N > 2 and S /N > 4. Emission-line maps are discussed in Sect. 3, but we anticipate that our kinematic analysis in Sect. 4 fits directly the 3D cubes, modeling the e↵ects of beam smearing on the emission-line profiles at each spatial location. Table 1 summarizes the properties of the ALMA data.
A106, page 3 of 21
The SINFONI data of zC-400569 were obtained from the archive of the SINS/zC-SINF Adapative Optics (AO) survey 1 (Förster Schreiber et al. 2018). We first add the WCS coordinate system to the cube and convert the third axis unit from µm to velocity assuming the H↵ redshift of 2.2405 (Förster Schreiber et al. 2018). To align the cube with the north direction, we perform a rotation of 30 degrees counterclockwise at the kinematic center of the galaxy and resample the data with a pixel size of 0.05 00 , the same as its original value. Finally, we checked the astrometry by overlaying the momentzero map from the SINFONI cube with images from the Hubble Space Telescope (HST), which are discussed in Sect. 2.3. The full-width half-maximum (FWHM) of the point spread function (PSF) of the SINFONI cube is ⇠0.32 00 , while the FWHM of the line spread function is V ' 87 km s 1 ( inst ' V /2.35 ' 37 km s 1 ) around the H↵ line (Förster Schreiber et al. 2018).
At the spatial location of the central galaxy, continuum emission is detected. The continuum emission is noisy and no stellar absorption line is detected. To estimate the continuum flux and subtract it from the SINFONI cube, we fit a first order polynomial to the line-free channels, excluding spectral regions with H↵ and [N II] emission. We explored the use of higher-order polynomials but they did not significantly improved the continuum subtraction. After continuum subtraction, we checked that the flux ratios of the two [N II] lines are consistent with the theoretical value within the errors. Finally, we produced a subcube that covers only the H↵ emission region. This continuumsubtracted H↵-only cube will be used in the rest of our analysis.
Both zC-400569 and zC-488879 have been detected by HST in various filters. Since we aim to model the stellar mass distribution of these galaxies, we analyzed the reddest HST images available in the Hubble Legacy Archive 2 , which were taken with the NIR F160W and F110W filters of the Wide Field Camera 3 (WFC3; Project ID 12578). The FWHM of the PSF in the F110W and F160 images are 0.13 00 and 0.15 00 , respectively. At the redshift of zC-400569 (z ' 2.24), the F110W filter is comparable to rest-frame Cousins U band, while the F160W filter is inbetween rest-frame V and B bands. At the redshift of zC-488879 (z ' 1.47), the F110W and F160W filters are comparable to rest-frame B and R bands, respectively.
The F160W image is shown in Fig. 1 overlaid with the CO moment-zero maps derived in Sect. 2.1. Overall, the HST astrometry is in good agreement with the ALMA one. For zC-400569 the peak emission in the HST image is slightly to the southwest of the centroid of the innermost CO contour, so one may have the visual impression that the HST astrometry needs to be rectified by 2-3 pixels to the top-left direction. On the other hand, for zC-488879 the peak emission in the HST image is slightly to the northeast of the centroid of the innermost CO contour, so one may have the visual impression that the HST astrometry needs to be rectified by 2-3 pixels in the opposite direction. Distortions in the HST camera cannot be so large, so the minor o↵-sets between stellar and CO distributions are not driven by technical issues and we apply no corrections to the HST astrometry.
1 https://www.mpe.mpg.de/ir/SINS/SINS-zcSINF-data 2 https://hla.stsci.edu/ We performed standard surface photometry using the Archangel softwarefoot_0 (Schombert 2011). Archangel has been widely used on local galaxies of all morphological types, from giant ellipticals to low-surface brightness dwarf galaxies, using images from various telescopes operating from NIR to ultra-violet bands (Tully et al. 2009(Tully et al. , 2013;;Schombert 2011Schombert , 2016Schombert , 2018;;Schombert & Smith 2012;Beygu et al. 2017;Greene et al. 2019). In short, after frame cleaning and sky determination, Archangel fits elliptical isophotes to the images and derives azimuthally averaged surface brightness profiles (see previous references for details). The resulting profiles are shown in Fig. 2 for the F160W filter. The profiles from the F110W filter show a similar shape but they have larger errors, so they will not be used in the rest of this paper.
For zC-400569, the surface brightness profile is well fitted by a classic bulge plus disk model, in which the bulge is described by a de Vaucouleurs (1948) profile and the disk by an exponential profile (Freeman 1970). The model has four free parameters: the e↵ective radius and e↵ective surface brightness of the bulge (R e, bul and µ e, bul ) and those of the disk (R e, disk and µ e, disk ). For zC-400569, the same model does not provide satisfactory results because the surface brightness profile displays a curvature at 1 00 < R < 2 00 that is reminiscent of "lenses" and "barlenses" in disk galaxies at z ' 0 (Laurikainen et al. 2011;Athanassoula et al. 2015). Thus, we add a lens-like component using a Sérsic (1963) profile with n = 0.2 and two further free parameters: R e, lens and µ e, lens . Table 2 provides the results of our fits, which were performed using the orthogonaldistance-regression algorithm in the SciPy package of Python (Virtanen et al. 2020). For both galaxies the bulge parameters are very uncertain (especially R e, bul ) because most of the bulge light is contained within the PSF (⇠0.15 00 or ⇠1.3 kpc for the adopted cosmology). To properly calculate the stellar gravitational potential, therefore, the bulge profile needs to be extrapolated below the HST spatial resolution (see Sect. 5).
3.1. The galaxy zC-400569 at z ' 2.4
Figure 3 compares the distributions of molecular and ionized gas in zC-400569 at z ' 2.24. The H↵ emission display a "clumpy" morphology with a northsouth orientation. The two brightest H↵ components to the north are encompassed by the CO emission, whereas the two weakest H↵ components to the south are undetected in any CO line. In these two components, which we name "A" and "B", the [N II] lines are undetected too (in agreement with Genzel et al. 2017). These two "clumps" are also visible in the HST images. Most likely, they are low-mass, low-metallicity galaxies that are interacting with the main one, as previously suggested by Genzel et al. (2017).
Figure 4 provides an overview of all emission line data, showing total intensity maps, line-of-sight velocity maps, and position-velocity (PV) diagrams along the galaxy major axis (see Sect. 4 for details). The CO distribution appears smoother than the H↵ distribution, but this may possibly be an e↵ect of the lower spatial resolution. Any giant clump in molecular gas must be significantly smaller than the ALMA beam (3 4 kpc). The spatial distributions of CO(4 3) and CO(3 2) emission are in good agreement.
In terms of kinematics, CO and H↵ emissions form regularly rotating disks, as shown by velocity fields and PV diagrams. This implies that stellar feedback does not destroy the overall kinematic regularity of the gas disk, similarly to the situation in local spirals where stellar feedback mostly drives a kpc-scale gas circulation via galactic fountains (e.g., Marasco et al. 2019;Li et al. 2021a). The possible interaction with components A and B appears to have a negligible e↵ect on the overall regularity of the main disk, albeit some minor kinematic disturbances are visible toward the southeast of the velocity fields. We stress that details in the velocity fields should be interpreted with caution because of beam smearing e↵ects; the gas kinematics is studied in Sect. 4 using 3D models.
The [C I](1 0) emission seems to trace the same rotating disk of the other lines but the major axis PV diagram is asymmetric and noisy. Moreover, the [C I] systematic velocity appears o↵set with respect to that from the H↵ and CO lines. We were not able to perform a detailed kinematic analysis of the [C I] data, so we use them solely as a total H 2 mass tracer in Sect. 5.3.
Figure 5 provides an overview of the emission-line data of zC-488879 at z ' 1.47. Similarly to zC-400569, the CO distribution is smooth and does not show the large clumps that are often seen in rest-frame UV and H↵ images of galaxies at cosmic noon (e.g., Zanella et al. 2019). The spatial distributions of CO(3 2) and CO(2 1) emission are in overall agreement, but the CO(2 1) emission peaks near the galaxy center while the CO(3 2) emission peaks toward the south. This may be possibly related to varying CO line excitation conditions in the galaxy.
A106, page 5 of 21 Fig. 4. Overview of zC-400569 at z ' 2.24. Left panels: total intensity maps for di↵erent emission lines, indicated in the top-left corner. Contours are at (3,6,12,24) map with the 3 map values given in Table 1. The white star shows the kinematic center. The gray ellipse to the bottom-left corner represents the PSF. Middle panels: velocity fields. The thick contour corresponds to the systematic velocity (set to zero) and the thin contours are in steps of ±60 km s 1 . The white star shows the kinematic center and the dashed line shows the disk major axis. Right panels: position-velocity diagrams along the disk major axis. The blue colorscale shows the observed gas emission. Black contours range from 2 cube to 8 cube in steps of 2 cube . Red contour (at 2 cube and 4 cube ) shows the best-fit 3D kinematic model from 3D Barolo. The red dots show the best-fit sky-projected rotation curve. The vertical and horizontal dashed lines correspond to the galaxy center and systemic velocity, respectively. Both CO(2 1) and CO(3 2) lines indicate regularly rotating disks, as shown by velocity maps and PV diagrams. The CO(3 2) emission has an extension to the southwest, which is kinematically connected with the main disk and display a small velocity gradient. This extension is reminiscent of the lopsided phenomenon seen in nearby H I disks at z ' 0 (Baldwin et al. 1980;Sancisi et al. 2008), which may be due to minor mergers and/or gas accretion events. Deeper CO observations are necessary to clarify the nature of this outer extension.
We model the gas kinematics in a nonparametric fashion without a-priori assumption on the gas density profile, rotation curve shape, and gravitational potential. This is important because beam-smearing e↵ects depend on several unknown quantities, such as the intrinsic gas distribution and rotation-curve shape (Warner et al. 1973;Begeman 1989). Parametric 3D models assume predefined functions for the rotation-curve shape and the gas density profile, so they necessarily make implicit assumptions on the intrinsic strength of beam-smearing e↵ects. For example, an intrinsically flat rotation curve down to small radii leads to more severe beam-smearing e↵ects than an intrinsically rising rotation curve. As a result, parametric models may only partially account for beam-smearing e↵ects. In addition, parametric models do not provide an actual derivation of the rotation curve because a smooth shape is imposed, neglecting possible real features (such as bumps and wiggles) that are often seen in rotation curves at z ' 0 (e.g., Sancisi 2004;Lelli et al. 2016a). On the other hand, parametric models have the advantage of hav-ing less free parameters than tilted-ring fits and may represent the only viable method to model rotating disks that are resolved with less than 2-3 resolution elements (e.g., Bouché et al. 2015).
The task 3Dfit of 3D Barolo (Di Teodoro & Fraternali 2015) performs a tilted-ring modeling (Warner et al. 1973;Rogstad et al. 1974;Begeman 1989) fitting the 3D cube rather than 2D moment maps, so beam-smearing e↵ects on each single line profile are taken into account. The disk is divided in a set of rings, where each ring is characterized by five geometric parameters center's coordinates (x 0 , y 0 ), systemic velocity (V sys ), position angle (PA), and inclination angle (i) and five physical parameters surface density (⌃ gas ), vertical thickness (z 0 ), rotation velocity (V rot ), radial velocity (V rad ), and velocity dispersion ( V ).
We use rings with a width of p a ⇥ b/2, where a and b are the major and minor axes of the synthesized beam. We adopt a fully axisymmetric disk, so the surface density of each ring is directly computed from the observed moment-zero map using azimuthal averages (option Norm=Azim). For the vertical density distribution, we assume an exponential law (option Ltype=3) with a fixed scale height of 100 pc: the precise value of z 0 is expected to have virtually no e↵ects on our results considering the ALMA beam size of 3 4 kpc. We also set V rad = 0 because there is no strong evidence of radial motions in the velocity maps, such as a clear nonorthogonality between the kinematic major and minor axis (see, e.g., Lelli et al. 2012a,b for galaxies where V rad , 0). Thus, we are left with seven free parameters in each ring. The systemic velocity is set to zero at the fiducial redshift of the galaxies (given in Table 1) but is kept as a free parameter to account for minor adjustments in units of km s 1 .
For data with relatively low S/N and resolution (as for most high-z galaxies), automated fitting codes such as 3D Barolo A106, page 7 of 21 Lelli, F., et al.: A&A 672, A106 (2023) 1. hV rot i is the mean rotation speed; its uncertainty is estimated using Eq. ( 3) in Lelli et al. (2016b).
must be used with caution. Firstly, it is important to define a proper mask within which the residuals are evaluated. After various trials, we find that this mask should be conservative to avoid that pixels with low S/N bias the fit. To build such a mask, we use the task Smooth & Search with parameters Factor=1.5 (factor for spatial smoothing), Snrcut=4 (primary S/N threshold), Growthcut=3 (secondary S/N threshold when growing the initial mask), and Minchannels=3 (minimum number of channels for a detection to be accepted, considering that the data have been Hanning smoothed over 3 channels). Next, our fitting strategy consists of three main steps: setting the disk geometry (Sect. 4.2), measuring the gas velocity dispersion (Sect. 4.3), and tracing the rotation curve (Sect. 4.4).
We run 3Dfit on each cube leaving all seven parameters free. In this first fit, all pixels are uniformly weighted (option Wfunc=0).
After various trials, we find that some best-fit parameters may slightly depend on some of the initial estimates, so we provide sensible initial guesses: PA = 35 and i = 50 for zC-400569 and PA = 20 and V sys = 50 km s 1 for zC-488879. Initial estimates for all the other parameters are automatically estimated by 3D Barolo. Next, we determine the geometric parameters taking the median of the best-fit values of di↵erent CO lines across all rings; uncertainties are estimated as the median absolute deviation. Thus, we fix the same disk geometry for all emission lines using the CO kinematics. We do not model warps because there is no strong indication for them in any of the data cubes at the available resolution.
The [C I](1 0) cube of zC-400569 has too low sensitivity to robustly determine the geometric parameters, while the H↵ kinematics may bias the results due to the H↵-bright companions to the south. In addition, the H↵ line could su↵er from varying dust extinction across the disk, unlike the CO lines that trace the gas kinematics unbiased by such e↵ects. The center and systemic velocity of the H↵ cube, however, are determined independently from CO lines to account for possible o↵sets in the astrometry and/or spectral calibration of the SINFONI and ALMA data.
Table 3 summarizes our results. The geometric parameters from CO kinematics are approximately consistent with the restframe optical morphology of the galaxies. The inclination of zC-400569, however, is more uncertain than the formal errors from the median absolute deviation: the CO(3 2) line points to an inclined disk (60 65 ), while the CO(3 4) and optical morphologies point to a more face-on disk (45 50 ). Thus, we assign an error of 5 that encompasses these inclinations within 2 .
Having fixed the disk geometry, we use the task SpacePar in 3D Barolo to look for a global minimum in the V rot V space. The parameter space is explored in steps of 1 km s 1 , considering a range in V from 1 to 60 km s 1 and a range in V rot of 200 km s 1 centered on the expected rotation velocity. At each step, SpacePar builds a 3D disk model with given V rot and V , then computes the average sums of residuals r = |M D|, where M and D are the flux values at each 3D pixel of model and data, respectively (default option Ftype=2). All pixels are uniformly weighted (option Wfunc=0).
For all emission lines, SpacePar returns very low values of V of just 1 2 km s 1 . For example, Fig. 6 shows the results for the H↵ cube of zC-400569, excluding the innermost ring in which beam smearing e↵ects are most severe and the outermost ring in which the S/N is low. Similar plots for the CO cubes are shown in Appendix A. Values of V ' 1 2 km s 1 are physically acceptable but the minimum in the V rot V space is very shallow; an intrinsic velocity dispersion of ⇠10 km s 1 increases the 2 by a mere 5%. Most likely, the available spectral resolution ( V ), spatial resolution, and sensitivity do not allow us to robustly measure the intrinsic velocity dispersion.
One may wonder whether the low velocity dispersions are driven by high surface brighntess pixels, which could dominate the average sum of residuals r. We reran SpacePar calculating r = |M D|/|M + D|, which normalizes the residuals at each pixel to the local flux density (option Ftype=3). Again, we found formal best-fit values of V of just 1-2 km s 1 , indicating that this result is not driven by how we weight the residuals.
To further explore the situation and set a fiducial upper limit on V , we built mock-observed cubes of rotating disks with known V rot and V . We smoothed the mock cubes to the same spatial and spectral resolutions of our observations, and add Gaussian noise with a peak S/N of 6 to 8, similarly to our data. We then ran SpacePar on each mock-observed cube. Independently of the input value of V rot , when the intrinsic V is equal or smaller than the instrumental dispersion inst ' V /2.35, SpacePar pushes the best-fit value of V to the lower boundary of the explored parameter space. This is the same behavior that we obtain for all emission-line cubes of zC-400569 and zC-488879. When the intrinsic V is slightly higher than inst , instead, SpacePar does not push V to the lower boundary of the explored parameter space, albeit it may not exactly recover the correct V due to the low S/N. Thus, we conclude that inst is a sensible upper limit for the intrinsic velocity dispersion.
Given that the ALMA cubes have a spectral resolution V ' 30 km s 1 , the formal upper limit on the CO velocity dispersion ( CO ) is about 13 km s 1 . To be conservative, we adopt the channel width of 15 km s 1 as our upper limit on CO . Given the SINFONI spectral resolution V ' 86 km s 1 around the H↵ line (Förster Schreiber et al. 2009), we adopt an upper limit of 37 km s 1 on the H↵ velocity dispersion ( H↵ ). We stress that these upper limits are themselves uncertain: while the low spectral and spatial resolutions increase the line broadening leading to an overestimate of the intrinsic V , the low S/N may lead to an underestimate of the intrinsic V .
To check the robustness of our upper limits, we created observed cubes with a channel width of 5 km s 1 and spectral resolution of 10 km s 1 (after Hanning smoothing). These cubes have significantly lower sensitivity than the cubes at 30 km s 1 A106, page 8 of 21 resolution. In fact 3D Barolo does not detect the source in most of them apart from the CO(4 3) cube of zC-400569 and the CO(3 2) cube of zC-488879. For these two CO cubes, SpacePar provides CO ' 4 12 km s 1 , which may be considered as proper measurements because inst ' 4 km s 1 . Given the low S/N of these cubes, however, we prefer to be conservative and base all our conclusions on the upper limit CO < 15 km s 1 .
Fixing the velocity dispersion to its formal upper limit, we run 3Dfit using only V rot as a free parameter in each ring. This provides the observed rotation curve for each emission line. In this last fit, pixels are weighted according to the function cos 2 (✓), where ✓ is the azimuthal angle around the major axis (option Wfunc=2). Thus, pixels around the major axis are maximally weighted to maximize information on rotational velocities.
The right panels of Figs. 4 and 5 compare major-axis PV diagrams from the observed cubes with those from our best-fit 3D disk models. The best-fit rotation curves are projected on the PV diagrams as V rot sin(i). Overall, the 3D disk models give a good description of the observations. In particular, the observed thickness of the PV diagrams (the line broadening at each radius) is well reproduced, indicating that the assumed velocity dispersion is sensible. The bulk of the observed line broadening, indeed, is driven by beam-smearing e↵ects, not by the intrinsic gas velocity dispersion. Another known e↵ect of severe beam smearing is that the best-fit rotation velocities are not necessarily near the peak of the emission line along the major axis, but they are closer to the uppermost and lowermost velocity edges of the PV diagrams. Similar rotation curves, indeed, could be obtained using the "envelope-tracing method", which considers the terminal velocity near the lowest density contour in a PV diagram (e.g., Sofue & Rubin 2001). The envelope-tracing method was originally developed for edge-on galaxies (Sancisi & Allen 1979) but can be applied in poorly resolved galaxies too because multiple line profiles are flux-averaged within the PSF (similarly to the line-of-sight integration in edge-on disks), so the true rotation velocity corresponds to the terminal velocity near the edge of the emission, not the velocity near the peak of the emission.
The rotation curves are shown in Fig. 7. Reassuringly, different emission lines provide consistent results, suggesting that the various gas phases are kinematically settled. The observed rotation velocities imply V rot / CO & 17 for zC-400569 and V rot / CO & 22 for zC-488879, so there is no need to apply corrections for pressure support (asymmetric drift). Overall, the different gas tracers point to the same result: flat rotation curves with no clear signs of a decline in the outer parts.
3D Barolo provides asymmetric errors ( + , ) that correspond to a variation of 5% of the residuals from the global minimum. We compute symmetric 1 errorbars as p + , which will be used as weights in subsequent Bayesian fits (Sect. 5.2). For zC-400569, rotation velocities from the CO(4 3) line have the smallest uncertainties, while those from the H↵ have the largest ones due to the relative S/N of the data. For zC-488879, rotation velocities from the CO(2 1) and CO(3 2) lines have comparable uncertainties.
The H↵ kinematics of zC-400569 has been previously studied by the SINS/zC-SINF team using the Dysmal parametric code (Genzel et al. 2017(Genzel et al. , 2020;;Price et al. 2021;Nestor Shachar et al. 2023). The previous H↵ results are mixed in terms of both the intrinsic H↵ velocity dispersion and the rotation curve shape.
Regarding the H↵ velocity dispersion, we give a formal upper limit of 37 km s 1 while the SINS/zC-SINF team reports actual measurements: 34 km s 1 (Genzel et al. 2017, consistent with our upper limit), 45 km s 1 (Genzel et al. 2020), 71.82 +5.97 8.40 km s 1 (Price et al. 2021), and 58 ± 5 km s 1 (Nestor Shachar et al. 2023). We do not fully understand why the SINS/zC-SINF team finds significantly di↵erent velocity dispersions from the same IFU data, but these varying values are in line with our basic conclusion: the H↵ velocity dispersion of zC-400569 cannot be reliably measured with the existing data. This demonstrates that fitting methodologies can play a major role in obtaining measurements and/or upper limits of V from poorly resolved data with low S/N.
Regarding the H↵ rotation curve, our measurements and previous ones di↵er in both the inner and outer parts. In the inner parts, the H↵ rotation curve from Dysmal fits (Genzel et al. 2017(Genzel et al. , 2020;;Price et al. 2021;Nestor Shachar et al. 2023) rises until R ' 0.4 00 whereas our rotation curve is already flat. As far as we understand, the previous H↵ works show the beam-smeared rotation curve, not the intrinsic one such as that provided by 3D Barolo (Fig. 7), so the di↵erences in A106, page 9 of 21 the inner rotation velocities can be ascribed to beam-smearing e↵ects (e.g., Lelli et al. 2010). In the outer parts, the H↵ rotation curve from Dysmal 1D fits (Genzel et al. 2017(Genzel et al. , 2020;;Nestor Shachar et al. 2023) shows a mild decline, whereas that from Dysmal 2D fits (Price et al. 2021) keeps rising until the last measured point at R ' 1 00 . Our 3D fits do not confirm any of these trends, but give a flat rotation curve out to R ' 1 00 . These di↵erences are likely driven by the H↵-bright companions to the south (Fig. 3) that may hinder a precise determination of the PA, as also pointed out by Price et al. (2021). The CO data do not su↵er from these complications and give a PA of 37.5 ; the resulting major-axis PV diagrams (Fig. 4) unambiguously show that the rotation curve is relatively flat. In fact, independently of fitting codes, one would obtain the same result by tracing the rotation curve by eye using a "human neural network" (e.g., Sancisi & Allen 1979;Verheijen & Sancisi 2001).
The di↵erent measurements of V and V rot from di↵erent studies imply di↵erent values of the ratio V rot / V , which determine whether pressure support is relevant or not. Di↵erently from previous studies, we conclude that the values of V rot / V are relatively high (&7 for H↵ and &17 for CO), so corrections for pressure support (asymmetric drift) are negligible. Thus, the observed rotation speeds trace the circular velocity of a test particle in the equilibrium gravitational potential and can be safely used to build mass models. In the case of zC-400569, the value of V plays the key role because the absolute di↵erences in V rot from di↵erent studies are relatively small, apart from the overall di↵erence in rotation-curve shape.
We fit the rotation curves with mass models that includes di↵erent gravitational contributions. We start with mass models that contain only baryonic components (Sect. 5.3), then mass models that also include a dark matter (DM) halo (Sect. 5.4), and finally mass models in the context of Milgromian dynamics (MOND; Sect. 5.5). In this section we describe the calculation of the Newtonian baryonic contributions: V bul from the stellar bulge, V disk from the stellar disk, and V gas from the cold gas disk.
The contribution of each baryonic component is computed using the task Rotmod in the Gipsy software (Vogelaar & Terlouw 2001). The total baryonic contribution is then given by
where ⌥ i = M i /(10 10 M ) are dimensionless factors. Basically, for numerical convenience, the gravitational contribution of each component is computed for an arbitrary mass of 10 10 M and rescaled using fitting parameters ⌥ i on the order of one. The total baryonic mass of the model is then given by M bar = (⌥ bul + ⌥ disk + ⌥ gas ) ⇥ 10 10 M .
(2)
The gravitational contribution of the stellar bulge is computed assuming spherical symmetry, deprojecting a given 2D surface density profile (Kent 1986). We adopt a de Vaucouleurs (1948) profile for which R e, bul is fixed from fits of the observed surface brightness profile (Sect. 2.3). Thus, the only free parameter is M bul or equivalently ⌥ bul . Importantly, it is necessary to extrapolate the bulge profile at radii smaller than the HST spatial resolution (⇠1.3 kpc) otherwise the innermost points of the rotation curve would not be correctly reproduced. Di↵erences in V bul between spherical and oblate geometries are on the order of 10% 20% (Noordermeer 2008) and degenerated with M bul , which is much more uncertain than that.
The gravitational contribution of the stellar disk is computed considering a disk of finite thickness with a given density profile ⇢(R, z) = ⌃(R)⇣(z) (Casertano 1983). For the radial density distribution ⌃(R), we use the best-fit model to the observed surface brightness profile after subtracting the bulge component (Sect. 2.3). Thus, for zC-400569 we adopt a pure exponential disk, while for zC-488879 the disk term includes an exponential plus a barlens-like component (see Fig. 2). Di↵erently from A106, page 10 of 21 mass modeling at z ' 0 (Lelli et al. 2016a), we use parameteric models for the radial density profile because they need to be extrapolated at R < 1.3 kpc due to limited spatial resolution of the available images. For the vertical density distribution ⇣(z), we assume an exponential function with constant scale height z disk . We estimate z disk using the scaling relation z disk = 0.196 (R e,disk /1.68) 0.633 that holds for edge-on disk galaxies at z ' 0 (Bershady et al. 2010). The vertical geometry has a small e↵ect on the resulting V disk : at fixed M disk a thick disk gives a smaller velocity contribution than a thin disk at small radii (Casertano 1983). The uncertainties on M disk , however, are significantly larger than plausible variations in z disk , so there is little value in leaving z disk free in the fit.
The gravitational contribution of the gas disk is computed in a similar fashion as the stellar disk. For the radial density distribution, we use the observed CO(3 2) profile for zC-400569 and the observed CO(2 1) profile for zC-488879 because they are the lowest CO transitions available, which are expected to best trace the H 2 surface density distribution irrespective of the ambient conditions (such as gas density and kinetic temperature). For the vertical density distribution, we assume an exponential profile with a constant scale height of 100 pc, which is reasonable for a cold gas disk with V rot / V & 20. Our computation of V gas neglects the contribution of atomic gas (H I), which may dominate the total gas budget of high-z galaxies (Chowdhury et al. 2022). In the inner regions of local spirals, however, the H 2 surface densities usually dominate over the H I surface densities, so the gas gravitational contribution is mostly due to H 2 at R . 5 10 kpc while H I prevails at larger radii (Martinsson et al. 2013;Frank et al. 2016). According to theoretical models of the interstellar medium (ISM), the transition from atomic to molecular gas depends mostly on local properties (such as pressure, metallicity, and far-UV radiation) and is expected to be spatially abrupt once the H I-to-H 2 phase transition criterion is satisfied, leading to a nearly fully molecular ISM (Elmegreen 1989(Elmegreen , 1993;;Papadopoulos et al. 2002;O↵ner et al. 2013). Given that we extract rotation curves using CO lines, it is safe to assume that these CO-bright regions are well within the H I-to-H 2 phase transition radius, thus dominated by molecular gas (Bisbas et al. 2021).
The parameters of the mass models are determined using a Markov-chain-Monte-Carlo (MCMC) method in a Bayesian context. We define the likelihood L = exp( 0.5 2 ) with
where V rot is the observed rotation velocity at radius R k , V rot is the associated error, and V mod is the model rotation velocity that depends on the fitting parameters p. As it is often the case in Astronomy, the errors V rot are "educated guesses" and do not truly represent a strict 1 deviation from a Gaussian distribution. However, scaling all V rot by an arbitrary factor would only a↵ect the width of the posterior probability distribution but not its overall shape. This implies that we can robustly determine the best-fit model parameters p, but the associated uncertainties should be taken with a grain of salt. For each galaxy, we consider the ensemble of rotation velocity measurements from di↵erent emission lines at the same time (Fig. 7). This increases the statistical significance of our data set, allowing us to more robustly infer the optimal model parameters and associated uncertainties.
The posterior probability distributions of the model parameters are mapped using emcee (Foreman-Mackey et al. 2013). The MCMC chains are initialized with 200 walkers. We run 1000 burn-in iterations, then the sampler is run for another 2000 iterations. The emcee parameter a, which controls the size of the stretch move, is set equal to 2. In general, this gives acceptance fractions larger than 50%.
Both V rot and V rot depend on the disk inclination i. The value of i is constrained by kinematic fits to the emission-line cubes with uncertainties of a few degrees (Sect. 4.2). Thus, we treat i as a nuisance parameter using a Gaussian prior with central value i 0 and standard deviation i 0 set equal to the kinematic estimates from 3D Barolo. Then, V rot and V rot transform as sin(i 0 )/ sin(i).
Fitting rotation curves can be a strongly degenerate problem (e.g., Li et al. 2019Li et al. , 2020Li et al. , 2021b)). For both galaxies, indeed, the velocity contributions V gas and V disk display a similar trend with radius, so the values of ⌥ gas and ⌥ disk are strongly degenerate and essentially unconstrained when left entirely free. To break this degeneracy, we impose two physically motivated priors: 1. A log-normal prior on ⌥ gas using order-of-magnitude estimates of M gas from CO luminosities. Adopting the average line ratios R 21 = L CO(2 1) /L CO(1 0) ' 0.85 and R 31 = L CO(3 2) /L CO(1 0) ' 0.60 (e.g., Bisbas et al. 2021), we find M gas ' 10 10 10 11 M assuming either a starburst-like conversion factor ↵ CO ' 0.4 M (K km s 1 pc 2 ) 1 or the Milky Way ↵ CO ' 4.3 M (K km s 1 pc 2 ) 1 . Thus, we center the prior at log(⌥ gas ) = 0.5 with a standard deviation of 0.5 dex. 2. A log-normal prior on µ gas = ⌥ gas /⌥ disk , the ratio between gas mass and stellar mass in the star-forming disk. Considering the results of Tacconi et al. (2018) and Liu et al. (2019b) for galaxies at z ' 1.5 2.5, we center the prior at log(µ gas ) = 0 with a standard deviation of 0.3 dex. These priors are used in all the mass models described in the following sections. The fitting results are summarized in Table 4.
Figure 8 (left panels) shows mass models considering only baryonic components, so V mod = V bar (⌥ bul , ⌥ disk , ⌥ gas ) in Eq. ( 3). Given the overwhelming evidence for the DM e↵ect in the Universe, these mass models are not entirely physical but provide hard upper limits on stellar and gas masses; they are analogous to "maximum disk" mass models of local galaxies (van Albada et al. 1985;Starkman et al. 2018). For both galaxies, the rotation curve is well reproduced with no need of DM. Similarly to local spiral galaxies, the bulge dominates the inner galaxy regions, while the stellar and gas disks become important at larger radii. Di↵erently from local galaxies, however, the molecular gas disk dominates over the stellar disk. This result is partially driven by our assumed priors (see Sect. 5.2): the relative contributions of V gas and V disk would be nearly unconstrained if left entirely free. The ratio M bul /M bar , instead, is a robust quantity because V bul declines with radius while V disk and V gas rise. For zC-400569 we find M bul /M bar ' 0.2, comparable to late-type spirals (Sc or Sb), while for zC-488879 we find M bul /M bar ' 0.5, comparable to early-type disks (Sa or S0).
We compare the stellar masses (M ? = M bul + M disk ) from our dynamical model with those independently measured from SED fitting using automated pipelines (Liu et al. 2019a). For zC-400569, we find M ? = (7.8 +2.3 2.8 ) ⇥ 10 10 M that is consistent with the SED value of (21.9 ± 6.5) ⇥ 10 10
M at the 2 level. For zC-488879, we find M ? = (21.1 +3.6 4.2 ) ⇥ 10 10 that is significantly Notes. For the uncertainties, see Appendix B because several posterior-probability distributions are not Gaussian, so formal ±1 errors can be misleading.
Fig. 8. Mass models for zC-400569 (top panels) and zC-488879 (bottom panels) using only baryonic components (left), baryons plus a NFW halo (middle) and baryons in a MOND context (right). In all panels, the observed rotation curve (dots with errorbars) is fitted considering the gravitational contributions of a gas disk (green dotted line), a stellar disk (blue dashed line), a stellar bulge (red dash-dotted line, if present), and a DM halo (magenta dash-dotted line, if present).
di↵erent from the SED value of (5.2 ± 1.2) ⇥ 10 10 M at the 3 level. The high stellar mass of this galaxy is dominated by the bulge component, which is unavoidably needed to explain the inner rotation speeds of ⇠300 km s 1 . This indicates that stellar masses from automated SED fitting may occasionally be uncertain up to a factor of 4.
In principle, stellar masses from SED fitting and/or galaxy colors may be improved by considering the bulge and stellar disk separately (e.g., Schombert et al. 2022). This requires spatially resolved images that are currently available only in the rest-frame UV and optical parts of the spectrum (see Sect. 2.3), which are most sensitive to the unavoidable assumptions on star-formation history, chemical enrichment, and dust extinction (e.g., Schombert et al. 2019). The situation may be improved with JWST images, probing the rest-frame NIR part of the spectrum, which is less sensitive to those assumptions and may therefore provide accurate stellar masses for the bulge and disk components separately.
As a consistency check, we compute an e↵ective ↵ CO = M gas /L CO(1 0) in units of M (K km s 1 pc 2 ) 1 , using the gas A106, page 12 of 21 mass from our dynamical model and the observed CO(2 1) and CO(3 2) luminosities with the assumption of average R 21 = 0.85 and R 32 = 0.60. We find ↵ CO = 2.1 +0.7 0.8 for zC-400569 and ↵ CO = 2.9 +0.7 0.9 M for zC-488879. Both values are consistent with the Milky-Way ↵ CO within the uncertainties and are within the range of values found in star-forming galaxies across cosmic time (see Table 14 in Dunne et al. 2022). For zC-400569, we make a similar consistency check for ↵ [C I] = M gas /L [C I](1 0) using the [C I](1 0) emission, which is too faint for rotationcurve measurements but is a good global tracer of the H 2 mass (Papadopoulos et al. 2004). We find ↵ [C I] = 6.9 +1.9
2.4 which is in tension at more than 4 with the latest calibration of 17.0 ± 0.3 from Dunne et al. ( 2022) using a sample of 407 galaxies across cosmic time. Our value, however, is very close to ↵ [C I] = 7.3 found by Crocker et al. (2019) in 18 nearby galaxies using spatially resolved, multitransition CO and [C I] observations. It is also consistent with theoretical expectations from 3D astrochemical simulations (Bisbas et al. 2021).
Figure 8 (middle panels) shows mass models including a DM halo. From a statistical perspective, adding a DM halo with additional free parameters is not necessary because the rotation curves are already well fitted by the baryon-only mass model. The following models, however, allow us to check whether the existing data are consistent with expectations from ⇤CDM cosmology.
We assume a spherical DM halo with a Navarro-Frenk-White (NFW) density profile (Navarro et al. 1996), which has two free parameters: the halo concentration C 200 and the halo mass M 200 (or equivalently the halo velocity V 200 ). These quantities are defined in the same way as in Li et al. (2020). Since the observed rotation curve can be fully explained by baryons, the halo parameters cannot be constrained using only the kinematic data. Thus, we impose two ⇤CDM scaling relations as Bayesian priors (following Li et al. 2020) M 200 relations derived by Legrand et al. (2019) at z = [1.1 1.5] and z = [2.0 2.5] on zC-488879 and zC-400569, respectively. Similarly, the mass-concentration relations from Dutton & Macciò (2014) at z ' 1 and z ' 2 are imposed on zC-488879 and zC-400569, respectively. We also consider the redshift-dependent scatters given by Legrand et al. (2019) on log(M 200 ) and by Dutton & Macciò (2014) on log(C 200 ). Both scatters are on the order of 0.1 dex.
With respect to the baryons-only model, the MCMC fits decrease the gravitational contributions of stellar and gas disks to leave room for the DM contribution. The bulge contribution, instead, is almost unchanged given its characteristic declining shape, so the value of M bul /M bar increases for both galaxies. For zC-400569, the best-fit masses decrease by a factor of ⇠1.7 in stars and ⇠2.1 in gas, while for zC-488879 they decrease by a factor of ⇠1.4 in stars and ⇠3.1 in gas. The actual uncertainties on stellar and gas masses of high-z galaxies are surely as large as a factor of 3, so we cannot categorically rule out the scenario in which the DM halo is dynamically important in the inner galaxy regions (R . 8 9 kpc). Moreover, the best-fit parameters lie on the M ? M 200 and M 200 C 200 relations (imposed as priors) within the uncertainties. Thus, we conclude that the observed rotation curves are consistent with "cuspy" ⇤CDM halos due to a severe disk-halo degeneracy, as we discuss in Sect. 6.5.
Milgromian dynamics (MOND;Milgrom 1983a,b,c) is the major alternative to particle DM. The MOND paradigm modifies the laws of gravity and/or inertia when accelerations are smaller than an acceleration scale a 0 ' 10 10 m s 2 , which is typical for the outer parts of spiral galaxies (see Famaey & McGaugh 2012;Milgrom 2014;Banik & Zhao 2022, for reviews). The rotation curves of the two galaxies studied here are limited to the inner high-acceleration regions where V 2 obs /R > 3 4a 0 , so they are expected to probe the Newtonian regime of the theory. Still, it is important to test MOND in these high-z galaxies for various reasons: (1) the MOND acceleration scale may possibly vary with cosmic time; (2) the ⇤CDM relation between angular distance (D A ) and redshift may not apply in a MOND cosmology, so rotation curve fits of high-z galaxies may provide empirical constraints on the D A (z) relation in MOND; (3) the MOND interpolation function, linking the Newtonian and Milgromian regimes, may not be constant with cosmic time.
Constraining the redshift evolution of a 0 is particularly important (Milgrom 2017). Empirically, it is known that a 0 ' c ⇥ H 0 and a 0 ' c 2 p ⇤, where c is the speed of light and ⇤ is the cosmological constant. Both equivalences may be mere numerical coincidences, or they may have a deeper physical meaning. For example, the first coincidence may suggest that a 0 (z) = c ⇥ H(z), while the second one that a 0 does not vary with z.
We perform MOND fits assuming (1) the empirical value a 0 = 1.2 ⇥ 10 10 m s 2 found at z = 0 (e.g., Lelli 2022), (2) the theoretical D A z relation from ⇤CDM, which is expected to hold in some relativistic MOND theories (Skordis & Złośnik 2021), and (3) the interpolation function that best fits the empirical radial acceleration relation (RAR) at z ' 0 (McGaugh et al. 2016;Lelli et al. 2017). Figure 8 (right panels) shows the mass models from MOND fits. With respect to the baryon-only Newtonian mass models, the stellar and gas masses are decreased by 10% 30% due to a moderate MOND e↵ect at large radii. The resulting values of M ? and M gas are physically acceptable as discussed in Sect. 5.3. This implies that rotation curves at z ' 1 2 are compatible with the value of a 0 and the interpolation function measured at z ' 0.
6.1. Rotation versus pressure support in high-z galaxies In Sect. 4.3 we find that the CO velocity dispersion of both zC-400569 and zC-488879 is surprisingly low. We cannot robustly measure CO because the observed line broadening is fully dominated by the e↵ects of spatial and spectral resolution, so we have a fiducial upper limit of 15 km s 1 set by the spectral resolution. The resulting lower limits in the V rot / CO ratio are on the order of 17-22, so the molecular gas disk is fully supported by rotation and there is no need of corrections for pressure support. These properties are comparable to those of spiral galaxies at z ' 0: local CO disks typically have CO ' 5 15 km s 1 (Mogotsi et al. 2016;Bacchini et al. 2020) and V rot / CO ' 5 50, with lower V rot / CO values at small radii and higher V rot / CO values in the outer parts (see Appendix in Bacchini et al. 2020). In addition, ALMA observations at exceptionally high spatial resolutions (10 30 pc) revealed that local galaxies can host nuclear CO disks with extremely low velocity dispersions of just 1 2 km s 1 and V rot / CO & 100 (Davis et al. 2017(Davis et al. , 2018)). This fact hints at the possibility that, A106, page 13 of 21 even at z ' 0, velocity dispersions from CO data with moderate spatial (⇠500 pc) and spectral (⇠5 km s 1 ) resolutions may be inflated by observational e↵ects. Evidently, measuring the intrinsic velocity dispersion of a rotating disk is not a trivial exercise.
Figure 9 (left panel) shows the redshift evolution of V inferred by Übler et al. (2019) using both IFU data of ionized gas and NOEMA data of molecular gas. We also show the prediction of a toy model based on the Toomre-instability criterion, which was built by Wisnioski et al. (2015) to reproduce ionized gas data at z . 3.5. For comparison, we consider a small sample of galaxies that satisfy two quality criteria:
(1) high-resolution ALMA data of cold gas tracers such as CO, [C I], and [C II] lines, and (2) detailed 3D kinematic modeling for consistency with the current work. In addition to the two main-sequence galaxies studied here, this sample include an AGN-host galaxy at z ' 2.6 with [C I](2 1) data (Lelli et al. 2018) and several submillimeter galaxies at z ' 4 5 with [C II] data, considering both lensed (Rizzo et al. 2020(Rizzo et al. , 2021) ) and unlensed (Sharda et al. 2019;Lelli et al. 2021;Fraternali et al. 2021) sources. The velocity dispersions from ALMA data appear to be systematically below the expectations from IFU data of ionized gas as well as NOEMA data of molecular gas. Figure 9 (right panel) compares the V rot / V ratios of this sample with the model of Wisnioski et al. (2015). The V rot / V of galaxies with ALMA data are higher than expected from the extrapolation of the model, but comparable to those of rotation-supported spiral galaxies at z ' 0.
In our view, it is possible that previous measurements from IFU and/or NOEMA data may have not been fully corrected for beam-smearing e↵ects, leading to systematic overestimates of V and underestimates of V rot / V (see also Di Teodoro et al. 2016). In fact, as the data quality increases, the measured velocity dispersion decreases (see Fig. 6 in Rizzo et al. 2021). In addition, ionized gas tracers may be more easily "contaminated" by large-scale wind components than H 2 -only tracers, leading to an artificial increase of the measured velocity dis-persion (see the discussion in Lelli et al. 2018). The H I-to-H 2 phase-transition criterion, indeed, is much more readily satisfied in the high-pressure regions inside the star-forming disks than in the wind regions outside them (though in some AGN and/or starburst galaxies, CO-rich H 2 winds do exist). Moreover, a CO disk would correspond to lower disk scale-heights z disk and thus smaller V because in the vertical direction one has V ' z disk p 2⇡Gh⇢ 0 i, where h⇢ 0 i is the average mid-plane mass density. In conclusion, it could well be that high-z disks are not as turbulent as generally thought, but more high-resolution and high-sensitivity ALMA observations are needed to confirm or refuse this possibility.
In Sect. 4.4, we find that the rotation curves of zC-400569 and zC-488879 are flat with no sign of decline out to R ' 8 kpc. Different emission lines, such as CO(2 1), CO(3 2), CO(4 3) and H↵, provide consistent results (Fig. 7). This fact suggests that the various gas phases are kinematically settled and that the resulting rotation curve probes the equilibrium gravitational potential.
Flat rotation curves have been previously found in galaxies at z ' 1 using H↵ kinematics from IFU data (Di Teodoro et al. 2016;Sharma et al. 2021). Specifically, Di Teodoro et al. (2016) analyzed a carefully selected sample of 18 galaxies with highquality H↵ data, while Sharma et al. (2021) studied a statistical sample of 344 galaxies with data of variable quality. Both studies used nonparametric tilted-ring modeling, so they provide the actual rotation curve on a ring by ring basis, not a functional parametrization thereof. Using parametric models, instead, Genzel et al. (2017) studied six galaxies at z ' 0.8 to 2.4 (including zC-400569) and found declining rotation curves at large radii. Subsequently, Genzel et al. (2020) increased their sample to 41 galaxies (including CO data from NOEMA for seven of them) and found a variety of rotationcurve shapes, from rising, to flat, to declining. A similar variety A106, page 14 of 21 Fig. 10. Rotation curves of zC-400569 at z ' 2.24 (red dots) and zC-488879 at z ' 1.47 (orange stars) are compared with those of disk galaxies at z ' 0 from the SPARC database (Lelli et al. 2016a). The latter ones are color-coded by the mean stellar surface density within the e↵ective radius. Rotation curves at z ' 0 are much more extended than those at high z thanks to the availability of H I observations. in rotation-curve shapes is found by Bouché et al. (2022) for nine galaxies at z ' 1, fitting parametric models to [O II] observations from MUSE.
These di↵erent results may be driven by the di↵erent modeling techniques, but they could also indicate a genuine variety among galaxies. At z ' 0 it is well established that the inner rotation curves (R . 2 3R e ) display a wide variety of shapes, depending on the galaxy morphology, surface brightness, and bulge-to-disk ratio (e.g., Corradi & Capaccioli 1990;Casertano & van Gorkom 1991;de Blok et al. 1996;Tully & Verheijen 1997;Noordermeer et al. 2007;Swaters et al. 2009;Lelli 2022). This evidence is nicely summarized by the so-called Renzo's rule: "For any feature in the luminosity profile of a galaxy there is a corresponding feature in the rotation curve, and vice versa" (Sancisi 2004). Quantitatively, the inner steepness of the rotation curve, proxy of the central dynamical surface density, closely correlates with the central surface brightness, proxy of the central baryonic surface density (Lelli et al. 2013(Lelli et al. , 2014a(Lelli et al. , 2016c). Thus, it may not be surprising that a similar diversity in the inner rotation curves exists at cosmic noon, given the observed diversity in the light (or stellar mass) distribution of galaxies.
The vast majority (if not all) of rotation curves at z ' 0 become flat at very large radii, outside the bright stellar component of galaxies, which can be probed with H I observations (Bosma 1981;van Albada et al. 1985;Begeman et al. 1991;Sanders & Verheijen 1998). To illustrate the situation, Fig. 10 compares the rotation curves of zC-400569 and zC-488879 with those of local galaxies from the SPARC sample (Lelli et al. 2016a). The H I rotation curves of massive spirals at z ' 0 reach radii out to 50 100 kpc, so they are much more extended than the CO or H↵ rotation curves of galaxies at cosmic noon, halting at R ' 8 kpc. Unfortunately, the H I line is weak and can-not be detected (nor spatially resolved) in individual galaxies at z ' 1 3, probably not even with the future Square Kilometre Array at full capacity.
We stress that zC-400569 is one of the galaxy with the best IFU data (Genzel et al. 2017), so represents the current state of the art in terms of H↵ kinematics at cosmic noon. In principle, to increase the sensitivity in the outer galaxy regions, one may stack the line emission of multiple galaxies after normalizing both the radial extent and the velocity extent in some way. Such stacking experiments have been performed (Lang et al. 2017;Tiley et al. 2019) but led to controversial results because the final stacked rotation curves depend on the adopted normalization, which is not surprising given the observed variety in the inner rotationcurve shapes of individual galaxies.
Generally speaking, H↵ and CO emissions are cospatial with the star-forming disk, so they cannot easily probe the ubiquitous flattening of rotation curves beyond the bright stellar component of galaxies. This problem is even more prominent for the high-J CO lines (J = 3 2 and higher) that are routinely imaged by ALMA at high z. In fact, high-J CO lines trace only dense and warm H 2 gas that is closely associated with active star-forming sites, unlike the lower-J CO lines that trace also less dense and cooler H 2 phases. A possible alternative is the S (0) : J u J l = 2 0 rotational line of H 2 at 28 µm that emanates from a warm molecular gas phase concomitant with the cold atomic phase in spirals (Papadopoulos et al. 2002), but there is currently no telescope to image this line in galaxies across cosmic time. The most promising line at high z is the bright [C II] emission at 158 µm, which traces a combination of atomic, molecular, and ionized gas, but this line is most easily observed by ALMA only at z > 4 (e.g., Lelli et al. 2021).
Despite all these limitations, Fig. 10 shows a remarkable fact: the inner rotation curves (R < 8 9 kpc) of zC-400569 and zC-488879 are very similar to those of local galaxies with similar stellar surface densities (⇠10 3 M pc 2 ). This result suggests that main-sequence galaxies at cosmic noon are dynamically evolved and may turn into local spiral galaxies without the need for a major redistribution of their total mass, at least within the inner 10 kpc. This fact is in line with ALMA kinematic studies of galaxies at z > 4 (Rizzo et al. 2020(Rizzo et al. , 2021;;Lelli et al. 2021;Tsukui & Iguchi 2021) and recent results at z > 3 (Ferreira et al. 2022) from the James Webb Space Telescope (JWST), pointing to speedy galaxy evolution in the first billions year of the Universe's lifetime.
In Sect. 5.3, we find that the rotation curves of both zC-400569 and zC-488879 can be fitted with a mass model that contains only baryonic components. This occurs even though the rotation curve is flat due to the sum of a declining bulge contribution at small radii and a rising disk contribution in the outer parts (see Fig. 8). Similar mass models are common for local spirals (Lelli et al. 2016a). Historically, indeed, the evidence for the DM e↵ect in galaxies did not come simply from "flat rotation curves" because most H↵ rotation curves can be fitted by mass models with only baryonic components (Kalnajs 1983;Kent 1986). The unambiguous evidence for the DM e↵ect came from rotation curves that remain flat outside the stellar component, beyond the peak of the disk contribution at 2 3R d , such as those from deep H I observations (van Albada et al. 1985;Kent 1987;Begeman 1989; see Sanders 2014 for an historical perspective).
To illustrate the situation, Fig. 11 shows the location of zC-400569 and zC-299569 on the empirical RAR at z ' 0 A106, page 15 of 21 Lelli, F., et al.: A&A 672, A106 (2023) . For both galaxies we assume the values of ⌥ gas , ⌥ disk , and ⌥ bul from the baryons-only mass model (Sect. 5.3). Both galaxies lie on the 1:1 relation (dotted line) at high accelerations, similarly to disk galaxies at z ' 0. To unambiguously probe the DM e↵ect (g obs > g bar ) one needs more extended rotation curves probing lower accelerations. (McGaugh et al. 2016;Lelli et al. 2017): the observed centripetal acceleration from rotation curves (g obs = V 2 rot /R) is plotted against that expected from the distribution of baryons (g bar = r bar ). For local galaxies, the value of g bar is well measured using Spitzer photometry at 3.6 µm, which provides a dust-free tracer of the stellar mass distribution and allows for the adoption of a common mass-to-light ratio for all galaxies (e.g., Lelli et al. 2017). For high-z galaxies, rest-frame NIR photometry is not available yet, so we need to make assumptions on stellar masses to locate them on the RAR. This is similar to the observational situation at z ' 0 about 20 yr ago (McGaugh 2004).
In Fig. 11, we assume the values of ⌥ disk , ⌥ bul , ⌥ disk from the baryon-only mass model (Sect. 5.3). As expected, both high-z galaxies lie on the unity line (g obs = g bar ) at accelerations higher than a few 10 9 m s 2 . Clearly, lower values of ⌥ disk , ⌥ bul , ⌥ disk (such as those found in the NFW model, Sect. 5.4) would shift g bar toward the left and move the data above the unity line. In local galaxies, the DM e↵ect (g obs > g bar ) appears at accelerations lower than ⇠10 10 m s 2 . If we assume that the rotation curves of our high-z galaxies remain flat at the observed values, then they should be traced out to 20 30 kpc to unambiguougly probe the DM-dominated regime. Unfortunately, these distances seem out of reach for H↵ and CO data from existing facilities.
In Sect. 5.4, we find that the rotation curves of both zC-400569 and zC-488879 can be fitted with a NFW halo in a ⇤CDM context. Adding an NFW halo comes at the expense of decreasing the disk mass by a factor of 2 3 with respect to the baryon-only model. As a result, for the NFW model, the DM contribution is already important in the inner galaxy regions. Given the large uncertainties on stellar and gas masses at z ' 1 3, a NFW halo cannot be categorically ruled out.
The current situation is reminiscent of the well-known "diskhalo degeneracy" (van Albada et al. 1985), which has plagued the mass modeling of galaxies for a few decades. Essentially, the halo and disk contributions are degenerate because they have a similar shape (see Fig. 8), so one needs to know the disk mass with high accuracy to determine the DM halo properties and inner DM fractions. One way to break this degeneracy is using rest-frame NIR photometry (e.g., Sanders & Verheijen 1998;Lelli et al. 2016a) because the stellar mass-to-light ratio is almost insensitive to the galaxy star-formation history, metallicity, and dust extinction (e.g., Schombert et al. 2019), so the stellar masses can be estimated with an accuracy of ⇠25%. This approach gives negligible DM fractions in the inner regions of massive galaxies at z ' 0 (see Fig. 11). For galaxies at cosmic noon, rest-frame NIR imaging will become available soon thanks to JWST, so we will be able to break the disk-halo degeneracy at z ' 1 3.
Our conclusions di↵er from those of Genzel et al. ( 2020) and Price et al. (2021), who report high baryonic fractions and cored DM distributions in galaxies at cosmic noon (see also Bouché et al. 2022). Possibly, their findings are driven by an underrating of the disk-halo degeneracy, which may occur when the DM halo concentration is fixed and/or strong priors on stellar and gas masses are imposed. The most conservative conclusion is that it is just impossible to draw firm statements on the amount and distribution of DM in galaxies at high z with the available data. However, the various similarities between our galaxies at z ' 1 2 and massive spirals at z ' 0 suggest that inner baryon dominance and core DM halos may, indeed, be the most likely scenario for galaxies at the cosmic noon. This scenario may soon be tested using rest-frame NIR images from JWST.
We presented high-resolution ALMA observations of multiple CO transitions for two main-sequence galaxies at cosmic noon: zC-400569 at z ' 2.24 and zC-488879 at z ' 1.47. We also reanalyzed H↵ data for zC-400569 from SINFONI at the VLT. Our main results can be summarized as follows: 1. Both galaxies have regularly rotating CO disks at R . 8 kpc and show possible hints of minor mergers and/or mild interactions at larger radii; 2. The rotation curves of both galaxies are flat out to R ' 8 kpc and show no sign of a decline in the outer parts; 3. The CO velocity dispersion seems low ( CO . 15 km s 1 ), smaller than the spectral resolution, so the gas disks have V rot / CO & 17 22 indicating a low turbulence environment; 4. Mass models from HST images reveal a severe disk-halo degeneracy: models with inner baryon dominance and models with a NFW halo can fit the observed rotation curves equally well with acceptable stellar and gas masses; 5. Mass models in Milgromian dynamics (MOND) can fit the observed rotation curves with the same acceleration scale and interpolation function measured at z ' 0. In the near future, rest-frame NIR imaging from JWST will provide accurate stellar masses (comparable to those at z ' 0) and allow us to distinguish between various mass models, breaking the disk-halo degeneracy. In particular, we may be able to discriminate between ⇤CDM models, which predict a significant DM e↵ect within 8 kpc, and MOND models, which do not.
In terms of galaxy evolution, the dynamical properties of these two galaxies (rotation-curve shapes, high V rot / V ratios, A106, page 16 of 21 and bulge-to-total ratios) are remarkably similar to those of spiral galaxies at z ' 0. This fact suggests that some massive galaxies at cosmic noon undergo weak dynamical changes over more than 10 Gyr of the Universe's lifetime, evolving into local spiral galaxies without the need of a severe mass redistribution in their inner regions.
http://abyss.uoregon.edu/~js/archangel/ A106, page
of 21