The ALMA-ALPAKA survey I. High-resolution CO and [CI] kinematics of star-forming galaxies at z = 0.5–3.5

Context. Spatially resolved studies of the kinematics of galaxies provide crucial insights into their assembly and evolution, enabling one to infer the properties of the dark matter halos, derive the impact of feedback on the interstellar medium (ISM), as well as measure and characterize the outﬂow motions. To date, most of the kinematic studies at z = 0 . 5 − 3 . 5 have been obtained using emission lines tracing the warm, ionized gas (e


Introduction
According to the current paradigm of galaxy formation and evolution, the assembly of galaxies is regulated by a variety of physical processes: the interplay between dark and baryonic matter, gas accretion, galaxy mergers, star formation, and stellar and active galactic nucleus (AGN) feedback (Mo et al. 2010;Cimatti et al. 2019;Vogelsberger et al. 2020).From the theoretical perspective, state-of-the-art cosmological simulations repro-duce most of the global properties of galaxies at different cosmic epochs (e.g., Schaye et al. 2015;Nelson et al. 2018;Pillepich et al. 2018;Roca-Fàbrega et al. 2021;Kannan et al. 2022;Pallottini et al. 2022).However, this success does not necessarily imply that we understand how galaxies form and evolve in detail, as the modeling of processes acting on scales below the resolution of simulations (e.g., stellar and AGN feedback, star formation) relies on assumptions and calibrations for the so-called sub-grid models (e.g., Naab & Ostriker 2017;Vogelsberger et al. 2020).
As a consequence, the kpc and sub-kpc spatial distributions of the baryonic and dark matter within galaxies strongly vary with the adopted models (Kim et al. 2016;Roca-Fàbrega et al. 2021).
From the observational perspective, the role played by various processes in driving galaxy evolution is still unclear, even in redshift ranges that have been widely studied in the last decades (e.g., z = 0.5 − 3).Understanding the assembly of galaxies requires simultaneous high-resolution, multi-wavelength studies of their morphologies and kinematics.The morphological analysis of the distribution of stars (e.g., Lang et al. 2014;van der Wel et al. 2014;Mowla et al. 2019;Ferreira et al. 2022;Kartaltepe et al. 2023) andgas/dust (e.g., Calistro Rivera et al. 2018;Gullberg et al. 2018;Hodge et al. 2019;Rujopakarn et al. 2019;Kaasinen et al. 2020;Tadaki et al. 2020;Puglisi et al. 2021;Ikeda et al. 2022) within galaxies allows one to trace the build up of structures.Galaxy kinematics, on the other hand, provide constraints on the impact of gas accretion, mergers, outflows, and environmental mechanisms on the growth of galaxies (e.g., Übler et al. 2018;Loiacono et al. 2019;Concas et al. 2022;Bacchini et al. 2023;Roman-Oliveira et al. 2023).For instance, the prevalence of rotating disks among the star-forming population at z ∼ 1 − 3 has been established as a key evidence that their assembly is mainly driven by smooth gas accretion as opposed to mergers (e.g., Wisnioski et al. 2015;Stott et al. 2016;Förster Schreiber & Wuyts 2020).In addition, for galaxy disks, measurements of gas rotation allow one to study the dynamics and infer the content and distribution of dark matter within galaxies (e.g., Straatman et al. 2017;Posti et al. 2019;Price et al. 2021;Sharma et al. 2022;Mancera Piña et al. 2022;Lelli et al. 2023).In contrast, measurements of gas velocity dispersion provide key insights into the mechanisms injecting turbulence into the interstellar medium (ISM; e.g., stellar feedback, release of gravitational energy through accretion from the cosmic web, radial flows within galaxies; Krumholz et al. 2018;Übler et al. 2019;Kohandel et al. 2020;Rizzo et al. 2021;Jiménez et al. 2023;Rathjen et al. 2023).
Three reasons could be responsible for the discrepancies in the two redshift ranges (i.e., z = 0.5 − 3.5 and z = 4 − 5): spectral resolution, beam-smearing corrections and the use of different tracers of the ISM phases.On the one hand, IFU observations are characterized by spectral resolutions that are typically a factor of 3 worse than those achievable with ALMA.Furthermore, discrepancies between different studies may arise due to the techniques employed for deriving beam-smearing corrected kinematic parameters.The beam-smearing effect causes, in fact, a degeneracy between the rotation velocity and velocity dispersion, mainly resulting in inflated values of σ (e.g., Begeman 1989).In the past decade, different techniques have been developed to correct for beam-smearing effects.They rely on either applying forward-modeling approaches to fit the data cubes (e.g., Bouché et al. 2015;Di Teodoro & Fraternali 2015) or correcting the velocity dispersion maps a-posteriori (e.g., Swinbank et al. 2012;Stott et al. 2016;Johnson et al. 2018).Using ALMA mock data at different angular resolutions, Rizzo et al. (2022) showed that the a-posteriori correction is suboptimal, resulting in overestimation of σ up to 200% when the galaxy is spatially resolved with only a few resolution elements.Instead, the typical errors of the velocity dispersion using forward-modeling techniques are of the order of 25% (Di Teodoro & Fraternali 2015;Rizzo et al. 2022).On the other hand, the distinct gas phases probed by different emission lines may also play a critical role.Kinematic studies at z ≈ 0 show that the velocity dispersion is typically higher when measured with warm ionized gas as opposed to molecular or atomic gas (Levy et al. 2018;Girard et al. 2021).At high-z, simultaneous studies of the cold and warm kinematics on the same galaxies are still in their infancy and limited in statistics (e.g., Übler et al. 2018;Molina et al. 2019;Lelli et al. 2023).Furthermore, the warm gas might be more affected by the presence of gas in outflows than the cold gas (Lelli et al. 2018;Ejdetjärn et al. 2022;Kretschmer et al. 2022).Should this be the case, the rotation velocity and velocity dispersion from warm gas might be uncertain as these two quantities are typically derived assuming that the gas moves in circular orbits, while radial or vertical motions could bias their measurements.
Systematic, high-resolution kinematic measurements obtained using cold gas tracers are needed not only to gain a comprehensive insight into the ISM properties of z ∼ 0.5 − 3.5 galaxies but also to make an accurate comparison with the dynamical properties of galaxies at z ≳ 4. Unfortunately, because of the sensitivity and frequency coverage of state-of-the-art facilities, high-resolution [CII] observations at z ≲ 3.5 are challenging or not feasible, due to the low atmospheric transmission (Carilli & Walter 2013).To study the z ∼ 0.5 − 3.5 kinematics using cold gas tracers, one can target carbon monoxide, CO, transitions or fine structure lines from atomic carbon, [CI].However, emission lines from CO and [CI] are typically fainter than [CII] (Carilli & Walter 2013;Bernal & Kovetz 2022).To date, the only instrument capable of achieving the combined sensitivity and angular resolution requirements of ≲ 0.3 ′′ needed for studying the CO or [CI] kinematics at z ≳ 0.5 is ALMA.However, even with ALMA, CO and [CI] observations are extremely time consuming.For this reason, studies of the kinematics using cold gas tracers at cosmic noon have been presented only for a handful of galaxies (e.g., Molina et al. 2019;Noble et al. 2019;Kaasinen et al. 2020;Xiao et al. 2022;Lelli et al. 2023).Finally, since the published kinematic measurements are derived using different algorithms and assumptions, a systematic compilation and comparison of the results is not straightforward.
Here we present the project "ALMA Archival Large Program to Advance Kinematic Analysis" (ALMA-ALPAKA), aimed at filling all the above gaps by collecting high data quality emission line observations of z = 0.5 − 3.5 galaxies from the ALMA public archive.We note that the use of archival data does not have an impact on the originality of the present work as, to date, the ALMA observations for ≈ 67% of the ALPAKA galaxies have not been published in any studies.With an increase in sample size by a factor of 3 compared to the total number of targets found in all previous works (e.g., Molina et al. 2019;Noble et al. 2019;Kaasinen et al. 2020;Xiao et al. 2022;Lelli et al. 2018Lelli et al. , 2023)), ALPAKA will be used to obtain the first high-resolution characterization of the morpho-kinematic properties of galaxies at z ∼ 0.5 − 3.5 using millimeter to optical observations.
In this paper, the first of a series, we present the ALPAKA sample, discuss the sample selection (Sect.2) and describe the ALMA and Hubble Space Telescope (HST) data analyzed in this work (Sect.3).In Sect.4, we derive the global properties of the ALPAKA targets using ancillary data.The kinematic modelling and assumptions, the identification of the disks and the description of their kinematic parameters is presented in Sect. 5.In Sect.6, we describe the kinematics of each ALPAKA target.In Sect.7, we discuss the potential bias due to selection effects and the dynamical properties of the subsample of ALPAKA disk galaxies.Finally, in Sect.8, we summarize the main findings and describe the main objectives of the ALPAKA project.
Throughout this paper, we assume a ΛCDM cosmology with parameters from Planck Collaboration et al. ( 2020) and a Chabrier Initial Mass Function (Chabrier 2003).

Selection criteria
The ALPAKA project is designed to study the kinematic and dynamical properties of galaxies at z = 0.5 − 3.5.To this end, we collected the sample by selecting data publicly available from the ALMA archive and with spectral setup covering CO and/or [CI] emission lines.We queried the database at the end of August 2022 to select galaxies with spectroscopic redshift in the range 0.5 − 4, with angular resolution ≲ 0.5 ′′ , spectral resolution ≲ 50 km s −1 and medium to high signal-to-noise ratio (SNR), i.e., data cubes with a SNR ≳ 3.5 per channel in at least 5 spectral channels.These requirements guarantee data quality sufficient to infer robust kinematic parameters (Rizzo et al. 2022).The resulting sample consists of 28 galaxies, whose ID, name, coordinates, and redshifts are reported in Table 1 (see Sect. 6 for details on each target).In Fig. 1 (upper panel), we show the histogram with the redshift distribution of the ALPAKA sample.We note that about half of the sample (12/28) covers the lowest redshift range (z ≲ 1.6) with two targets at z ≈ 0.6, while the remaining 16 galaxies are at z ≳ 2. The peak of 8 galaxies in the redshift bin centred at z = 1.5 is due to targets belonging to the same cluster (see Sect. 2.2).

