# The ALPINE-ALMA [C ii] Survey: kinematic diversity and rotation in massive star-forming galaxies at z 4.4–5.9

G.C. Jones, D. Vergani, M. Romano, M. Ginolfi, Y. Fudamoto, M. Bethermin, S. Fujimoto, B.C. Lemaux, L. Morselli, P. Capak, et al.

► To cite this version:

G.C. Jones, D. Vergani, M. Romano, M. Ginolfi, Y. Fudamoto, et al. The ALPINE-ALMA [C ii] Survey: kinematic diversity and rotation in massive star-forming galaxies at z 4.4–5.9. Mon.Not.Roy.Astron.Soc., 2021, 507 (3), pp.3540-3563. 10.1093/mnras/stab2226. hal-03210320

HAL Id: hal-03210320

https://hal.science/hal-03210320

Submitted on 4 May 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# The ALPINE-ALMA [C II] Survey: kinematic diversity and rotation in massive star-forming galaxies at $z \sim 4.4–5.9$

G. C. Jones<sup>1,2*</sup> D. Vergani,<sup>3</sup> M. Romano,<sup>4,5</sup> M. Ginolfi<sup>6</sup> Y. Fudamoto<sup>7,8</sup>

M. Béthermin,<sup>9</sup> S. Fujimoto<sup>10,11</sup> B. C. Lemaux<sup>12</sup> L. Morselli,<sup>4,5</sup> P. Capak,<sup>13</sup>

P. Cassata,<sup>4</sup> A. Faisst,<sup>13</sup> O. Le Fèvre,<sup>9</sup> D. Schaerer,<sup>14</sup> J. D. Silverman,<sup>15,16</sup> Lin Yan,<sup>17</sup>

M. Boquien,<sup>18</sup> A. Cimatti,<sup>19,20</sup> M. Dessauges-Zavadsky,<sup>14</sup> E. Ibar,<sup>21</sup> R. Maiolino,<sup>1,2</sup>

F. Rizzo<sup>10,11</sup> M. Talia<sup>3,19</sup> and G. Zamorani<sup>3</sup>

*Affiliations are listed at the end of the paper*

Accepted 2021 July 29. Received 2021 July 29; in original form 2021 April 7

## ABSTRACT

