We present highresolution inlineformulatexmath idTM0002 notationLaTeXsim 0.14texmathinlineformula arcsec  710 pc Atacama Large Millimetresubmillimetre Array C II 158 inlineformulatexmath idTM0003 notationLaTeXmutexmathinlineformulam and dust continuum followup observations of REBELS25 a C IIluminous inlineformulatexmath idTM0004 notationLaTeXLmathrmCII1.7pm 0.2times 109 mathrmLodot texmathinlineformula galaxy at redshift inlineformulatexmath idTM0005 notationLaTeXz7.3065pm 0.0001texmathinlineformula. These highresolution high signaltonoise observations allow us to study the subkpc morphology and kinematics of this massive inlineformulatexmath idTM0006 notationLaTeXM  842 times 109 mathrmrm Modot texmathinlineformula starforming SFRinlineformulatexmath idTM0007 notationLaTeXmathrmUVIR  19910163 mathrmrm Modot  mathrmyr1texmathinlineformula galaxy in the Epoch of Reionization. By modelling the kinematics with inlineformulatexmath idTM0008 notationLaTeXmathrm3DtexmathinlineformulaBAROLO we find it has a lowvelocity dispersion inlineformulatexmath idTM0009 notationLaTeXbarsigma   3397texmathinlineformula km sinlineformulatexmath idTM0010 notationLaTeX1texmathinlineformula and a high ratio of orderedtorandom motion inlineformulatexmath idTM0011 notationLaTeXVmathrmrot maxbarsigma   11 65texmathinlineformula indicating that REBELS25 is a dynamically cold disc. Additionally we find that the C II distribution is well fit by a nearexponential disc model with a Srsic index n of inlineformulatexmath idTM0013 notationLaTeX1.3 pm 0.2texmathinlineformula and we see tentative evidence of more complex nonaxisymmetric structures suggestive of a bar in the C II and dust continuum emission. By comparing to other high spatial resolution cold gas kinematic studies we find that dynamically cold discs seem to be more common in the highredshift Universe than expected based on prevailing galaxy formation theories which typically predict more turbulent and dispersiondominated galaxies in the early Universe as an outcome of merger activity gas accretion and more intense feedback. This higher degree of rotational support seems instead to be consistent with recent cosmological simulations that have highlighted the contrast between cold and warm ionized gas tracers particularly for massive galaxies. We therefore show that dynamically settled disc galaxies can form as early as 700 Myr after the big bang
['Astrophysics - Astrophysics of Galaxies']
REBELS25 discovery of a dynamically cold disc galaxy at z  7.31
{'REBELS-25: Discovery of a dynamically cold disc galaxy at 𝑧 = 7 . 31': "Lucie E. Rowland 1 ★ , Jacqueline Hodge 1 , Rychard Bouwens 1 , Pavel Mancera Piña 1 , Alexander Hygate 1 , Hiddo Algera 2 , 3 , Manuel Aravena 4 , Rebecca Bowler 5 , Elisabete da Cunha 6 , 7 , Pratika Dayal 8 , Andrea Ferrara 9 , Thomas Herard-Demanche 1 , Hanae Inami 2 , Ivana van Leeuwen 1 , Ilse de Looze 10 , Pascal Oesch 11 , 12 , Andrea Pallottini 9 , Siân Phillips 13 , Matus Rybak 14 , 1 , 15 , Sander Schouws 1 , Renske Smit 13 , Laura Sommovigo 16 , Mauro Stefanon 17 , 18 , Paul van der Werf 1  \n- 1 Leiden Observatory, Leiden University, P.O. Box 9513, 2300 RA Leiden, The Netherlands\n- 2 Hiroshima Astrophysical Science Center, Hiroshima University, 1-3-1 Kagamiyama, Higashi-Hiroshima, Hiroshima 739-8526, Japan\n- 3 National Astronomical Observatory of Japan, 2-21-1, Osawa, Mitaka, Tokyo, Japan\n- 4 Núcleo de Astronomía, Facultad de Ingeniería y Ciencias, Universidad Diego Portales, Av. Ejército Libertador 441, Santiago, Chile\n- 5 Jodrell Bank Centre for Astrophysics, Department of Physics and Astronomy, School of Natural Sciences, The University of Manchester, Manchester M13 9PL, UK\n- 6 International Centre for Radio Astronomy Research, University of Western Australia, 35 Stirling Hwy, Crawley, WA 6009, Australia\n- 7 ARC Centre of Excellence for All Sky Astrophysics in 3 Dimensions (ASTRO 3D)\n- 8 Kapteyn Astronomical Institute, University of Groningen, 9700 AV Groningen, The Netherlands\n- 9 Scuola Normale Superiore, Piazza dei Cavalieri 7, 56126, Pisa, Italy\n- 10 Sterrenkundig Observatorium, Ghent University, Krĳgslaan281-S9, B-9000 Ghent, Belgium\n- 11 Department of Astronomy, University of Geneva, Chemin Pegasi 51, 1290 Versoix, Switzerland\n- 12 Cosmic Dawn Center (DAWN), Copenhagen, Denmark\n- 13 Astrophysics Research Institute, Liverpool John Moores University, 146 Brownlow Hill, Liverpool L3 5RF, United Kingdom\n- 14 Faculty of Electrical Engineering, Mathematics and Computer Science, Delft University of Technology, Mekelweg 4, 2628 CD Delft, The Netherlands\n- 15 SRON - Netherlands Institute for Space Research, Niels Bohrweg 4, 2333 CA Leiden, The Netherlands\n- 16 Center for Computational Astrophysics, Flatiron Institute, 162 5th Avenue, New York, NY 10010, USA\n- 17 Departament d'Astronomia i Astrofísica, Universitat de València, C. Dr. Moliner 50, E-46100 Burjassot, València, Spain\n- 18 Unidad Asociada CSIC 'Grupo de Astrofísica Extragaláctica y Cosmología' (Instituto de Física de Cantabria - Universitat de València)  \nAccepted XXX. Received YYY; in original form ZZZ", 'ABSTRACT': 'We present high resolution ( ∼ 0 . 14" = 710 pc) ALMA [CII] 158 𝜇 mand dust continuum follow-up observations of REBELS-25, a [CII]-luminous ( 𝐿 [ CII ] = ( 1 . 7 ± 0 . 2 ) × 10 9 L ⊙ ) galaxy at redshift 𝑧 = 7 . 3065 ± 0 . 0001. These high resolution, high signal-tonoise observations allow us to study the sub-kpc morphology and kinematics of this massive ( 𝑀 ∗ = 8 + 4 -2 × 10 9 M ⊙ ) star-forming (SFRUV + IR = 199 + 101 -63 M ⊙ yr -1 ) galaxy in the Epoch of Reionisation. By modelling the kinematics with 3D BAROLO , we find it has a low velocity dispersion ( ¯ 𝜎 = 33 ± 9 km s -1 ) and a high ratio of ordered-to-random motion ( 𝑉 rot , max / ¯ 𝜎 = 11 + 8 -4 ), indicating that REBELS-25 is a dynamically cold disc. Additionally, we find that the [CII] distribution is well fit by a near-exponential disc model, with a Sérsic index, 𝑛 , of 1 . 3 ± 0 . 2, and we see tentative evidence of more complex non-axisymmetric structures suggestive of a bar in the [CII] and dust continuum emission. By comparing to other high spatial resolution cold gas kinematic studies, we find that dynamically cold discs seem to be more common in the high redshift Universe than expected based on prevailing galaxy formation theories, which typically predict more turbulent and dispersion-dominated galaxies in the early Universe as an outcome of merger activity, gas accretion and more intense feedback. This higher degree of rotational support seems instead to be consistent with recent cosmological simulations that have highlighted the contrast between cold and warm ionised gas tracers, particularly for massive galaxies. We therefore show that dynamically settled disc galaxies can form as early as 700 Myr after the Big Bang.  \nKey words: galaxies: evolution - galaxies: high-redshift - galaxies: kinematics and dynamics', '1 INTRODUCTION': 'In the current picture of galaxy formation, galaxies in the early Universe should be more dominated by turbulent motion than their lower redshift counterparts, as a result of an increase in merger activity and  \nmore violent gas accretion (Conselice et al. 2008; Dekel et al. 2009), as well as more intense active galactic nuclei (AGN) and stellar feedback (e.g., Hayward & Hopkins 2017; Nelson et al. 2019). With the advent of the Hubble Space Telescope (HST), high angular resolution rest-frame UV images of redshift 𝑧 > 1 galaxies indeed revealed more clumpy and irregular morphologies (e.g., Conselice et al. 2003; Papovich et al. 2005; Elmegreen & Elmegreen 2005; Conselice et al.  \n2008), in comparison to the rotationally supported disc galaxies that make up 50-80% of the 𝑧 ∼ 0-1 Universe (e.g., Kassin et al. 2012; Wisnioski et al. 2015; Swinbank et al. 2017). In addition, multiple integral field (IFU) surveys of warm ionised gas kinematics have found that the turbulence (as measured by the velocity dispersion, 𝜎 , of the gas) within galaxies increases from 𝑧 = 0 to 𝑧 ∼ 3, and that the degree of rotational support (rotational velocity over velocity dispersion, 𝑉 rot / 𝜎 ) decreases over the same redshift range (e.g., Wisnioski et al. 2015; Di Teodoro et al. 2016; Johnson et al. 2018; Wisnioski et al. 2019; Übler et al. 2019).  \nHowever, these findings have recently been challenged on multiple fronts. JWST, with its unprecedented high-angular resolution and sensitivity in the near- and mid-infrared, has now revealed more massive ( 𝑀 ∗ ≳ 10 9 -10 10 M ⊙ ) disc-like galaxies at 𝑧 > 6 (Adams et al. 2022; Akins et al. 2023; Labbé et al. 2023; Rodighiero et al. 2023; Casey et al. 2023; Atek et al. 2023) than predicted from HST observations, hinting that significant disc formation may have occurred much earlier than expected. Similarly, both JWST and the Atacama Large Millimetre/submillimetre Array (ALMA) have observed an increasing number of dynamically cold ( 𝑉 rot / 𝜎 > 2) disc galaxies at 𝑧 > 3 (Smit et al. 2018; Sharda et al. 2019; Neeleman et al. 2020; Rizzo et al. 2020; Neeleman et al. 2021; Lelli et al. 2021; Rizzo et al. 2021; Fraternali et al. 2021; Tsukui & Iguchi 2021; Roman-Oliveira et al. 2023; Posses et al. 2023; Pope et al. 2023; Neeleman et al. 2023; de Graaff et al. 2023; Parlanti et al. 2023a; Übler et al. 2024; Xu et al. 2024). Cosmological simulations often struggle to reproduce these findings, with simulated rotationally supported gas discs typically only forming at 𝑧 ≲ 2 (Simons et al. 2019; Pillepich et al. 2019).  \nThese discrepancies and recent discoveries have made the early build-up of galaxies a forefront challenge in extragalactic astronomy today. To address the tension between observations and theory, some studies over the past few years have suggested that the \'disciness\' of galaxies is more dependent on the mass of the halo, rather than on the redshift (Dekel et al. 2020; Gurvich et al. 2022; Tiley et al. 2021; Kohandel et al. 2023), with gas discs only surviving in 𝑀 ∗ ≳ 10 9 M ⊙ haloes. However, discrepancies remain between the previous findings of an increase in turbulence with redshift, in comparison to the growing number of discoveries of dynamically cold discs at 𝑧 > 3. Likely explanations could include the spatial resolution, signal-to-noise ratio (SNR), sample selection, and kinematic tracer used for observations (e.g. Rizzo et al. 2022; Kohandel et al. 2023). It has therefore become clear that deep, high spatial resolution, multi-wavelength observations of a range of galaxy types are critical to fully investigate the morphological and kinematic evolution of galaxies, and thus inform how galaxies build up their mass over cosmic time. ALMA has been fundamental in such investigations, in particular thanks to detections of the [CII] (157.7 𝜇 m) far-infrared (FIR) emission line (e.g. Stacey et al. 1991; Malhotra et al. 1997, and see review of Hodge & da Cunha 2020). This line is one of the brightest interstellar medium (ISM) cooling lines, tracing the cold ( 𝑇 ∼ 100K) neutral, molecular and ionised ISM (Shibai et al. 1996; Stacey et al. 2010; Pineda et al. 2013; Vallini et al. 2015).  \nThe Reionization Era Bright Emission Line Survey (REBELS, Bouwens et al. 2022) is an ALMA large program (LP) that has carried out spectral scans for the [CII] line in ALMA\'s Band 6 for ∼ 40 UV-bright star-forming galaxies at 𝑧 ∼ 6 -8. As well as multiple, significant [CII] detections, this LP has also enabled simultaneous dust continuum detections of these galaxies (Inami et al. 2022). This sample therefore presents an opportunity to study mass assembly during the Universe\'s last major phase transition; the Epoch of Reionization (EoR). This program has highlighted ALMA\'s capabilities as a \'redshift machine\' (see also Smit et al. 2018; Schouws et al. 2022), as  \nTable 1. Table showing the properties of REBELS-25 from previous studies.  \nReferences: [1] Schouws et al. in prep, [2] Bouwens et al. (2022), [3] Stefanon et al. (in prep), [4] Hygate et al. (2023), [5] Inami et al. (2022), [6] Algera et al. (2023).  \nwell as enabling studies of, for example, dust-obscured star formation (Fudamoto et al. 2021), specific star formation rates (sSFR, Topping et al. 2022), Ly 𝛼 emission (Endsley et al. 2022), and dust and ISM properties (Sommovigo et al. 2022; Dayal et al. 2022; Ferrara et al. 2022) at 𝑧 > 6. Some of these REBELS galaxies show evidence of velocity gradients in the LP data (Schouws et al. in prep), including REBELS-25, which is the most infrared luminous and most promising rotating disc candidate from the REBELS LP (Hygate et al. 2023). However, at the relatively low ( ≳ 1") resolution of the LP observations, morphological and kinematic classification of galaxies remains challenging (e.g., Gonçalves et al. 2010; Simons et al. 2019; Rizzo et al. 2021). Indeed, in Hygate et al. (2023), the angular resolution of the LP observations of REBELS-25 ( ∼ 1 . 3" = 6 . 7 kpc) was found to be insufficient to distinguish between a rotating disc and merger scenario. In addition, analysis of HST rest-frame UV observations by Stefanon et al. (2019) revealed a clumpy UV morphology for this galaxy, which could be indicative of a merger. For these reasons, follow-up high resolution [CII] and dust continuum observations have been obtained for REBELS-25 (ID 2021.1.01603.S, PI J. Hodge) and planned/executed for a subsample of other [CII]-luminous REBELS galaxies (ID 2022.1.01131.S, PI R. Smit).  \nIn this work, we present the morphological and kinematic analysis at sub-kpc resolution for REBELS-25. Analysis of the remaining high resolution REBELS subsample will be the focus of subsequent works (Phillips et al in prep). In Section 2, we discuss the observations and the data reduction used in this high resolution follow-up analysis of REBELS-25. In Section 3, we present the morphological analysis of REBELS-25, and then in Section 4.1 we present the methodology used to investigate its kinematics, with the best-fit models and tests used to distinguish between the merger and rotating disc scenario discussed in Section 4.2. In Section 5, we discuss and interpret our findings and make comparisons to other similar studies, and finally in Section 6 we summarise our main results and their implications. Throughout the paper, we adopt a standard Λ CDM cosmology with Hubble constant 𝐻 0 = 70 km s -1 Mpc -1 , matter density Ω 𝑚 = 0 . 3 and vacuum energy Ω Λ = 0 . 7. At the redshift of REBELS-25, this corresponds to a luminosity distance of 72519 Mpc.', '2 OBSERVATIONS': 'REBELS-25, also known as UVISTA-Y-003 or UVISTA-Y3 (Stefanon et al. 2019; Schouws et al. 2022), has J2000 RA, Dec coordinates of 10 h 00 m 32 . 32 s , + 01 · 44 \' 31 . 3" from rest-frame UV imaging (Stefanon et al. 2019) and a [CII] spectroscopic redshift of 𝑧 = 7 . 3065 ± 0 . 0001 (Bouwens et al. 2022). Analysis of the REBELS LP data of this galaxy has been detailed in Hygate et al. (2023). Here, we present the sub-kpc resolution observations used to further investigate the nature of this source.', '2.1 ALMA Data reduction and imaging': 'The high resolution Band 6 observations of REBELS-25 were obtained in ALMA Cycle 8 with ID 2021.1.01603.S (PI J. Hodge). These observations were carried out in November 2021 in the C-6 configuration, under good weather conditions (typical precipitable water vapour of 0.4 mm). The baseline lengths spanned between 41.4 m and 3.6 km, resulting in a maximum recoverable scale (MRS) of ∼ 1 . 7" and an angular resolution of ∼ 0 . 1". In total, 256.67 minutes were spent on source. The calibration and flagging were done using the standard pipeline running on CASA (Common Astronomy Software Applications for Radio Astronomy, The Casa Team et al. 2022) version 6.2.1.7, and we use the same pipeline version for all subsequent imaging. Inspection of the pipeline-calibrated data tables revealed data of high quality, and the 𝑢𝑣 -data were therefore used without further modification to the calibration scheme or flagging. This Cycle 8 data set contains a spectral window covering the [CII] line with a bandwidth of 1875 MHz, and three other spectral windows for detecting the FIR dust continuum surrounding the line ( ∼ 150 𝜇 m), each with a bandwidth of 2000 MHz.  \nTo obtain a [CII] line-emission cube, we first created a continuumfree measurement set using the task uvcontsub with a zeroth-order fit to line-free channels, excluding a region that is two times the full width at half maximum (FWHM) of the [CII] line. We created a dirty line-emission cube using tclean in the \'cube\' mode to find the root mean square (RMS) noise level per velocity channel and then produced a final line-emission cube cleaned down to two times this RMS noise. This cleaning was done iteratively using CASA \'s automask mode with the sub-parameters initially determined by the recommendations of the ALMA automasking guide 1 for long baselines. On inspection of the mask created for the cleaning and the residuals produced, we lowered the \'noisethreshold\' parameter to 4 and the \'sidelobethreshold\' to 1.5 (both in units of the RMS) to ensure cleaning of source emission in channels not masked in the preliminary clean. The rest of the sub-parameters are left at the recommended values (\'lownoisethreshold\' = 1.5, \'negativethreshold\' = 7, \'minbeamfrac\' = 0.3). The native channel width of the data is ∼ 5.1 km s -1 . In the final line-cube, we bin the channels by a factor of 3 to a channel width of ∼ 15 . 4 km s -1 , as a compromise between sensitivity and velocity resolution.  \nAll imaging was carried out with a multi-scale CLEAN ; a modification of the classical CLEAN algorithm that assumes sources in the sky are extended structures of different scales. We use a small-scale bias of 0.9 and deconvolution scales of 0, 1 × the beam FWHM, and 3 × the beam FWHM. We find that the use of this multi-scale CLEAN , in comparison to the default \'hogbom\' deconvolver, is more successful at masking both compact and extended emission during the cleaning process.  \nA ∼ 150 𝜇 m continuum map was similarly produced from the channels not containing line or atmospheric emission in the \'mfs\' mode. The cleaning was also done iteratively with automasking and cleaned down to two times the RMS noise of the dirty continuum image. We will refer to this continuum map as the dust continuum hereafter.  \nWith the cleaning process described above, we produced both natural-weighted and Briggs-weighted line cubes and continuum maps. The Briggs weighted images were obtained with a robust parameter of 0.5, which is a compromise between sensitivity and spatial resolution. For the [CII] imaging, the natural-weighted beam size is 0.14" × 0.13" (710 × 660 pc), and the Briggs-weighted beam size is 0.12" × "0.10 (610 × 510 pc). For the dust continuum, the natural-weighted beam size is 0.14" × 0.13" (710 × 660 pc), and the Briggs-weighted beam size is 0.11" × 0.09" (560 × 460 pc). Although the Briggs-weighted imaging achieves a higher spatial resolution, the increased sensitivity of the natural weighted data allows for a higher signal-to-noise, particularly in the outer regions of the galaxy, and results in more resolution elements across the source, whilst still providing sub-kpc resolution. For the subsequent analysis, we therefore focus on the naturally weighted imaging. For the [CII] line-cube, the final RMS per 15 km s -1 channel is ∼ 110 𝜇 Jy beam -1 . For the dust continuum map, the final RMS is 𝜎 ∼ 5.7 𝜇 Jy beam -1 . We find consistent [CII] and dust continuum flux densities, within the uncertainties, when comparing the high resolution data with the low resolution REBELS LP data.  \nWhilst it would be possible to increase the sensitivity of these data by concatenating the high resolution data set with the low resolution LP data, the combination of very different array configurations means we would be impacted by the "JvM effect" (see Jorsater & Van Moorsel 1995; Czekala et al. 2021; Posses et al. 2024). Whilst it is possible to correct for this effect (Czekala et al. 2021), this correction means that the residuals have been rescaled, and so the RMS is no longer a true representation of the sensitivity of the observations. The findings of Casassus & Carcamo (2022), in particular, caution against the use of this JvM correction. Since the focus of this paper is a morphological and kinematic analysis, with tools dependent on some SNR masking, these rescaled residuals would impact all subsequent analysis. For these reasons, we choose not to combine the low and high resolution imaging, although we discuss tests using a JvM-corrected, concatenated [CII] line cube (as well as line cubes with different weighting schemes and different channel binnings) in Appendix A. Overall, we find that the main results and conclusions of this paper are robust against the data set and reduction process used.', '2.2 Rest-frame UV imaging': 'In order to compare the [CII] and dust continuum morphology with the rest-frame UV emission, we also make use of the HST Wide Field Camera 3 (WFC3) F160W imaging from the COSMOS-DASH survey (PID 13868, PI D. Kocevski). Details of the HST imaging extracted from the COSMOS-DASH mosaic (Mowla et al. 2019) and the astrometry correction are given in Hygate et al. (2023).', '3 MORPHOLOGY': 'In Figure 1, we present the [CII], dust and rest-frame UV morphology of REBELS-25 at sub-kpc resolution. The [CII] emission reveals a clear extended disc, whereas the rest-frame UV morphology can be  \nseparated into clumps that are offset from the [CII] and dust continuum. Both the [CII] and dust maps show an inner, bright region that is misaligned from the extended gas disc. In the following sections, we detail the fitting tools used to investigate these morphological features.', '3.1 [CII] morphology': 'To investigate the morphology of the [CII] emission in REBELS-25, we use the CASA task immoments to create a moment-0 map by integrating a spectral region covering 2 × FWHM of the [CII] line emission for the naturally weighted line cube (Figure 1). From this moment-0 map, we find a peak SNR of 14.3. We then fit a 2D Sérsic profile to the moment-0 map using the Astropy Sersic2D modelling class and the task TRFLSQFitter with the models convolved with the beam of the observations using PetroFit (Geda et al. 2022). From this fitting, we find a position angle of 224 ± 28 · , an effective radius, 𝑟 e, of 0.42" ± 0.06" (2 . 1 ± 0 . 2 kpc), and a Sérsic index, 𝑛 , of 1 . 3 ± 0 . 2, indicating that the [CII] emission is well-fit by a near-exponential disc. The best-fit ellipticity is 0.13 ± 0.09, which corresponds to an inclination, 𝑖 , of 30 + 10 -15 · , assuming the disc is razor-thin. However, previous studies have found that high𝑧 galaxies are likely to have thicker discs due to increased turbulence. The assumption of a thin disc would therefore result in an underestimation of the inclination, which subsequently impacts the inclination-corrected velocities.  \nWe therefore use the CANNUBI 2 package to also fit for the thickness and inclination of the disc. As introduced in Mancera Piña et al. (2020, 2022a); Fraternali et al. (2021) and Roman-Oliveira et al. (2023), CANNUBI fits the inclination ( 𝑖 ), the disc thickness ( 𝑍 0 ), the morphological centre ( 𝑥 0 , morph , 𝑦 0 , morph ), and morphological position angle ( 𝑃𝐴 morph ) of galactic discs. A more detailed description of CANNUBI is given in Roman-Oliveira et al. (2023), but in summary, we use CANNUBI to fit the observed [CII] flux map with resolutionmatched maps of 3D tilted-ring model disc galaxies. This means that CANNUBI models the galaxy morphology without a prior assumed parametric description of the surface brightness distribution. We set a lower limit of 15 · to the fitting for the inclination, which is a reasonable lower limit based on the inclination derived above from the ellipticity fitted using Sersic2D , assuming a razor-thin disc. The best-fit parameters to the disc geometry are returned via a Markov Chain Monte Carlo (MCMC) routine with 30 walkers and run until convergence, as estimated according to the auto-correlation times of the parameters (Goodman & Weare 2010). The posterior distributions from the CANNUBI fitting to REBELS-25 are shown in Figure B1, and we list the median values and their uncertainties (from the 16th and 84th percentiles) in Table 2.  \nWe find that, overall, the [CII] morphological parameters are well constrained via the CANNUBI fitting procedure. One exception is the fitting for 𝑍 0 , from which we can only infer that the best-fit models have a thickness lower than the ∼ 0 . 14\' (710 pc) resolution of our observations. From tests carried out on mock galaxies in RomanOliveira et al. (2023), we anticipate that the poorly-fitted disc thickness will not significantly impact the derived inclination, since this effect is found to be more important for thicker discs. Having confirmed the thin disk assumption, we therefore set 𝑍 0 = 0 for the subsequent modelling. The posterior for the inclination also appears to have a tail extending up to the 15 · lower limit, however this is likely due to the CANNUBI models appearing very similar for low inclination galaxies (Mancera Piña et al. 2022a). We note that if the inclination is  \nin fact lower, this would only increase the derived intrinsic rotational velocities and further strengthen our findings, as discussed in Section 4.2. The 𝑃𝐴 morph posterior also shows a secondary, lower peak at ∼ 160 · , which is potentially caused by a bright, central, misaligned component which we discuss in Section 3.3.  \nThe morphological centre and 𝑃𝐴 morph returned via the Sérsic and CANNUBI models are consistent within errors. CANNUBI also returns the surface density within each ring, or resolution element. Wefind that this surface brightness profile is best fit with a 1D Sérsic profile with 𝑛 = 1 . 13 ± 0 . 06, which is also consistent with the 2D Sérsic fit. The inclination derived from CANNUBI , 𝑖 = 25 ± 6 · , is lower but within error of the 𝑖 = 56 ± 29 · value determined from the ratio of the major and minor axes of the low resolution LP data in Hygate et al. (2023), where the [CII] emission is only marginally resolved and therefore the axis ratio is strongly dependent on the beam. As mentioned above, the inclination has a significant impact on the derived kinematic properties, highlighting the necessity of high spatial resolution observations for both morphological and kinematic analyses.  \nIn both the CANNUBI and Sersic2D fitting, we see some residual ( ∼ 6 𝜎 RMS in the CANNUBI fitting) emission at the centre of the [CII] map. This feature is discussed in more detail below, in Section 3.3. Despite these residuals, the CANNUBI model produces a better reduced chi squared value ( 𝜒 2 ∼ 0 . 9 for CANNUBI ) than the Sérsic model ( 𝜒 2 ∼ 0 . 18 for Sérsic2D ). Additionally, CANNUBI uses the same software as will be used for our kinematic analysis in Section 4.1. For these reasons, we adopt the morphological parameters for the [CII] emission returned by CANNUBI as our fiducial values, and list these in Table 2.', '3.2 Dust continuum morphology': 'In Figure 1, we also show the dust continuum map of REBELS-25, where we obtain a peak SNR of 10.1. As with the [CII] morphology, we also fit a 2D Sérsic model to the dust continuum and find 𝑟 𝑒 = 0 . 50" ± 0 . 08" (2 . 5 ± 0 . 7 kpc) and 𝑛 = 2 . 2 ± 0 . 4, indicating emission that is more centrally peaked than would be expected for an exponential disc. We find ∼ 3 𝜎 residuals in clumps to the North and South of the centre from the Sersic2D fitting. These clumps are around 0.5" ( ∼ 2 . 5 kpc) from the fitted centre, and each contain around 10% of the flux of the central component.', '3.3 Comparing the [CII], dust and rest-frame UV morphology': 'From the lower resolution REBELS LP data, a potential offset of 0 . 17\' ± 0 . 04\' (around 0 . 9 ± 0 . 2 kpc) between the [CII] emission and the dust continuum is identified in Hygate et al. (2023). With the follow-up, sub-kpc resolution observations analysed here, however, we find that the centres of the dust and [CII] emission are in good agreement, again highlighting the necessity of observations with sufficient angular resolution.  \nThe [CII] morphology is best described by an exponential disc, whereas the dust morphology is not as well fit by a single component Sérsic model, likely due to the two faint clumps to the North and South of the peak of emission. These clumps drive the 𝑃𝐴 morph of the [CII] and dust to differ by ∼ 60 · at the outer extent of emission.  \nWe also find that the effective radii of the [CII] and dust emission are comparable, with the 𝑟 𝑒 of the dust even slightly larger (although they are consistent within the uncertainties). This is in contrast with other studies of high redshift star forming galaxies where the dust emission is found to be around 2 × more compact than the [CII]  \nFigure 1. Left : [CII] moment-0 map from the naturally weighted data cube, with contours showing 2, 3, ...14 𝜎 RMS emission (where 𝜎 RMS = 11 mJy beam -1 km s -1 ). The size of the beam is indicated by the turquoise ellipse in the bottom left corner. Centre : Naturally weighted dust continuum map, with orange contours showing 2, 3, ...10 𝜎 RMS emission (where 𝜎 RMS = 5 . 5 𝜇 Jy beam -1 ). The size of the beam is indicated by the orange ellipse in the bottom left corner. Right : HST WFC3 F160W image from the COSMOS-DASH mosaic (Mowla et al. 2019) with the [CII] emission and dust continuum shown by the turquoise and orange contours, respectively.  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \nFigure 2. Left : We show the best-fit CANNUBI model, with the observed [CII] emission from Figure 1 overplotted in black contours. The dashed magenta line, magenta cross and magenta ellipse mark the morphological position angle, the centre and the radial extent of the best-fit model, respectively. Right : The residuals from the best-fit CANNUBI model, with 2, 3, ..., 6 𝜎 RMS contours shown by the solid black lines, and -2, -3 𝜎 RMS contours show by the dashed black lines, where 𝜎 RMS is the local RMS in the [CII] moment-0 map.  \n<!-- image -->  \nemission (e.g., Fudamoto et al. 2022; Fujimoto et al. 2020). However, we note that the Sérsic model is a poor fit to the dust emission, and this large 𝑟 𝑒 is likely driven by the clumpy structures to the North and South, and we therefore caution against comparing the effective radii of these Sérsic fits.  \nNotably, both the [CII] and dust show a central region in the inner ∼ 1 kpc (around 1.4-2 beams wide) that appears misaligned with the 𝑃𝐴 morph of the extended gas disc. We see also see this misaligned, central region in the residuals around the fitted centre of the galaxy in both the [CII] and dust continuum 2D fits, as well as a tail in the posterior distribution for PA morph with CANNUBI . To illustrate this, we fit elliptical isophotes to the [CII] moment-0 and dust map using the \'isophote.Ellipse.fit\\_image\' in PHOTUTILS from the Astropy package (Bradley et al. 2020) with a sigma clip of 2 𝜎 RMS . For the [CII] fit, we follow the methodology outlined in Amvrosiadis et al. (2024), that is, we first leave all parameters free, and then for the final iteration we fix the centre of all ellipses to the median value. However, the fit for the dust fails to converge when the centre is fixed, and we therefore  \nFigure 3. Top left : The elliptical isophote fit to the [CII] emission is shown by the black contours. Top right : The elliptical isophote fit to the dust continuum map. Bottom : The position angle of the fitted ellipses as a function of radius for the [CII] (turquoise) and dust (orange) emission.  \n<!-- image -->  \nleave all parameters free for the dust fitting. The resulting fits are shown in Figure 3.  \nWithin the central 1 kpc, the ellipses fitted to the [CII] and dust have an average PA of ∼ 120 · . At roughly the same radii ( ∼ 1 kpc), the PA for the [CII] suddenly increases to an average PA of ∼ 220 · at ≳ 1 kpc (consistent with 𝑃𝐴 morph from CANNUBI and Sersic2D ), and for the dust it changes more gradually to ∼ 190 · (consistent with 𝑃𝐴 morph from Sersic2D for the dust). In addition, the ellipticity in the central 1 kpc in both the [CII] and dust is found to be ∼ 0 . 3, whereas for the extended [CII] disc this is found to be ∼ 0 . 1 (nearcircular).  \nSimilar changes in both the ellipticity and position angle from elliptical isophote fits have been used to identify and analyse bars within nearby galaxies (e.g., Abraham et al. 1999; Laine et al. 2002; Jogee et al. 2004; Menendez-Delmestre et al. 2007; Consolandi 2016), and now even for some higher redshift galaxies with JWST (e.g., Huang et al. 2023 at 𝑧 = 2 . 467, Amvrosiadis et al. 2024 at 𝑧 = 3 . 762, Le Conte et al. 2024 at 1 ≤ 𝑧 ≤ 3). However, we note that this tentative bar-like central component in REBELS-25 is barely resolved with these observations.  \nThe rest-frame UV morphology differs significantly from the [CII] anddust. As previously identified in Hygate et al. (2023) and Schouws et al. (2022), the HST imaging reveals 3-4 clumps that are ∼ 0 . 6" offset from the centre of both the [CII] and dust continuum emission. We see in Figure 1 that these UV clumps lie close to and outside the outer edge of the [CII] and dust emission. As discussed in Hygate et al. (2023), these UV clumps could have a number of explanations, including differential dust obscuration, regions of intense star formation or potential merger activity. We further discuss these possible interpretations in Section 5.1.', '4 KINEMATICS': 'We next focus on the [CII] kinematics of REBELS-25. In Section 4.1, we first provide an overview of the techniques used to investigate the kinematic properties. We then present the results of our kinematic modelling in Section 4.2, along with some investigation of the potential mechanisms that could explain the observed properties.', '4.1 Methods': 'For the kinematic modelling of REBELS-25, we use 3D BAROLO (Di Teodoro & Fraternali 2015, v1.7). 3D BAROLO is well-tested at a range of redshifts and from low to high resolution observations (e.g. at low resolution: Di Teodoro & Fraternali 2015; Mancera Piña et al. 2019; Gray et al. 2023, at mid resolution: Bacchini et al. 2020; Mancera Piña et al. 2022b; Roman-Oliveira et al. 2023, and at high resolution: Iorio et al. 2016; Di Teodoro & Peek 2021). Its algorithm fits 3D tilted ring models to emission line data cubes on a channelby-channel basis, and uses Monte Carlo methods to return the best-fit kinematic and/or morphological parameters. The tilted ring models are convolved to the same resolution as the input observations, and it therefore accounts for the effects of beam smearing.  \nBefore the modelling, we extract a 3.5" × 3.5" sub-cube centred on the [CII] emission that covers ± 4 × the FWHM of the [CII] emission in the spectral dimension from our full data cube. We then use the SEARCH algorithm to carry out a source scan on the inputted sub cube, with a primary SNR cut of 3 𝜎 RMS . To grow a mask around the detected source, a range of GROWTHCUT parameters were tested. Based on these preliminary tests, we find a threshold of 2.2 × 10 -4 Jy beam -1 , which is is equivalent to 2 𝜎 RMS , is sufficient to produce a mask that encloses all of the [CII] emission without introducing too much noise. We set the radial separation of the rings to 0.11" ( ∼ 80% of the beam FWHM), resulting in five rings within the mask produced by 3D BAROLO . This maximises the number of rings whilst ensuring each ring is nearly comparable to an independent resolution element. For all fits with 3D BAROLO , we use an azimuthal normalisation of the surface density, we fit the kinematics using both sides of the rotation curve, and we minimise the reduced chi squared statistic.  \nWith 3D BAROLO , there are nine possible parameters to fit: the rotation velocity ( 𝑉 rot), velocity dispersion ( 𝜎 ), radial velocity ( 𝑉 rad ), systemic velocity ( 𝑉 sys), kinematic position angle ( 𝑃𝐴 kin ), kinematic  \ncentre ( 𝑥 0 , kin , 𝑦 0 , kin ), 𝑖 , and 𝑍 0 . To limit the number of free parameters, we leave only 𝑉 rot, 𝜎 , and 𝑃𝐴 kin of each ring free. For the remaining parameters, we fix the centre and inclination to the values returned by CANNUBI , and we fix 𝑍 0 to zero following our finding that the thickness of the disc is unresolved ( < 710 pc). We also fix 𝑉 rad to zero, meaning we assume no radial motion, following preliminary tests described below. For the systemic velocity, ( 𝑉 sys), we allow 3D BAROLO to estimate this from the data cube, which results in 𝑉 sys = 35 . 72 km s -1 from the [CII] line at a redshift 𝑧 = 7 . 3065, which is consistent with 𝑉 sys derived from a single component Gaussian fit to the [CII] spectrum and consistent with the reported uncertainty on the redshift from Table 1. The results from the 3D BAROLO fitting are given in Table 2 and shown in Figure 4, 5, and 6. We note that, whilst the 𝑃𝐴 kin was allowed to vary for each ring, the best-fit model has a constant 𝑃𝐴 kin across the disc.  \nTo further test the robustness of our modelling, we also run 3D BAROLO with all parameters left free. The resulting inclination ( 𝑖 = 32 ± 5) is consistent within the uncertainties with the 𝑖 fitted with CANNUBI , as is the fitted centre (within one beam FWHM of 𝑥 0 , morph , 𝑦 0 , morph ). The fitted values for 𝑉 rad are consistent with zero, within the uncertainties, for all five rings. The fitted 𝑍 0 is also found to be consistent with zero. Whilst we adopt as our fiducial model the fit described above with only three free parameters ( 𝑉 rot, 𝜎 and 𝑃𝐴 kin ), we note that the overall results of this paper are consistent when the fit with all parameters left free is instead used.', '4.2 Results': 'We show the final velocity maps from our fitting with 3D BAROLO in Figure 4 and the position velocity diagrams (PVDs) in Figure 5. We see a clear velocity gradient in the velocity field maps of Figure 4 and a near-S-shape curve in the major axis PVD in Figure 5, which are characteristic features of rotating discs.  \nIn Table 2, we report the final kinematic parameters derived by 3D BAROLO , including the maximum rotation velocity ( 𝑉 rot , max), average velocity dispersion ( ¯ 𝜎 ) and the 𝑉 rot , max / ¯ 𝜎 ratio for REBELS-25. The ¯ 𝜎 gives a measure of the overall disc turbulence, and 𝑉 rot , max / ¯ 𝜎 measures the ratio of ordered to turbulent motion. For this analysis, we assume that the 𝑉 rot , max value is comparable to the rotational velocity at which the velocity curve flattens, since we see some flattening beyond 2 kpc in the rotation curve and major-axis PVD.  \nWe also plot the 𝑉 rot and 𝜎 values of each ring as a function of radius in Figure 6. We note that the rotation velocities reported are not corrected for pressure support, i.e., we have assumed that the rotation velocity we measure is an optimal tracer of the circular velocity, 𝑉 circ . This is a reasonable assumption for REBELS-25 thanks to the high 𝑉 rot / 𝜎 values found, meaning that the pressure term in the Virial theorem is likely negligible (see e.g., Equations 6 and 7 in Iorio et al. 2016). We therefore take the 𝑉 rot , max value as the maximum circular velocity of the system in subsequent calculations.  \nThe rotation velocity curve of REBELS-25 in Figure 6 shows the velocity rising quickly in the centre and then flattening to the maximum value of ∼ 370 km s -1 at around 2 kpc. Similar rotation curves are observed in disk galaxies at lower redshifts (e.g., Noordermeer et al. 2007; Lelli et al. 2016, 2023; Di Teodoro et al. 2016). We also find a potential \'bump\' feature in the rotational velocity curve within the first ∼ 1 kpc. Whilst the uncertainties are such that this feature is not certain, this bump persists also in the fits with 3D BAROLO where all the parameters are left free, and could indicate a potential bulge/bar component, as found in the Milky Way (e.g. Portail et al. 2017), local galaxies (e.g. Lelli et al. 2016), and in some recent stud-  \nFigure 4. 3D BAROLO fitting for REBELS-25. Emission is masked at 2 𝜎 RMS. The first column on the left shows the observed data, the middle column the model and the column on the right shows the residuals. The first row is for the intensity map, the second row for the velocity field map and the bottom row for the velocity dispersion map. In the first row, the black cross, ellipse and dashed line show the centre, radial extent and position angle of the 3D BAROLO model, respectively. In the second row, the grey dots give an indication of the separation of each ring along the velocity field (0.11"). In the velocity field map of the data and model, we also plot the iso-contours from -180 to 180 km s -1 in 45 km s -1 increments. In all maps, the beam size is indicated by the turquoise ellipse in the bottom left corner.  \n<!-- image -->  \nFigure 5. Major- (upper) and minor- (lower) axis position-velocity diagrams (PVDs) of REBELS-25, with black contours showing 2 𝜎 RMS and 4 𝜎 RMS emission, and the dashed grey contours -2 𝜎 RMS, where 𝜎 RMS = 168 𝜇 Jy/beam and is equal to one standard deviation above the median value as calculated by 3D BAROLO . The red contours illustrate the best-fit model from 3D BAROLO . In the major-axis PVD, the yellow markers show the rotation velocities of each ring projected along the line of sight.  \n<!-- image -->  \nies of galaxies at 𝑧 > 4 (e.g. Lelli et al. 2021; Rizzo et al. 2021; Lelli et al. 2023; Tripodi et al. 2023; Roman-Oliveira et al. 2024).  \nA decrease of 𝜎 with radius is also observed, as found in some local disk galaxies (e.g., Bacchini et al. 2020; Fraternali et al. 2002; Boomsma et al. 2008; Iorio et al. 2016). This decrease is usually ascribed to the radial change of the mechanisms sustaining the turbulence within galaxies (e.g., supernova feedback, accretion and gravitational instabilities).', '4.2.1 Non-circular features in the 3D BAROLO models': 'Before drawing further conclusions based on the above modelling, we first investigate any features that are not well fit by 3D BAROLO , which could indicate that the observed velocity gradient might be caused by outflows, interactions, or mergers, rather than solely the presence of a rotating disc (as discussed in e.g., Loiacono et al. 2019; Simons et al. 2019). Firstly, we see some evidence of twisted iso-velocity contours in the moment-1 map of REBELS-25 that are not well-fit by the rotating disc model with a constant 𝑃𝐴 kin and 𝑖 . These features could indicate a potential warped disc or distortions due to radial motions (e.g., Di Teodoro & Peek 2021). In some cases, twisted iso-velocity contours are also attributed to bars or other nonaxisymmetric structures within galaxies (e.g., Wu et al. 2021; Buta 1987; de Naray et al. 2009). Secondly, in the velocity field map and the residual map of Figure 4 we can see a region on the approaching side (to the West) of the galaxy with high velocities which is not reproduced within the model. This feature can also be seen in the  \nchannel maps at -228 and -274 km s -1 (Figure C1). This could have a number of causes, such as recently accreted gas, streaming motions along spiral arms, or due to a warped gas disc. These features may also be consistent with noise in the data. In an attempt to investigate these features further, we tested 3D BAROLO fits that incorporate radial motions. As mentioned in Section 4.1, we find 𝑉 rad consistent with zero given the uncertainties, and we also find that these models were not able to reproduce similar high velocity regions.  \nAdditionally, we see some faint emission beyond the rotating disc model in the PVDs. For example, in the lower half of the minoraxis PVD, and also in the lower half of the major-axis PVD, we see some faint high receding velocity regions. As with the features described above, these could signify non-circular motions, such as vertical/radial motions due to inflows, outflows or mergers, the presence of a bar or the presence of spiral arms. Similar deviations due to non-circular motions are common in rotating galaxies at 𝑧 = 0 (e.g., Fraternali et al. 2001; Jorsater & Van Moorsel 1995; Zurita et al. 2004; Trachternach et al. 2008; Erroz-Ferrer et al. 2015), and similar faint emission outside ideal rotating disc models can also be seen in mock data of simulated rotating discs (e.g., see Figure 6 of Rizzo et al. 2022). Distinguishing between the different processes potentially responsible for these features is not feasible at the current resolution and SNR, and we note that this emission is only detected at ∼ 2 𝜎 RMS .  \nOverall, we find that the majority of the emission is well described by the rotating disc model, although we cannot exclude some disturbed kinematics due to, for example, outflows or a warped disc. Follow-up data on the stellar morphology would be needed to confirm the presence of any non-axisymmetric features that could also be introducing non-circular motions.', '4.2.2 Secondary [CII] component': 'In Hygate et al. (2023), a potential secondary [CII] component is identified in the range ∼ + 250 to + 650 km s -1 in the [CII] spectrum of REBELS-25 from the lower resolution LP data. Within Hygate et al. (2023), it is hypothesised that this spectral component could be the result of an outflow (where it is modelled as a broad component with width ∼ 380 km s -1 ) or a minor merger (modelled as a narrow component with width ∼ 100 km s -1 ). We have searched for this secondary [CII] component in our data set following a similar methodology as in Hygate et al. (2023), whereby we extracted a spectrum for REBELS-25 using a circular aperture of 1.75" and attempted to fit a main and secondary [CII] component, where we limit the centre of the secondary peak to the same limits used in Hygate et al. (2023) for both the outflow and minor merger model. However, we find no convincing secondary component in the [CII] spectrum from the high resolution data. We also extracted spectra from beam sized apertures across the source in case the component has been diluted in the global spectrum, as found for a different galaxy by Herrera-Camus et al. (2021), but again find no significant secondary component. Indeed, Hygate et al. (2023) suggested the component might instead be extended on ∼ 9 kpc scales. This would make the large-scale outflow scenario the more likely option. Investigating the nature of this potential component further would therefore require deeper data with increased surface brightness sensitivity to such extended structure.', '4.2.3 Merger or rotating disc?': "In some cases at low spatial resolution, the observed velocity map of a merging system can be similar to that of a smooth rotating disk, as  \nTable 2. Table showing the [CII] morphological and kinematic parameters derived from a 2D Sérsic model, CANNUBI and 3D BAROLO , as described in the text.  \nFigure 6. The rotation velocity ( 𝑉 rot, dark blue) and velocity dispersion ( 𝜎 , red) curves for REBELS-25. The size of the band indicates the estimated uncertainty from the 3D BAROLO fit for 𝑉 rot and 𝜎 . The dashed black line indicates the effective radius from Sersic2D fitting.  \n<!-- image -->  \nobservations can smooth out the irregularities and asymmetries of a merger (e.g., Simons et al. 2019; Kohandel et al. 2020; Rizzo et al. 2022). Within the literature, there exist multiple methods and criteria used for differentiating between rotating discs and mergers. Here we describe the results from applying the rotating disc criteria laid out in Wisnioski et al. (2015), and then adapted for [CII] ALMA-ALPINE data in Jones et al. (2021), and also the PV Split method described in Rizzo et al. (2022).  \nWe use the same five rotating disc criteria listed in Jones et al. (2021) with the 2D intensity, rotational velocity and velocity dispersion maps outputted from 3D BAROLO . We also make use of the morphological parameters fitted with CANNUBI , and we use 3D BAROLO 's 2D fit capabilities to an inputted velocity field for the third and fifth criteria.  \nIn Jones et al. (2021), typically at least 3 of these criteria are met for the ALPINE galaxies classified as rotating discs. For REBELS-25, we find that all 5 criteria are met:  \n- (i) We find that a slice along the kinematic major axis reveals a  \nslope in position-velocity space with > 3 𝜎 significance, implying a smooth velocity gradient.  \n- (ii) For each of the five fitted rings, we find 𝑉 rot > 𝜎 (see Figure 6).\n- (iii) We find that the kinematic centre derived by a 2D fit with 3D BAROLO and the peak in the velocity dispersion map are consistent within the uncertainties.\n- (iv) The 𝑃𝐴 morph from CANNUBI and the average 𝑃𝐴 kin across all rings are consistent within the uncertainties, as expected for an ideal rotating disc.\n- (v) We find that the kinematic centre fitted by 3D BAROLO and the morphological centre fitted by CANNUBI are consistent within error.  \nFrom these tests, we may conclude that REBELS-25 is a rotating disc. However, Rizzo et al. (2022) find that these criteria can mistakenly classify galaxy kinematics, and we therefore also use the PV Split method described in the same paper.  \nIn Rizzo et al. (2022), mock ALMA data cubes of SERRA simulated rotating discs, disturbed discs, and merging galaxies (Pallottini et al. 2022) were used to show that the morphology and symmetry of the major and minor axis PVDs can be used to kinematically classify galaxies, with rotating discs typically showing a symmetric, S-shaped position-velocity profile. Three parameters are defined to measure the morphology and symmetry of the PVDs, and a plane that divides these two kinematic classes in the 3D PV Split parameter space is defined and tested using local disc galaxies and mergers from the WHISP survey (van der Hulst et al. 2001) in Roman-Oliveira et al. (2023). Using the PVDs returned by 3D BAROLO for REBELS-25, we find that it lies in the rotating disc parameter space of the 3D PV Split diagram. However, we again cannot exclude a disturbed disc or minor merger scenario with this method.  \nAs described above, we have found some evidence of non-circular features and therefore cannot rule out the existence of a minor merging component or outflows/inflows with the current data. However, following visual inspection of the 3D BAROLO fits, tests using the rotating disc criteria of Wisnioski et al. (2015), and tests using the PV Split method of Rizzo et al. (2022), we find that REBELS-25 is best described as a rotation-dominated disc galaxy.", '4.2.4 Dynamical mass and gas mass': 'For a system in virial equilibrium, the dynamical mass, 𝑀 dyn , enclosed within a radius, 𝑅 , is given by:  \n𝑀 dyn ( < 𝑅 ) = 𝑘 ( 𝑅 ) 𝑅𝑣 2 circ ( 𝑅 ) 𝐺 , (1)  \nwhere 𝑘 ( 𝑅 ) is the virial coefficient at radius 𝑅 . Price et al. (2022) use Sérsic models to derive the virial coefficients for a range of 𝑛 and intrinsic axis ratios ( 𝑞 ). As we have some evidence that REBELS-25 is a thin exponential disc from CANNUBI and 2D Sérsic fits, we adopt a total virial coefficient (which relates the total dynamical mass of the system to the circular velocity) of 𝑘 tot = 1 . 8, which is the virial coefficient for 𝑛 = 1 and 𝑞 = 0 . 2 (typically quoted for thin discs, e.g. van der Wel et al. 2014). We take the 𝑉 rot in the 4th ring ( = 366 + 30 -31 km s -1 ), since this is the closest to the 𝑟 𝑒 as derived from the 2D Sérsic fitting. This results in 𝑀 dyn , tot = ( 1 . 2 ± 0 . 3 ) × 10 11 M ⊙ . This is larger but within the uncertainties of the value derived in Hygate et al. (2023), which used the low resolution LP data.  \nAssuming that the contribution of dark matter to the mass budget within 𝑟 𝑒 is negligible (as found at high𝑧 in e.g., van Dokkum et al. 2015; Price et al. 2016; Wuyts et al. 2016; Genzel et al. 2017; Gelli  \net al. 2020), we can estimate the total gas mass from 𝑀 gas = 𝑀 dyn -𝑀 ∗ . Adopting a stellar mass of 𝑀 ∗ = 8 + 4 -2 × 10 9 M ⊙ results in 𝑀 gas = ( 1 . 1 + 0 . 4 -0 . 3 )× 10 11 M ⊙ , with a gas fraction of 𝑓 gas = 𝑀 gas , tot /( 𝑀 gas , tot + 𝑀 ∗ ) = 0 . 93 + 0 . 03 -0 . 07 . However, we note that if we adopt the stellar mass derived in Topping et al. (2022) using a non-parametric star formation history ( 𝑀 ∗ = 19 + 5 -8 × 10 9 M ⊙ ), 𝑀 gas , tot = ( 9 . 8 + 4 . 3 -3 . 5 ) × 10 10 M ⊙ and 𝑓 gas = 0 . 84 0 . 09 -0 . 11 .  \nWith these constraints on the gas mass, we can also derive a conversion factor between the [CII] luminosity and 𝑀 gas, 𝛼 [ CII ] , assuming that the [CII] emission is a good tracer of the total gas mass (Swinbank et al. 2012; Gullberg et al. 2018; Zanella et al. 2018). With a double Gaussian fit to the [CII] spectrum, we find 𝐿 [ CII ] = ( 1 . 8 ± 0 . 2 ) × 10 9 𝐿 ⊙ , which is consistent within the uncertainties with the value derived from the low resolution LP data in Hygate et al. (2023). This results in 𝛼 [ CII ] = 62 + 32 -22 M ⊙ / L ⊙ (or 55 + 35 -23 if the Topping et al. 2022 𝑀 ∗ is used). These values are a factor of ∼ 10 times higher than the median value found for dusty star-forming galaxies at 𝑧 4 . 5 in Rizzo et al. (2021), and a factor of ∼ 2 times higher than the value of 30 M ⊙ / L ⊙ found for main sequence galaxies in Zanella et al. (2018) (although we note that this conversion is derived for the molecular gas mass, and not total gas mass). However, we note that the uncertainties are high, and a more statistically significant sample would be necessary to investigate 𝛼 [ CII ] values for star forming galaxies at 𝑧 > 6.  \nFrom the properties derived by 3D BAROLO , and an estimate for the gas mass, it is also possible to determine the Toomre parameter (Toomre 1964), 𝑄 , of REBELS-25. This 𝑄 parameter is used to measure the local instability of a disc, with 𝑄 < 1 indicating that it is susceptible to local gravitational perturbations (LGIs) (e.g., Leung et al. 2020). Using Equation (1) from Neeleman et al. (2020), we find an average 𝑄 of 0 . 4 + 0 . 5 -0 . 2 over the radial extent of the [CII] emission, suggesting it may be locally unstable. However, some studies at 𝑧 ∼ 0 have found that 𝑄 does not show any correlation with fragmentation or star formation (e.g Leroy et al. 2008; Romeo & Falstad 2013; Elmegreen & Hunter 2015; Romeo et al. 2023), and a recent study also indicates that high𝑧 galaxies are not as locally unstable as expected from the 𝑄 parameter (Bacchini et al. 2024).', '5 DISCUSSION': "Through our kinematic analysis of REBELS-25, we have found a dynamically cold disc with a high degree of rotational support ( 𝑉 rot , max / ¯ 𝜎 = 11 8 -4 ) and low overall turbulence ( ¯ 𝜎 = 33 ± 9 km s -1 ). In Kohandel et al. (2023), galaxy discs are classed according to their 𝑉 / 𝜎 ratios, with 𝑉 / 𝜎 > 10 being defined as 'super cold' discs. Following this classification, we find that REBELS-25 is the most distant super cold disc galaxy observed to date. In the following, we first discuss our findings on this galaxy in more detail (Section 5.1) before placing REBELS-25 in context with other high resolution gas kinematic studies.", '5.1 An evolved, dynamically cold rotating disc at 𝑧 = 7 . 31': "In Section 3.1, we found that the [CII] emission of REBELS-25 is well-fit by an exponential disc model, similar to what is found for galaxies in the local Universe and for the dust continuum in some submillimetre galaxies at 𝑧 ∼ 1 -3 (e.g., Hodge et al. 2016; Barro et al. 2016; Rivera et al. 2018; Fujimoto et al. 2018; Gullberg et al. 2019). This is in contrast to the clumpy and irregular morphologies observed and anticipated at 𝑧 > 6 (see Section 1). We also find  \nthat the gas disc of REBELS-25 is likely thinner than the 𝑍 0 ∼ 1 kpc derived from galaxies at 𝑧 ∼ 4 . 5 in Roman-Oliveira et al. (2023). This contradicts the general expectation that galaxies at higher redshifts will tend to have thicker discs due to increased turbulence, as found in Λ CDMgalaxy formation simulations (e.g., Brook et al. 2012; Bird et al. 2013; Park et al. 2019; Bird et al. 2021; Renaud et al. 2021).  \nWe also see some tentative evidence that REBELS-25 is a barred disc (Section 3.3), with a sudden change in both the position angle and ellipticity of the [CII] and dust continuum emission at ∼ 1 kpc. With the current data, it is not feasible to search for kinematic signatures of a bar, and we therefore cannot rule out alternative scenarios, such as inflows/outflows or a late-stage minor merger. However, a bar may be feasible since barred structures are thought to form relatively quickly (order of a hundred million years) in massive dynamically cold disc galaxies (e.g., Hohl 1971; Sellwood & Wilkinson 1993; Athanassoula 2002; Rosas-Guevara et al. 2022; Bland-Hawthorn et al. 2023). Barred galaxies are common in the nearby Universe (e.g., Knapen et al. 2000; Whyte et al. 2002; Laurikainen et al. 2004; Buta et al. 2015), but the fraction of barred galaxies has been found to decrease between 𝑧 = 0 to 1 (e.g., Sheth et al. 2008). However, some recent observations have identified barred galaxies at 𝑧 = 1 -4 . 4 (Simmons et al. 2014; Hodge et al. 2019; Tsukui et al. 2023; Huang et al. 2023; Smail et al. 2023; Guo et al. 2023; Amvrosiadis et al. 2024; Le Conte et al. 2024), further indicating that such structures can form in galaxies within the first few Gyrs after the Big Bang.  \nWe additionally see some very tentative evidence of arm-like or clump-like extended features at ∼ 3 kpc to the North and South of the peak of emission in the dust continuum map (Figure 1). Similar features have been observed in the dust emission of SMGs (e.g., Gullberg et al. 2019; Hodge et al. 2019). However, these features are very faint. Three approved JWST programs (GO-1626, GO-6036 and GO-6480) targeting this galaxy and providing both low resolution spectroscopy and high spatial resolution imaging with grism spectroscopy of the stellar and ionised gas content may shed more light on whether or not any bar, spiral arm and/or clump structures are present.  \nSome clumpy structure has already been identified in REBELS-25 from the rest-frame UV imaging with HST. As discussed in Hygate et al. (2023), these UV clumps could have a number of explanations, including differential dust obscuration, regions of intense star formation or potential merger activity. From the analysis described thus far, the high resolution [CII] observations disfavour the hypothesis that these clumps are the result of an ongoing merger, since we see no strong kinematic evidence of merging activity. Our analysis instead indicates that these clumps may be due to differential dust obscuration, particularly at the centre of the disc where the dust emission peaks, or due to star forming clumps. Similar clumpy morphology in the rest-frame UV embedded within a rotating gas disc has also been found in, for example, a 𝑧 = 6 . 072 galaxy (Fujimoto et al. 2024a).  \nWe also find that REBELS-25's dynamical mass is a factor of ∼ 14 times greater than the adopted stellar mass from Table 1, where this stellar mass is obtained with the BEAGLE SED modelling code assuming a constant star formation history (SFH) and a Chabrier (2003) IMF, with full details described in Stefanon et al. (in prep). If we instead adopt the stellar mass derived in Topping et al. (2022)( 𝑀 ∗ = 19 + 5 -8 × 10 9 / 𝑀 ⊙ ), 𝑀 dyn / 𝑀 ∗ becomes ∼ 6. de Graaff et al. (2023) have recently reported 𝑀 dyn / 𝑀 ∗ values as high as 30 in high𝑧 galaxies, whereas studies of 𝑧 ∼ 2 galaxies find 𝑀 dyn / 𝑀 ∗ ∼ 4 (Maseda et al. 2013). For the case of REBELS-25, we have seen that the gas fraction is likely very high, however we cannot yet rule out contributions from dark matter components to the mass budget.", '5.2 Selecting a comparison sample': 'To place REBELS-25 into context with the evolution of the gas kinematic properties of galaxies across cosmic time, a comparison sample across a broad redshift range with comparable observations is necessary. One important point to consider when comparing kinematic studies at a range of redshifts is the ISM tracer used. At 𝑧 ∼ 0, gas kinematics of galaxies are typically traced by warm ionised gas tracers, such as H 𝛼 . These warm ionised gas tracers have been found to result in higher velocity dispersion compared to values derived from atomic or molecular gas (e.g., Levy et al. 2018; Übler et al. 2019; Girard et al. 2021; Kohandel et al. 2023). We therefore mainly compare with studies using cold gas tracers, and we discuss the implications of any comparisons made with different tracers in subsequent sections.  \nAsdiscussed in Section 1, the [CII] 157.7 𝜇 memissionline is a useful, strong-line tracer of cold gas kinematics. However, observations of this line become increasingly difficult at 𝑧 ≲ 3 . 5 from groundbased instruments due to the poor transmission through Earth\'s atmosphere. For this reason, we compile literature results of both high resolution [CII] kinematics studies, typically at 𝑧 > 3 . 5, and high resolution [CI] and CO kinematics studies, typically at 0 . 5 ≲ 𝑧 ≲ 3 . 5. For this selection from the literature, we take into consideration the criteria recommended in Rizzo et al. (2021) for distinguishing between mergers and rotating discs. We therefore select only literature studies that use observations with ≳ 3 independent resolution elements across the semi-major axis based on the reported radial extent of emission of the galaxies and the beam size of their corresponding observations. Our final literature compilation consists only of galaxies that are classified as rotating discs based on these high-quality kinematic data. We discuss the literature sample in more detail below.  \nRizzo et al. (2023) presented the ALMA-ALPAKA survey; a collection of 28 star forming galaxies with high resolution, high SNR ALMA archival observations of CO and [CI] emission at 𝑧 = 0 . 5 -3 . 5. The physical resolution of these sources varies from 1 to 4 kpc. Whilst sub-kpc resolution observations are preferable to robustly analyse kinematics, this resolution is found to be sufficient to sample the velocity curves with ≳ 3 resolution elements for 13 discs (classified as rotating discs using the PV Split method) within the sample. Of these 13 discs within the ALMA-ALPAKA survey, seven are classified as starbursts, and six as main sequence galaxies (the remaining galaxy has no constraint on 𝑀 ∗ or SFR). Additionally, five of these are classified as AGN galaxies or AGN candidates by Rizzo et al. (2023). Many of these 13 discs were also the focus of previous kinematic studies, as described in Rizzo et al. (2023).  \nParlanti et al. (2023b) collected and analysed all available Band 4, 5, and 6 ALMA archival observations of [CII] (and [OIII]88 𝜇 m) emission in galaxies at 𝑧 > 3 . 5 with angular resolution < 1.5" and SNR > 7. From this sample, we initially select seven galaxies that meet our criteria. Of these seven galaxies, five were studied using the same data set in other works. These five galaxies are DLA0817g, ALESS 073.1, HZ4, HZ7, and COS-29, analysed previously by Neeleman et al. (2020) (and Roman-Oliveira et al. 2023), Lelli et al. (2021); Herrera-Camus et al. (2022); Lambert et al. (2022), and Posses et al. 2023, respectively. Lelli et al. (2021), Posses et al. (2023), and Roman-Oliveira et al. (2023) use 3D BAROLO for their kinematic analyses, and Neeleman et al. (2020) use QUBEFIT but find that 3D BAROLO produces consistent results. To be more comparable to our analysis of REBELS-25 and the kinematic analyses of the ALMA-ALPAKA galaxies in Rizzo et al. (2023), who also use 3D BAROLO , we therefore use the kinematic properties derived from Roman-Oliveira et al. (2024), Lelli et al. (2021), and Posses et al. (2023) for DLA0817g, ALESS 073.1, and COS-29, respectively. The  \ngalaxy HZ7 is found to be a potential merger by the kinematic analysis of Lambert et al. (2022), and we therefore exclude it from our selection. In Herrera-Camus et al. (2022), the 3D parametric code DYSMAL is used to model the kinematics of HZ4. The resulting kinematic properties derived in Herrera-Camus et al. (2022) and Parlanti et al. (2023b) are consistent for HZ4, and therefore for HZ4 and the remaining two galaxies (HZ9, J1211, and VR7) selected from Parlanti et al. (2023b), we use the kinematic properties derived therein.  \nTo increase the size of our comparison sample, we also consider starburst, quasar host, sub-millimetre, and gravitationally lensed galaxies that do not show signs of merging activity from RomanOliveira et al. (2023); Rizzo et al. (2020); Liu et al. (2023); Fraternali et al. (2021); Hodge et al. (2012); Fujimoto et al. (2024b); Neeleman et al. (2023); Lelli et al. (2023); and Amvrosiadis et al. (2024). There exist additional studies of cold gas kinematics (including HI kinematics) with observations of sufficient angular resolution and sensitivity, particularly at lower redshifts, however we choose to focus on higher redshifts ( 𝑧 > 0 . 5) where observations are more comparable to those of REBELS-25. Additional galaxies considered but not included in the comparison sample include AZTEC-1, a sub-millimetre galaxy at 𝑧 = 4 . 6 with high resolution [CII] ALMA observations initially analysed by Tadaki et al. (2018), who found a rotating disc with 𝑉 rot , max ∼ 227 km s -1 and ¯ 𝜎 ∼ 74 km s -1 . However, when the same dataset was analysed in Roman-Oliveira et al. (2023), they found evidence of merging activity according to the disc criteria of Wisnioski et al. (2015), the PV Split method, and from a visual inspection of the velocity field and PVDs. Roman-Oliveira et al. (2023) therefore conclude that AZTEC-1 is likely a merger. Other notable high redshift cold gas kinematic studies include Rivera et al. (2018); Rybak et al. (2019); Tamura et al. (2023); Pope et al. (2023); Smit et al. (2018), and Neeleman et al. (2021). However, these studies do not match our observation selection criteria, increasing the uncertainty in the kinematic properties derived when compared to the analysis of REBELS-25.  \nOverall, our selection results in a total of 37 objects with [CII], [CI], or CO high resolution kinematic studies from the literature. We plot the stellar mass and SFR for the selected galaxies in Figure 7 at three different redshift ranges (0 . 5 ≤ 𝑧 < 3, 3 ≤ 𝑧 < 4 . 5, 4 . 5 ≤ 𝑧 < 7). Due to the limited sensitivities and resolution of observations at high redshift, we see that our selection is biased towards massive ( 𝑀 ∗ ≳ 10 9 M ⊙ ) main-sequence or starburst (SFR ≳ 10M ⊙ yr -1 ) galaxies. We note that, as discussed in Hygate et al. (2023), the definition of a main-sequence and a starburst galaxy differs in the literature. In Figure 7, we adopt the Schreiber et al. (2015) main sequence relation, and define a starburst threshold at 4 times above this line (e.g., Rodighiero et al. 2011). Here, we see that REBELS-25 sits close to the starburst threshold when using our adopted stellar mass from Bouwens et al. (2022). However, as discussed in Hygate et al. (2023), the classification of REBELS-25 can change depending on the SED fitting applied and stellar mass value adopted.', '5.3 The cosmic evolution of turbulence within galaxies': 'In Figure 8, we plot the average velocity dispersion, ¯ 𝜎 , for REBELS25 and the literature sources as a function of redshift. We plot also the analytical fit to the large IFU survey KMOS3D at 𝑧 ∼ 1 -3 (Übler et al. 2019) for both ionised and atomic/molecular gas (solid black and grey lines, respectively) and extrapolate these fits to higher redshifts (dashed black and grey lines). As discussed in Section 1, this Übler et al. (2019) study of 175 star-forming disc galaxies finds  \n<!-- image -->  \n<!-- image -->  \nFigure 7. Distribution of REBELS-25 and a few comparison samples in the 𝑀 ∗ - SFR plane. We split the sample into three redshift bins and show the main-sequence line from Schreiber et al. (2015) for 𝑧 ∼ 1.25, 3.75 and 6.25 in panels (a), (b) and (c), respectively, by solid black lines. The grey shaded areas show the 1 𝜎 uncertainty on the main-sequence lines, and the green shaded areas show the region occupied by starbursts according to a threshold of 4 × the main-sequence lines (e.g. Rodighiero et al. 2011). Galaxies observed with [CII] emission are plotted with orange markers (apart from REBELS-25, which is plotted with a large red star-shaped marker in panel c), whilst galaxies observed with either [CI] or CO are shown with blue markers. Empty markers indicate galaxies where the kinematic properties are given as limits only.  \n<!-- image -->  \nFigure 8. Velocity dispersion as a function of redshift for REBELS-25 (red star-shaped marker) and other well-resolved galaxies from the literature with CO or [CI] (blue) and [CII] (orange) observations, and The markers are the same as in Figure 7. Empty markers show where one or more of the kinematic parameters ( 𝑉 rot or 𝜎 ) is given as an upper or lower limit. The grey shaded region shows the predicted evolution of 𝜎 from the semi-analytic model of Wisnioski et al. (2015). The yellow-to-purple solid lines show predictions from Illustris TNG50 (Pillepich et al. 2019) in increasing mass bins, whilst the yellow-to-red dashed lines and shaded areas show the predictions from the SERRA simulations (Kohandel et al. 2023) and their 1 𝜎 uncertainties. The empty magenta box shows the prediction from the Krumholz et al. (2018) model that only considers feedback effects as the driver of turbulence, whilst the hatched magenta box shows the prediction from the model that incorporates both feedback effects and gravitational perturbations, as described in Section 5.3.  \n<!-- image -->  \nFigure 9. As with Figure 8, but instead the ratio of ordered to random motion (the ratio of the maximum rotational velocity to the average velocity dispersion) is plotted as a function of redshift.  \n<!-- image -->  \nthat turbulence increases with redshift for both warm ionised gas and cold neutral gas tracers. Wisnioski et al. (2015) find that this evolution is consistent with being driven by a redshift evolution of the gas fraction according to the Toomre stability criterion (Toomre 1964). We plot the predictions from this semi-analytical Wisnioski et al. (2015) model for the H 𝛼 emission of a disc galaxy of mass log ( 𝑀 ∗ / 𝑀 ⊙ ) ∼ 10 with 𝑄 = 1 by the shaded grey area in Figure 8. Additionally, we plot the predictions from IllustrisTNG50 (Nelson et al. 2019) simulations for H 𝛼 emission (Pillepich et al. 2019) and from SERRA simulations (Pallottini et al. 2022) for [CII] emission (Kohandel et al. 2023) for different mass bins. The predictions from the SERRAsimulations are specifically for [CII] emission from mock ALMA observations at different inclinations of 3218 galaxies at 8 ≤ log ( 𝑀 ∗ / 𝑀 ⊙ ) ≤ 10 . 3 from 4 < 𝑧 < 9, with full details of their kinematic analyses given in Kohandel et al. (2023).  \nWefindthat the average velocity dispersion of REBELS-25 agrees better with the SERRA simulations than the other predictions at high𝑧 . This likely follows from the fact that the Wisnioski et al. (2015) and TNG50 predictions are for warm ionised gas rather than cold gas tracers. Übler et al. (2019) find an average difference of 10 -15 km s -1 between velocity dispersions measured from molecular or atomic gas and corresponding measurements from ionised gas at 𝑧 ∼ 1 -3, although this offset may evolve with redshift. Similarly, Kohandel et al. (2023) find that 𝜎 measured from H 𝛼 emission is ∼ 2 × larger than 𝜎 as measured from [CII] emission. However, following from the expectation that higher redshift galaxies should be dominated by turbulent motions, we would still expect the same general trend with the turbulence increasing with redshift. Instead, we see that the average velocity dispersion for REBELS-25 (¯ 𝜎 = 33 ± 9kms -1 ) is comparable to 𝑧 < 3 galaxies.  \nTo investigate why this may be the case, we consider the processes responsible for driving the turbulence, and hence the gas  \nvelocity dispersion, within galaxies. Firstly, for REBELS-25, the derived 𝑓 gas = 0 . 93 + 0 . 03 -0 . 07 is consistent with the predictions of Wisnioski et al. (2015) according to its redshift and sSFR, and one may therefore expect it to have a similar velocity dispersion as predicted by the Wisnioski et al. (2015) model (albeit lower for the cold gas tracer used), where the dispersion is given by:  \n¯ 𝜎 = 1 √ 2 𝑉 rot , max 𝑓 gas ( 𝑧 ) . (2)  \nThis would predict a ¯ 𝜎 of ∼ 240 km s -1 (for a warm ionised gas tracer, so ∼ 120 km s -1 for the [CII] emission according to Kohandel et al. 2023). However, as discussed in Wisnioski et al. (2015), the redshift evolution of the necessary parameters to define 𝑓 gas ( 𝑧 ) , such as the depletion time and the sSFR, are uncertain at 𝑧 ≳ 3, and may therefore not be applicable at the high redshift of REBELS-25.  \nSecondly, we consider the models defined in Krumholz et al. (2018), which find that the turbulence within galaxies is driven by gravitational instabilities and/or feedback effects. For the galaxies in Parlanti et al. (2023b), from which we have selected seven mainsequence star forming galaxies for comparison with REBELS-25, the models that include both feedback effects and gravitational instabilities from Krumholz et al. (2018) were found to agree well with their data. However, in Rizzo et al. (2021), from which we have plotted the kinematic properties of four starburst galaxies and one main-sequence galaxy, models with only energy injected by feedback effects were found to be sufficient. The same is also true for the four galaxies in Roman-Oliveira et al. (2024), which are a mix of main sequence and starburst galaxies. For REBELS-25, we have plotted the predicted velocity dispersion according to both of these Krumholz et al. (2018) models in Figure 8 (for cool atomic or molecular gas). Here we see that, as with the galaxies in the Rizzo et al. (2021)  \nand Roman-Oliveira et al. (2024) samples, the models considering only feedback effects better agree with the velocity dispersion obtained for REBELS-25. This is perhaps due to the lower mass of the galaxies in the Parlanti et al. (2023b) sample in comparison to the Rizzo et al. (2021) and Roman-Oliveira et al. (2024) samples, which could mean that gravitational instabilities have more of an effect on the bulk motion of the galaxies. For REBELS-25, the stellar mass quoted in Table 1 is comparable to those in the Parlanti et al. (2023b) sample, but this value is likely to be an underestimate if there is a large amount of dust obscuration. Improved constraints on the stellar mass will be possible with the upcoming JWST data on this galaxy.  \nFinally, we also consider that, for this sample of galaxies that is biased towards the most active and massive star-forming galaxies, the difference between the velocity dispersion traced by the warm ionised and cold gas emission may be more extreme than in other studies such as Übler et al. (2019), with potentially more intense outflows and other radial motions that could increase the turbulence in the warm ionised gas. In high𝑧 studies, comparing both cold and warm galaxy kinematics is still limited due to scarce data, although some progress is now being made with the advent of JWST (Parlanti et al. 2023a; Fujimoto et al. 2024a; Übler et al. 2024). We note that upcoming analysis of ALMA Cycle 9 Band 8 high resolution observations of the [OIII]88 𝜇 m line for REBELS-25 (Algera et al. in prep, ID:2022.1.01324.S) will help to shed light on the effect of different gas tracers on the fitted kinematics of this galaxy.', '5.4 The cosmic evolution of rotational support': "Another kinematic property found to evolve with redshift is the ratio of maximum rotational velocity to average velocity dispersion, 𝑉 rot , max / ¯ 𝜎 . As with Figure 8, we plot this ratio as a function of redshift for REBELS-25, the literature sources and various predictions in Figure 9. Here, we see that REBELS-25, and most of the galaxies in our comparison sample, have 𝑉 rot , max / ¯ 𝜎 > 2. This places the majority of these sources above the predictions from the Wisnioski et al. (2015) model and from the TNG50 simulations, although we again note that these predictions are for warm ionised gas tracers, which are likely to result in higher velocity dispersions and therefore lower 𝑉 rot , max / ¯ 𝜎 .  \nIn Figure 9, we again plot the predictions for the evolution of the degree of rotational support from the SERRA simulations (Kohandel et al. 2023). Here we see that REBELS-25 is consistent within the uncertainties with the predictions for the log ( 𝑀 ∗ / 𝑀 ⊙ ) > 9 sample. This is also the case for most of the galaxies in our comparison sample at 𝑧 > 4, although we note that the galaxies from Parlanti et al. (2023b) have large uncertainties such that they are consistent with most of the predictions shown here.  \nThe discrepancies between the predictions from SERRA and TNG50 could be due to a combination of reasons. The first possible cause is, as discussed above, the SERRA predictions shown in Figure 9 are for kinematics as traced by [CII] emission, whereas the TNG50 predictions are for H 𝛼 emission. Another possible cause is that SERRA uses mock ALMA observations (Kohandel et al. 2023) whereas TNG50 takes the intrinsic kinematic properties from a face-on projection of the simulated galaxies (Pillepich et al. 2019). However, these reasons do not explain why the SERRA predictions show almost no redshift evolution. Instead, the SERRA predictions mainly indicate that the stellar mass of the galaxy has a more significant impact on its degree of rotational support. Similar emphasis on this mass dependence is also made in Dekel et al. (2020), where both the VELA simulations and SAMI observations indicate a significant increase in the degree of rotational support at log ( 𝑀 ∗ / 𝑀 ⊙ ) ∼ 9. Sim-  \nt al. (2021) find that stellar mass, rather than redshift, most strongly correlates with the disc fraction amongst star-forming galaxies at 𝑧 ∼ 1 . 5, observing only a modest increase in the prevalence of discs between 𝑧 ∼ 1 . 5 and 𝑧 ∼ 0 . 04 at fixed stellar mass. Additionally, recent FIRE-2 simulations have implied that disc galaxies can exist at any redshift in sufficiently massive haloes (Gurvich et al. 2022). However, as can be seen from the mass ranges of our comparison sample in Figure 7, there is a lack of high resolution observations of cold gas kinematics in low-mass ( 𝑀 ∗ < 10 9 M ⊙ ) galaxies at high𝑧 to confidently validate these hypotheses.  \nTo investigate this mass dependence further, we therefore also make comparisons to a recent study of the ionised gas kinematics at high𝑧 with JWST NIRSpec Multi-Object Spectroscopy (MOS) observations (de Graaff et al. 2023), and also studies of individual galaxies with NIRSpec IFU observations, including a lensed galaxy at 𝑧 = 6 . 072 (Cosmic Grapes, Fujimoto et al. 2024a), ALESS073.1 at 𝑧 = 4 . 76, and GN-z11 at 𝑧 = 10 . 60 (Xu et al. 2024). We note that due to the compactness of GN-z11, it is not currently possible to rule out a merger scenario from the current spatial resolution (Xu et al. 2024). Additionally, the observations from de Graaff et al. (2023) use the NIRSpec MOS, and so it is also not possible to differentiate between merging and rotating disc systems since the objects are resolved by only a few resolution elements along a single spatial direction. However, the galaxies studied in de Graaff et al. (2023) are amongst the lowest mass galaxies studied kinematically at 𝑧 > 1 so far, and therefore make for an interesting additional comparison sample. We also add that the kinematics of ALESS073.1 have been studied with both ALMA ([CII], Lelli et al. 2021) and with JWST (H 𝛼 , Parlanti et al. 2023a) observations at high resolution, as has the 'Cosmic Grapes' 𝑧 = 6 . 072 galaxy from Fujimoto et al. (2024a). In Figure 10, we plot the 𝑉 rot , max / ¯ 𝜎 against stellar mass for the same galaxies as in Figure 9, along with the seven galaxies from de Graaff et al. (2023) in magenta diamond markers.  \nFigure 10 appears to show a tentative positive correlation between 𝑉 rot , max / ¯ 𝜎 and 𝑀 ∗ . For the warm ionised gas tracers, it is likely that the velocity dispersions are higher than if they were traced by cold gas (as seen for Cosmic Grapes and ALESS073.1). We could therefore assume that the cold gas velocity dispersion would be lower, which would in turn increase 𝑉 rot , max / ¯ 𝜎 . It is therefore not possible to confidently assess this trend. It should also be noted that there are large uncertainties in deriving the stellar mass at high redshift, which may also affect any observed relationship. Future studies with spatially resolved SED fitting may help to improve stellar mass estimates at high𝑧 (e.g., Abdurro'uf et al. 2023; Giménez-Arteaga et al. 2023).", '6 SUMMARY AND CONCLUSIONS': 'In this paper, we have presented follow-up high-resolution ( ∼ 710 pc) ALMA [CII] and dust ( ∼ 150 𝜇 m) continuum observations of REBELS-25; a massive star-forming galaxy at 𝑧 = 7 . 31, originally targeted as part of the ALMA REBELS LP (Bouwens et al. 2022). In previous studies (Stefanon et al. 2019; Hygate et al. 2023), this galaxy was found to have multiple UV clumps, indicating a potential disturbed, merging or clumpy morphology. However, the follow-up high resolution observations presented here indicate that the [CII] emission in this galaxy is well fit by a near-exponential disc profile, with a Sérsic index of 1 . 3 ± 0 . 2. These UV clumps that are offset from the [CII] and dust are therefore likely to be the result of dust obscuration, or are potentially star forming clumps embedded into the gas disc. As well as 2D Sérsic fitting, we have used CANNUBI to fit additional morphological parameters of this galaxy, and find  \nFigure 10. As with Figures 8 and 9, but instead the ratio of ordered-torandom motion is plotted as a function of stellar mass. The magenta markers are from recent JWST studies of ionised gas kinematics (de de Graaff et al. 2023; Fujimoto et al. 2024a; Parlanti et al. 2023a; Xu et al. 2024). The black vertical dashed line indicates a cut-off stellar mass at ∼ 10 9 M ⊙ described in Gurvich et al. (2022) as a threshold mass where rotationally supported discs can form.  \n<!-- image -->  \nthat it has a low inclination ( 𝑖 = 25 ± 6 · ) and potentially a very thin disc ( 𝑍 0 < 710 pc). In addition, we see some evidence of more complex morphological features, including tentative evidence of a bar, identified by fitting elliptical isophotes to both the [CII] and dust continuum emission.  \nWe find the [CII] kinematics of REBELS-25 are well explained by a rotation-dominated disc using the 3D tilted ring fitting tool, 3D BAROLO , and several independent criteria to distinguish between a disc and a major merger, including the five-disc criteria from Wisnioski et al. (2015) and Jones et al. (2021), and also the PV Split method from Rizzo et al. (2021). The best-fit rotating disc model with 3D BAROLO reveals a low overall turbulence ( ¯ 𝜎 = 33 ± 9 km s -1 ) and a high ratio of 𝑉 rot , max / ¯ 𝜎 = 11 + 8 -4 . This low average dispersion velocity obtained for REBELS-25 is consistent with stellar feedback as the main driver of turbulence within this galaxy, as has also been found for other dynamically cold discs in the high redshift Universe (Rizzo et al. 2021; Roman-Oliveira et al. 2024). However, we do see some evidence of non-circular motions, which could be due to inflows/outflows, a minor merging component, a central bar and/or spiral arms.  \nUsing the rotational velocity at 𝑟 𝑒 (equivalently, this is where the curve begins to flatten), we find a total dynamical mass of (1 . 2 ± 0 . 3 )× 10 11 M ⊙ . This results in a gas mass of 𝑀 gas , tot = ( 1 . 1 + 0 . 4 -0 . 3 ) × 10 11 M ⊙ and 𝛼 [ CII ] = 62 + 32 -22 M ⊙ / L ⊙ for 𝑀 ∗ = 8 + 4 -2 × 10 9 M ⊙ , however, the total stellar mass of this galaxy is likely very uncertain thanks to dust obscuration, therefore we find that estimates of 𝑀 gas and 𝛼 [ CII ] are currently also highly uncertain. Upcoming JWST observations (GO-1626, GO-6036 and GO-6480) will help improve estimates of the stellar mass and stellar morphology of the galaxy, enabling future work on dynamical modelling and rotation curve decomposition of REBELS-25 with the kinematic information obtained in this paper.  \nAn increasing number of dynamically cold discs have been identified in the high𝑧 Universe, although these observations are often not sufficiently resolved to confidently distinguish between mergers and discs. With these sub-kpc resolution observations, REBELS-25  \nis amongst the most distant robustly confirmed cold discs observed to date. This finding of a very distant ( 𝑧 = 7 . 31), very dynamically cold ( 𝑉 rot , max / ¯ 𝜎 = 11 + 8 -4 ) disc challenges the predictions from many state-of-the-art models and cosmological simulations, which tend to predict very turbulent and dispersion-dominated discs at 𝑧 > 3. By comparing to other cold gas kinematics studies of 𝑧 > 0 . 5 galaxies with observational data of similar quality, we find that dynamically cold discs seem to be more common in the high𝑧 Universe than predictions based on warm ionised gas tracers (e.g., Wisnioski et al. 2015; Pillepich et al. 2019; Übler et al. 2019), although we note an observational bias towards massive star forming galaxies. However, these observations are more consistent with recent predictions from the SERRAcosmological simulations (Pallottini et al. 2022) based on mock observations of [CII] emission (Kohandel et al. 2023), suggesting that the kinematic tracer used significantly impacts the derived velocity dispersion, and therefore degree of rotational support. For the case of REBELS-25, both high resolution ALMA observations of [OIII]88 𝜇 memission (Algera et al. in prep, ID:2022.1.01324.S) and JWSTobservations of [OIII]5007 kinematics (GO-6036 PI J. Hodge, and GO-6480 PI S. Schouws) will enable a comparison between cold and warm ionised gas tracers.  \nREBELS-25is a massive ( 𝑀 ∗ = 8 + 4 -2 × 10 9 M ⊙ ) galaxy, and several studies have suggested a stronger dependence on the mass of the galaxy, rather than its redshift, in setting the evolution of the gas kinematics (e.g. Dekel et al. 2020; Gurvich et al. 2022; Kohandel et al. 2023). A comparison to lower-mass galaxies tentatively provides observational evidence to support this (Figure 10), although further work, particularly at lower stellar masses, is necessary to study this mass dependence. JWST will be a powerful tool for such studies, with ionised gas kinematic studies of fainter, lower mass galaxies at high𝑧 now becoming more feasible (de Graaff et al. 2023).  \nOverall, in this work we have shown that dynamically settled rotating disc galaxies, such as REBELS-25, can form as early as just 700 Myr after the Big Bang. We therefore expect that future, high resolution studies of cold gas kinematics at high𝑧 will reveal even more cold, massive discs. In particular, ongoing ALMA observations of other REBELS galaxies will enable robust kinematic modelling of additional rotating disc candidates at 𝑧 ∼ 6 -8 (Phillips et al. in prep).', 'ACKNOWLEDGEMENTS': 'The authors thank Jorge González-López for their useful discussions about the JvM effect and correction. The authors also acknowledge Mahsa Kohandel for useful discussions relating to the SERRA simulations. JH acknowledges support from the ERC Consolidator Grant 101088676 (VOYAJ). M. R. is supported by the NWO Veni project \'Under the lens\' (VI.Veni.202.225). PD acknowledges support from the NWO grant 016.VIDI.189.162 (\'ODIN") and from the European Commission\'s and University of Groningen\'s CO-FUND Rosalind Franklin program. AF acknowledges support from the ERC Advanced Grant INTERSTELLAR H2020/740120. IDL acknowledges funding support from ERC starting grant 851622 DustOrigin. M.A. acknowledges support from FONDECYT grant 1211951, and ANID BASAL project FB210003. RAAB acknowledges support from an STFC Ernest Rutherford Fellowship [grant number ST/T003596/1]. PEMP acknowledges support from the Dutch Research Council (NWO) through the Veni grant VI.Veni.222.364.', 'DATA AVAILABILITY': 'The ALMA observations used in this article are available in the ALMA archive https://almascience.eso.org/aq/ and can be accessed with their project codes: 2019.1.01634.L and 2021.1.01603.S. The COSMOS-DASH HST image used available from the Mikulski Archive for Space Telescopes on the COSMOS-DASH page: https://archive.stsci.edu/hlsp/cosmos-dash. The data underlying this article will be shared on reasonable request to the corresponding author.', 'REFERENCES': '- Dayal P., et al., 2022, Monthly Notices of the Royal Astronomical Society, 512, 989\n- Hohl F., 1971, The Astrophysical Journal, 168, 343  \n- Neeleman M., Prochaska J. X., Kanekar N., Rafelski M., 2020, Nature, 581, 269  \n- 487, 4305  \nvan der Wel A., et al., 2014, The Astrophysical Journal, 792, L6 Übler H. D. N., et al., 2019, The Astrophysical Journal, 880, 48 Übler H., et al., 2024, GA-NIFS: NIRSpec reveals evidence for non-circular motions and AGN feedback in GN20, http://arxiv.org/abs/2403.  \n- 03192', 'APPENDIX A: TESTS USING DIFFERENT [CII] LINE CUBES': 'In Table A1, we list the properties of the ALMA [CII] line cubes used to test our main findings. The first three rows correspond to Cycle 8 only line cubes created with a range of channel widths (5, 10 and 15 km s -1 ) and using both natural and Briggs-weighted data sets. For the measurement set labelled "LP + Cy 8, JvM-corr", we concatenate the low resolution REBELS LP data with the high resolution Cycle 8 observations. We follow the same reduction steps as outlined in Section 2.1 for each data set, however for the combined cube we add additional steps, which are described below, to account for the different array configurations used in the two observations.  \nFor the combined cube, we are impacted by the JvM effect, as mentioned in Section 2.1. This effect occurs when the volume of the CLEAN beam differs from the volume of the dirty beam, resulting in an incorrect flux scaling of the CLEAN residual map. Consequently, if no correction is applied, we derive a [CII] flux from the concatenated dataset that is ∼ 3 . 4 times higher than the [CII] flux from the LP data (Hygate et al. 2023) alone. Recently, a correction for this effect has been suggested by Czekala et al. (2021). In summary, a correction factor, 𝜖 , is determined from the ratio of the dirty to CLEAN beam volume (taken at the first null of the dirty beam), and the residuals are then scaled by this factor before being added to the convolved CLEAN model to produce the final imaging.  \nAs in Posses et al. (2024), we find that cleaning down to a threshold of 1 𝜎 RMS minimises this effect, and we find that binning to a higher channel width (20 km s -1 ) also reduces this effect. Since these observations were taken over two years apart and with multiple array configurations, we also apply the statwt task to the continuum-subtracted measurement set to re-assign weights based on the measured noise. For this statwt task, we again fit only line-free channels, since emission could contaminate the weight calculation. Wethen apply the JvM correction with 𝜖 = 0 . 18, as determined from the ratio of the CLEAN and dirty beam volumes. Applying this correction results in a consistent total [CII] flux (in a circular aperture with a radius of 1.75") between the combined cube, the Cycle 8 data only cubes, and the LP data. However, as discussed in, for example, Casassus & Carcamo (2022); Posses et al. (2024), the application of the JvM correction means that the residuals have been rescaled, and so the RMS is no longer a true representation of the sensitivity of the observations. This can therefore significantly exaggerate the SNR(Casassus & Carcamo 2022). When fitting the morphology and kinematics on the concatenated cube with CANNUBI and 3D BAROLO , we therefore increase the GROWTHCUT parameter to 6. We note that the uncertainties on the fitted parameters for the concatenated cube are likely underestimated.  \nWe run some initial morphology fitting tests for each data set with Sersic2D and CANNUBI . Notably, we find a Sérsic index of ∼ 1 for all cubes tested, and we find the bright, misaligned component discussed in Section 3.3 is present in all of the different cubes tested. We overall find morphological parameters that are consistent with those in Section 3.1 and Table 2, and we therefore fix all the morphological parameters to those in Table 2 when fitting with 3D BAROLO . For these kinematic fits, we leave only 𝑉 rot, 𝜎 , and PA kin as free parameters,  \nTable A1. Table showing some of the observation details of the different [CII] line cubes tested, and the resulting ratio of ordered to random motion ( 𝑉 rot / ¯ 𝜎 ).  \nas done with the final kinematic fitting in Section 4.1. For the cube binned to 5 km s -1 and for the cube with Briggs weighting, we fit only 4 rings due to the lower SNR. The 𝑉 rot , max / ¯ 𝜎 and their uncertainties are then derived using the same method as in Section 4.2, and are listed in Table A1, where we see that all results are consistent with the main findings of this paper. We also show in Figure A1 the PVDs for the kinematic fit to the JvM-corrected cube. In all rotation curves produced from these tests, the bump feature, where there is a higher 𝑉 rot value in the first ring, is present. From these tests, we therefore conclude that the kinematic parameters we have derived in the main body of the paper are robust.', 'APPENDIX B: CANNUBI CORNER PLOTS': 'In Figure B1, we show the posterior distributions obtained by applying CANNUBI to our [CII] observations of REBELS-25. The inclination (incl) and position angle (pa) are given in degrees, the disc thickness (Z0) in arcseconds and the geometric centre (x0 and y0) in arbitrary pixel units that are equivalent to J2000 RA, Dec of 10:00:32.3397, +1:44:31.167.', 'APPENDIX C: CHANNEL MAPS': 'In Figure C1, we show representative channel maps of the [CII] emission of REBELS-25 (top panels) and of the best-fit kinematic model (bottom panels) obtained with 3D BAROLO . We show the data in the upper panels with blue contours for the positive emission and grey, dashed contours for the negative emission. The model is shown with red contours in the lower panels. The contours follow the emission intensity according to 2 𝜎 RMS and 4 𝜎 RMS , with the negative contours at -2 𝜎 RMS . The green cross shows the kinematic centre fitted by 3D BAROLO .  \nThis paper has been typeset from a T E X/L A T E X file prepared by the author.  \nFigure A1. The major-axis (top) and minor-axis (bottom) PVDs for the bestfit 3D BAROLO model to the JvM corrected cube. The contours for the data (in black) and the model (in red) are at 6, 12 and 18 𝜎 RMS as calculated from 3D BAROLO . Negative contours for the data at -12 and -6 𝜎 RMS are indicated by the dashed grey contours. Note, however, that this 𝜎 RMS is not a true representation of the sensitivity of the data due to the scaling of the residuals in the JvM correction. The yellow markers in the top panel indicate the separation of the six rings fitted by 3D BAROLO and their fitted rotation velocities along the line-of-sight. The asymmetry in the model of the minoraxis PVD is the result of a changing PA across the disc which is preferred for this fit.  \n<!-- image -->  \n200  \n100  \n200  \n100  \nFigure B1. A triangle plot of the posterior distributions for the morphological parameters fitted by CANNUBI .  \n<!-- image -->', 'Beam': '<!-- image -->  \n<!-- image -->  \nFigure C1. Representative channel maps of REBELS-25, with contours at 2 and 4 𝜎 RMS for the data in blue (upper panels) and the model in red (lower panels). Negative contours at -2 𝜎 RMS are indicated by the grey dashed contours. As with Figure 5, 𝜎 RMS = 168 𝜇 Jy beam -1 and is equal to one standard deviation above the median value as calculated by 3D BAROLO . The centre of the galaxy as determined from CANNUBI and fixed in the kinematic fit is marked by a green cross.  \n<!-- image -->'}