Target characterization
High-resolution ALMA observations of high-z galaxies require long integration times.Therefore, to optimize the observing strategy, such observations are usually performed on sources for which the fluxes of the target emission lines are known through previous low angular resolution and shallow campaigns.Despite the fact that our targets were selected using an agnostic approach regarding the multi-wavelength coverage, we find that all AL-PAKA galaxies lie in broadly studied survey areas.In Table 1, we list the main features of the targets as reported in the literature: the fields in which they lie, and information regarding the environment and the presence of AGN (see Sect. 6 for details).The latter are based on previous diagnostic features (e.g., rest- frame UV or mid-infrared spectra, X-ray luminosities) depending on spectroscopic coverage of the galaxies in different fields.
Given that the selection criteria we adopt are only based on the data quality of archive data, the ALPAKA sample is heterogeneous: 68% (19 out of 28) of the targets lie in overdense regions (clusters, groups and protocluster); 7% (2 out of 28) were previously considered as merging systems and 1 of the two is in a group; 25% (7 out of 28) host an AGN.Among the AGN hosts, 3 sources are protocluster members.A discussion on the potential impact of AGN feedback on the ISM kinematics properties of the host galaxies is provided in Sect.7.

ALMA data
In Table 2, we list the main properties of the ALMA datasets: ALMA project ID, frequency coverage and emission line for each target.When multiple observations at similar angular resolutions are available (e.g., ID10 and 14), we combine them to increase the SNR.Due to the selection criteria, and to the ALMA sensitivity and frequency coverage, the ALPAKA galaxies have observations of various emission lines: 18 targets have low-J CO transition observations -CO(2-1) or CO(3-2) -, while for 10 targets only high-J CO transitions or [CI] emission lines are available (see bottom panel in Fig. 1 and Table 2).The on-source integration times for each target are listed in Table 2 and range  from 43 minutes to 14 hours, for a total of ≈ 90 hours.The total integration time, including overheads, is of 147 hours, corresponding to the duration of an ALMA Large Program.
In this paper, we make use of the calibrated measurement sets provided by the European ALMA Regional Centre (Hatziminaoglou et al. 2015), that calibrated the raw visibility data using the standard pipeline script delivered with the raw observation sets.All of the post-processing steps were handled using the Common Astronomy Software Applications (Casa) suite (Mc-Mullin et al. 2007), version 6.The calibrated data were first inspected to confirm the quality of the pipeline calibration and that no further flagging was required.For data sets containing one single target, we then subtracted the continuum from the line spectral windows using UVCONTSUB.Most of the data are averaged into groups of between two to six channels in order to obtain an average SNR of at least 3.5 per velocity channel in at least 5 spectral channels of the data cube1 .This procedure results in channels with typical velocity widths ranging from 16 to ≈ 39 km s −1 .The continuum and spectral lines were imaged using the TCLEAN routine in the Casa package, assuming a natural weighting of the visibilities to maximise the SNR.Targets with high SNR are imaged using a Briggs weighting of the visibilities (robust parameter set equal to 0.5; Briggs 1995), in order to enhance the angular resolution of the output images, without significantly degrading their effective sensitivity and meet the requirements on the SNR per channel defined above.The CLEAN algorithm is run down to a flux threshold of 2 × RMS, where RMS is the root mean square of the data measured within the dirty data cubes.
For datasets containing multiple targets, to better account for the source-to-source variation of the continuum signal, we perform the continuum subtraction using the IMCONTSUB task.In Table 2, we present the main properties of the ALPAKA data cubes, their beam size, channel widths, and RMS per spectral channel.The beam of the observations of the ALPAKA targets ranges from 0.1 ′′ to 0.5 ′′ (median value of 0.25 ′′ ) and the corresponding resolution in physical units varies between 1 kpc and 4 kpc (median value of 2 kpc).

HST data
Spatially resolved observations of the stellar continuum of galaxies allow one to derive not only the structural parameters defining the light distribution (e.g., effective radius, Sérsic index) but also to constrain the geometrical parameters (e.g., inclination angle) of the sources.As discussed in detail in Sect.5.1.1,accurate measurements of the inclination and position angles provide robust velocity estimates and kinematic characterizations.The comparison between the stellar and gas morphology and the gas kinematics of a galaxy helps to identify any merger or outflow features (see Sect. 5).Further, combined measurements of the stellar and gas distributions and the rotation curves provide a unique means to infer the dark matter content within galaxies.
For 23 ALPAKA galaxies, HST observations are publicly available and are taken from the Complete Hubble Archive for Galaxy Evolution (CHArGE, Brammer et al., in prep.).The latter performs uniform processing of all archival HST imaging and slitless spectroscopy observations of high-z galaxies.
In CHArGE, the data were processed with the Grizli pipeline, which creates mosaics for all filter exposures that cover a given area of the sky (Brammer & Matharu 2021).All exposures are aligned to each other using different techniques (see Kokorev et al. 2022, for details) resulting in a typical astrometric precision ≲ 100 mas.
For each source we select the reddest HST filter (see Table 3) in which the source is detected in order to minimize biases in the determination of the structural parameters due either to the irregular morphology of the galaxies that host UV-bright star-forming regions (Guo et al. 2018;Zanella et al. 2019) or the presence of dust attenuation and its patchy distribution (Cibinel et al. 2017).
In Table 3 we report the central rest-frame wavelength covered by the selected HST filter.For 17 galaxies the HST data cover the rest-frame optical emission (≳ 4000Å), while 6 galaxies are covered only in the near-UV range.Fig. 2 shows the HST images for the 23 ALPAKA targets, while the white contours show the integrated CO or [CI] emission lines from ALMA data.By visually inspecting the overlap of the HST and ALMA data, it is evident that the bulk of the stellar and gas emissions largely overlap.However, for some galaxies (e.g., ID1, 2, 12, 15) the gas emission seems more compact than the stellar counterpart (see Sects. 6 and 7.1).A detailed comparison between the sizes and morphologies of the stellar and gas emission is contained in the MSc thesis by D. Frickmann 2 and it will be extended and published as part of the ALPAKA series (Sect.8).

Global properties
Being in well characterized survey areas, extensive studies of the global, unresolved properties (e.g., stellar mass, M ⋆ , and starformation rate, SFR) of the ALPAKA galaxies are already available.However, these parameters are derived using different algorithms and assumptions.For this reason, we collect UV-to-radio photometric data using public multi-wavelength catalogues (Fu et al. 2013;Ivison et al. 2013;Magnelli et al. 2013;Jin et al. 2018;Hayashi et al. 2018;Liu et al. 2019;Weaver et al. 2022) with the aim of fitting the spectral-energy distribution (SED) in a consistent way.The only exceptions are ID11 and ID12, for which photometric data are not publicly available.Throughout the rest of the paper, for these two galaxies, we will refer to the stellar masses and SFR reported in Noble et al. (2017) and obtained after making assumptions consistent with those used for 2 https://nbi.ku.dk/english/theses/masters-theses/ditlev-frickmann/ the rest of our sample (see Kokorev et al. 2021, for a detailed discussion on the M ⋆ and SFR obtained with different tools).
To derive the stellar masses and SFR for the ALPAKA galaxies, we perform the SED fitting using Stardust (Kokorev et al. 2021).Stardust performs a multi-component fit that linearly combines three classes of templates: stellar libraries from an updated version of the Stellar Population Synthesis models described in Brammer et al. (2008); AGN torus templates from Mullaney et al. (2011); infrared models of dust emission arising from star formation from Draine et al. (2007Draine et al. ( , 2014)).These three components are fitted simultaneously but independently from each other, i.e., without assuming an energy balance between the absorbed UV/optical radiation and the infrared emission.This approach allows one to account for spatial offsets between the stellar and dust distributions within a galaxy (see discussion in Kokorev et al. 2021).For fitting the ALPAKA targets, we include the AGN templates only for galaxies that are previously identified as AGN hosts (Table 3).
In Table 3, we show the best-fit stellar masses, SFR, and the star-formation infrared luminosity (L IR ) from Stardust.We note that, for AGN-host galaxies, in addition to star-formation, the dusty tori contribute to the total infrared luminosity.However, the values of L IR listed in Table 3 are obtained after integrating the best-fit star-formation dust model employed by Stardust in the range 8 -1000 µm.For all ALPAKA galaxies, we consider only the dust-obscured SFR, traced by L IR , and we neglect the contribution from UV emission.For three galaxies -ID16, 17, 24 -due to the lack of good coverage in the optical/nearinfrared, the uncertainties on the stellar masses are much larger than the parameter values itself.Therefore, the respective estimates of the stellar mass are not stastically meaningful.In addition, due to the lack of a good far-infrared coverage, measuring the SFR for ID24 is challenging.In Fig. 3, we show the distribution of the ALPAKA sample in the SFR-M ⋆ plane for the 25 galaxies with a reliable estimate (i.e., uncertainties smaller than the value) of M ⋆ .We divide our targets in three redshift bins and compare them with the main-sequence relation at the corresponding redshift.For the latter, we use the parametrization of Schreiber et al. (2015, solid lines) using a Chabrier IMF and show the 0.3 dex scatter at the average redshifts of the galaxies in the three bins, z = 1.3, 2.2, 3 (dashed lines).The dot-dashed lines in Fig. 3 show the empirical lines, located 4 times above the SFR of main sequence that is usually used to divide mainsequence from starburst galaxies (Rodighiero et al. 2011).According to this definition, a large fraction of ALPAKA galaxies are starbursts (12/25 or 48%; see column 5 in Table 3).In the low redshift bin (z =0.5 -1.5), 70% of ALPAKA galaxies lie within the ±1σ scatter of the main sequence relation, while this fraction falls to 37% at z ∼ 3.
The ALPAKA sample covers high stellar mass galaxies, ≳ 10 10 M ⊙ with SFR ranging from 8 to ∼3000 M ⊙ yr −1 .This can be ascribed to a selection effect: being discovered as bright sources in the infrared or sub-mm wavelength, ALPAKA galaxies have high SFRs and gas fractions.To visualize this, in Fig. 4, we show the distribution of the ALPAKA sample in the emission line-infrared luminosity planes, an observational proxy of the Kennicutt-Schmidt relation (Schmidt 1959;Kennicutt 1998), with respect to compilations of local and highz galaxies from the literature (e.g., Liu et al. 2015;Silverman et al. 2018;Valentino et al. 2020;Boogaard et al. 2020;Birkin et al. 2021;Valentino et al. 2021).The CO fluxes and respective luminosities for the ALPAKA sample, as listed in Table 3, are derived by summing the flux above 3 × RMS in the highresolution flux-integrated maps presented in Sect.3.1.We check Notes.For targets ID4 -9, the observation consists of three mosaic pointings; the corresponding on-source integration time refers to the total value summed over all pointings.For targets with multiple observations (i.e., ID10 and 14), we image the CO after combining the corresponding measurement sets and we list the resulting properties of the cubes (beam, channel width, RMS).
Here, we provide the list of principal investigators (PIs) for the ALMA projects employed in this work: 2017.that these values are within ±20% from the ones obtained by fitting the total-flux maps with the IMFIT tool within Casa that adopts a 2D Gaussian.In all cases, these fluxes are consistent with previous estimates, mostly obtained through unresolved observations, showing that we are not missing any significant emission on large scales.The samples of galaxies from the literature comprise: the compilation in Valentino et al. (2020) consisting of 30 main-sequence galaxies at z ∼ 1 and 65 submillimeter galaxies (SMGs) and quasars at z ∼ 2.5 (Walter et al. 2011;Alaghband-Zadeh et al. 2013;Cañameras et al. 2018;Yang et al. 2017;Andreani et al. 2018) and 146 local starbursts from Liu et al. (2015); 12 starbursts at z ∼ 1.6 (Silverman et al. 2018); 22 main-sequence and starburst galaxies at z ∼ 1.4 from the ALMA Spectroscopic Survey in the Hubble Ultra Deep Field survey (Boogaard et al. 2020); 47 SMGs at a median z of 2.5 from Birkin et al. (2021).Considering the substantial effort in obtaining high SNR spatially resolved ALMA observations at z ∼ 0.5 − 3.5, the ALPAKA sample typically covers the brightest part of the distributions in each panel of Fig. 4, with a median luminosity of 2 × 10 10 K km s −1 pc 2 .ALPAKA galaxies are, therefore, somewhat biased towards the most actively star-forming galaxy populations (see discussion in Sect.7).