While the kinematics of galaxies up to $z \sim 3$ have been characterized in detail, only a handful of galaxies at high redshift ($z > 4$) have been examined in such a way. The Atacama Large Millimeter/submillimeter Array (ALMA) Large Program to INvestigate [C II] at Early times (ALPINE) survey observed a statistically significant sample of 118 star-forming main-sequence galaxies at $z = 4.4–5.9$ in [C II] 158 µm emission, increasing the number of such observations by nearly 10×. A preliminary qualitative classification of these sources revealed a diversity of kinematic types (i.e. rotators, mergers, and dispersion-dominated systems). In this work, we supplement the initial classification by applying quantitative analyses to the ALPINE data: a tilted ring model (TRM) fitting code ($^{3D\,BAROLO$), a morphological classification (Gini-$M_{20}$), and a set of disc identification criteria. Of the 75 [C II]-detected ALPINE galaxies, 29 are detected at sufficient significance and spatial resolution to allow for TRM fitting and the derivation of morphological and kinematic parameters. These 29 sources constitute a high-mass subset of the ALPINE sample ($M_* > 10^{9.5} M_\odot$). We robustly classify 14 of these sources (six rotators, five mergers, and three dispersion-dominated systems); the remaining sources showing complex behaviour. By exploring the $G-M_{20}$ of $z > 4$ rest-frame far-infrared and [C II] data for the first time, we find that our 1 arcsec $\sim 6$ kpc resolution data alone are insufficient to separate galaxy types. We compare the rotation curves and dynamical mass profiles of the six ALPINE rotators to the two previously detected $z \sim 4–6$ unlensed main-sequence rotators, finding high rotational velocities ($\sim 50–250$ km s$^{-1}$) and a diversity of rotation curve shapes.

**Key words:** galaxies: evolution – galaxies: high-redshift – galaxies: kinematics and dynamics.

## 1 INTRODUCTION

The past century has seen a massive broadening of scope in the study of galaxy kinematics. Early works focused on spectral line emission from single nearby galaxies (e.g. Pease 1918; Burbidge, Burbidge & Prendergast 1959), while current studies are able to probe the kinematics of many local objects (e.g. Garrido et al. 2002; Conselice et al. 2005; de Blok et al. 2008; Shapiro et al. 2008; Puech 2010; Gómez-López et al. 2019; Korsaga et al. 2019; den Brok et al. 2020) or single objects in much greater detail (e.g. Carignan et al. 2006; Chemin, Carignan & Foster 2009; Cramer et al. 2019; North et al. 2019; Braine et al. 2020). This detailed analysis has also been extended to galaxies at intermediate redshift (i.e. $z \sim 1–3$; $t_H \sim 6–2$ Gyr) with integral field spectroscopy, mostly in the rest-frame optical and near-infrared (NIR; e.g. Förster Schreiber et al. 2006; Epinat et al. 2012; Burkert et al. 2016; Harrison et al. 2017; Molina et al. 2017; Swinbank et al. 2017; Turner et al. 2017; Loiacono et al. 2019; Wisnioski et al. 2019). These studies have revealed galaxies with a variety of kinematic types (i.e. rotators, mergers, and dispersion-dominated sources), as well as inflows and outflows.

One major finding from studies up to $z \sim 2$ is that most of these galaxies show evidence for significant dark matter haloes. By combining KMOS and MUSE data for $\sim 1500$ star-forming galaxies

<footnote>* E-mail: gj283@cam.ac.uk</footnote>
(SFGs) at $z \sim 0.6-2.2$, Tiley et al. (2019) find that galaxies in this epoch mainly feature flat or rising rotation curves at large radii. This is quite similar to the rotation curves of local galaxies, which show nearly constant or only slightly declining velocities at large radii (e.g. de Blok et al. 2008), indicating the additional gravitational force from an underlying dark matter halo (e.g. Rubin & Ford 1970).

While a multitude of galaxies at $z < 3$ have been well characterized kinematically, it is only recently, with the advent of the Atacama Large Millimeter/submillimeter Array (ALMA) and the Jansky Very Large Array (JVLA), that we are able to characterize the kinematics of high-redshift ($z > 4$) galaxies, using rest-frame far-infrared (FIR) emission lines. This includes the clumpy nature of $z \sim 5-7$ galaxies observed with ALMA in [C II] emission (e.g. Carniani et al. 2018), rotational model fitting of $z \sim 4-6$ galaxies (e.g. De Breuck et al. 2014; Jones et al. 2017; Shao et al. 2017; Rizzo et al. 2020; Fraternali et al. 2021), and even evidence for ordered rotation of galaxies at $z > 6$ (Smit et al. 2018; Bakx et al. 2020). However, the number of galaxies at $z > 4$ that are kinematically characterizable (i.e. observed with great enough S/N and resolution) is still much less than at $z < 4$. These high-redshift galaxies represent some of the first galaxies to form in the Universe and sit on the ramp-up of star formation rate density (e.g. Madau & Dickinson 2014), so their morphological and kinematic (i.e. morpho-kinematic) state is crucial for informing theories of early galaxy formation evolution, as well as constraining ongoing cosmological simulations (e.g. Pallottini et al. 2017; Kohandel et al. 2019).

Addressing this lack of kinematically characterized galaxies at $z > 4$ was one of the main driving goals of the ALMA Large Program to INvestigate C II at Early Times (ALPINE; Béthermin et al. 2020; Faisst et al. 2020; Le Fèvre et al. 2020). This survey observed 118 galaxies at $z = 4.4-5.9$ in [C II] 158 µm and the surrounding rest-frame FIR emission, following the success of the pilot program (Capak et al. 2015) and increasing the number of such observations by an order of magnitude. Source selection was based on UV luminosity ($L_{UV} > 0.6$ L*), pre-existing spectroscopic redshifts, and lack of type 1 active galactic nucleus (AGN). ALPINE sources lie on the star-forming main sequence (e.g. Noeske et al. 2007; Faisst et al. 2020) and are thus broadly representative of the underlying population of galaxies at these redshifts.

As the first large spectroscopic survey of $z > 4$ SFGs, ALPINE probes a unique era of cosmic history and provides the first measure of the morpho-kinematic diversity in the $z \sim 4-6$ epoch. Using a custom line search algorithm, Béthermin et al. (2020) found that 75 of the 118 targets were detected in [C II] emitting emission. The morpho-kinematic diversity of this sample was then assessed by a team within the ALPINE collaboration (Le Fèvre et al. 2020) who examined the [C II] channel maps, integrated intensity (moment 0) maps, velocity fields (moment 1), position-velocity diagrams (PVDs) along the major and minor axes, integrated spectra, and ancillary photometry (Faisst et al. 2020). Based on this information, each member independently classified each galaxy as rotating (class 1), merging (class 2), extended dispersion-dominated (class 3), compact dispersion-dominated (class 4), or too weak to characterize (class 5), and the class for each galaxy was agreed upon.

This initial qualitative classification was an invaluable first step towards characterization of the sample, and we may add further measures to contextualize and interpret these classes. This addition is crucial, as the [C II] observations of this sample are relatively limited. That is, the spatial resolution of this data set is low, with at best only $\sim 2-3$ beams across each source. In addition, while these sources are 'normal' for this redshift range (i.e. on the star-forming main sequence), the modest integration times (i.e. $< 1$ h) result in

limited sensitivity. These effects cause the true nature of the galaxy to be partially obfuscated by random noise peaks, and to be spatially smoothed.

However, by investigating the morpho-kinematics of each source using different methods, it is possible to refine their classification, adding evidence for their intrinsic behaviour. To this end, we conduct a quantitative analysis of the morpho-kinematics of the sample by fitting each galaxy with a tilted ring fitting code and applying well-tested local morpho-kinematic classification criteria.

While tilted ring model (TRM) fitting codes have mainly been implemented to derive the physical parameters of rotating disc galaxies at low and high redshift (e.g. Fan et al. 2019; Sharma et al. 2021), they may also be used to classify the morpho-kinematic class of galaxies. That is, if the galaxy is truly a rotating disc, then the models will be well fit and yield rotation curves, morphological parameters, and kinematic parameters. On the other hand, early stage mergers will be poorly fit (but see discussions of Simons et al. 2019). Dispersion-dominated galaxies will be well fit with a TRM, but since their intrinsic morphologies are not necessarily flat discs, the best-fitting physical parameters may not be physical. So by fitting a TRM to each of the [C II] data cubes in ALPINE, we may add additional evidence to each kinematic classification.

In this work, we examine the 75 [C II]-detected galaxies in the ALPINE survey using TRMs and morpho-kinematic classification criteria, with the goals of studying the kinematic diversity of this unique sample, testing the applicability of low-redshift classification criteria to high-redshift observations, and characterizing the properties (e.g. rotation curves, morphological parameters, velocity dispersion) of rotation-dominated ALPINE galaxies. We begin by describing the ALPINE survey and the creation of the [C II] data cubes in Section 2. The $^{3D}$BAROLO modelling procedure and results are then presented in Section 3. We test the applicability of two quantitative morpho-kinematic classifiers from low-redshift studies in Section 4. Following a discussion of these results in Section 5, we conclude in Section 6. We assume a standard concordance cosmology ($\Omega_{\Lambda}$, $\Omega_{m}$, $h$) = (0.7, 0.3, 0.7) throughout.

## 2 OBSERVATIONS AND DATA CUBE CREATION

The overview of the ALPINE survey is presented in Le Fèvre et al. (2020), while the observations, data reduction, and the creation of the public catalogue<fnref n="1" /> are detailed in Béthermin et al. (2020), and the multiwavelength photometry analysis of the ALPINE sources is described in Faisst et al. (2020). The [C II] 158 µm and the surrounding rest-frame FIR continuum emission of each source was observed with ALMA in cycles 5 and 6. Each observation was calibrated using the standard heuristic-based CASA (Common Astronomy Software Applications; McMullin et al. 2007) pipeline, and the results were inspected for quality. After some minor additional flagging, these data were used to create preliminary data cubes, which are now publicly available.

Each cube was continuum-subtracted in the uv-plane using the CASA task UVCONTSUB, resulting in a line-only, continuum-free data cube. Since the restoring beam of each observation was comparable (average of 1.13 arcsec × 0.85 arcsec; Béthermin et al. 2020), each image was created using a uniform cell size of 0.15 arcsec. The line cubes were constructed with channels of 25 km s⁻<fnref n="1" />, or $\sim 30$ MHz.

<footnote n="1">¹Catalogue and data products are available at https://cesam.lam.fr/a2c2s/data_release.php.</footnote>
All cubes were cleaned down to $3\sigma$ (CASA TCLEAN) and were created using natural weighting.

# 3 TILTED RING MODEL FITTING

While there is a multitude of codes that extract the physical parameters and kinematic details of a galaxy from observational data, they may be divided into two broad classes: codes that fit models to the velocity field (e.g. GIPSY ROTCUR, van der Hulst et al. 1992; NEMO ROTCURSHAPE, Teuben 1995; RESWRI, Schoenmakers 1999; RINGFIT, Simon et al. 2003; KINEMETRY, Krajnović et al. 2006; DISKFIT, Sellwood & Spekkens 2015) and codes that use the entire data cube (e.g. TiRiFiC, Józsa et al. 2007; DYSMAL, Cresci et al. 2009; KINMS, Davis et al. 2013; $^3$DBAROLO, Di Teodoro & Fraternali 2015; GALPAK 3D, Bouché et al. 2015; BLOBBY3D, Varidel et al. 2019).

We choose to apply the 3D TRM fitting code $^3$DBAROLO to the [C II] emission of the ALPINE sample. This code has already been well-tested on data cubes with relatively low resolution and sensitivity (Di Teodoro & Fraternali 2015) and has been successfully applied to high-redshift ($z > 2$) observations (e.g. Talia et al. 2018; Fan et al. 2019; Neeleman et al. 2020; Fraternali et al. 2021).

## 3.1 Signal isolation

To fit our data with a TRM, we must first mask the [C II] emission and exclude sources that are detected at too low resolution or sensitivity for kinematic modelling. First, the data cubes for each [C II]-detected ALPINE source were trimmed to only include the sideband (i.e. two spectral windows; SPWs) containing [C II] emission, and the rest frequency of each was set to that of [C II] at the redshift derived by Béthermin et al. (2020). In some cases, the spectral setup created a gap between the two SPWs, which we masked.

Next, the line emission in each continuum-subtracted cube was identified using the SEARCH algorithm in $^3$DBAROLO, which is based on the code DUCHAMP (Whiting 2012). After automatically determining the median noise level of each channel, this algorithm searches for pixels with intensities above a user-provided signal-to-noise ratio (SNR<sub>upper</sub> = 3.0). Here, SNR is defined as the value of a given pixel divided by the noise level of its channel. A three-dimensional (i.e. RA, Dec., velocity) search is then conducted around these peaks for emission above a second user-provided SNR (SNR<sub>lower</sub> = 2.5), creating a three-dimensional mask of the signal in the cube. Decreasing SNR<sub>lower</sub> (e.g. to 2.0) results in the inclusion of more pixels for every source, but in many cases also encloses noise peaks that complicate kinematic characterization. On the other hand, SNR<sub>upper</sub> is simply used to find the bright pixels at the centre of each mask, so slightly varying this value (e.g. from 2.5 to 3.0) has little effect for most of our [C II]-detected sources (>3.5σ; Béthermin et al. 2020).

To avoid a contribution from non-target sources (i.e. noise peaks and serendipitous sources distant from the phase centre; see Loiacono et al. 2021 for analysis of this latter type of objects), we remove all sources that are ≥5 beams from the centre from our signal mask. In addition, we exclude all unresolved sources by removing sources that feature a masked area (i.e. $2.5\sigma$ contour) smaller than the half-power area of the synthesized beam of the input [C II] data cube.

### 3.1.1 SEARCH results

Of the 118 galaxies observed in ALPINE, Béthermin et al. (2020) found [C II] emission (>3.5σ) in 75 at the expected frequency from the UV spectroscopy results (Faisst et al. 2020). By applying $^3$DBAROLO SEARCH to the continuum subtracted cubes of these 75 [C II]-detected ALPINE galaxies, we identify emission in 40 galaxies, given the constraints outlined above.

The disparity in the number of detected galaxies is mainly due to the fact that the line search algorithm used by Béthermin et al. (2020) is more sensitive to low-level, broad emission than SEARCH. That is, the prior algorithm is able to search the cube using kernels of different velocity widths, while SEARCH searches for pixels above a single value and then merges them. Since our morphokinematic analysis will examine each spectral channel of the data cube (see Section 3.2), we require significant (i.e. >2σ) emission in all spectral channels and thus must exclude broad, low-level emission.

An additional important difference is our criterion that the identified signal must have a $2.5\sigma$ contour that is larger than the half power contour of the restoring beam. This limit on compactness, which is not present in the Béthermin et al. (2020) search, is necessary to avoid fitting marginally resolved or unresolved galaxies. It excludes all compact dispersion-dominated sources (morphokinematic class 4 in Le Fèvre et al. 2020).

The five morphokinematic classes of Le Fèvre et al. (2020) contain 9, 31, 15, 8, and 12 galaxies (75 total) for classes 1–5, respectively. Of these sources, we recover 6, 21, 13, 0, and 0 galaxies, respectively (40 total). The total lack of compact or weak (class 4 or 5, respectively) galaxies is explained by the algorithm criteria noted above. However, we do recover ~67–87 per cent of the sources originally classified into classes 1–3.

## 3.2 3DFIT details

Once the [C II] signal from resolved sources has been isolated using $^3$DBAROLO SEARCH, the main function of $^3$DBAROLO (i.e. 3DFIT) is then used to fit a TRM to the line emission. In general, TRMs assume that a galaxy may be approximated as a set of concentric rings, each with its own geometric and kinematic parameters (e.g. Rogstad, Lockhart & Wright 1974). Three-dimensional model fitting considers some variables that are included in two-dimensional fits: inclination (i), position angle (PA), rotational velocity ($v_{rot}$), systemic redshift ($z$ or $v_{sys}$), and the spatial centroid ($x_0, y_0$). By expanding to the spectral dimension, we may also consider velocity dispersion ($\sigma_v$), radial brightness profile ($I(r)$), and scale height ($Z_o$). While this task is able to automatically predict many parameters and is well-tested for low-SNR and low-resolution observations, its performance is dependent on the initial parameter estimates and the overall model geometry, including the radius of the model and the width of each ring. Before we detail the fitting process, we will discuss the methods used to derive physically motivated estimates for each parameter.

### 3.2.1 Parameter estimation

To begin, each cube is collapsed over all channels containing line emission (as identified by $^3$DBAROLO SEARCH) using the CASA toolkit task IMAGE.MOMENTS. The resulting integrated intensity (i.e. moment 0) map is fit with a two-dimensional Gaussian using the CASA toolkit task IMAGE.FITCOMPONENTS, yielding a best-fitting peak intensity, integrated flux density, beam-deconvolved full width at half-maximum (FWHM) of the major and minor axes, position
angle, and central position. In addition, a velocity field (moment 1) and velocity dispersion (moment 2) map are created using the same task, including only pixels identified through the SEARCH algorithm.

To avoid under- or oversampling the data, the ring geometry is determined in an automatic fashion. If the emission in the moment zero map is sufficiently resolved to allow CASA to fit a 2D Gaussian and deconvolve the intrinsic emission from the synthesized beam, then we set the maximum model radius to $0.8 \times$ (the deconvolved FWHM of the major axis of the 2D Gaussian), based on the high-resolution kinematic analysis of Neeleman et al. (2020). If the source is not sufficiently resolved, we exclude it from further analysis.

The minimum radius of each ring is fixed to be the FWHM of the minor axis of the restoring beam divided by 2.5. This factor of 2.5 is comparable to other beam/ring width ratios (e.g. 2.2, Shao et al. 2017; 2.5, Talia et al. 2018; Fan et al. 2019) and is chosen as a middle-ground between the number of pixels per beam FWHM (~5) and the lower factors adopted by other studies of galaxy morphokinematics (e.g. ≤1, Mancera Piña et al. 2020; Salak et al. 2020; Shelest & Lelli 2020; ~2, de Blok et al. 2008). Note that while this high ratio may oversample the data in the case of 2D (i.e. velocity field) modelling, our use of 3D modelling reduces the effects of beam smearing, allowing us to sample spatial scales slightly smaller than the synthesized beam (e.g. Di Teodoro & Fraternali 2015).

The number of rings is then determined by dividing the maximum model radius by the minimum ring radius and rounding down to the closest integer. If this number is ≥2, then the width of each ring is set to the maximum model radius divided by the number of rings. We do not consider models containing only one ring.

While $^3\text{BAROLO}$ allows the user to fit the central position of each ring at the same time as the inclination and position angle, we fix the central position of all model rings to the morphological central position, as found through a 2D Gaussian fit to the moment zero map, in order to reduce the free parameters of the fit. Since we wish to test how well the data are fit by a rotating disc model, this is tantamount to assuming that the morphological and kinematic centres of the galaxy are coincident.

The initial guess for the velocity dispersion is based on the moment two map, and a thin disc is assumed ($Z_c = 0.01$ arcsec ~60 pc at $z \sim 5$). Each galaxy is given an initial systemic velocity estimate of $0 \text{ km s}^{-1}$ with respect to the assumed redshift. The inclination is allowed to vary between 10° and 80°, with an initial guess of 45°. $^3\text{BAROLO}$ provides an initial guess for the rotational velocity based on the data cube, while we estimate the position angle by examining the high-velocity channel maps.

### 3.2.2 Fitting procedure

Using these initial morphological and kinematic estimates and a list of ring radii and widths, $^3\text{BAROLO}$ first creates a model of the innermost ring by populating a physical volume with discrete clouds, such that all estimated parameters are recreated. This model is then converted into an observational cube (i.e. with axes of RA, Dec., velocity), convolved with the synthesized beam of the input data cube, normalized by setting the integrated flux per spaxel in the model and input data cubes equal, and compared to the input cube. Each parameter (i.e. rotational velocity, velocity dispersion, systemic velocity, inclination, and position angle) is varied, until the absolute residual (i.e. $|\text{model} - \text{observation}|$ over all pixels) is minimized. When the residual is minimized, the fitting stops, and the next ring is analysed. After all rings have been fit, the best-fitting spatial parameters ($i$, PA) and $v_{\text{sys}}$ across all rings are averaged and

fixed. The fitting is then repeated for the remaining parameters (i.e. rotational velocity and velocity dispersion).

The final outputs of this process are morphological parameters from a 2D Gaussian fit to the moment zero map (i.e. inclination, position angle, central position, radius), morpho-kinematic parameters from the $^3\text{BAROLO}$ fit (i.e. inclination, position angle, velocity dispersion profile, rotation curve, systemic velocity/redshift), and a best-fitting model data cube that may be directly compared to the data. For further discussion of this procedure, see Appendix A.

### 3.2.3 Tilted ring model results

When $^3\text{BAROLO}$ 3DFIT is run on the 40 recovered sources (see Section 3.1.1 for details of signal isolation), 15 per cent (six sources) are not successfully modelled due to being only marginally resolved (i.e. it was not possible to deconvolve the moment zero map from the synthesized beam), while another ~13 per cent (five sources) are resolved but are excluded due to only featuring one ring. For the 29 galaxies that are successfully fit, we present comparisons of moment maps, PVDs, and spectra for the data and models in Figs 1 and B1–B6, list the best-fitting parameters in Table 1, and show the best-fitting rotation curves and velocity dispersion profiles in Fig. 2. Each source is also discussed in Appendix B.

These figures show that the morpho-kinematic diversity noted by Le Fèvre et al. (2020) is indeed present. Some galaxies depict perfect velocity gradients and are well fitted with TRMs, while others show significant residuals, and others seem dispersion-dominated.

We apply $^3\text{BAROLO}$ to each source with the intent of examining the residuals in each moment map, PVD, and the integrated spectrum, in order to refine the morpho-kinematic classification for the 40 ALPINE galaxies whose emission was identified though $^3\text{BAROLO}$ SEARCH (see Section 3.1). That is, we do not pre-select rotators from the ALPINE sample before fitting the data with TRMs. Due to its optimization for low-SNR and low-resolution observations and large number of free parameters, it easily fits mergers and dispersion-dominated sources with rotating disc models. Because of this, the kinematic values presented in Table 1 and curves in Fig. 2 should be interpreted with care. Only the results for galaxies that are robustly classified as rotators (see Section 5.3) are fully reliable.

## 3.3 Subsample properties

The previous subsections have detailed the methods used to isolate [C II] line emission in the ALPINE data cubes and fit this signal with the TRM fitting code $^3\text{BAROLO}$. Before proceeding, it is worth briefly discussing the features of our 29 $^3\text{BAROLO}$-fit galaxies.

Our signal isolation and model fitting routine has successfully fit the majority of the bright [C II]-detected ALPINE sources. Of the 21 sources detected at $\text{SNR}_{[C II]} \ge 10$ by Béthermin et al. (2020), we only exclude four sources due to weak (i.e. spectrally broad) and/or marginally resolved [C II] emission. While this does result in the exclusion of some strong but compact [C II] emitters, such as DC488399 ($\text{SNR}_{[C II]} = 26.2$) and DC630954 ($\text{SNR}_{[C II]} = 11.2$, Béthermin et al. 2020), we are exploring the three-dimensional (i.e. RA, Dec., velocity) distribution of [C II] emission in these galaxies, and so resolved data are required.

To further examine the subset of ALPINE sources studied in this work, we compare their location on the SFR-$M_*$ plane (see top panel of Fig. 3, blue squares) to the ALPINE sources excluded by this analysis but detected by the line search algorithm of Béthermin
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Moment maps, PVDs, and spectra for observed data, model, and residual for an example rotator (upper left), merger (upper right), dispersion dominated (lower left), and uncertain class (lower right) galaxy (see Section 5.1 for details of each kinematic class). In each of these figures, the first three rows show (from top to bottom) the moment 0 (integrated intensity), moment 1 (velocity field), and moment 2 (velocity dispersion field; see Section 3.2.1 for details of moment map creation). For these rows, the three columns denote (from left to right) the observed data cubes, model cubes, and the corresponding residuals. The white crosses in the lower right corner of each panel in the first row show a 5 kpc × 5 kpc physical scale. Solid lines in the second row represent the kinematic major axis, while the dashed lines represent the minor axis. The restoring beam is shown as a red ellipse. The bottom row shows (from left to right) the major axis PVD, minor axis PVD, and integrated spectra. For each PVD, the observed data are shown by the background colour, while the contours represent the model at 20 per cent, 50 per cent, and 80 per cent of its maximum value. The data (D), model (M), and residual (R) spectra are depicted by the green, purple, and orange lines, respectively. The data spectrum is extracted from the continuum-subtracted [CII] cube over all spaxels where line emission was identified by 3D BAROLO SEARCH. The $1\sigma$ uncertainty, calculated as (average rms noise level per channel) × √number of beams in source, is shown by the shaded grey area. See Appendix B for plots of other sources.</figcaption></figure>
<table><caption>**Table 1.** Morphological and kinematic parameter values for ALPINE galaxies successfully fit by $^3D_{\text{BAROLO}}$, as well as the two previously detected $4 < z < 6$ rotating un lensed SFGs (see Section 5.3). Reported uncertainties on redshift are the greater of the fit uncertainty and the width of a velocity channel. The morphological position angles and inclinations are taken from two-dimensional Gaussian fits to moment zero maps created from SEARCH-identified [C II] emission, while the kinematic values are the best-fitting values from $^3D_{\text{BAROLO}}$. The W15 criteria correspond to the five tests of disc-like behaviour suggested by Wisnioski et al. (2015), as explored in Section 4.2. We list both the morpho-kinematic class derived by Le Fèvre et al. (2020) (L20) and the class derived by combining the results of our analyses (J21, see Section 5.1).</caption><thead><tr><th rowspan="2">Name</th><th rowspan="2">z</th><th rowspan="2">PA<sub>M</sub><br>(°)</th><th rowspan="2">PA<sub>K</sub><br>(°)</th><th rowspan="2">i<sub>M</sub><br>(°)</th><th rowspan="2">i<sub>K</sub><br>(°)</th><th colspan="5">W15</th><th rowspan="2">KC1<br>L20</th><th rowspan="2">KC2<br>J21</th></tr><tr><th>1</th><th>2</th><th>3</th><th>4</th><th>5</th></tr></thead><tbody><tr><td>CG32</td><td>4.4105 ± 0.0005</td><td>90 ± 60</td><td>238 ± 14</td><td>50 ± 20</td><td>55 ± 19</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>2</td><td>ROT</td></tr><tr><td>DC396844</td><td>4.5424 ± 0.0005</td><td>161 ± 28</td><td>197 ± 14</td><td>60 ± 19</td><td>57 ± 31</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>1</td><td>ROT</td></tr><tr><td>DC417567</td><td>5.6700 ± 0.0006</td><td>160 ± 80</td><td>49 ± 4</td><td>40 ± 50</td><td>68 ± 13</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>2</td><td>UNC</td></tr><tr><td>DC432340</td><td>4.4045 ± 0.0005</td><td>11 ± 18</td><td>187 ± 8</td><td>59 ± 19</td><td>68 ± 12</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>2</td><td>UNC</td></tr><tr><td>DC434239</td><td>4.4883 ± 0.0005</td><td>28 ± 8</td><td>217 ± 3</td><td>63 ± 8</td><td>73 ± 13</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>2</td><td>MER</td></tr><tr><td>DC454608</td><td>4.5835 ± 0.0005</td><td>140 ± 40</td><td>193 ± 9</td><td>67 ± 25</td><td>57 ± 8</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>2</td><td>UNC</td></tr><tr><td>DC494057</td><td>5.5446 ± 0.0005</td><td>60 ± 90</td><td>185 ± 13</td><td>31 ± 25</td><td>46 ± 4</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>1</td><td>ROT</td></tr><tr><td>DC519281</td><td>5.58 ± 0.02</td><td>100 ± 80</td><td>70 ± 10</td><td>60 ± 30</td><td>53 ± 9</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>2</td><td>UNC</td></tr><tr><td>DC552206</td><td>5.5016 ± 0.0005</td><td>138 ± 23</td><td>322 ± 7</td><td>52 ± 13</td><td>61 ± 17</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>2</td><td>ROT</td></tr><tr><td>DC627939</td><td>4.5333 ± 0.0006</td><td>170 ± 90</td><td>278 ± 4</td><td>50 ± 30</td><td>66 ± 8</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>2</td><td>UNC</td></tr><tr><td>DC683613</td><td>5.5421 ± 0.0005</td><td>170 ± 40</td><td>353 ± 2</td><td>58 ± 30</td><td>49 ± 19</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>3</td><td>UNC</td></tr><tr><td>DC733857</td><td>4.5445 ± 0.0005</td><td>87 ± 27</td><td>84 ± 16</td><td>84 ± 19</td><td>74 ± 1</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>3</td><td>UNC</td></tr><tr><td>DC773957</td><td>5.6773 ± 0.0023</td><td>101 ± 11</td><td>35 ± 5</td><td>76 ± 11</td><td>66 ± 19</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>2</td><td>UNC</td></tr><tr><td>DC818760</td><td>4.5611 ± 0.0005</td><td>94 ± 2</td><td>90 ± 5</td><td>74 ± 2</td><td>62 ± 11</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>2</td><td>MER</td></tr><tr><td>DC848185</td><td>5.2930 ± 0.0005</td><td>149 ± 11</td><td>320 ± 9</td><td>53 ± 9</td><td>52 ± 0</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>3</td><td>DIS</td></tr><tr><td>DC873321</td><td>5.1544 ± 0.0005</td><td>118 ± 5</td><td>281 ± 10</td><td>73 ± 6</td><td>65 ± 11</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>2</td><td>MER</td></tr><tr><td>DC873756</td><td>4.5457 ± 0.0005</td><td>121 ± 19</td><td>298 ± 0</td><td>42 ± 10</td><td>39 ± 28</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>2</td><td>DIS</td></tr><tr><td>DC881725</td><td>4.5778 ± 0.0005</td><td>140 ± 40</td><td>323 ± 9</td><td>76 ± 33</td><td>48 ± 11</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>1</td><td>ROT</td></tr><tr><td>VC5100537582</td><td>4.5502 ± 0.0005</td><td>60 ± 90</td><td>78 ± 6</td><td>40 ± 50</td><td>49 ± 3</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>3</td><td>UNC</td></tr><tr><td>VC5100541407</td><td>4.5632 ± 0.0005</td><td>69 ± 15</td><td>69 ± 7</td><td>59 ± 11</td><td>66 ± 9</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>2</td><td>UNC</td></tr><tr><td>VC5100559223</td><td>4.5626 ± 0.0005</td><td>167 ± 28</td><td>12 ± 6</td><td>75 ± 26</td><td>49 ± 4</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>3</td><td>UNC</td></tr><tr><td>VC5100822662</td><td>4.5205 ± 0.0005</td><td>0 ± 80</td><td>11 ± 15</td><td>40 ± 50</td><td>53 ± 7</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>2</td><td>MER</td></tr><tr><td>VC5100994794</td><td>4.5800 ± 0.0005</td><td>55 ± 27</td><td>243 ± 9</td><td>69 ± 16</td><td>49 ± 17</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>3</td><td>UNC</td></tr><tr><td>VC5101209780</td><td>4.5701 ± 0.0016</td><td>23 ± 16</td><td>31 ± 13</td><td>57 ± 12</td><td>72 ± 8</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>2</td><td>MER</td></tr><tr><td>VC5101218326</td><td>4.5741 ± 0.0006</td><td>51 ± 25</td><td>9 ± 2</td><td>54 ± 27</td><td>44 ± 4</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>3</td><td>DIS</td></tr><tr><td>VC510786441</td><td>4.463 ± 0.004</td><td>4 ± 6</td><td>1 ± 2</td><td>65 ± 7</td><td>66 ± 10</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>2</td><td>UNC</td></tr><tr><td>VC5110377875</td><td>4.5506 ± 0.0005</td><td>146 ± 18</td><td>147 ± 7</td><td>56 ± 15</td><td>50 ± 6</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>1</td><td>ROT</td></tr><tr><td>VC5180966608</td><td>4.5296 ± 0.0005</td><td>125 ± 34</td><td>311 ± 9</td><td>45 ± 22</td><td>52 ± 5</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>2</td><td>UNC</td></tr><tr><td>VE530029038</td><td>4.4297 ± 0.0005</td><td>130 ± 60</td><td>306 ± 7</td><td>50 ± 50</td><td>47 ± 15</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>1</td><td>UNC</td></tr><tr><td>J0817, low-res</td><td>4.2605 ± 0.0009</td><td>65 ± 24</td><td>93 ± 24</td><td>58 ± 23</td><td>48 ± 5</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>—</td><td>ROT</td></tr><tr><td>J0817, high-res</td><td>4.2600 ± 0.0005</td><td>109 ± 9</td><td>108 ± 7</td><td>46 ± 6</td><td>64 ± 7</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>—</td><td>ROT</td></tr><tr><td>HZ9</td><td>5.5413 ± 0.0004</td><td>83 ± 22</td><td>17 ± 10</td><td>57 ± 14</td><td>58 ± 7</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>—</td><td>ROT</td></tr></tbody></table>

et al. (2020; green diamonds), and the [C II]-undetected sources (red circles). Note that the ALPINE sample was chosen to contain SFGs, so the vast majority of these sources are within the $1\sigma$ scatter of the main-sequence relation. This plot shows that the subset fit here contains massive systems (i.e. $M_* \gtrsim 10^{9.5} M_\odot$) that exhibit a wide range of star formation rate (SFR), including a few sources above and below the star-forming main sequence.

We may also examine the redshift distribution of the $^3D_{\text{BAROLO}}$-fit galaxies compared to the 75 [C II]-detected ALPINE sources (Béthermin et al. 2020; see bottom panel of Fig. 3). While each redshift bin contains $\gtrsim 2 \times$ as many [C II]-detected sources as $^3D_{\text{BAROLO}}$-fit sources, we can see that the $^3D_{\text{BAROLO}}$-fit galaxies cover nearly the same redshift range as the full sample. Thus, our sample of 29 $^3D_{\text{BAROLO}}$-fit ALPINE galaxies contains massive SFGs at a range of SFR/$M_*$ and redshifts.

## 4 ADDITIONAL CLASSIFICATION METHODS

Thus far, we have presented the best-fitting TRMs for a subset of the ALPINE sample. To further examine the morpho-kinematics of these sources, we apply two objective, quantitative classification

methods tested at $z < 4$ to the sample: the Gini-$M_{20}$ method of Lotz, Primack & Madau (2004) and the five disc-like criteria of Wisnioski et al. (2015). Each of these methods has not been applied to rest-frame FIR and [C II]158 µm observations of $z > 4$ SFGs, so this section represents an investigation of their applicability to our data.

### 4.1 Gini-$M_{20}$ analysis

First, we test whether the ALPINE sample may be classified based solely on its observed morphology by examining the distribution of our sources in the 'G-$M_{20}$' plot (e.g. Lotz et al. 2004). This approach has been used to distinguish different morphological types in the local Universe.

The Gini coefficient (or G) is a measure of how uniformly the brightness is distributed. This value approaches 0 when all pixels have the same value (e.g. diffuse emission), and approaches 1 when all brightness is focused in a few pixels (e.g. central nuclei). On the other hand, the moment of light ($M_{20}$) is the normalized second-order moment of the pixels that make up the brightest 20 per cent of the galaxy, and is high when non-axisymmetric features (e.g.
<figure><img src="image_3.png" /><figcaption>**Figure 2.** The best-fitting rotation velocity (blue) and velocity dispersion (orange) curves for each galaxy. Vertical extent of each point shows the fit uncertainty output by $^{3D_{BAROLO}}$ (i.e. not including systematic uncertainty on the inclination), while the horizontal extent shows the width of each ring. The letters after each galaxy name denote the morpho-kinematic classification we assign in Section 5.1: ROT (R), MER (M), DIS (D), or UNC (U).</figcaption></figure>

bars, clumps, spiral arms) are present. If both of these measures
are calculated for each galaxy, then it is possible that mergers and
single sources will occupy separate parameter spaces (Lotz et al.
2008).

4.1.1 *Implementation*

This analysis is applicable to two-dimensional images (i.e. not
data cubes), so we apply it to both rest-frame FIR continuum (i.e.
~160 µm) maps and [C II] moment zero maps from the ALPINE
observations. Both sets of maps are taken from the ALPINE Data
Release 1 (DR1), and the details of their creation are available
<figure><img src="image_4.png" /><figcaption>Figure 3. Properties of the 29 $^{3D}$BAROLO-fit galaxies, compared to the 75 [C II]-detected galaxies and full sample. Top: SFR and $M_*$ for all 118 ALPINE galaxies, based on modelling of ancillary photometry with LePhare (Faisst et al. 2020). Sources undetected (i.e. <3.5σ) by the line search algorithm of Béthermin et al. (2020) are shown as red circles. The sources identified by $^{3D}$BAROLO SEARCH are shown by blue squares, while the Béthermin et al. (2020) detected sources that are not identified by SEARCH are shown by green diamonds. The star-forming main-sequence relation of Speagle et al. (2014) for a representative ALPINE redshift ($z = 5.15$) is shown by the solid black line, while the grey shaded region shows the scatter induced by varying each variable in the Speagle relation by $1\sigma$, as well as varying $z$ over the redshift range of the ALPINE sample ($z = 4.4–5.9$). Bottom: Redshift distribution for the 29 $^{3D}$BAROLO-fit sources (blue) compared to that of the parent sample of 75 sources detected in [C II] emission by Béthermin et al. (2020; green).</figcaption></figure>

in Béthermin et al. (2020). This set of maps contains 75 [C II] moment zero maps for the [C II]-detected set of sources and 23 continuum images for the FIR-detected set. Since our goal is to compare the G and $M_{20}$ distribution of these maps to their [C II] morpho-kinematic diversity, we exclude two sources that are detected in FIR emission but not [C II] (CANDELS_GOODSS_19 and DEIMOS_COSMOS_460378), bringing the number of analysed continuum maps to 21.

To isolate the signal in each map, we follow a procedure similar to that of $^{3D}$BAROLO SEARCH. The noise level of each map is determined by sigma clipping above and below $3\sigma$ and taking the standard deviation of the resulting pixels (see ASTROPY task sigma_clipped_stats). All pixels above $3 \times$ this standard deviation are identified. All such peaks outside of 20 pixels (i.e. ~4 beams) from the centre are rejected as noise or serendipitous sources. A mask is then created by extrapolating out from the central peaks to their surrounding $1\sigma$ contours.<fnref n="2" /> By applying this mask to the original map, we are able to isolate the pixels that are associated with the target source. We choose to mask the data in this way, rather than analysing the full image or using an aperture, to ensure that each masked map contains only signal.

To calculate G, it is first necessary to sort the value of each pixel in the map ($f_i$) into increasing order (i.e. $f_1$ for the smallest value, $f_n$ for the largest). G may then be calculated as

$$G = \frac{\sum_{i=1}^{n} (2i - n - 1) f_i}{\bar{f} n(n - 1)}, \quad (1)$$

where n is the number of pixels and $\bar{f}$ is the mean pixel value.

$M_{20}$ requires the values of all pixels sorted in decreasing order (i.e. $f_1$ for the largest value, $f_n$ for the smallest). We then calculate the second spatial moment of each pixel ($M_i$) and their sum ($M_{tot}$):

$$M_i = f_i \left[ (x_i - x_c)^2 + (y_i - y_c)^2 \right] \quad (2)$$

$$M_{tot} = \sum_{i=1}^{n} M_i, \quad (3)$$

where $(x_c, y_c)$ is the source centre computed by minimizing $M_{tot}$. Next, we find the index $i_{20}$ such that $\sum_{i=1}^{i_{20}} f_i = 0.2 \sum_{i=1}^{n} f_i = 0.2 f_{tot}$, and find $M_{20}$ using

$$M_{20} = \log_{10} \left( \frac{\sum_{i=1}^{i_{20}} M_i}{M_{tot}} \right). \quad (4)$$

Our approach does not allow for the calculation of uncertainties on G or $M_{20}$, but we do note that these values are affected by imaging parameters (e.g. image weighting, cell size, cleaning threshold), as well as the choice of 3D mask.

### 4.1.2 Results

The bottom panel of Fig. 4 shows that the [C II] moment zero maps feature a distinct distribution compared to the continuum maps. All of the ‘compact dispersion-dominated’ and ‘weak’ objects lie at low G (i.e. ~0.2–0.3). While the class 1–3 objects spread to higher G (i.e. ~0.2–0.55), there again is no separation between the kinematic classes.

<footnote n="2">²A higher threshold was tested ($2\sigma$) which shifted each source to a higher G value but was not used because it excluded low-level emission.</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Gini-$M_{20}$ for the continuum (top) and [CII] moment zero (bottom) maps. Points are coloured by morpho-kinematic classification of Le Fèvre et al. (2020).</figcaption></figure>

Since this is the first application of this technique to ALMA FIR continuum and [CII] emission of $z > 4$ SFGs, there is no comparison sample to which we may compare. One successful application of this procedure was that of Lotz et al. (2008), which analysed the stellar distribution of lower redshift ($0.2 < z < 1.2$) galaxies as observed with HST and found that mergers featured high $G$ values (i.e. clumpy emission), E/S0/Sa had low $M_{20}$ values (i.e. symmetric morphology), and Sb-Irr galaxies featured low $G$ values (i.e. significant diffuse emission). This analysis of stellar morphology is not directly applicable to FIR and [CII] emission at $4 \lesssim z \lesssim 6$, since the latter trace warm dust, star formation regions, and molecular gas.

This lack of separation between morphological types of galaxies was also seen in a sample of 494 SFGs at $2.5 \lesssim z \lesssim 3.0$, as observed with HST H160 (Talia et al. 2014). Indeed, this HST study found that ellipticals, compact, disc-like, and irregular galaxies were focused between $-2 < M_{20} < -1$, with significant overlap between ellipticals, compact galaxies, disc-like galaxies, and irregulars.

<figure><img src="image_6.png" /><figcaption>**Figure 5.** G-$M_{20}$ for a set of model point sources convolved with a 2D Gaussian beam. For each of the 118 ALPINE synthesized beams, we explore five upper SNR thresholds ($3, 5, 10, 30, 50$) for $\text{SNR}_{\text{lower}} = 1.0$.</figcaption></figure>

Since the galaxies in the ALPINE sample are marginally resolved (i.e. <5 beams across the source) and feature a variety of SNR, it is conceivable that the observed $G-M_{20}$ scatter could be affected by observational properties. To determine the effects of the degree of resolution on the recovered $G$ and $M_{20}$, we assume the extreme case of each source being completely unresolved, and convolve a point source with the set of ALPINE synthesized beams (Béthermin et al. 2020). Using these point spread functions, we explore the effects of a range of SNR values on the recovered quantities. In practice, we take each of the 118 restoring beam axial ratios and explore five upper SNR values ($3, 5, 10, 30, 50$), assuming that the lower SNR threshold is 1. We construct a spatial model grid identical to that of the data (i.e. 0.15 arcsec cells) and calculate $G$ and $M_{20}$ for each $\text{SNR}_{\text{upper}}$.

As seen in Fig. 5, the 118 beams of ALPINE are similar, and thus generate only a small scatter in the plane for a given $\text{SNR}_{\text{upper}}$. The SNR threshold generates a significant effect, as stronger point sources result in higher $G$ values. The ALPINE continuum and [CII] moment zero maps feature peak SNRs of 3.6–22.3 and 3.5–32.6, respectively. This suggests that the scatter in $G$ values of the ALPINE sample may be explained by a point source convolved with a Gaussian. On the other hand, the point source models in Fig. 5 occupy a narrow range of $M_{20}$, so the observed range of ALPINE $M_{20}$ sources is likely dominated by intrinsic source morphology. This simple test demonstrates that the observed data are not sufficient for classification on morphology alone. However, this is among the first applications of the Gini-$M_{20}$ morphological classification technique to a statistical sample of FIR observations of high-redshift ($z > 4$) SFGs and will be useful for future investigations of galaxy morphology in the early Universe.

Future high-resolution spectral and continuum observations will likely enable this classification to be performed unambiguously. These include lines observable with ALMA (e.g. [CII] 158 µm, [OIII] 88 µm) or future observatories such as JWST (e.g. H$\alpha$, H$\beta$, [OIII] 5007 Å, [OII] 3727 Å), as well as the underlying rest-frame FIR and IR continuum emission.
## 4.2 Disc identification

One of the strengths of this survey is the three-dimensional nature of our data (i.e. RA, Dec., velocity). This allows us to move beyond the two-dimensional morphological classification described above to examine the morpho-kinematics of each [C II] data cube. This was done in Section 3 using $^3\text{BAROLO}$, and here we add modified versions of the five disc identification criteria of Wisnioski et al. (2015):

(i) If a slice of the velocity field is taken along the kinematic major axis, does it show a slope that is significant (i.e. $>3\sigma$)?

(ii) Is the average rotational velocity greater than the average velocity dispersion, across all rings?

(iii) Is the average location of the two extreme values of the velocity field (moment 1 map) within one beamwidth of the peak of velocity dispersion (moment 2)?

(iv) Is the difference between the morphological and kinematic position angles less than 30°?

(v) Is the average location of the two extreme values of the velocity field (moment 1 map) within one beamwidth of the peak of intensity (moment 0)?

If some or all of these criteria are met for a given source, then it is likely a rotating disc galaxy. Conversely, if none are met, the source is likely disturbed or a merger. Each criterion may be assessed using the three moment maps (i.e. integrated intensity, line of sight velocity, velocity dispersion) and the results of a $^3\text{BAROLO}$ TRM fit. The results for each galaxy are presented in Table 1.

The first criterion, which originally was the presence of a continuous velocity gradient, was made quantitative by requiring a significant slope in position-velocity space along the kinematic major axis. This replicates the standard 'rotating disc' criterion used by many high-redshift investigations (e.g. Smit et al. 2018; Bakx et al. 2020), and agrees with a by-eye inspection that some sources show obvious gradients (e.g. VE.9038) and others do not (e.g. DC519281).

The second criterion requires the best-fitting model to be strongly rotation-dominated. Perhaps surprisingly, this is met by almost all sources.

Both the third and fifth criteria relate to how well ordered the galaxy is. For an ideal disc, the extreme values of velocity would be at the extreme edges of the major axis, while both the integrated intensity and velocity dispersion would peak at the galaxy centre. On the other hand, the intensity peak in a merging system may be located in one of the component galaxies, resulting in a failure of criterion three. Similarly, the velocity dispersion map may be more complex, due to tidal features and shocks introduced by ongoing merging activity, resulting in a failure of criterion five.

The fourth criterion is a measure of how well the morphology and rotation pattern agree. Because the morphological position angle ($\text{PA}_M$) is derived by fitting a 2D Gaussian to the [C II] moment zero map, a marginally resolved source will return a $\pm\text{PA}_M$ that is strongly influenced by the synthesized beam. The most illustrative case is DC396844, where the restoring beam is perpendicular to the kinematic axis, resulting in a failure of this criterion. On the other hand, CG32 is barely resolved (as seen by the best-fitting $\text{PA}_M = 90 \pm 60^\circ$), but the kinematic angle happens to agree ($\text{PA}_K = 238 \pm 14^\circ$), so this criterion returns a positive value.

Only two sources in our subset meet all five criteria: DC432340 (classified as 'uncertain') and DC818760 (classified as 'merger', see Section 5.1). Since neither of these are robustly determined to be rotating disc galaxies, we find that the resolution and sensitivity of our data are insufficient for this specific analysis alone.

While the W15 criteria only highlight one measure of the data each (e.g. agreement between moment maps, the presence of a velocity gradient), they are of use when combined with the detailed best-fitting source characteristics from $^3\text{BAROLO}$ (Figs 1 and B1–B6), as discussed in the next section.

# 5 DISCUSSION

## 5.1 Synthesis of ALPINE morpho-kinematic classification

By combining our new $^3\text{BAROLO}$ fits and Wisnioski et al. (2015) disc criteria, we are now able to describe the morphology and kinematics of a massive subset of the galaxies in the ALPINE sample in unprecedented detail. Individual descriptions of each galaxy fitted with $^3\text{BAROLO}$ are included in Appendix B.

We classify each source into one of three robust classes or a fourth 'uncertain' classification:

(i) ROT: Rotators

(ii) MER: Mergers

(iii) DIS: Dispersion-dominated

(iv) UNC: Uncertain

The ROT and MER classes are identical to classes 1 and 2 of Le Fèvre et al. (2020), respectively, while DIS is essentially a combination of their classes 3 and 4 (i.e. ‘extended dispersion-dominated’ and ‘compact dispersion-dominated’, respectively).<fnref n="3" /> The ‘UNC’ class contains sources that we are unable to classify due to low SNR, low spectral resolution, and/or conflicting evidence from our combination of analyses. Note that the UNC class is quite different from class 5 (i.e. ‘weak’) of Le Fèvre et al. (2020), as it also contains strongly detected lines that we are unable to classify into ROT, MER, or DIS due to complex morpho-kinematics.

When classifying galaxies, we choose to be conservative. That is, a galaxy must clearly show multiple spatial and/or spectral components to be a merger, while it must exhibit a well-ordered velocity gradient and a single resolved source to be a rotator. Similarly, dispersion-dominated sources must show identical PVDs along orthogonal slices and be well-fit by a dispersion-dominated $^3\text{BAROLO}$ model. If a source does not fit into these classes for any reason, it is classified as ‘UNC’, regardless of its strength.

Of course, this approach is not perfect, as it is still possible that close-separation mergers may be misidentified as rotators (Simons et al. 2019), and face-on rotators may be interpreted as dispersion-dominated galaxies (e.g. Kohandel et al. 2019). More importantly, these classes are a vast oversimplification of the true morpho-kinematics of each source, as each galaxy may contain ordered rotation, minor and major mergers, star formation and AGN-driven outflows, filamentary accretion, and a number of other evolutionary processes that affect the kinematics, gas content, stellar mass, metallicity, and other properties of the galaxy and its environment (e.g. Law et al. 2009; Wisnioski et al. 2011, see also discussion in Valore, Dainotti & Kopczyński 2020). Indeed, cosmological simulations show that high-redshift galaxies are expected to interact with their environments through merging and accretion (e.g. Pallottini et al. 2017; Kohandel et al. 2019; Graziani et al. 2020; Zanella et al. 2021). So while these classes are worthwhile indicators of the general morpho-kinematic state of a galaxy, they may always be improved with higher resolution, higher sensitivity observations.

<footnote n="3">³Since our signal isolation procedure (see Section 3.1) has excluded all compact sources, no class 4 objects are included in our analysis.</footnote>
Of the 29 galaxies that were fit with ≥2 rings by $^{3D}$BAROLO, 6, 5, 3, and 15 were found to be ROT, MER, DIS, and UNC, respectively. As discussed in Section 3.3, these 29 galaxies constitute a representative subset (i.e. in redshift and SFR) of the high-$M_*$ portion of the ALPINE sample. The significance of this diversity is explored in the next section.

## 5.2 Morpho-kinematic diversity comparison

To examine the evolution of morpho-kinematic diversity with redshift, we may compare the fractions of galaxies in various kinematic classes in the ALPINE sample found in this work to those of lower redshift samples. In addition, we will place the results of Le Fèvre et al. (2020) in this framework for the first time.

In the relatively low-redshift Universe (i.e. $z = 0.4–0.75$), the IMAGES survey observed 68 intermediate-mass galaxies with VLT/FLAMES-GIRAFFE, resulting in $[\text{OII}]\lambda\lambda3726, 3729~\AA$ data cubes for each (Yang et al. 2008). Using the resulting moment 1 and 2 maps, they separate 39 galaxies into 9 (23 per cent) rotating discs, 11 (28 per cent) perturbed rotators, 16 (41 per cent) having complex kinematics, and 3 (8 per cent) unclassifiable. Note that the ‘perturbed rotator’ class contains sources that show disc-like rotation but have a velocity dispersion peak that is not coincident with the galactic centre. While this classification technique contains no possibility of dispersion-dominated sources or mergers, they may be hidden within the ‘complex kinematics’ class, due to the low spatial resolution and field of view of the utilized instrument.

Using VLT/SINFONI, Epinat et al. (2012) examined the $H\alpha$ velocity fields and I-band images of 45 galaxies at $z = 0.9–1.6$ from MASSIV (Mass Assembly Survey with SINFONI in VVDS; Contini et al. 2012). This sample is chosen to be representative of SFR $>5~M_{\odot}~\text{yr}^{-1}$ galaxies at $z \sim 1.3$ with $\log(M_*/M_{\odot}) \sim 9–11$. Using a blind group classification method, they find that 19 (42 per cent) are rotating discs, 16 (36 per cent) are non-rotating discs, and 10 (22 per cent) are not robustly classified. Separately, they carried out a search for companions, finding that of the 41 galaxies where this search was able to be performed robustly, 28 (68 per cent) were isolated and 13 (32 per cent) were interacting or merging.

A tremendous amount of work has been done to characterize the morpho-kinematics of galaxies around cosmic noon (i.e. $z \sim 2$). By taking deep NIR adaptive optics-assisted (AO) observations from the Spectroscopic Imaging survey in the Near-infrared with SINFONI (SINS; Förster Schreiber et al. 2009) and zCOSMOS (Lilly et al. 2007, 2009) samples, Förster Schreiber et al. (2018) classified a representative subset of 35 sources from $z = 1.35$ to 2.38 (i.e. SINS/zC-SINF). This sample contains 70 per cent rotation-dominated discs and 30 per cent other classes. This value agrees with the results of KMOS$^{3D}$ (Wisnioski et al. 2019), who found that 77 per cent of their $0.6 < z < 2.7$ galaxies were rotation-dominated.

At higher redshift ($z \sim 2.6–3.8$), Gnerucci et al. (2011) observed a sample of galaxies from the AMAZE (Maiolino et al. 2008) and LSD (Mannucci et al. 2009) programs with VLT/SINFONI, resulting in a set of $[\text{OII}]\lambda\lambda7007, 4959~\AA$ and $H\alpha$ data cubes. By fitting each of their 32 velocity fields with a simple plane model, they find that 11 (34 per cent) are rotators, 14 (44 per cent) are non-rotators, and 7 (22 per cent) are unclassifiable.<fnref n="4" /> This percentage of rotators agrees well with that of the KMOS Deep Survey (KDS;

Turner et al. 2017), which examined 77 $z \sim 3.5$ SFGs and found that $34 \pm 8$ per cent were rotation dominated.

Le Fèvre et al. (2020) used a similar blind group classification approach as Epinat et al. (2012): a group of collaborators examined moment maps, multiwavelength photometry, spectra, and channel maps of the 75 ALPINE galaxies detected in $[C\text{II}]$ at $>3.5\sigma$ by Béthermin et al. (2020). They find 9 (12 per cent) rotators, 31 (41 per cent) mergers, 15 (20 per cent) extended dispersion-dominated galaxies, 8 (11 per cent) compact dispersion-dominated galaxies, and 12 (16 per cent) galaxies that were too weak to classify. There are two notes of interest here: the fraction of sources classified as mergers is higher than at lower redshifts (49 per cent, excluding the unclassifiable galaxies), and the ratio of rotators to dispersion-dominated systems is low (0.39). When compared to the $z \sim 1$ kinematic diversity of Epinat et al. (2012), this suggests that mergers were more common for SFGs at high redshift. A more in-depth study of these merging systems in ALPINE is presented in Romano et al. (2021).

By combining our analyses (see Section 5.1), we are able to robustly classify 14 high-$M_*$ ALPINE galaxies, finding 6 (43 per cent) rotators, 5 (36 per cent) mergers, and 3 (21 per cent) dispersion-dominated sources. Surprisingly, this massive subset of the ALPINE sample contains a higher number of rotators than mergers. The differences between these fractions and those of Le Fèvre et al. (2020) are likely primarily due to two reasons: our exclusion of unresolved sources, and our stricter criteria for classification (see Section 3). The latter results in different classifications for some galaxies compared to Le Fèvre et al. (2020). As an example, Le Fèvre et al. (2020) classified DC432340 as a merger (class 2) due to the presence of a northern extension and an asymmetric PVD. In this work, we classify it as UNC, as it is too asymmetric to be a ROT or DIS, and only shows a single source, excluding it from MER.

In general, we find that for the ALPINE sample, both the full set of $[C\text{II}]$-detected galaxies and the massive subset studied in this work feature significant morpho-kinematic diversity, containing rotators, mergers, and dispersion-dominated systems. This first statistical test of the diversity of star-forming main-sequence galaxies only 1–1.5 Gyr after the big bang suggests a higher number of mergers compared to lower redshift sources, and a number of massive rotators. This agrees with cosmological simulations, which show that while high-redshift galaxies are expected to interact with their environments through merging and accretion, the central galaxy may still feature strong rotation (e.g. Kohandel et al. 2019; Kretschmer, Dekel & Teyssier 2021).

## 5.3 Main-sequence rotators at $z > 4$

A vast sample of rotating galaxies has been well studied at $z < 4$ (e.g. Epinat, Amram & Marcelin 2008; Förster Schreiber et al. 2009; Gnerucci et al. 2011), revealing dynamical masses and empirical scaling relations. While a number of clumpy or merging galaxies have been detected at $z > 4$ (e.g. Carniani et al. 2018; Díaz-Santos et al. 2018; Pavesi et al. 2018), and a handful of rotating starburst or quasar host galaxies at $z \sim 4–6$ have been well-modelled (e.g. Jones et al. 2017; Pensabene et al. 2020; Tadaki et al. 2020; Fraternali et al. 2021; Lelli et al. 2021), only two unlensed rotating main-sequence galaxies have been observed at $z \sim 4–6$: HZ9 (Capak et al. 2015)

<footnote n="4">⁴Based on the results reported in their table 2.</footnote>
and J0817 (Neeleman et al. 2017, 2020).<fnref n="5" /> In this section, we add six new $z \sim 4-6$ rotators from the ALPINE sample to this class: CG32, DC396844, DC494057, DC552206, DC881725, and VC.7875.

### 5.3.1 Previous data

To compare each source on equal footing, we re-analyse the [C II] kinematics of the two previous rotators by fitting TRMs using the ³D BAROLO analysis detailed in Section 3. Both of these sources were previously fit with the 2D tilted ring fitting code GIPSY ROTCUR (Jones et al. 2017), resulting in decent reproductions of each velocity field. However, the two-dimensional nature of this previous work made it impossible to account for beam smearing or velocity dispersion, resulting in fits that were not quite physical.

The data for HZ9 come from Capak et al. (2015), and are the same data cube analysed by Jones et al. (2017). The galaxy J0817 has been observed in [C II] emission with ALMA at two different resolutions: ~1 arcsec (Neeleman et al. 2017) and ~0.2 arcsec (Neeleman et al. 2020). To have similar spatial resolution to the ALPINE data, we choose to analyse the 1 arcsec J0817 [C II] data cube of Neeleman et al. (2017), and will examine the ~0.2 arcsec observation of Neeleman et al. (2020) in Section 5.4. Each of these cubes was passed through the same ³D BAROLO analysis as the ALPINE sources, resulting in the best-fitting parameters listed in Table 1 (see also Fig. 6).

Both of these data cubes are well-fit by ³D BAROLO, with small residuals in each moment map and agreement between the model and data PVDs.

### 5.3.2 Addition of ALPINE sources

We compare each of the best-fitting rotation curves and dynamical mass profiles of the eight un lensed, main-sequence $z \sim 4-6$ rotators in Fig. 7. The dynamical mass is estimated using the radius and circular velocity of each ring, assuming perfectly circular rotation:

$$M_{dyn} = \frac{v^2 r}{G}, \qquad (5)$$

where G is the gravitational constant (e.g. Wang et al. 2010; De Breuck et al. 2014; Lang et al. 2020). The values used to create the rotation curves, velocity dispersion profiles, and dynamical mass profiles are listed in Appendix C (Table C1).

Due to the highly non-linear nature of the kinematic fitting method used in ³D BAROLO, the inclusion of this inclination uncertainty in the rotation curves or dynamical mass profiles is non-trivial. Since ³D BAROLO fits each ring individually, a different inclination will result in both a different spatial distribution and kinematic properties of each ring, vastly changing the model. Due to this complexity, we simply note that the uncertainty in the best-fitting inclination induces an additional source of uncertainty in each rotation curve and dynamical mass profile, which is not depicted in Fig. 7.

The dynamical profiles of these sources may be compared to those derived through 2D modelling by Jones et al. (2017), which includes both main-sequence and starburst galaxies. We find smaller spatial extents (i.e. maximum model ring radii), due to our accounting for beam smearing, resulting in smaller dynamical masses. This may be

seen in HZ9, where the 2D approach yielded a maximum dynamical mass of ~4 × 10¹⁰ M☉ and a spatial extent of 3.8 kpc, while we find $M_{dyn} \sim 2 \times 10^{10} M☉$ and a spatial extent of 2.7 kpc. This confirms the well-known fact that beam smearing has severe impacts on the recovered kinematic parameters of marginally resolved objects, and must be taken into account (e.g. Begeman 1989; de Blok et al. 2008; Molina et al. 2019).

When examining the rotation curves in Fig. 7, it is evident that the rotators in this sample show high inclination-corrected rotation speeds (~50–250 km s⁻¹) and modest spatial extents ($r < 7$ kpc). The dynamical mass profiles for these sources show some variety, with some rising rotation velocity profiles (i.e. DC552206, J0817, HZ9), approximately flat profiles (i.e. CG32, DC494057, VC.7875), and declining profiles (i.e. DC396844, DC881725).

At first glance, these declining rotation curves could be interpreted as evidence for high baryon fractions (or a weak dark matter halo), as suggested in some analyses of $z \sim 1–3$ galaxies (e.g. Genzel et al. 2017; Lang et al. 2017). However, there are several main points that do not allow such a strong conclusion. Primarily, the [C II] observations of these sources have synthesized beam sizes of ~7 kpc, and thus smooth over large swathes of each galaxy (i.e. beam smearing). This effect of low resolution is also seen in the small number of rings in each model. Similarly, while the cores of these galaxies are very well detected, the weaker outskirts may be confused by noise contributions. Finally, each of these TRMs includes a non-zero velocity dispersion for each ring, which is not accounted for in these dynamical masses.

Interestingly, half of the ALPINE rotators also show significant [C II] haloes (DC396844, DC881725, VC.7875; Fujimoto et al. 2020). The coincidence of enriched (i.e. not pristine) haloes and ordered rotation may be caused by a number of phenomena, including star formation-driven outflows, AGN activity (unlikely in the ALPINE sample, which was chosen to exclude type 1 AGN), and past merger activity. Indeed, cosmological simulations (e.g. Kohandel et al. 2019) have shown that galaxies may undergo major mergers and revert to a rotation-dominated state a short time (i.e. tens of Myr) later.

Here, we are able to show that our current data strongly argues for order, disc-like rotation in each of these high-redshift galaxies, and that their rotation curves show a variety of shapes over $r \sim 5$ kpc. The detailed comparison of how these rotators compare to low-redshift rotating galaxies (e.g. $v/\sigma$, $M_{dyn}$ versus $M_*$, resolved gas fractions) will be presented in a future work (Rizzo et al., in preparation). In addition, high-resolution observations will allow us to delineate the small-scale kinematics of these rotating galaxies, and low-resolution observations will reveal the kinematics of their surrounding gaseous haloes.

## 5.4 Effects of higher spatial resolution

Throughout this work, we have focused on [C II] data sets from the ALPINE survey, which feature spatial resolutions of ~1 arcsec. For the ALPINE sources (i.e. $4.4 < z < 5.9$), this spatial scale corresponds to ~6–7 kpc. This resolution is sufficient to identify a number of rotators, mergers, and dispersion-dominated objects, and to estimate the dynamical mass of each. This morpho-kinematic characterization may be made even more precise and informative by increasing the spatial resolution of observations.

To illustrate the effects of increased spatial resolution on the recovered kinematics, we apply our ³D BAROLO analysis to the recent high-resolution ALMA [C II] observation of J0817 ($z = 4.2603$; 0.2 arcsec ~ 1.4 kpc; Neeleman et al. 2020). The 1 arcsec ~

<footnote n="5">⁵Note that we exclude the galaxy HZ10 (likely a close-separation merger, Pavesi et al. 2016) and the strongly gravitationally lensed rotators of Rizzo et al. (2020, 2021).</footnote>
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Moment maps, PVDs, and spectra for observed data, model, and residual for the two previously detected $z > 4$ main-sequence, unlensed rotators not in the ALPINE sample. In each of these figures, the first three rows show (from top to bottom) the moment 0 (integrated intensity), moment 1 (velocity field), and moment 2 (velocity dispersion field; see Section 3.2.1 for details of moment map creation). For these rows, the three columns denote (from left to right) the observed data cubes, model cubes, and the corresponding residuals. The white crosses in the lower right corner of each panel in the first row show a 5 kpc × 5 kpc physical scale. The solid lines in the second row represent the kinematic major axis, while the dashed lines represent the minor axis. The bottom row shows (from left to right) the major axis PVD, minor axis PVD, and extracted spectra. For each PVD, the observed data are shown by the background colour, while the contours represent the model. The data (D), model (M), and residual (R) spectra are depicted by the green, purple, and orange lines, respectively. The $1\sigma$ uncertainty, calculated as (average rms noise level per channel)×√number of beams in source, is shown by the shaded grey area.</figcaption></figure>

6.8 kpc [C II] observations (Neeleman et al. 2017) of this source were analysed in Section 5.3.

We present the results of applying the 3D BAROLO analysis (Section 3) to the high-resolution J0817 [C II] data in Fig. 8. The ~1.4 kpc resolution of these data enables us to resolve the galaxy over many beams (i.e. >5), revealing the presence of a central bulge (see moment zero map), a complex velocity field, and an non-axisymmetric velocity dispersion map. 3D BAROLO is able to reproduce the moment 1 and 2 maps, as well as the PVDs. The spectrum features a significant red residual, which suggests that these high-resolution data expose non-rotating components of the J0817 system.

We may compare the observed morpho-kinematics of J0817 at 0.2 arcsec (Fig. 8) and 1 arcsec (Fig. 6) spatial resolution. The effects of beam smearing are drastically reduced in the high-resolution data, as seen by the more compact moment zero map, complex velocity field, and sharp pair of PVDs. The rotator nature of this source is discernible from the low-resolution data but high-resolution observations allow for more detailed characterization.

Finally, we may compare the recovered rotation curves and dynamical mass profiles for this galaxy at each resolution (Fig. 7). Using the 0.2 arcsec data, we are able to explore the rotation curve at higher spatial resolution. The observed curve is in agreement with a steep initial rise in circular velocity, followed by a flattening at large radii (as seen by Neeleman et al. 2020). The recovered velocities are slightly lower than those recovered from the 1 arcsec resolution data, suggesting a different kinematic inclination angle. However, the two

dynamical mass profiles are in complete agreement, suggesting that the 1 arcsec resolution data are sufficient for this analysis.

## 6 SUMMARY AND CONCLUSIONS

In this work, we have investigated the morpho-kinematics of the star-forming main-sequence galaxies in the ALMA ALPINE survey using three quantitative analyses: a TRM fitting code (3D BAROLO), a morphological classifier ($G-M_{20}$), and the disc-like criteria of Wisnioski et al. (2015).

(i) For each galaxy that we fit with 3D BAROLO 3DFIT (29 sources), we derive the best-fitting morphological and kinematic parameters (e.g. position angle, inclination, rotation curve, velocity dispersion profile, and dynamical mass) and present moment maps, PVDs, and spectra, yielding a more detailed view into their kinematic and morphological nature. By examining the placement of this subset of 29 galaxies on the $M_*$-SFR plane, we find that it contains high-$M_*$ (>10$^{9.5}$ $M_\odot$) objects that are representative of the star-forming main sequence at this redshift. The morpho-kinematic diversity of these $z \sim 4-6$ galaxies is further supported, with six rotators, five mergers, and three dispersion-dominated galaxies in our high-$M_*$ sample. The morpho-kinematic diversity of the ALPINE sample is placed in context with low-redshift studies, suggesting an increase in number of mergers at $z > 4$. This diversity implies for a number of pathways to building up the mass of these sources (e.g. mergers, filamentary accretion).
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Comparison of the best-fitting rotation curves (top) and dynamical mass profiles (bottom) for the six new un lensed, main-sequence $z \sim 4-6$ rotators from ALPINE and the two previous such objects. For J0817, we show the results for both the 1 arcsec (solid grey lines) and 0.2 arcsec data (dashed grey lines; see Section 5.4). Vertical extent of each point in the top panel shows the fit uncertainty output by $^3\text{BAROLO}$ (i.e. not including systematic uncertainty on the inclination), while the vertical extent of each bar in the lower panel is the dynamical mass uncertainty (see equation 5) accounting for the fit uncertainty in rotational velocity output by $^3\text{BAROLO}$ and a radius uncertainty of half a ring width. The horizontal extent of each point shows the width of a ring.</figcaption></figure>

(ii) To examine the morphology of the ALPINE sample in an additional quantitative manner, we perform one of the first investigations of the distribution of high-redshift ($z > 4$) galaxies on the Gini-$M_{20}$ plane (e.g. Lotz et al. 2004) using maps of the integrated [C II] and rest-frame FIR emission. We find that our sources do not separate by morpho-kinematic class and feature low G (signifying a large presence of diffuse emission). To test whether this scatter is due to the low resolution of our observations (i.e. 1–3 beams per source), we simulate the scatter in this plane due to point sources of different SNR, finding that the scatter in G is primarily due to the low resolution of our observations (i.e. $\sim 1$ kpc $px^{-1}$). This groundbreaking test suggests that future, high-resolution ALMA observations of high-$z$ galaxies will enable the use of this diagnostic plot.

<figure><img src="image_9.png" /><figcaption>**Figure 8.** Moment maps, PVDs, and spectra for observed data, model, and residual for the high-resolution (0.2 arcsec) data of J0817 from Neeleman et al. (2020). In each of these figures, the first three rows show (from top to bottom) the moment 0 (integrated intensity), moment 1 (velocity field), and moment 2 (velocity dispersion field; see Section 3.2.1 for details of moment map creation). For these rows, the three columns denote (from left to right) the observed data cubes, model cubes, and the corresponding residuals. The white crosses in the lower right corner of each panel in the first row show a 5 kpc $\times$ 5 kpc physical scale. The solid lines in the second row represent the kinematic major axis, while the dashed lines represent the minor axis. The bottom row shows (from left to right) the major axis PVD, minor axis PVD, and extracted spectra. For each PVD, the observed data are shown by the background colour, while the contours represent the model. The data (D), model (M), and residual (R) spectra are depicted by the green, purple, and orange lines, respectively. The $1\sigma$ uncertainty, calculated as (average rms noise level per channel) $\times \sqrt{\text{number of beams in source}}$, is shown by the shaded grey area.</figcaption></figure>

(iii) We also apply the disc criteria of Wisnioski et al. (2015) to the 29 galaxies where $^3\text{BAROLO}$ 3DFIT modelling was possible due to sufficient resolution and S/N. We find that the first two Wisnioski et al. (2015) criteria are useful for identifying rotating galaxies, but these criteria are not suitable for the current data when used alone.

(iv) We compare the six confirmed main sequence $z \sim 4-6$ rotators in the ALPINE sample to the two previously detected un lensed sources, finding high inclination-corrected rotational velocities, little evidence for declining rotation curves, and comparable dynamical mass profiles. While future high-sensitivity observations will be required to place stronger constraints on the behaviour of these rotation curves at large radii, these current results show strong evidence for ordered rotation in a sample of $z \sim 4-6$ SFGs.

(v) To test the reliability of our recovered morpho-kinematics, we apply the same TRM fitting procedure as above to both a low-resolution ($\sim 1$ arcsec) and high-resolution ($\sim 0.2$ arcsec) [C II] observation of the $z \sim 4.26$ SFG J0817. The low-resolution observation is well-fit by such a model, revealing the intrinsic rotation of the galaxy. Thus, our $\sim 1$ arcsec resolution ALPINE data are appropriate for morpho-kinematic characterization and dynamical mass derivation. By using higher resolution data, we may
confirm these, and derive more detailed kinematics of the system (e.g. deviations from circular rotation).

By applying these analyses to the statistical sample of ALPINE, we may gain insight into the morpho-kinematics of main-sequence galaxies at $z > 4$. Future high-resolution [C II] observations will allow us to precisely determine the merger rate, rotation properties, turbulence support, and dynamical mass of each source. In parallel, additional detailed investigations of single targeted (e.g. Jones et al. 2020; Ginolfi et al. 2020b) and serendipitous (e.g. Romano et al. 2020) galaxies will reveal the morpho-kinematics of early galaxies. Together, these statistical and single observations will inform current cosmological simulations (e.g. Kohandel et al. 2019) and allow us to precisely determine the kinematics (e.g. merger rate, relative populations of kinematic classes, baryon cycle) of galaxies in the early Universe.

## ACKNOWLEDGEMENTS

We thank the anonymous referee for useful feedback that resulted in an improved manuscript. This paper is based on data obtained with the ALMA Observatory, under Large Program 2017.1.00428.L. We also analyse archival data from ALMA programs 2012.1.00523.S and 2017.1.01052.S. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), MOST and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO, and NAOJ. Based on data products from observations made with ESO Telescopes at the La Silla Paranal Observatory under ESO programme ID 179.A-2005 and on data products produced by TERAPIX and the Cambridge Astronomy Survey Unit on behalf of the UltraVISTA consortium. This program is supported by the national program Cosmology and Galaxies from the CNRS in France. GCJ and RM acknowledge ERC Advanced Grant 695671 'QUENCH' and support by the Science and Technology Facilities Council (STFC). JDS was supported by the JSPS KAKENHI Grant Number \*JP18H04346\*, and the World Premier International Research Center Initiative (WPI Initiative), MEXT, Japan. AC and MT acknowledge the support from grant PRIN MIUR 2017 - 20173ML3WW_001. EI acknowledges partial support from FONDECYT through grant N° 1171710. This paper is dedicated to the memory of Olivier Le Fèvre, PI of the ALPINE survey.