Kinematic analysis
In this section, we describe how we analyze the kinematics of the ALPAKA galaxies by fitting the data using rotating disk models (Sec.5.1).This allows us to provide a first-order description of the gas motions within the ALPAKA sample.In Sect.5.2, we show how we identify any deviations from pure circular orbits, likely due to radial and vertical motions driven by outflows and interactions.The kinematic properties of each ALPAKA target are described in Sect.6.In Sect.7.1, we present the kinematic parameters of the disk subsample.2.5 < z ≤ 3.6 Fig. 3. Distribution of the ALPAKA galaxies in the stellar mass (M ⋆ ) -SFR plane, divided in three redshift bins.The dark green squares show ALPAKA galaxies in the field, while the green circles show galaxies in overdense regions (e.g., clusters, protoclusters, groups).The solid line in the three panels show the empirical main-sequence relations from Schreiber et al. (2015) at z = 1.3, 2.2, 3, which are the average redshifts of ALPAKA targets in the three bins.The dashed and dot-dashed lines show the ±1σ scatter and the line dividing the main-sequence and starburst galaxies, respectively (Rodighiero et al. 2011).We note that only the 25 ALPAKA targets with good estimates of both M ⋆ and SFR are shown.

Disk modeling
We fitted the kinematics of the ALPAKA galaxies using the software 3D Barolo (Di Teodoro & Fraternali 2015).3D Barolo produces three dimensional (two spatial, one spectral axis) realizations of a so-called tilted-ring model (Rogstad et al. 1974).The latter consists of a disk divided into a series of concentric rings, each with its kinematic (i.e., systemic velocity V sys , rotation velocity V rot and velocity dispersion σ) and geometric properties (i.e., center, inclination angle i and position angle PA3 ).For a thin disk model, the line-of-sight velocity V los at a radius R is given by where ϕ is the azimuthal angle in the disk plane.
The best-fit model is obtained by means of a least-square minimization.At each step of the model optimization and before calculating the residuals between the data and the model, 3D Barolo convolves the model disk with a Gaussian kernel with sizes and position angle equal to the beam of the corresponding observation.In the case of the ALMA observations, this is set to be equal to the synthesized beam of the cleaned image.This approach allows for a robust recovery of the rotation velocity and velocity dispersion profiles, since it largely mitigates the effects of beam smearing also in the case of data with relatively low angular resolution (Di Teodoro & Fraternali 2015;Di Teodoro et al. 2016;Rizzo et al. 2022).3D Barolo is a suitable tool to model the gas kinematics at the typical spatial resolution and SNR of the galaxies in ALPAKA.This tool has been extensively tested using mock and real data over a wide range of data quality (e.  ) Fig. 4. Distribution of the ALPAKA galaxies and samples of local and high-z galaxies from the literature (e.g., Liu et al. 2015;Silverman et al. 2018;Valentino et al. 2020;Boogaard et al. 2020;Birkin et al. 2021) in the emission line-infrared luminosity planes.The ALPAKA galaxies typically cover the bright part of the distributions, especially for high-J transitions, as spatially resolved ALMA observations are feasible only for very luminous galaxies.

Geometrical parameters from morphological fitting
3D Barolo can estimate the geometrical parameters, namely the center, inclination, and position angle of each tilted ring component.However, due to the relatively small number of resolution elements covering the CO/[CI] line emission, we prefer to reduce the number of free parameters by fixing the center and inclination.In particular, estimating the inclination of the galaxies is crucial.Correcting for it can in fact account for a large fraction of the rotation velocity if the galaxies are seen at low inclinations, due to the sin i dependence of Eq. ( 1).When dealing with low-resolution observations of low-z and high-z galaxies, the inclination is usually fixed to the one estimated from the optical images (de Blok et al. 1996;Lelli et al. 2016;Wisnioski et al. 2019;Kaasinen et al. 2020).However, since only a fraction of ALPAKA galaxies have HST data covering the rest-frame optical emission, we use two methods for estimating their inclinations: -Galfit on HST data.For the 23 galaxies with HST data, we used Galfit (Peng et al. 2002) to model their 2D surface brightness using one Sérsic component (Sersic 1968, see details in Appendix A).Galfit fits the center, the three parameters describing the Sersic profile (total magnitude, Sérsic index, effective radius), the position angle PA HST and the axis ratio b/a between the projected major and minor axis.The latter allows for computing the inclination angles, i HST = arccos(b/a).In Carlo algorithm that models the geometry of galaxies without assuming parametric descriptions of the surface brightness distribution.Cannubi uses 3D Barolo to fit either the total-flux map or the entire cube using resolution-matched 3D tilted-ring models of rotating disks (see details in Roman-Oliveira et al. 2023;Mancera Piña et al. 2020).The free parameters of the fit are the center of the disk, its radial extent, the position and inclination angles (PA ALMA , i ALMA ).
To perform a brief check of the independence of our results on the specific method used to model each data set, we repeat the analysis of the ALMA data using Galfit.Since Galfit is optimized for dealing with optical/near-infrared images (e.g., units in counts, magnitude zero-points, Poissonian error in each pixel), we applied some arbitrary conversion factors to fit the total-flux maps obtained from ALMA data.The resulting best-fit inclinations are consistent within 5% with the ones obtained with Cannubi.Fig. 5, left panel, shows the distribution of the difference between the inclination angles found with the two methods (upper left panel).For 21 out of the 23 ALPAKA galaxies with HST data, i HST and i ALMA values are consistent within the ±1.5σ and the difference between these two angles are within 20 deg.Since the rotation velocity has a dependence on the inclination angle that goes as 1/ sin i, Eq.( 1), an uncertainty of 20 deg results in relative uncertainties on the rotation velocity of 7% and 18% for a nearly edge-on (e.g., 75 deg) and face-on (e.g., 26 deg) galaxy, respectively4 .The only galaxies with a difference between the two inclination angles ≳ 20 deg are ID3 and ID22, with i HST − i ALMA of 54 and 42 deg, respectively.However, the HST data for both ID3 and ID22 cover only their rest-frame near UV/blue optical emission and they are likely not representative of the bulk of the stellar population.In both cases, there is a difference between the CO and HST morphologies.For instance, ID3 has, two UV bright clumps, clearly visible in Fig. A.1, while the CO emission has a smooth distribution and its center is located between them (Fig. 2).The CO(5-4) emission from ID22 cover only the innermost regions of the corresponding HST data.To summarize, for most of the galaxies with HST data, the value of i HST are consistent with i ALMA , despite that these two values are derived by fitting two different components of the galaxy (i.e., stellar continuum and gas emission line) and using different tools and assumptions.
For the rest of this subsection, we attribute to each ALPAKA galaxy the value of i HST if HST images are available (with the exception of ID3 and 22, see Sect.6), and i ALMA otherwise.In Appendix B, we discuss potential bias towards low inclinations for the ALPAKA sample.

Assumptions for the kinematic fitting
In this section, we describe the assumptions we made to run 3D Barolo and fit the kinematics of the ALPAKA sample: -Mask.Before fitting the data, 3D Barolo uses the source finder derived from Duchamp (Whiting 2012) to build a mask around the regions that are identified as containing the emission from the target.Within 3D Barolo, different parameters can be used to define how to build the mask.For our sample, we selected the parameters SNRCUT and GROWTHCUT that define the primary and secondary SNR cuts applied to the data.Once the emission pixels with flux above a threshold defined by SNRCUT are identified, the algorithm increases the detection area by adding nearby pixels that are above some secondary threshold defined by GROWTHCUT and not already part of the detected object.For ALPAKA galaxies, we used SNRCUT values of 2.5 -4 and GROWTHCUT of 2 -3 depending on the quality of each data cube.We note that the best-fit parameters are robust against the shape of the mask.Performing the fitting with masks obtained with the SMOOTH&SEARCH task results in values of V and σ consistent within the uncertainties with the values reported in this paper.-Radial separation.To keep the number of modeled rings as close as possible to the number of resolution elements, we used a radial separation between rings close to 0.5 -1 times the angular resolution (see Sect.7.1).This assumption ensures that the rotation curves and the velocity dispersion profiles are sampled with almost independent points, with the number of fitted rings ranging between 2 and 5 per galaxy.-Center.We fixed the galactic centers to the values obtained with Cannubi.For 18 out of the 23 galaxies with HST data, the center from Cannubi are consistent within ≲ 0.1 ′′ with the ones found by Galfit on the HST data.The only exceptions are ID3, 14, 17, and 22 (see Sects.5.1.1 and 6 for further details on these targets).In some cases, we changed the coordinates of the center by less than 2 pixels5 after visu-

±5 D
Notes.Columns 2 -5: values of the position and inclination angles, as derived by Galfit using the HST data (when available) and Cannubi using the ALMA total-flux map.Column 6: kinematic position angles derived by 3D Barolo.Column 7: kinematic class (KC) of each ALPAKA source (D: disk, U: uncertain, M: merger).Galaxies for which the HST filter covers the rest-frame wavelength at λ rest,eff ≲ 4000Å (Table 3) are marked with an asterisk * .In these cases, the geometric parameters (PA HST and i HST ) may be not representative of the bulk of stellar emission.
ally inspecting the position-velocity diagrams (PVDs).The latter are cuts of the cubes along the major and minor axis (see Sect.5.2 for details).The coordinates of the final adopted centers are listed in Table 1.-Inclination angle.As discussed in Sect.5.1.1,for 21 galaxies with HST data, we fixed the inclination angles to the ones found by Galfit fitting, i HST , while for the remaining 7 galaxies, we fixed them to the values i ALMA found by Cannubi.
Using the assumptions described above, we run 3D Barolo and we fit the rotation velocity V rot , the velocity dispersion σ and the position angles.For each galaxy, the number of free parameters is equal to N × V rot + N × σ + PA kin = × N + 1, where N is the number of rings over which the galaxy disk is divided.With the angular resolution of the ALPAKA galaxies, N range from 2 to 5. We note that, despite some galaxies show evidence of a warp, we assume a constant PA kin across the galaxy disks as the quality (SNR and angular resolution) of the data is not sufficient for constraining the change of PA as a function of radius.

Outputs
For each galaxy, we show in Fig. 6 and Figs. 8 -10, the total-flux and velocity field maps and the PVDs along the major and minor axis for the data (dashed black contours) and the model (red solid contours).In Appendix C, we show 7 representative channel maps of the data, model, and residuals for each ALPAKA galaxy.The 2D kinematic maps (i.e., velocity field and velocity dispersion field) of the models are not shown here because, as extensively discussed in Rizzo et al. (2022), these are not as informative as the PVDs and channel maps.The velocity field and velocity dispersion fields are, in fact, strongly affected by the angular resolution and SNR of the data.
The best-fit PA kin as fit by 3D Barolo are listed in Table 4.In Fig. 5, we compare these values with PA HST and PA ALMA .As shown in the scatter plots in the central panel, the difference between the morphological PA and the kinematic ones are consistent, within 1.5-σ, with ± 30 deg for 20 of the 23 galaxies with HST data, indicating a general remarkable regularity of the sources.The only exception are ID13, 14, and 28, with PA HST − PA kin in the range 40 -53 deg.Similarly, the absolute difference between PA ALMA and PA kin is consistent, within 1-σ, with ±30 deg for 27 out of the 28 ALPAKA galaxies, while ID27 has a difference of ∼ 45 deg.In Sect.5.2, we discuss potential explanations of the discrepancies between PA kin , PA ALMA and PA HST for this subsample of ALPAKA galaxies.A detailed description of the best-fit rotation velocity and velocity dispersion values is, instead, given in Sect.7.1.

Kinematic classification
The velocity fields of the ALPAKA galaxies are characterized by a smooth gradient.Therefore, we could, in principle, conclude that all ALPAKA galaxies are smooth rotating disks.However, the velocity map of a merging system can be similar to that of a smooth rotating disk as, due to the angular resolution of observations, the irregularity and asymmetries are smoothed out (e.g., Simons et al. 2019;Kohandel et al. 2020;Rizzo et al. 2022).Further, even outflows could result in velocity gradients that could be erroneously interpreted as rotation (Loiacono et al. 2019).In this Section, we will discuss how we identify any non-circular motions in the ALPAKA targets and build a subsample of galaxies where the presence of a rotating disk can be considered robust.

Visual inspection
For rotating disks with non solid-body rotation curves, the PVDs have specific features: the major-axis PVD has an S-shape profile, and the minor-axis PVD has a diamond shape, symmetric with respect to the axes defining the center and the systemic velocity.Using mock ALMA data of simulated galaxies, Rizzo et al. (2022) show that these features are imprinted in the data even for barely resolved observations.At the typical resolutions of high-z observations and for galaxies with flat rotation curves, the flux distribution along the major-axis PVD has two symmetric brightest emission regions in the approaching and receding sides, along the horizontal parts of the S-shape.In addition, using geometrical arguments, it can be shown that for an ax-isymmetric rotating disk, the kinematic position angle should be aligned with its projected morphological major axis.Differences between these two angles of ≳ 30 deg can be ascribed to a variety of reasons -e.g., presence of outflows (e.g., Lelli et al. 2018;Hogarth et al. 2021) or non-axisymmetric structures (e.g., bar or interaction features; Krajnović et al. 2011).By visually inspecting the PVDs of the data and models, the channel maps and comparing the morphological and kinematic position angles, we identify three classes of galaxies: -Disk These are galaxies with PVDs and channel maps typical of rotating disks, and alignment between the kinematic and morphological position angles, with PA kin − PA ALMA ≲ 30 deg.Some ALPAKA galaxies host disks rotating with remarkable regularity (e.g., ID1, 13, 18, 23, 24, see Sect.6 for details), while others have a few features indicating the presence of kinematic anomalies mostly identified at low SNR: asymmetries along the minor axis (e.g., ID3, 6, 7, 11), excess emission at high velocities in the inner regions (e.g., ID 28, see Sect.6.15 for details).The former can be ascribed to disturbances driven by environmental effects (e.g., ram pressure stripping) or gravitational perturbations; the latter are likely due to emission from outflows.Overall, 19 out of the 28 ALPAKA galaxies are disks.-Merger These are systems where either there are two interacting galaxies that can be clearly identified in the PVDs, channel maps, and integrated-flux maps (ID17) or galaxies for which a rotating disk model does not reproduce the emission in the PVDs and channel maps, and the velocity fields are strongly irregular (ID27).The latter case likely indicates the presence of a late-stage interaction.-Uncertain We choose to be conservative in this classification.Therefore, if a galaxy does not fit into the above classes for any reason, it is classified as uncertain.The 7 ALPAKA galaxies in this class have distorted iso-velocity contours and asymmetries in the PVDs.In most cases, the angular resolution of the observations is not sufficient for definitely identifying their nature (e.g., ID5,16) and discriminating between irregularities due to the presence of two interacting/merging galaxies or kinematic anomalies driven by non-circular motions within a disk structure.
The kinematic class for each ALPAKA galaxy is listed in Table 4 and shown in Fig. 6 and Figs. 8 -10 with a letter: D for disk, M for merger or interacting system, U for uncertain.Details regarding the motivation of the kinematic class assigned to each target are provided in Sect.6.

PVsplit analysis
Most of the kinematic classification methods rely on the analysis of the moment maps.However, recent studies show that at the typical resolution and SNR of current observations, the success rate of these methods can be as low as 10% (Rizzo et al. 2022).On the other hand, Rizzo et al. (2022) present a new classification method -"PVsplit"-that relies on the morphological and symmetric properties of the major-axis PVDs, quantitatively defined by three parameters: P major , P V , and P R .The parameter P major quantifies the symmetry of the PVD with respect to the systemic velocity: a rotating disk with an S-shape profile should have a completely symmetric PVD and a values of P major = 0.The parameters P V and P R define the position of the centroid of the brightest regions in the major-axis PVDs with respect to the line-of-sight velocity and center position, respectively.The PVsplit method has been tested using both low-z systems and ALMA mock data of simulated galaxies (Kohandel et al. 2020;Pallottini et al. 2022), known to be disks, disturbed disks, and interacting systems.Rizzo et al. (2022) find that disks and mergers occupy different locations in the 3D space defined by the PVsplit parameters (gray circles and squares in Fig. 7) and Roman-Oliveira et al. ( 2023) define a plane that divides these two kinematic classes (red plane in Fig. 7).Simulated disks with disturbances driven by interactions are located both in the disk and merger PVsplit sections.
Despite the high quality of the ALPAKA data allowing for the accurate derivation of kinematic properties and identification of merger features, we applied the PVsplit method to get a further quantitative confirmation that the three classes described in Sect.5.2.1 are reliable.In Fig. 7, we show the location of the galaxies that we classified as disks (green circles), mergers (yellow squares), and uncertain (blue diamond) according to the visual inspection of the spectral channels, PVDs, and residuals.Overall, ALPAKA disks and non disks (mergers and uncertain) lie in the two distinct regions of the PVsplit diagram.Some AL-PAKA disks with kinematic anomalies (e.g., ID 3, 28) lie in the merger section or close to the dividing plane, similarly to the simulated perturbed disks analyzed in Rizzo et al. (2022).

ALPAKA in detail
In this Section, we summarize the main physical properties of each ALPAKA target based on previous results from the literature, along with a description of the kinematic fitting and properties.The 2D maps (total flux and velocity fields), the PVDs of the data and models are shown in Figs. 6, 8-13.

ID1
ID1 is part of the IMAGES survey that investigated the [OII]3726, 3729 Å kinematics using the FLAMES-GIRAFFE multi-object IFU (Yang et al. 2008).Based on their kinematic analysis, (Yang et al. 2008) classified ID1 as a regular disk.
The HST data for this galaxy clearly show that ID1 is a spiral galaxy (Fig.
A.1).The residuals in the HST fitting are due, indeed, to the presence of spiral arms (Fig. A.1).The channel maps and PVDs from the CO(2-1) ALMA data have the typical features of a rotating disk (e.g., see the S-shape profile in Fig. 6).With our analysis, we find a rotation velocity for ID1 of 200 km s −1 , consistent within the 1σ uncertainties with the value of 230 ± 33 km s −1 found by the IMAGES collaboration (Puech et al. 2008).

ID2
This galaxy is part of the PHIBBS2 survey with ID: XV53.
Based on the morphology derived by visually inspecting the HST I-band (HST/ACS F814W) images, ID2 is classified as a disk with asymmetric features (Freundlich et al. 2019).The HST image shows, in fact, a smooth disk and a clump in the south east direction (Fig. A.1). Lackner et al. (2014) interpreted the presence of these two clumps as due to two interacting systems and classify this galaxy as a late-stage minor merger.Here, we model the HST image using two off-center Sérsic components.The CO(3-2) emission is compact and does not extend to the south-east peak.Overall, the channel maps and the PVD are well reproduced by a rotating disk model.However, in the velocity field, the iso-velocity contours on the receding side of the velocity fields are more distorted than the approaching side.This is also visible in the PVDs: along the major-axis, there are two bright peaks at positive velocities, one symmetric with respect to the negative side at ∼ 250 km s −1 and the other located at ∼ 50 km s −1 , possibly corresponding to the interacting companion (black arrow in Fig. 6); the minor axis PVD is asymmetric.Considering the disturbances on one side of ID2, we fitted V rot and σ only using the approaching side of the galaxy.To perform such a fitting, we fixed the PA kin and the systemic velocity to the values obtained by using both sides of the galaxy.