## DATA AVAILABILITY

The data analysed in this work are available from the ALMA data archive (https://almascience.nrao.edu/asax/) under project code 2017.1.00428.L. Catalogues and data products are available from the ALPINE Data Release page (https://cesam.lam.fr/a2c2s/data_release.php).

## REFERENCES

* Aretxaga I. et al., 2011, MNRAS, 415, 3831

* Bakx T. J. L. C. et al., 2020, MNRAS, 493, 4294

* Begeman K. G., 1989, A&A, 223, 47

* Béthermin M. et al., 2020, A&A, 643, A2

* Bouché N., Carfantan H., Schroetter I., Michel-Dansac L., Contini T., 2015, AJ, 150, 92

* Braine J., Hughes A., Rosolowsky E., Gratier P., Colombo D., Meidt S., Schinnerer E., 2020, A&A, 633, A17

* Burbidge E. M., Burbidge G. R., Prendergast K. H., 1959, ApJ, 130, 739

* Burkert A. et al., 2016, ApJ, 826, 214

* Capak P. et al., 2008, ApJ, 681, L53

* Capak P. L. et al., 2015, Nature, 522, 455

* Carignan C., Chemin L., Huchtmeier W. K., Lockman F. J., 2006, ApJ, 641, L109

* Carilli C. L. et al., 2008, ApJ, 689, 883

* Carniani S. et al., 2018, MNRAS, 478, 1170

* Cassata P. et al., 2020, A&A, 643, A6

* Chemin L., Carignan C., Foster T., 2009, ApJ, 705, 1395

* Conselice C. J., Bundy K., Ellis R. S., Brichmann J., Vogt N. P., Phillips A. C., 2005, ApJ, 628, 160

* Contini T. et al., 2012, A&A, 539, A91

* Cramer W. J., Kenney J. D. P., Sun M., Crowl H., Yagi M., Jáchym P., Roediger E., Waldron W., 2019, ApJ, 870, 63

* Cresci G. et al., 2009, ApJ, 697, 115

* Davis T. A. et al., 2013, MNRAS, 429, 534

* de Blok W. J. G., Walter F., Brinks E., Trachternach C., Oh S.-H., Kennicutt R. C., Jr, 2008, AJ, 136, 2648

* De Breuck C. et al., 2014, A&A, 565, A59

* den Brok M. et al., 2020, MNRAS, 491, 4089

* Di Teodoro E. M., Fraternali F., 2015, MNRAS, 451, 3021

* Díaz-Santos T. et al., 2018, Science, 362, 1034

* Epinat B., Amram P., Marcelin M., 2008, MNRAS, 390, 466

* Epinat B. et al., 2012, A&A, 539, A92

* Faisst A. L. et al., 2020, ApJS, 247, 61

* Fan L., Knudsen K. K., Han Y., Tan Q.-h., 2019, ApJ, 887, 74

* Förster Schreiber N. M. et al., 2006, ApJ, 645, 1062

* Förster Schreiber N. M. et al., 2009, ApJ, 706, 1364

* Förster Schreiber N. M. et al., 2018, ApJS, 238, 21

* Fraternali F., Karim A., Magnelli B., Gómez-Guijarro C., Jiménez-Andrade E. F., Posses A. C., 2021, A&A, 647, A194

* Fujimoto S. et al., 2019, ApJ, 887, 107

* Fujimoto S. et al., 2020, ApJ, 900, 1

* Garrido O., Marcelin M., Amram P., Boulesteix J., 2002, A&A, 387, 821

* Genzel R. et al., 2017, Nature, 543, 397

* Ginolfi M. et al., 2020a, A&A, 633, A90

* Ginolfi M. et al., 2020b, A&A, 643, A7

* Gnerucci A. et al., 2011, A&A, 528, A88

* Gómez-López J. A. et al., 2019, A&A, 631, A71

* Graziani L., Schneider R., Ginolfi M., Hunt L. K., Maio U., Glatzle M., Ciardi B., 2020, MNRAS, 494, 1071

* Harrison C. M. et al., 2017, MNRAS, 467, 1965

* Jones G. C. et al., 2017, ApJ, 850, 180

* Jones G. C. et al., 2020, MNRAS, 491, L18

* Józsa G. I. G., Kenn F., Klein U., Oosterloo T. A., 2007, A&A, 468, 731

* Kohandel M., Pallottini A., Ferrara A., Zanella A., Behrens C., Carniani S., Gallerani S., Vallini L., 2019, MNRAS, 487, 3007

* Korsaga M., Epinat B., Amram P., Carignan C., Adamczyk P., Sorgho A., 2019, MNRAS, 490, 2977

* Krajnović D., Cappellari M., de Zeeuw P. T., Copin Y., 2006, MNRAS, 366, 787

* Kretschmer M., Dekel A., Teyssier R., 2021, preprint (arXiv:2103.06882)

* Lang P. et al., 2017, ApJ, 840, 92

* Lang P. et al., 2020, ApJ, 897, 122

* Law D. R., Steidel C. C., Erb D. K., Larkin J. E., Pettini M., Shapley A. E., Wright S. A., 2009, ApJ, 697, 2057

* Le Fèvre O. et al., 2020, A&A, 643, A1

* Lelli F., Di Teodoro E. M., Fraternali F., Man A. W. S., Zhang Z.-Y., De Breuck C., Davis T. A., Maiolino R., 2021, Science, 371, 713

* Lilly S. J. et al., 2007, ApJS, 172, 70

* Lilly S. J. et al., 2009, ApJS, 184, 218

* Loiacono F., Talia M., Fraternali F., Cimatti A., Di Teodoro E. M., Caminha G. B., 2019, MNRAS, 489, 681

* Loiacono F. et al., 2021, A&A, 646, A76

* Lotz J. M., Primack J., Madau P., 2004, AJ, 128, 163

* Lotz J. M. et al., 2008, ApJ, 672, 177

* McMullin J. P., Waters B., Schiebel D., Young W., Golap K., 2007, in Shaw R. A., Hill F., Bell D. J., eds, ASP Conf. Ser. Vol. 376, Astronomical Data
Analysis Software and Systems XVI. Astron. Soc. Pac., San Francisco, p. 127

Madau P., Dickinson M., 2014, ARA&A, 52, 415

Maiolino R. et al., 2008, A&A, 488, 463

Mancera Piña P. E. et al., 2020, MNRAS, 495, 3636

Mannucci F. et al., 2009, MNRAS, 398, 1915

Molina J., Ibar E., Swinbank A. M., Sobral D., Best P. N., Smail I., Escala A., Cirasuolo M., 2017, MNRAS, 466, 892

Molina J. et al., 2019, MNRAS, 482, 1499

Neeleman M., Kanekar N., Prochaska J. X., Rafelski M., Carilli C. L., Wolfe A. M., 2017, Science, 355, 1285

Neeleman M., Prochaska J. X., Kanekar N., Rafelski M., 2020, Nature, 581, 269

Noeske K. G. et al., 2007, ApJ, 660, L43

North E. V. et al., 2019, MNRAS, 490, 319

Pallottini A., Ferrara A., Gallerani S., Vallini L., Maiolino R., Salvadori S., 2017, MNRAS, 465, 2540

Pavesi R. et al., 2016, ApJ, 832, 151

Pavesi R. et al., 2018, ApJ, 861, 43

Pavesi R., Riechers D. A., Faisst A. L., Stacey G. J., Capak P. L., 2019, ApJ, 882, 168

Pease F. G., 1918, Proc. Natl. Acad. Sci., 4, 21

Pensabene A., Carniani S., Perna M., Cresi G., Decarli R., Maiolino R., Marconi A., 2020, A&A, 637, A84

Puech M., 2010, MNRAS, 406, 535

Rizzo F., Vegetti S., Powell D., Fraternali F., McKean J. P., Stacey H. R., White S. D. M., 2020, Nature, 584, 201

Rizzo F., Vegetti S., Fraternali F., Stacey H., Powell D., 2021, MNRAS, in press

Rogstad D. H., Lockhart I. A., Wright M. C. H., 1974, ApJ, 193, 309

Romano M. et al., 2020, MNRAS, 496, 875

Romano M. et al., 2021, A&A, preprint (arXiv:2107.10856)

Rubin V. C., Ford W. K., Jr, 1970, ApJ, 159, 379

Salak D., Nakai N., Sorai K., Miyamoto Y., 2020, ApJ, 901, 151

Schinnerer E. et al., 2008, ApJ, 689, L5

Schoenmakers R. H. M., 1999, PhD thesis, Univ. Groningen

Sellwood J. A., Spekkens K., 2015, preprint (arXiv:1509.07120)

Shao Y. et al., 2017, ApJ, 845, 138

Shapiro K. L. et al., 2008, ApJ, 682, 231

Sharma G., Salucci P., Harrison C. M., van de Ven G., Lapi A., 2021, MNRAS, 503, 1753

Shelest A., Lelli F., 2020, A&A, 641, A31

Simon J. D., Bolatto A. D., Leroy A., Blitz L., 2003, ApJ, 596, 957

Simons R. C. et al., 2019, ApJ, 874, 59

Smit R. et al., 2018, Nature, 553, 178

Speagle J. S., Steinhardt C. L., Capak P. L., Silverman J. D., 2014, ApJS, 214, 15

Swinbank A. M. et al., 2017, MNRAS, 467, 3140

Tadaki K.-i. et al., 2020, ApJ, 889, 141

Talia M., Cimatti A., Mignoli M., Pozzetti L., Renzini A., Kurk J., Halliday C., 2014, A&A, 562, A113

Talia M. et al., 2018, MNRAS, 476, 3956

Teuben P., 1995, in Shaw R. A., Payne H. E., Hayes J. J. E., eds, ASP Conf. Ser. Vol. 77, Astronomical Data Analysis Software and Systems IV. Astron. Soc. Pac., San Francisco, p. 398

Tiley A. L. et al., 2019, MNRAS, 485, 934

Turner O. J. et al., 2017, MNRAS, 471, 1280

Valore P., Dainotti M. G., Kopczyński O., 2020, preprint (arXiv:2008.04746)

van der Hulst J. M., Terlouw J. P., Begeman K. G., Zwitser W., Roelfsema P. R., 1992, in Worrall D. M., Biemesderfer C., Barnes J., eds, ASP Conf. Ser. Vol. 25, Astronomical Data Analysis Software and Systems I. Astron. Soc. Pac., San Francisco, p. 131

Varidel M. R. et al., 2019, MNRAS, 485, 4024

Wang R. et al., 2010, ApJ, 714, 699

Whiting M. T., 2012, MNRAS, 421, 3242

Wisnioski E. et al., 2011, MNRAS, 417, 2601

Wisnioski E. et al., 2015, ApJ, 799, 209

Wisnioski E. et al., 2019, ApJ, 886, 124

Yang Y. et al., 2008, A&A, 477, 789

Zanella A., Pallottini A., Ferrara A., Gallerani S., Carniani S., Kohandel M., Behrens C., 2021, MNRAS, 500, 118

## APPENDIX A: ADDITIONAL NOTES ON TILTED RING MODELLING

In Section 3, we present an overview of our signal isolation and TRM fitting procedure, as well as the results and properties of the analysed subsample. Here, we discuss several assumptions that we have made in developing this process.

We have excluded all serendipitous (i.e. non-target) sources, which are analysed in Loiacono et al. (2021). By definition, each of these sources is not at the phase centre of our observations, and thus require a separate primary beam correction for accurate morpho-kinematic analysis. Since they were not selected as part of the ALPINE sample, they do not necessarily represent SFGs at $z \sim 4.4-5.9$. For example, they include sources at $z \sim 1$ (S5110377875, S460378) as well as a well-studied SMG at $z = 4.54$ (S842313, also known as AzTEC/C17; Aretxaga et al. 2011, COSMOS J100054+02343; Carilli et al. 2008; Capak et al. 2008; Schinnerer et al. 2008). Because of this variety, they will not be used in our current goal of probing the kinematic diversity of this specific class of galaxies. However, a number of these sources were successfully fit with our TRM analysis (e.g. S842313), so their morpho-kinematics may be analysed in a future work.

As part of the fitting procedure, we only allow inclination to vary from $10^\circ$ to $80^\circ$, excluding face-on ($i < 10^\circ$) and edge-on ($i > 80^\circ$) discs. Face-on discs feature large spatial extents, but low line-of-sight velocities, resulting in dispersion-dominated, Gaussian line profiles (e.g. Kohandel et al. 2019). These sources are indistinguishable from true dispersion-dominated galaxies. On the other hand, edge-on discs are dominated by velocities along the line of sight, but have small axial ratios, and are likely to be poorly resolved along their minor axis. Because of these features, $^3\text{BAROLO}$ has been found to return poor fits for extreme inclinations (e.g. Di Teodoro & Fraternali 2015).

Our analysis is based on signal isolation using $^3\text{BAROLO}$ SEARCH, rather than a threshold-based 2D or 3D mask. This has the benefit of excluding low-level noise within the mask, but the relatively high lower SNR threshold ($\text{SNR}_{\text{lower}} = 2.5$) causes the $^3\text{BAROLO}$ fits to miss diffuse, low-level emission, as seen in the non-zero residuals in the spectra and moment zero maps. Indeed, recent ALMA studies (e.g. Fujimoto et al. 2019, 2020; Ginolfi et al. 2020a) have revealed extended [C\,II] emission around $z > 4$ galaxies, implying the presence of a diffuse, enriched circumgalactic medium (CGM). Because this CGM emission is not captured with our signal isolation technique, the scope of the kinematic analysis in this paper is focused on the relatively compact interstellar medium (ISM) of the galaxy.

The $^3\text{BAROLO}$ 3DFIT routine allows the user to choose from two methods of normalizing the model. The first ('local') rescales the flux density of each pixel in the model data cube so that the moment zero map of the model is equal to the moment map of the masked input data cube. The second ('azimuthal') performs a similar rescaling procedure, but requires the model to have azimuthal symmetry, resulting in an elliptical model. To account for the complex morphologies of the ALPINE sources, we have assumed 'local' normalization. This choice has no effect on the recovered moment 1 or moment 2 maps, but allows our model to fit models more complex than a simple symmetric disc. In addition, objects that are obviously not simple rotators (such as the merger
VC9780 in Fig. 1) return best-fitting moment 0 maps that have multiple components or disturbed morphologies, adding additional evidence to their merger classifications.

The majority of our sources have only 2–3 rings per model (Fig. 2). This is much fewer than low-redshift studies (e.g. >20 for de Blok et al. 2008; Lang et al. 2020), but is comparable to other $z > 4$ studies (e.g. Jones et al. 2017; Fraternali et al. 2021). This low number of rings is a direct result of the relatively low spatial resolution of these observations (i.e. ~1 arcsec ~ 7 kpc), as this dictates the spatial scale that we may examine. Each of these TRMs is able to reproduce the morphology and kinematics of each galaxy (see Figs 1 and B1–B6). So while our current analysis is not adversely affected by a limited number of free parameters (or ring number), future high-resolution observations will allow for the examination of small-scale morpho-kinematics (e.g. Neeleman et al. 2020).

## APPENDIX B: GALAXY DETAILS

The diversity of this sample makes it worthwhile to comment on each source successfully fitted with $^3\text{DBAROLO}$ individually. Here, we present moment maps, PVDs, and spectra created from the data, $^3\text{DBAROLO}$ model, and residual for each source. In addition, we discuss the reasoning behind each morpho-kinematic classification (see Section 5.1) and compare them to the classifications of Le Fèvre et al. (2020). The rotation curves and velocity dispersion profiles of each source are depicted in Fig. 2, while the best-fitting parameters and W15 criteria are listed in Table 1.

**CANDELS GOODS 32 (ROT):** Originally classified as a merger, we can see that the observed velocity field shows a strong but not perfectly symmetric, gradient (Fig. B1). The spectral residuals are quite low, the best-fitting model is rotation-dominated, and four of the W15 criteria are met, so there is strong evidence that this source is a rotator.

**DEIMOS COSMOS 396844 (ROT):** Since the major axis of the restoring beam is nearly orthogonal to the observed velocity gradient of this source (Fig. B1), its morphological and kinematic position angles are not in agreement (see failure of W15 criterion 4). The model is mostly well fit, adding additional credence to its original classification as a rotator. However, a broad region of residual emission is apparent at ~−250 km s$^{-1}$. At first glance, this residual emission may be interpreted as evidence for a broad blueshifted outflow or low-level satellite. But since it is so weak, we suggest that it is only noise. Due to the low residuals of a rotation-dominated model, we agree with the Le Fèvre et al. (2020) classification of a rotator.

**DEIMOS COSMOS 417567 (UNC):** Due to its complex PVD behaviour (i.e. multiple peaks, no overall trends), this source was originally classified as a merger. This is supported by the poor $^3\text{DBAROLO}$ fit, where significant residuals exist in every plot (Fig. B1). However, we note that this source is one of the weakest in our present sample, and we may simply be experiencing the negative effects of poor SNR. While all W15 criteria besides PA agreement are met, the large residuals and odd appearance of each PVD earn this source a ‘uncertain’ classification.

**DEIMOS COSMOS 432340 (UNC):** The curved blue half of this galaxy was originally interpreted as evidence for an ongoing merger. $^3\text{DBAROLO}$ returns a reasonable rotation-dominated fit, and all five W15 criteria are met. However, the morphology of this source is quite disturbed, and the minor axis PVD shows a ‘C’-shape (Fig. B1). The SNR of this source is too low to definitely classify it as an ongoing merger, a galaxy with an outflow, or a

low-SNR well-ordered rotating disc, so we classify this source as ‘uncertain’.

**DEIMOS COSMOS 434239 (MER):** With a double intensity peak and unique PVD, this source was originally classified as a merger. We note that this source returns the same W15 criteria as DC552206 (a rotator), while the two galaxies show drastically different morpho-kinematics, highlighting the need to examine the morpho-kinematics of each source in multiple ways. The large residuals in the velocity field and spectrum, as well as the poor match in the major axis PVD (Fig. B2), cause us to classify this source as a merger.

**DEIMOS COSMOS 454608 (UNC):** The double-peaked major axis PVD originally earned this galaxy the classification of merger, but our $^3\text{DBAROLO}$ fits are inconclusive (Fig. B2). The TRM returns large moment one residuals, and a positive-velocity clump is included by SEARCH, making interpretation difficult. This clump may be excluded by adopting a higher SNR$_{lower}$ threshold, but then this source is too compact to be fit by $\ge 2$ rings, and so the fit fails. Due to the low-SNR and low spatial resolution, we classify this source as ‘uncertain’.

**DEIMOS COSMOS 494057 (ROT):** This source is nearly a perfect example of a rotator, with only small residuals in the moment zero and two maps and spectrum, and well-fit PVDs (Fig. B2). We do note that the moment one map shows residuals around the perimeter of the source, but this is likely due to the low-SNR of these regions. In addition, the morphological and kinematic position angles disagree (resulting in a failure of W15 criterion four), but this is likely due to the near-unity axial ratio (or low inclination) of the source.

**DEIMOS COSMOS 519281 (UNC):** While most sources in the present sample show monotonic velocity gradients, DC519281 features a central blueshifted region flanked by two redshifted components (Fig. 1). This may be interpreted as evidence for a counterrotating merger (see Q1623-BX528 in Förster Schreiber et al. 2006). However, the small spatial extent of this source makes physical interpretation impossible, and $^3\text{DBAROLO}$ returns a poor fit. Of the five Wisnioski et al. (2015) criteria, this source only meets the second ($v_{rot} > \sigma_v$), agreeing that this is not a rotator. We also note that when applying SEARCH, the standard lower SNR of 2.5 returned a noise peak to the south-west, so we applied an SNR$_{lower}$ of 2.7 for this source. While the complex velocity field of this source could be interpreted as evidence for ongoing merging activity, the low resolution and SNR of this source cause us to classify it as uncertain.

**DEIMOS COSMOS 552206 (ROT):** The redshift of each source in the ALPINE sample was taken from previous UV spectroscopic observations. However, the velocity difference between these redshifts and those of [C II] may be substantial (i.e. up to ~500 km s$^{-1}$, see Cassata et al. 2020 for details and physical interpretation). In the case of DC552206, this velocity shift caused the [C II] line to fall on the edge of an ALMA sideband, resulting in an incomplete line profile (Fig. B2). However, since the observed profile features a prominent velocity gradient (seen in the moment 1 map and major axis PVD), $^3\text{DBAROLO}$ is still able to fit an agreeable model, suggesting that this source is truly a rotator, rather than a merger, as suggested previously. Of course, a follow-up observation containing the blue edge of this line should be conducted to confirm this.

**DEIMOS COSMOS 627939 (UNC):** This source was originally classified as a merger, mainly due to its complex PVDs. $^3\text{DBAROLO}$ returns a reasonable fit, which reproduces the primary velocity gradient, as well as the overall shape of each PVD (Fig. B3).
<figure><img src="image_10.png" /><figcaption>**Figure B1.** Moment maps, PVDs, and spectra for observed data, model, and residual. In each of these figures, the first three rows show (from top to bottom) the moment 0 (integrated intensity), moment 1 (velocity field), and moment 2 (velocity dispersion field; see Section 3.2.1 for details of moment map creation). For these rows, the three columns denote (from left to right) the observed data cubes, model cubes, and the corresponding residuals. The white crosses in the lower right corner of each panel in the first row show a 5 kpc × 5 kpc physical scale. The solid lines in the second row represent the kinematic major axis, while the dashed lines represent the minor axis. The bottom row shows (from left to right) the major axis PVD, minor axis PVD, and extracted spectra. For each PVD, the observed data are shown by the background colour, while the contours represent the model. The data (D), model (M), and residual (R) spectra are depicted by the green, purple, and orange lines, respectively. The $1\sigma$ uncertainty, calculated as (average rms noise level per channel)×√number of beams in source, is shown by the shaded grey area.</figcaption></figure>

The returned W15 criteria are identical to those of DC396844 and DC494057, which are both classified as rotators. However, substantial residuals are present in both the velocity map and spectrum, hinting at ongoing merging below our spatial resolution.

Since we can neither fit a rotating disc model nor distinguish multiple components, we classify this source as ‘uncertain’.

**DEIMOS COSMOS 683613 (UNC):** While this source was originally classified as an extended dispersion-dominated galaxy,
<figure><img src="image_11.png" /><figcaption>Figure B2. Moment maps, PVDs, and spectra for observed data, model, and residual. See Fig. B1 for details.</figcaption></figure>

further examination of the [C II] moment maps and PVDs presented in Fig. B3 show that the [C II] appears to mainly emanate from a central source, with a weaker, blue, compact source to the south. The dominant source is moderately well fit by a TRM, and appears to be a weaker version of the dispersion-dominated source DC848185. However, since it is not clear if the southern component is evidence of merging or simply noise, we classify this as uncertain.

**DEIMOS COSMOS 733857 (UNC):** With no obvious velocity gradient (resulting in a failure of W15 criterion one) and a poorly fit rotation-dominated model (Fig. B3), this source is unlikely to be

a rotator. It was originally classified as an extended, dispersion-dominated galaxy, but since it is only marginally resolved and features a complex major axis PVD, we classify it as uncertain.

**DEIMOS COSMOS 773957 (UNC):** As a relatively weak, spatially compact source, this source is difficult to characterize. The integrated spectrum is moderately well fit, but each of the PVDs show that the 3D BAROLO fit is obviously not sufficient (Fig. B3). The present observations are insufficient to determine whether this morpho-kinematic complexity is noise or truly signal. Therefore, we classify this as 'uncertain'.
<figure><img src="image_12.png" /><figcaption>Figure B3. Moment maps, PVDs, and spectra for observed data, model, and residual. See Fig. B1 for details.</figcaption></figure>

**DEIMOS COSMOS 818760 (MER):** As previously explored in detail (Jones et al. 2020), DC818760 is most likely a system of three galaxies undergoing a merger. The two dominant galaxies are closely related in both space and velocity, and an east-west monotonotic velocity gradient is evident (Fig. B4). Due to our relatively low resolution, a TRM is able to replicate the flux density distribution, overall velocity gradient, and velocity dispersion of these two merging sources. However, the spectrum shows a significant red residual due to the third galaxy in this system not being captured by the fit, the model velocity dispersion and rotational velocity are nearly identical at all radii (Fig. 2), and the spectral width of the

eastern galaxy is poorly captured. So while this good fit highlights the ambiguity of close mergers and rotating discs, we conclude that this source is a triple merger.

**DEIMOS COSMOS 848185 (DIS):** Also known as HZ6 (Capak et al. 2015) or LBG-1 (Pavesi et al. 2019), this source exhibits strong [C II] emission but little evidence for rotation (Fig. B4). However, 3D BAROLO still returns an excellent fit, with a large velocity dispersion (see also Fig. 2). Thus, we agree that this source is extended and dispersion-dominated.

**DEIMOS COSMOS 873321 (MER):** Originally described as a merger in Capak et al. (2015, called HZ8), DC873321 features two
<figure><img src="image_13.png" /><figcaption>Figure B4. Moment maps, PVDs, and spectra for observed data, model, and residual. See Fig. B1 for details.</figcaption></figure>

bright peaks of emission at the same velocity (Fig. B4). $^{3D}$BAROLO returns a poor fit with comparable rotation velocity and velocity dispersion, strengthening this classification of ‘merger’.

**DEIMOS COSMOS 873756 (DIS):** Even with the brightest [C II] luminosity in the ALPINE sample, this source shows very little evidence for ordered rotation (Fig. B4). Indeed, $^{3D}$BAROLO is able to fit a dispersion-dominated model, resulting in only small residuals. The two PVDs are identical, with no evident gradient. There is some evidence for extended emission to the west at $v \sim 0$ km s$^{-1}$, which is also present as blue residuals in the spectrum. It is not clear if this

represents a minor merger, an outflow perpendicular to the line of sight, or another feature. Since the main galaxy dominates the field, we conclude that this source is extended and dispersion-dominated.

**DEIMOS COSMOS 881725 (ROT):** With a strong velocity gradient and central velocity dispersion peak, this source is well fit by $^{3D}$BAROLO, resulting in low spectral residuals (Fig. B5). The moment one residuals are complex but are mostly non-zero around the low-SNR perimeter. Thus, we agree with the Le Fèvre et al. (2020) classification of a rotator.

**VUDS COSMOS 5100537582 (UNC):** This source shows little to no velocity gradient, and PVDs along both axes are nearly
<figure><img src="image_14.png" /><figcaption>Figure B5. Moment maps, PVDs, and spectra for observed data, model, and residual. See Fig. B1 for details.</figcaption></figure>

identical (Fig. B5), earning it a Le Fèvre et al. (2020) classification of 3 (extended dispersion-dominated). However, $^3D_{\text{BAROLO}}$ fits it with a rotation-dominated model, and its W15 values are the same as those of DC733857, a UNC source. Visual inspection of the moment maps and PVDs shows that we lack the sensitivity and velocity resolution to determine the true morpho-kinematics of this source. Due to this ambiguity, we classify this source as ‘uncertain’.

**VUDS COSMOS 5100541407 (UNC):** With a significant extension to the north-east (Fig. B5) and a number of W15 criteria failures, VC.1407 is difficult to characterize. $^3D_{\text{BAROLO}}$ fits a moderately dispersion-dominated model, which results in significant

residuals. Its extended morphology earned it a merger classification by Le Fèvre et al. (2020), but we classify it as ‘uncertain’.

**VUDS COSMOS 5100559223 (UNC):** This source is quite weak, as our SEARCH technique includes two single-channel clumps (Fig. B5, similarly to DC454608). It was originally classified as an extended dispersion-dominated source, but it is poorly fit with a dispersion-dominated $^3D_{\text{BAROLO}}$ model and shows a possible velocity gradient in the major axis PVD. Until higher sensitivity observations are taken, we classify this source as ‘uncertain’.
<figure><img src="image_15.png" /><figcaption>Figure B6. Moment maps, PVDs, and spectra for observed data, model, and residual. See Fig. B1 for details.</figcaption></figure>

**VUDS COSMOS 5100822662 (MER)**: Similarly to DC873321, VC.2662 features two components at the same velocity (Fig. B6). But unlike this other source, one component (which is the targeted galaxy) dominates the moment zero map. 3D BAROLO returns a poor fit, adding further evidence that this is a merging system.

**VUDS COSMOS 5100994794 (UNC)**: 3D BAROLO returns a rotation-dominated model (see passed W15 criterion two), and the residuals are low, so it is possible that this source is a rotator, rather than the previous classification of an extended dispersion-dominated galaxy. In addition, the set of W15 criteria are identical to those of DC552206, a rotator. However, its major axis PVD shows

a broad region of red emission that dominates the velocity gradient (Fig. B6), and while the model is rotation-dominated at small radii, it is dispersion-dominated at large radii. Thus, we classify it as ‘uncertain’.

**VUDS COSMOS 5101209780 (MER)**: The moment zero map of this source is qualitatively similar to that of VC.2662: a dominant source to the south-west, with a secondary source to the north-east (Fig. 1). However, here the targeted ALPINE source is the subdominant [C II] emitter. The [C II] emission is relatively well-fit with a rotation-dominated TRM, but the discrepancies in each PVD reveal the underlying merger nature of this source. Its kinematics
VE530029038

<figure><img src="image_16.png" /><figcaption>Figure B7. Moment maps, PVDs, and spectra for observed data, model, and residual. See Fig. B1 for details.</figcaption></figure>

and other properties are investigated in depth in Ginolfi et al. (2020b).

**VUDS COSMOS 5101218326 (DIS):** Similarly to DC873756, this source is quite luminous in [C II], and the ³D BAROLO fit strongly suggests an extended dispersion-dominated galaxy (Fig. 1). In addition, nearly all W15 criteria fail, strongly arguing against rotation.

**VUDS COSMOS 510786441 (UNC):** With an elongated morphology and disturbed major axis PVD, VC.6441 was originally classified as a merger. However, ³D BAROLO fits a marginally rotation-dominated model to this source (Fig. B6), and the W15 criteria are similar to rotators in the sample (e.g. DC552206). Examination of the PVDs shows that the underlying morphokinematics are too complex to be classified as a rotator, but our observations are too low resolution to unambiguously classify this as a merger. We thus classify it as ‘uncertain’.

**VUDS COSMOS 5110377875 (ROT):** This source is perhaps the best case of an ordered rotator in the current sample but it is not quite ideal (Fig. 1). While the spectral and moment zero residuals are low, the velocity dispersion peak is not coincident with the intensity peak, the major axis PVD is not symmetric across $v = 0$ km s⁻<fnref n="1" />, and the velocity field residuals are significant. As suggested by Fujimoto et al. (2020), these discrepancies may be caused by interaction with an extended halo or an extended outflow. Despite these imperfections, we classify this as a rotator.

**VUDS COSMOS 5180966608 (UNC):** VC.6608 is one of the most mysterious sources in the ALPINE sample. On one hand, its spectrum is well fit by a high-velocity dispersion TRM, suggesting a single, extended source (Fig. B6). However, the two intensity peaks evident in the major axis PVD are nearly coincident with point sources in HST /ACS F814W, Subaru i+, and Subaru r+ images (Faiset et al. 2020). These point sources may indicate the bright cores of closely interacting galaxies or may be caused by a dense screen of dust blocking the central rest-frame optical continuum emission. The FIR continuum emission peaks between these two point sources but extends over both of them. High-resolution spectral line imaging

of this source is required to determine its nature. Until then, we classify this as ‘uncertain’.

**VUDS ECDFS 530029038 (UNC):** Similarly to other sources, VE.9038 features a significant velocity offset between its [C II] and UV spectral line redshifts (Cassata et al. 2020). Instead of shifting the [C II] profile to the edge of a sideband, this emission is shifted to the intersection of the two SPWs that compose the sideband (Fig. B7). Due to the bandpass shape of ALMA SPWs, a few edge channels at the edge of each SPW must be flagged. This extra flagging occasionally creates gaps in frequency coverage between bands, as seen here, where the central few channels of the [C II] profile are flagged. These missing channels at $v \sim 0$ km s⁻<fnref n="1" /> have an interesting effect, as the resulting moment 1 map is less smoothed in frequency, and the velocity gradient therefore appears more pronounced than it would with full velocity coverage. Despite this mostly cosmetic feature, ³D BAROLO is able to fit a rotating disc model to the [C II] emission. However, the large velocity residual on the red half, and the incomplete velocity coverage mean that additional observations must be made before this source is characterized.

## APPENDIX C: BEST-FITTING PROFILES OF ROTATORS

<table><caption>**Table C1.** Best-fitting rotational velocities, velocity dispersions, and dynamical masses for each model ring, as output by ³D BAROLO. In addition to the six confirmed rotators in the ALPINE sample, we also present the results of fitting the two previously detected $4 < z < 6$ main-sequence unensled rotators. Radius values are assumed ring radii (see Section 3.2.1). Dynamical masses are estimated using equation (5).</caption><thead><tr><th>Source</th><th><i>R</i><br>(kpc)</th><th><i>v</i><sub>rot</sub><br>(km s<sup>-1</sup>)</th><th><i>σ</i><sub>v</sub><br>(km s<sup>-1</sup>)</th><th><i>M</i><sub>dyn</sub><br>(M<sub>☉</sub>)</th></tr></thead><tbody><tr><td rowspan="2">CG32</td><td>1.17</td><td>117.43 ± 20.71</td><td>39.06 ± 14.95</td><td>(3.8 ± 2.3) × 10<sup>9</sup></td></tr><tr><td>3.5</td><td>115.04 ± 26.96</td><td>19.14 ± 12.47</td><td>(1.1 ± 0.5) × 10<sup>10</sup></td></tr><tr><td rowspan="2">DC396844</td><td>1.25</td><td>120.88 ± 10.74</td><td>44.42 ± 7.75</td><td>(4.2 ± 2.3) × 10<sup>9</sup></td></tr><tr><td>3.75</td><td>80.42 ± 17.66</td><td>19.84 ± 11.75</td><td>(5.6 ± 2.6) × 10<sup>9</sup></td></tr><tr><td rowspan="2">DC494057</td><td>1.04</td><td>78.44 ± 10.25</td><td>53.24 ± 6.21</td><td>(1.5 ± 0.8) × 10<sup>9</sup></td></tr><tr><td>3.13</td><td>80.31 ± 12.75</td><td>44.98 ± 7.1</td><td>(4.7 ± 1.7) × 10<sup>9</sup></td></tr><tr><td rowspan="3">DC552206</td><td>1.11</td><td>71.28 ± 21.74</td><td>116.65 ± 11.85</td><td>(1.3 ± 1.0) × 10<sup>9</sup></td></tr><tr><td>3.32</td><td>164.98 ± 20.23</td><td>67.4 ± 13.24</td><td>(2.1 ± 0.6) × 10<sup>10</sup></td></tr><tr><td>5.53</td><td>172.84 ± 27.63</td><td>65.34 ± 14.61</td><td>(3.8 ± 1.3) × 10<sup>10</sup></td></tr><tr><td rowspan="2">DC881725</td><td>1.15</td><td>98.37 ± 9.54</td><td>40.45 ± 6.94</td><td>(2.6 ± 1.4) × 10<sup>9</sup></td></tr><tr><td>3.44</td><td>62.07 ± 12.54</td><td>48.4 ± 7.77</td><td>(3.1 ± 1.3) × 10<sup>9</sup></td></tr><tr><td rowspan="2">VC.7875</td><td>1.28</td><td>118.97 ± 8.1</td><td>47.04 ± 5.69</td><td>(4.2 ± 2.2) × 10<sup>9</sup></td></tr><tr><td>3.84</td><td>102.85 ± 19.84</td><td>60.84 ± 12.74</td><td>(9.4 ± 4.0) × 10<sup>9</sup></td></tr><tr><td rowspan="2">J0817</td><td>0.95</td><td>247.46 ± 8.87</td><td>32.17 ± 11.36</td><td>(1.4 ± 0.7) × 10<sup>10</sup></td></tr><tr><td>2.84</td><td>252.09 ± 14.94</td><td>35.98 ± 10.47</td><td>(4.2 ± 0.9) × 10<sup>10</sup></td></tr><tr><td rowspan="3">HZ9</td><td>0.54</td><td>155.85 ± 16.84</td><td>71.1 ± 8.38</td><td>(3.0 ± 1.7) × 10<sup>9</sup></td></tr><tr><td>1.61</td><td>156.77 ± 19.02</td><td>75.12 ± 9.22</td><td>(9.2 ± 2.7) × 10<sup>9</sup></td></tr><tr><td>2.68</td><td>176.63 ± 25.45</td><td>4.82 ± 8.0</td><td>(1.9 ± 0.6) × 10<sup>10</sup></td></tr></tbody></table>

<footnote n="1">¹Cavendish Laboratory, University of Cambridge, 19 J. J. Thomson Ave., Cambridge CB3 0HE, UK</footnote>

<footnote n="2">²Kavli Institute for Cosmology, University of Cambridge, Madingley Road, Cambridge CB3 0HA, UK</footnote>

<footnote n="3">³INAF – Osservatorio di Astrofisica e Scienza dello Spazio, via Gobetti 93/3, I-40129 Bologna, Italy</footnote>
⁴Dipartimento di Fisica e Astronomia, Università di Padova, vicolo dell'Osservatorio 3, I-35122 Padova, Italy

⁵INAF – Osservatorio Astrofisico di Padova, vicolo dell'Osservatorio 5, I-35122 Padova, Italy

⁶European Southern Observatory, Karl-Schwarzschild-Str 2, D-85748 Garching, Germany

⁷Research Institute for Science and Engineering, Waseda University, 3-4-1 Okubo, Shinjuku, Tokyo 169-8555, Japan

⁸National Astronomical Observatory of Japan, 2-21-1 Osawa, Mitaka, Tokyo 181-8588, Japan

⁹Aix-Marseille Université, CNRS, LAM (Laboratoire d'Astrophysique de Marseille), UMR 7326, F-13388, Marseille, France

¹⁰Cosmic Dawn Center (DAWN), Jagtvej 128, DK-2200 Copenhagen N, Denmark

¹¹Niels Bohr Institute, University of Copenhagen, Lyngbyvej 2, DK-2100 Copenhagen Ø, Denmark

¹²Department of Physics & Astronomy, University of California, Davis, One Shields Ave., Davis, CA 95616, USA

¹³Infrared Processing and Analysis Center, California Institute of Technology, Pasadena, CA 91125, USA

¹⁴Observatoire de Genève, Université de Genève, 51 Ch. des Maillettes, CH-1290 Versoix, Switzerland

¹⁵Kavli Institute for the Physics and Mathematics of the Universe, The University of Tokyo, Kashiwa, Chiba 277-8583, Japan

¹⁶Department of Astronomy, School of Science, The University of Tokyo, 7-3-1 Hongo, Bunkyo, Tokyo 113-0033, Japan

¹⁷The Caltech Optical Observatories, California Institute of Technology, Pasadena, CA 91125, USA

¹⁸Centro de Astronomía (CITEVA), Universidad de Antofagasta, Avenida Angamos 601, Antofagasta 1270300, Chile

¹⁹Department of Physics and Astronomy (DIFA), University of Bologna, Via Gobetti 93/2, I-40129 Bologna, Italy

²⁰INAF – Osservatorio Astrofisico di Arcetri, Largo E. Fermi 5, I-50125 Firenze, Italy

²¹Instituto de Física y Astronomía, Universidad de Valparaíso, Avda. Gran Bretaña 1111, Valparaíso 2340000, Chile

This paper has been typeset from a TEX/LATEX file prepared by the author.