ID3
This galaxy, also know as PACS819 is a Hershel detected galaxy (Rodighiero et al. 2011) whose global properties are extensively studied in the far-infrared wavelength range (Silverman et al. 2015(Silverman et al. , 2018;;Chang et al. 2020).In addition, by examining the BPT diagram (Baldwin et al. 1981), Silverman et al. (2015) show that PACS-819 is close to the line separating star-forming and AGN galaxies due to the strong [NII] emission line (Kewley et al. 2013).In the rest-frame UV HST image, ID3 shows two bright clumps, while the center of the CO(5-4) emission is located between them.Due to the different morphology between the CO and UV data, for ID3 we fixed the inclination angle to the one obtained with Cannubi.
The iso-velocity contours and the PVDs show that ID3 is kinematically lopsided, meaning that the velocity gradient in its approaching and receding sides are different from each other (Swaters et al. 1999;Bacchini et al. 2023).We therefore fit the approaching and the receding sides separately (note that in Fig. 6 we show the receding-side model).The rotating disk models reproduce the bulk of the emission from ID3.However, the disturbances of the two external contours at 2.5 and 5σ along the minor and major axis PVDs (see arrows in Fig. 6), as well as the residuals in the channel maps (Fig. C.1), indicate the presence of non-circular motions, in particular at negative velocities.In Fig. 6, we show with a green arrow the gas that is moving at lower rotation than those predicted by the model.The black % majo r 3.0 2.5 2.0 1.5 1.0 0.5 0.0 % majo r 3.0 2.5 2.0 1.5 1.0 0.5 0.0

ALPAKA:
% majo r 3.0 2.5 2.0 1.5 1.0 0.5 0.0  arrow shows, instead, gas that is rotating faster than predicted by the model, in particular in the inner regions.Such kinematic anomalies are usually attributed to extra-planar gas that is inflowing or outflowing from the disk (e.g., Fraternali et al. 2001;Heald et al. 2011;Marasco et al. 2019).The minor-axis PVD is not well reproduced by the model: the contours of the data and the model do not overlap and the shape of the model appears rounder than the data.This is because 3D Barolo tries to fit the features at anomalous velocities, visible in the major-axis PVD, by increasing the velocity dispersion values.Therefore, for ID3, we consider the σ values obtained by 3D Barolo as upper limits.
The velocity gradient in the velocity field of ID4 is likely due to the presence of rotation, but the distorted iso-velocity contours may be ascribed to the presence of a bar (e.g., Hogarth et al. 2021).However, due to the angular resolution of these data discriminating between radial motions driven by a bar or other mechanisms is not feasible.The best-fit Sérsic index obtained by modeling the HST data is consistent with a disk (n ∼ 1) rather than a bar structure (n ≲ 0.5).In addition, the HST data of ID4 show a clump in the north-west direction that is not detected in CO and it is an [OII] emitter at z ∼ 1.46, detected using a narrow-band filter (Hayashi et al. 2014;Ikeda et al. 2022).
ID5 has strong asymmetric features in the PVDs.The brightest regions of the major axis PVD are located at 0 and 100 km s −1 , respectively; they may be the nuclei of two unresolved interacting galaxies (see black arrows in Fig. 8).However, no clumpy structures are identified in the HST data.The minor axis PVD of ID5 is asymmetric and not reproduced by the rotating disk model.Another potential explanation of the irregularity of the ID5 PVDs is that this galaxy is a lopsided disk (Noordermeer et al. 2001).Follow-up observations at higher angular resolution are needed to robustly identify its nature.
Galaxies 6 -9 have, instead, the typical features of rotating disks.However, they show some asymmetries, especially along the minor-axis PVD (see arrows in Fig. 8), and distorted isovelocity contours in the velocity fields that may be driven by environmental effects (e.g., ram pressure stripping Lee et al. 2017;Zabel et al. 2019;Cramer et al. 2020;Bacchini et al. 2023), gravitational interactions (Boselli et al. 2022) or gas accretion (de Blok et al. 2014).Due to the quality of the current data, discriminating between these mechanisms is challenging.For instance, at low-z, a clear indication of ram-pressure stripping is a distortion of the gas distribution with respect to the stellar component (Lee et al. 2017;Boselli et al. 2022).On the contrary, it is expected that gravitational interactions have an impact on both the gas and stellar components.For ID6 -9, the total-flux maps appear smooth in all cases, likely because of the resolution of the observations, and aligned with the stellar distribution (Fig. 2).However, the HST data of ID7 show a close-by source (Fig. 2) that is not visible in the ALMA data.For ID7, the asymmetric features visible in the velocity field and in the minor-axis PVD (see black arrow in Fig. 8) are in the same direction of the nearby source, indicating that the asymmetries of the ID7 disk could be due to gravitational interactions.3D Barolo tries to reproduce these kinematic perturbations by inflating the velocity dispersion (see the comparison between the models and the data in the minor-axis PVDs).For this reason, in Sect.7, we consider the σ values derived by 3D Barolo for ID7 as upper limits.
ID8 has, instead, a close-by source that is visible both in HST and in CO(2-1) at low significance (yellow controur in Fig. 2; see also Ikeda et al. 2022).To recover this source using the 3D source finder described in Sect.5.1, we assumed SNRCUT of 2 and GROWTHCUT of 1.8.This object is located at an angular distance of about ∼ 1.0 ′′ (8 kpc) in projection from the main galaxy while its spectral distance from the systemic redshift of ID8 is of ∼ 800 km s −1 .We note that the SNR of the ID8's companion is not sufficient to study its kinematics.

ID10
This galaxy is part of the SHiZELS survey (Swinbank et al. 2012;Molina et al. 2017;Gillman et al. 2019) that map the Hα emission line with VLT/SINFONI.The ALMA observations of the CO(2-1) line used here were previously analyzed by Molina et al. (2019), who classify ID10 as a disturbed dispersiondominated disk.Here, we confirm the presence of a smooth gradient in the velocity field, but the strong disturbances and asymmetries in both PVDs, the systematic residuals in the channel maps, and the presence of bright emission only on one side of the major-axis (see arrow in Fig. 9) lead us to classify ID10 as uncertain.

ID11 -12
ID11 and 12 are members of two galaxy clusters that are recently studied in Cramer et al. (2023) (see also Noble et al. (2019) for ID11).Both galaxies have the typical features of a rotating disk (e.g., S-shape major-axis PVD) but they show some asymmetries in the minor-axis PVDs -despite at low significance, see arrows in Fig. 9 -and residuals in the channel maps, likely due to tidal disturbances or ram-pressure stripping.

ID13
This galaxy is a protocluster member.Its global properties (stellar masses, SFR, dust content) are studied in numerous papers (e.g., Hung et al. 2016;Zavala et al. 2019).In the HST cutout, two additional sources are visible, one located in the north east and the other in the south east directions.However, the redshifts of these sources are unknown and no dust or CO emission are detected with the ALMA data employed in this paper.The PVDs of ID13 show the typical features of a rotating disk.Despite this, ID13 has a relative difference between PA kin and PA HST of 53 deg.Two reasons may be responsible for such a large difference: 1. being almost face-on, a robust estimation of the morphological position angle is not straightforward; 2. the light contamination from the two sources in the HST cutout may bias the estimation of the Galfit parameters.

ID14
ID14, also known as BX610, has been extensively studied both in optical/near-infrared (Förster Schreiber et al. 2009, 2014;Tacchella et al. 2018) and far-infrared wavelength (Aravena et al. 2014;Bolatto et al. 2015;Brisbin et al. 2019).Its kinematics was previously analyzed using Hα data from VLT/SINFONI survey (Förster Schreiber et al. 2018).Based on Hα, ID14 is described as a rotating disk containing massive star-forming clumps (Förster Schreiber et al. 2011, 2018).It also exhibits signatures of gas outflows driven by a weak or obscured AGN in the center and by star formation at the location of the bright southern clump visible both in Hα and UV (Förster Schreiber et al. 2014). Recently, Brisbin et al. (2019) suggest that the high CO(7-6)/L IR ratio measured in BX610 can be ascribed to shock excitation caused by a recent merger event.
The data cube for ID14 contains both CO(4-3) and [CI]( 3 P 1 − 3 P 0 ).However, we show here only the model obtained by fitting CO(4-3) which has higher SNR than [CI]( 3 P 1 − 3 P 0 ).The CO(4-3) data show the presence of strong kinematic anomalies: the major axis PVD has the brightest emission only in the central regions; the minor-axis PVD is strongly axisymmetric, the iso-velocity contours are distorted in the south-east direction where a peak in the CO(4-3) distribution is visible (see arrows in Fig. 10).The latter may be due to an interacting companion.In addition, ID14 has a large difference of 65 deg between PA HST and PA kin .This misalignment is a further indication that ID14 is not a regularly rotating disk.However, the CO(4-3) kinematics is not sufficient to discriminate between the presence of a merging system and a rotating disk where the CO emission is contaminated by non-circular motions, likely driven by outflows.Therefore, in this paper, we put ID14 into the uncertain class.The combination of Hα and CO kinematics and the multi-wavelength morphological analysis will allow us to constrain the nature of ID14 (Deveraux et al. in prep.).

ID15
ID15 is known as the prototypical example of a compact starforming galaxy that is rapidly consuming its gas reservoir and is expected to evolve into a quiescent galaxy (Popping et al. 2017).The kinematics of the CO(6-5) emission line is studied in Talia et al. (2018) that identify the presence of a rotating disk with a velocity V = 320 +92 −53 km s −1 .Further, Loiacono et al. (2019) analyzed the Hα and [OIII] kinematics, finding evidence of rotation motions, with V ∼ 200 +17 −20 km s −1 for an inclination of 75 deg.Here, we assume an inclination of 60 deg derived from modeling the HST data.The CO(3-2) data clearly show the presence of a compact rotating disk.However, similarly to previous studies (Talia et al. 2018;Loiacono et al. 2019), the velocity dispersion values of ID15 are poorly constrained, having uncertainties of 50-60% at all radii because of the low angular resolution of ALMA data.We note that ID15 has distorted iso-velocity contours (Fig. 10) and residuals in the channel maps that may indicate the presence of a warp or non-circular motions.However, higher angular resolution data are needed to robustly analyze them and constrain the driving mechanisms.

ID16 -17
ID16 and 17 are two galaxies, separated by 3 ′′ and connected by a bridge of gas and dust (Fu et al. 2013).These sources are mildly magnified by two foreground galaxies, with magnification factors of 1.5.Due to their close proximity, ID16 and 17 were first identified as a unique bright SMG, HXMM01, in the Hershel Multi-tiered Extragalctic survey (Oliver et al. 2012).Subsequently, higher-angular resolution observations resolve HXMM01 into three sources, ID16, 17 and its companion (Fu et al. 2013), also visible in the total-flux maps in Fig. 2. The bright clump visible in the HST image at the southern end of ID17 has an SED that is consistent with either a less obscured galaxy at z = 2.3 or a physically unrelated contaminating source (Fu et al. 2013).The ALMA cube employed here contains two emission lines, CO(7-6) and [CI]( 3 P 2 − 3 P 1 ) that slightly overlap along the frequency axis.In Figs. 10 and  CO(7-6), characterized by higher SNR than [CI]( 3 P 2 − 3 P 1 ).A kinematic analysis of this complex system was already presented in Xue et al. (2018) who interpreted ID16 and 17 as two rotating disks.Despite ID16 showing some symmetric features along the minor-axis PVD (Fig. 10), the major-axis PVD is strongly asymmetric, being bright only at negative velocities.Similarly, ID17 is strongly asymmetric as it is interacting with a companion, clearly visible in the total-flux map (Fig. 2) and major-axis PVD (Fig. 11) located at 1 ′′ in the north direction with respect to the ID17 center.Due to the presence of strong disturbances, we classify ID16 as uncertain and 17 as an interacting system.

ID18
ID18 lies in a well-studied protocluster.This structure was originally found in the largest extragalactic Hershel survey (H-ATLAS).Ivison et al. (2013Ivison et al. ( , 2019) ) find that HATLAS J084933.4+021443hosts an exceptionally high SFR ∼ 7000 M ⊙ yr −1 , spread over 5 galaxies (W, T, C, M, E), confined within a scale of ∼100 kpc .The ALMA data set used here maps galaxy W and T at high resolution and SNR, while the data quality is not sufficient for performing a kinematic analysis on galaxies C, M, and E. Further, galaxy T is not included in the ALPAKA sample as it is a strongly gravitationally lensed source, while galaxy W, that is, ID18 is not lensed (Ivison et al. 2013).
The CO(4-3) data analyzed in this paper indicates that ID18 is a regularly rotating disk.
6.12.ID19 -21 These galaxies are members of the forming cluster CLJ1001 at z = 2.51 (Wang et al. 2016;Gómez-Guijarro et al. 2019;Champagne et al. 2021).The CO(3-2) kinematics of ID19 -21, along with another galaxy member of the same cluster (130901), are recently studied by Xiao et al. (2022) using the same data employed here.We exclude 130901 from our analysis due to the low SNR.The comparison between the 3D Barolo model and the data indicates that the bulk of the motions in ID19 and 20 are reproduced by a rotating disk.The value of the velocity dispersion for ID19 is consistent within 1.5-σ with the one found by Xiao et al. (2022), while we find a value of rotation velocity which is a factor of 1.5 smaller than the value found by Xiao et al. (2022).This discrepancy is due to the different inclination angles used to recover the intrinsic rotation velocity.On the contrary, for ID20 we find rotation velocities consistent with those found in Xiao et al. (2022), while the velocity dispersion values are not well constrained, having relative uncertainties ≳ 60%, due to the low angular resolution and SNR of the data.ID21 has, instead, asymmetries in the minor-axis PVD, while its major-axis PVD lacks the bright features typical of rotating disks.One potential explanation for this latter feature is that ID21 has a rotation curve that is slowly rising in the inner regions.However, since the angular resolution of the data employed in this paper does not allow us to robustly test this sce-nario, we classify this target as an uncertain galaxy.Xiao et al. (2022) interpreted, instead, the gradient in the velocity map of ID21, also visible in Fig. 12, as due to the presence of a rotating disk.

ID22, 26 -27
These galaxies are recently studied by Cassata et al. (2020) using ALMA observations with angular resolutions of 0.6 ′′ .The 0.17 ′′ observations of the CO(5-4) used in our work allow us to resolve their kinematic structures.A rotating disk is a good description of the data only for ID22, despite the strong misalignment between the morphological position angle derived from HST data and PA kin .However, we note that, similarly to ID3, the HST data for ID22 probe the rest-frame UV and, therefore, the morphological parameters may be biased due to the contamination from the light of young stars and/or by dust attenuation.
ID26 has irregularities in the minor-axis PVD (see arrow in Fig. 13), systematic residuals in the channel maps and distorted iso-velocity contours in the approaching side of the galaxy.These features may be due to the presence of a warp in a disk.However, the angular resolution of the observations do not allow us to verify this scenario.Therefore, to be conservative, we classify ID26 as uncertain.On the contrary, we classify ID27 as a merger.This target has a very irregular PVDs and velocity field and a difference between PA ALMA and PA kin of 45deg, another indication that it is, likely, a non-resolved interacting system.

ID23 -25
These galaxies are members of the SSA 22 protocluster (Steidel et al. 1998) that was extensively studied with ALMA (Umehata et al. 2015(Umehata et al. , 2017(Umehata et al. , 2018)).Lehmer et al. (2009) identified X-ray luminous AGNs in ID 23 and 25 using observations from the Chandra Space telescope.Extended Lyman-α emission from multiple filaments between galaxies within SSA 22 are identified using MUSE observations (Umehata et al. 2019) and are thought to be responsible for the accretion of gas within the protocluster and the growth of galaxies and their supermassive black holes.The high-angular resolution observations employed here allow us to clearly see the typical features of edge-on rotating disks in ID 23, 24, and 25.The three galaxies have very extended S-shape major-axis PVDs (Fig. 10).ID24 has gas emission at −400 km s −1 along the minor-axis PVD that is not reproduced by the model (see arrow in Fig. 10).However, due to the high inclination of the galaxy, these asymmetries are difficult to interpret.
The kinematics of ID28 is peculiar as it shows the typical features of a rotating disk (e.g., S-shape along the major-axis PV, diamond shape along the minor axis PV) but also strong emission in the inner 1 kpc regions (see arrows in the PVDs, Fig. 13) that indicates the presence of gas moving at a line-of-sight velocity of 900 km s −1 relative to the systemic velocity.This emission, not reproduced by the symmetric rotating disk model (see The first possibility is that the CO distribution is asymmetric between the approaching and the receding sides and there is an inner rise of the rotation curve caused by the presence of a compact bulge.The second possibility is that this emission is due to non-circular motions driven by outflows.We note that strong outflow motions were already identified in ID28 from the analysis of rest-frame UV spectrum (Ginolfi et al. 2022).
The V values in the inner regions for this galaxy should thus be taken with caution as the emission from the disk may be strongly contaminated by the one from the outflow.Furthermore, to reproduce the emission at high velocities in the inner regions, 3D Barolo inflates the velocity dispersion (e.g., the inner contours of the model are rounder than the data contour in the minor-axis PVD).For this reason, we consider the σ values as upper limits.By analyzing the same CO(6-5) data presented in this paper, Ginolfi et al. (2022) find velocity dispersion and rotation velocity profiles consistent with those found here.

Kinematics of the disk subsample
In this subsection, we discuss the kinematic properties of the 19 ALPAKA secure disks.The rotation velocities and velocity dispersions of mergers are, instead, unreliable as they are derived under the assumption that the observed kinematics is dominated by circular motions.As such, we decide not to show nor discuss them.For the same reason, the kinematic properties of the 7 galaxies in the uncertain class are not shown here.In Fig. 14, we show the profiles of the rotation velocity and velocity dispersion as derived by 3D Barolo (Sect.5.1).The angular resolution of the ALPAKA data allows us to sample the kinematic profiles with only 2 independent resolution elements for 6 AL-PAKA galaxies and ≳ 3 for the remaining 13 targets.For the 16 ALPAKA disks with HST data, we quantified the relative extension of the kinematic profiles with respect to the rest-frame optical/UV effective radius R e obtained with Galfit.The extension of the kinematic profile is estimated by adding half of the beam size to the outermost radius used for the kinematic fitting.The comparison between the CO/[CI] and rest-frame optical/UV extension is not straightforward as differences may depend not only on the heterogeneous sensitivities of ALPAKA data but also on the different emission lines that we are using to trace the kinematics.Low-J CO transitions may trace, in fact, more diffuse, extended molecular gas than high-J CO transitions (e.g., Lagos et al. 2012).In addition, the values of R e may vary with the rest-frame wavelength (e.g., Vulcani et al. 2014;Kennedy et al. 2015).Considering these caveats, we find that only 7 sources -4 with CO(2-1), 1 with CO(3-2), 1 with [CI]( 3 P 2 − 3 P 1 ), 1 with CO(6-5) -have kinematic profiles that extend up to radii ≳ R e (see gray dotted lines in Fig. 14), while for the others the ALMA data trace only the innermost regions.
Both galaxies with compact and extended kinematic profiles show a variety of shapes: flat (e.g., ID12, 19), slowly increasing and then flattening (e.g., ID1, 2), declining in the inner regions and then flattening (e.g., ID13, 28).We caution the reader that the velocity profiles in Fig. 14 show the rotation and not the circular speed.The latter is a direct proxy of the gravitational potential and can be derived from V rot after applying the asymmetric drift correction which is a function of the velocity dispersion profile and gas distribution (e.g., Iorio et al. 2018).Such estimates, as well as the derivation of the gravitational potential, will be part of a future paper in the ALPAKA series.
For most of the ALPAKA disks, there is an indication of velocity dispersion profiles with a declining trend, where the σ values are higher in the inner regions than at the outermost radius.To quantify this trend, we fitted the velocity dispersion profiles using a linear function and we find negative slopes for 14/19 disks.However, due to the large uncertainties of the velocity dispersion values, for 9/14 galaxies these negative slopes are consistent with 0 within the 1-σ uncertainties.The decrease of σ with radius is also observed in the CO profiles of local disk galaxies (Bacchini et al. 2020), and it is ascribed to the radial decline of the SFR surface density.The analysis of the ALPAKA galaxies shows, therefore, that caution must be taken when fitting the velocity dispersion profile assuming a constant value across the disks.This assumption, often used at high-z (e.g., Bouché et al. 2015;Neeleman et al. 2020), may result in the estimation of velocity dispersion σ that is mainly influenced by the emission from the bright regions at the center of a galaxy, in turn resulting in larger σ measurements compared to the average values.On the other hand, another usually used assumption is to measure the velocity dispersion from the outer regions of the observed velocity dispersion field, where the contribution from the beamsmearing effect is minimized (e.g., Wisnioski et al. 2015;Harrison et al. 2017).In this latter case, it is not straightforward to evaluate whether the resulting σ values overestimate or underestimate the average σ across the galaxy disks, as it may strongly depend on the combination of galaxy size, angular resolution, shape of the rotation velocity and velocity dispersion profiles (e.g., Di Teodoro & Fraternali 2015; Burkert et al. 2016).
To facilitate the comparison with previous works, in Tab. 5, we refactor our result in terms of the kinematic parameters mostly common in the literature: the maximum values of the rotation velocity V max ; the radial average values of the velocity dispersion σ m ; the values of velocities and dispersions computed by averaging the two outermost values in their profiles, V ext and σ ext ; the ratios V max /σ m and V ext /σ ext .The latter are used to define the rotational support of galaxies: V/σ ≳ 2 is the threshold used to indicate rotation-dominated systems (Förster Schreiber et al. 2018;Wisnioski et al. 2019).In Fig. 15, we show σ m and V max /σ m as a function of redshift for the ALPAKA disks.In this figure, galaxies hosting an AGN are indicated with black dots.As discussed in Appendix 5.2, some ALPAKA disks have kinematic anomalies that 3D Barolo tries to reproduce by inflating the velocity dispersion.Therefore, for at least three galaxies (i.e., ID3, 7, 28), the velocity dispersion values can be considered upper limits (and thus lower limits for V/σ).
As discussed in Sect.4, ALPAKA sample is biased towards galaxies hosting energetic mechanisms (e.g., stellar and AGN feedback, interactions due to the environment) that are expected to either boost the velocity dispersion values or completely destroy the disk structures (e.g., Dubois et al. 2016;Penoyre et al. 2017;Gurvich et al. 2023;Kretschmer et al. 2022).Despite this, we find that σ m ranges between 10 and 52 km s −1 for all galaxies, except for ID 15 and 28, two AGN hosts with values of σ m of 99 km s −1 and an upper limit of 167 km s −1 , respectively.The median value of σ m , obtained after excluding only the three upper limits is 35 +11 −9 km s −1 (dashed gray line in Fig. 15, upper panel).The V max /σ m values range between 5 and 21, with a median value of 9 +7 −2 (dashed gray line in Fig. 15, bottom panel).The circles show the location of the rings used for fitting the data with 3D Barolo and they can be considered independent from each other.The dashed vertical lines show the location of the optical/UV effective radius R e for ALPAKA galaxies with HST data and with R e comparable to or smaller than the extent of the CO/[CI] emission.
In the second paper of this series, we will compare the distribution of σ and V/σ of the ALPAKA sample with redshift-matched samples of galaxies with warm gas kinematics, and we will study the evolution of the cold gas kinematics across cosmic time.

Selection effects and the impact of energetic mechanisms on the disk properties
The selection criteria used to build the ALPAKA sample are not based on the global physical properties of the sources, but on the quality of the available data in the ALMA archive.This, combined with the intrinsic faintness of CO transitions and the limited sensitivities of ALMA, results in a sample that is biased towards massive (i.e., ≳ 10 10 M ⊙ ) main-sequence or starburst galaxies, mostly in overdense regions.Numerous studies show that the gas-to-stellar mass ratio are typically higher in cluster than in field galaxies at z ≳ 1 (Noble et al. 2017;Hayashi et al. 2018;Tadaki et al. 2019).Further, at least 7 out of the 28 AL-PAKA galaxies are known to host an AGN.According to current models of galaxy formation and evolution, the growth of high-z massive star-forming galaxies is mainly driven by galaxy mergers and intense gas accretion which drive large amounts of gas towards their centers, boosting the SFR and the growth of supermassive black holes.In the following subsections, we discuss the impact of different astrophysical mechanisms on the dynamics of the ALPAKA disks.
Table 5. Global kinematic parameters of the ALPAKA galaxies classified as disks.

Starbursts
Among the 12 starbursts in ALPAKA, we find that 8 are rotating disks.The dynamical time scale of these disks is ≈ 10 Myr6 , a factor of 10 smaller than the typical depletion timescales (∼ 100s Myr) of starburst galaxies at these redshifts (Scoville et al. 2017;Liu et al. 2019).Despite the small number statistics, the presence of disks among the ALPAKA starbursts suggests that after a merger or an intense accretion event, galaxies quickly transition into a stable dynamical stage and form a disk with V/σ ∼ 10.These high V/σ values are consistent with recent findings of dynamically cold disks among dusty starburst galaxies at z ≳ As discussed in Sect. 1, the gas tracer may impact the derived rotation velocity and velocity dispersion.A detailed comparison with the results from the literature is beyond the scope of the paper and is part of a future study of the ALPAKA series.

AGN feedback
Among the subsample of 7 AGN-host galaxies, only one belongs to the uncertain class while the others are rotating disks.This result might indicate that AGN feedback does not prevent the formation of a rotating disk and it may have also a limited impact on the ISM properties of the host galaxies.With the exception of ID28 (see Sect.7.1), whose kinematic properties are likely contaminated by outflows, all AGN hosts have velocity dispersions consistent the with rest of the disk subsample.Interestingly, ID28 is the only ALPAKA galaxy belonging to the rare, extreme population of hyper-luminous dust-obscured AGN (Eisenhardt et al. 2012;Wu et al. 2012).The latter are thought to represent a short-lived phase during which there is an evolutionary transition from dusty starburst galaxies to UV-bright quasars (e.g., Wu et al. 2018;Díaz-Santos et al. 2021).

Overdense environments
Previous results from the literature (e.g., Ivison et al. 2013;Umehata et al. 2015, see Sect. 6 for details) indicate that most ALPAKA galaxies are in overdense regions, clusters or groups (13/28) or protoclusters (6/28).In these dense environments, galaxy mergers and tidal interactions could be efficiently enhanced with respect to the field (e.g., Merritt 1983;Moore et al. 1996;Kronberger et al. 2006;Cortese et al. 2021).Recent studies show that environmental effects, like ram-pressure stripping, are responsible for the kinematic asymmetries of the molecular gas distribution and kinematics in cluster galaxies at z ∼ 0 (e.g., Lee et al. 2017;Cramer et al. 2020;Bacchini et al. 2023).Among the 13 ALPAKA galaxies in clusters and groups, we classify 1 as an interacting system, 4 as uncertain, and the remaining 8 as disks.With the exception of ID8 and 9, the remaining cluster disks have kinematic anomalies likely driven by environmental mechanisms.On the contrary, among the 6 protocluster galaxies, 5 are disks and 1 is a disk with outflow contamination (ID28).This finding suggests that environmental mechanisms may not have a relevant impact on the ISM of galaxies during the early stages of cluster formation.However, the statistics is too low to draw robust conclusions.

CO and [CI] transitions
The ALPAKA sample comprises a large variety of CO and [CI] transitions.Overall, these lines trace the cold molecular gas (T < 100 K), a phase of the ISM different from the photoionized gas observed through forbidden (e.g., [OIII], [OII]) or recombination (e.g., Hα) lines.However, the different critical densities of the CO rotational transitions (e.g., 10 4 cm −3 for the CO(2-1) to ≈ 5×10 5 for the CO(7-6)) make them sensitive to distinct ISM conditions.Low-J CO transitions (J ≲ 3) trace the diffuse molecular component at T ≈ 5 -33 K, while higher J transitions trace increasingly dense gas at T ≈ 100 K (Carilli & Walter 2013).
In addition, high rotational transition levels (e.g., CO(5-4) and higher) if excited, may trace the presence of gas heated by various mechanisms (e.g, X-rays, van der Werf et al. (2010); shocks induced by mergers, radio jets, and supernova-or AGN-driven outflows, Kamenetzky et al. (2016); Vallini et al. (2019)).Sim-ilarly, [CI]( 3 P 1 − 3 P 0 ) is more sensitive to the cold, diffuse gas, than [CI]( 3 P 2 − 3 P 1 ) (Valentino et al. 2020(Valentino et al. , 2021)).Among the ALPAKA targets, 20 have CO(2-1), CO(3-2), CO(4-3) or [CI]( 3 P 1 − 3 P 0 ) observations, with a large fraction (15/20) containing a rotating disk.Instead, among the 8 AL-PAKA galaxies with CO(5-4), CO(6-5), CO(7-6) and [CI]( 3 P 2 − 3 P 1 ) observations, four are interacting or uncertain galaxies with strong asymmetries in their kinematics and 2/4 of the disks (i.e., ID3 and 28) have kinematic anomalies potentially driven by outflows.Despite the small statistics of the high-J subsample not allowing us to draw significant conclusions, the different fractions of regularly rotating disks in the low and high-J subsamples suggest that the high-J transitions tend to trace non-circular motions.On the other hand, the comparison between the distributions of σ m and V/σ for the low and high-J subsamples (green circles and orange squares in Fig. 15) suggest that gas components traced by these transitions are characterized by similar levels of turbulence.Follow-up observations of the two ALPAKA subsamples at high and low-J transitions, respectively, would be needed to gain insights into the impact of the specific transitions on the turbulence and the presence of non-circular motions.

Summary and conclusion
Studying the kinematics of galaxies using cold gas tracers is crucial for gaining insights into the formation and evolution of structures across cosmic time.Nevertheless, before the advent of the Next Generation Very Large Array (ngVLA; Carilli & Shao 2018) in 2030s and the improvement of the ALMA sensitivity thanks to the correlator upgrade (Carpenter et al. 2020), highresolution observations of CO transitions at z ≳ 0.5 are going to be feasible only for a few bright targets, due to the long integration times required to achieve sensitivities sufficient to recover robust kinematic measurements.For this reason, systematic investigations of the cold gas kinematics at z ∼ 0.5 − 3.5 are still lacking.In this context, the ALPAKA project aims at alleviating this limitation by providing a systematic derivation of the kinematic properties of a sample of 28 star-forming galaxies at z ∼ 0.5 − 3.6 (median z of 1.8).Data for the ALPAKA galaxies are obtained by collecting ALMA archive high-data quality observations (median angular resolution of 0.25 ′′ ) of CO and [CI] transitions.
In this paper, we present the ALPAKA galaxies and derive their global (M ⋆ , SFR), morphological and kinematic properties.We find that ALPAKA galaxies have high stellar masses (M ⋆ ≳ 10 10 M ⊙ ) and SFRs that range from 10 to 3000 M ⊙ yr −1 .A large fraction of ALPAKA are starbursts (12 out of the 25 galaxies with good estimates of M ⋆ and SFR), while the remaining are on the Main sequence of star-forming galaxies.Further, after combining a heterogeneous set of works from the literature, we conclude that 19/28 ALPAKA galaxies lie in overdense regions (clusters, protoclusters, groups).We exploited the ALMA data cubes and modeled the kinematics of the ALPAKA targets using 3D Barolo, a forward-modeling tool which fits galaxy kinematics assuming a rotating disk.Our analysis and the main results can be summarized as follow.
-We divide the 28 ALPAKA galaxies into three kinematic classes: rotating disks, mergers and uncertain.We base this classification on the visual inspection of the PVDs, channel maps, velocity fields and the comparison with the pure rotating disk model obtained with 3D Barolo.This kinematic classification is confirmed by PVsplit, a tool that discerns between disks and mergers based on the measurements of asymmetries of the major-axis PVD.We find that 19 targets are rotating disks, 2 are interacting systems, and 7 are uncertain.For the latter, the data quality is not sufficient to reliably determine their kinematic states.-We present the rotation curves of the disk subsample; their semi-major axes are sampled by at least 3 independent resolution elements for 13 out of the 19 disks.In these cases, we can trace the shape of the rotation curves, and we find that they flatten.-We present the velocity dispersion profiles of the 19 disks.
For 14 of them, there is a marginal indication that the profiles are not constant and the σ values are larger in the inner regions and decline up to a factor of 3 in the outskirts, indicating that there may be a radial change in the intensity of the mechanisms driving turbulence (e.g., decline in the SFR surface density).
-We present values for the global kinematics of the disk galaxies in ALPAKA: the maximum rotation velocity range from 204 to 509 km s −1 ; the velocity dispersion averaged across the radii, σ m , ranges from 10 to 100 km s −1 ; the V/σ range from 5 to 21.The ALPAKA disk sample has a median σ m of 35 +11 −9 km s −1 and an overall large V/σ with a median value of 9 +7 −2 .This work is the first of a series that will allow us to fully exploit the ALMA data presented here to characterize the ALPAKA galaxies.In particular, we plan to use the kinematic analysis described in this manuscript, in combination with the HST data presented in this paper and JWST observations that have been taken for some ALPAKA sources, to: systematically compare the dynamics from warm and cold gas tracers; investigate the evolution of cold-gas velocity dispersion with redshift and study the energy sources of turbulence; infer the dark matter halo properties and content within the ALPAKA sample; derive the dynamic scaling relations and study their evolution; constrain the impact of cold outflows and AGN feedback on the gas disks; compare the sizes and morphologies of the stellar, gas, and dust distributions.
Considering the limitations of current telescopes in the mm and sub-mm wavelength range, the ALPAKA sample -although biased towards massive, actively star-forming galaxies in overdense environment -will constitute a legacy for high-angular resolution studies in the next decade.We note, indeed, that highangular resolution observations of the cold gas kinematics of typical star-forming galaxies at z ∼ 2 are not feasible with ALMA, even considering the current 50 hours limit for normal observing programs.For instance, a Milky-Way progenitor with stellar and gas masses of 1 × 10 10 M ⊙ at z ∼ 2.2 (van Dokkum et al. 2013) is expected to have a CO(3-2) luminosity of 10 9 K km s −1 pc 2 and a flux of 0.08 Jy km s −1 , that is a factor of ≈ 10 smaller than the CO(3-2) fluxes of ALPAKA galaxies at similar redshifts (e.g., ID19-21).To observe such a Milky-Way progenitor at the same level of detail we have for ID19 (i.e., similar SNR and angular resolution) an on-source integration time of 300 hours -a factor of 100 times the on source time of ID19 -would be required.On the other hand, ngVLA, with a factor of 6 increase in sensitivity over ALMA (Carilli & Neeleman 2022), will enable the study of the dynamics of Milky-like progenitors at z ∼ 2 in 8 hours of integration time 7 .
We make the data cubes of the ALPAKA targets and the kinematic profiles derived in this manuscript publicly available, 7 ngVLA will cover the CO(2-1) line at z ∼ 2. Therefore, we compute the comparison of the required integration times assuming I CO(3−2) ≈ I CO(2−1) .
trusting that the community will exploit this sample well beyond what we had originally envisioned.

Fig. 1 .
Fig. 1.Distribution of the redshift (upper panel) and the rotational transitions of CO and atomic fine structure line transitions of [CI] (bottom panel) for the ALPAKA data set.

Fig. 2 .
Fig. 2. HST images with contours of the CO or [CI] total-flux maps from ALMA for the 23 ALPAKA galaxies with HST observations.The lowest white contour is a "pseudo-contour" -see Appendix B in Roman-Oliveira et al. (2023) for details -at 4 RMS.The additional yellow contour in the ID8 panel shows a second source detected in CO(2-1) (see Sect. 6.4 for details).The CO transitions and HST filters shown here are listed in Tables2 and 3, respectively.The white bar in the bottom left of each panel shows a scale of 1 ′′ .

Fig. 5 .
Fig. 5. Left panel: comparison between the inclination angles derived from the morphological fitting of HST and ALMA data using Galfit and Cannubi respectively.Central panel: comparison between the position angles derived from the morphological fitting of HST data with Galfit and kinematic fitting of ALMA data with 3D Barolo.Right panel: comparison between the position angles derived from the morphological fitting of ALMA data with Cannubi and the corresponding kinematic fitting.The gray line shows the 1:1 relation and the gray dotted lines show deviations at ±20 deg (left panel) and ±30 deg (central and right panels).

Fig. 6 .
Fig.6.For each target with ID in the upper left, we show from top to bottom: the total-flux map and the velocity field and the major and minor-axis PVDs.In the total-flux maps, the first external contour is a "pseudo-contour" (seeAppendix B in (Roman-Oliveira et al. 2023)  for details) at 4 RMS.In the velocity field, the black lines show the iso-velocity contours, with the thickest one indicating the systemic velocity.The black dashed and gray dotted lines are the kinematic and morphological position angles, respectively.For the morphological position angle, we show the one obtained from fitting HST data when available and the total-flux map otherwise.The beam is shown in the bottom left.In the PVDs, the y-axis shows the line-of-sight velocities centred on the systemic velocity or redshift and the x-axis shows the distance with respect to the center of the kinematic model.The contours for the data (solid black) and the model (red) are at[1, 2, 4, 8, 16, 32] × 2.5 RMS.The gray dotted contours are at -2.5 RMS.The white circles are the best-fit line-of-sight rotation velocities.For each target, we indicate in the velocity field whether it is classified as a disk (D), merger (M) or uncertain (U).The arrows show kinematic anomalies that are described and discussed in detail in Sect.6.

Fig. 7 .
Fig. 7. Distribution of the ALPAKA galaxies in the PVsplit parameter space.The gray circles and squares show simulated disks and mergers, respectively fromRizzo et al. (2022).The red plane divides the regions occupied by the two kinematic classes.The colored markers show ALPAKA galaxies classified as disks, merger and uncertain (right panel) based on the visual inspection of their PVDs and channel maps.Two slightly different projections are shown for better visualizing the position of all ALPAKA galaxies.

Fig. 14 .
Fig.14.Rotation velocity (upper panels, pink) and velocity dispersion profiles (bottom panels, green) for the ALPAKA galaxies classified as disks.The circles show the location of the rings used for fitting the data with 3D Barolo and they can be considered independent from each other.The dashed vertical lines show the location of the optical/UV effective radius R e for ALPAKA galaxies with HST data and with R e comparable to or smaller than the extent of the CO/[CI] emission.

Fig. 15 .
Fig. 15.Distribution of the ALPAKA disks in the velocity dispersionredshift plane (upper panel) and rotation-to-velocity dispersion ratioredshift plane (bottom panel).The values of σ m and V max /σ m (Table5) are plotted here.The redshifts of ID3, ID6 -9, and ID19 are shifted by |∆z| ≲ 0.25 for a better visualization of all the data points.Galaxies with kinematics from CO(2-1), CO(3-2) or [CI]( 3 P 1 − 3 P 0 ) are shown with green circles, while galaxies with CO(5-4), CO(6-5) or [CI]( 3 P 2 − 3 P 1 ) are shown with orange squares.The dashed gray line (and the gray area) show the median (and the 16th and 84th percentiles) σ m and V max /σ m values for the ALPAKA disks.
(Rizzo et al. 2020;Lelli et al. 2021;Rizzo et al. 2021;Roman-Oliveira et al. 2023) and they are in contrast with previous studies reporting values of V/σ ≲ 3 for this galaxy population(Swinbank et al. 2011;Alaghband-Zadeh et al. 2013;Olivares et al. 2016;Birkin et al. 2023).However, we note that the latter were mostly obtained with low-angular resolution observations and with no(Alaghband-Zadeh et al. 2013;Olivares et al. 2016) or suboptimal beam-smearing correction(Swinbank et al. 2011;Birkin et al. 2023, see Sect. 1).Another potential reason for the systematic difference with previous studies of z ≲ 4 starbursts is the emission lines employed to trace the galaxy kinematics.With the exception of the target studied in(Swinbank et al. 2011) with CO transitions, the others works employ restframe optical emission lines(Alaghband-Zadeh et al. 2013;Olivares et al. 2016;Birkin et al. 2023) tracing warm, ionized gas.

Table 2 .
Description of the ALMA observations and datasets of the ALPAKA sample.

Table 3 .
Properties of the ALPAKA sample and description of the HST dataset.
Noble et al. (2017)015)e distance of the ALPAKA galaxy with respect to the main-sequence relation, ∆ MS = log(SFR/SFR MS ), where SFR MS is the SFR defined bySchreiber et al. (2015)at fixed redshift and stellar mass.Column 5 indicates whether the galaxy is a main-sequence (MS) or a starburst (SB).To be conservative, galaxies with ID4, 11, 12, 25 are classified as starburst as their ∆ MS values are ≳ log(4) = 0.6 within the 1-σ uncertainties.Column 7 shows the flux for the CO or [CI] transitions listed in Table1.The last column shows the rest-frame effective wavelength probed by the HST data.*Fortargets 11 and 12, the M ⋆ and SFR values are taken fromNoble et al. (2017).

Table 4
, we report the best-fit geometrical parameters, PA HST , i HST .In Figs.A.1 and A.2, we show the HST data and the corresponding Galfit models and residuals.-Cannubion ALMA data.Cannubi is a Markov Chain Monte

Table 4 .
Geometric parameters of the ALPAKA galaxies.