REBELS-25: discovery of a dynamically cold disc galaxy at z = 7.31

ABSTRACT

We present high-resolution (⁠|$\sim 0.14$| arcsec = 710 pc) Atacama Large Millimetre/submillimetre Array [C ii] 158 |$\mu$|m and dust continuum follow-up observations of REBELS-25, a [C ii]-luminous (⁠|$L_{\mathrm{[CII]}}=(1.7\pm 0.2)\times 10^9\, \mathrm{L_{\odot }}$|⁠) galaxy at redshift |$z=7.3065\pm 0.0001$|⁠. These high-resolution, high signal-to-noise observations allow us to study the sub-kpc morphology and kinematics of this massive (⁠|$M_* = 8^{+4}_{-2} \times 10^9 \mathrm{{\rm M}_{\odot }}$|⁠) star-forming (SFR|$_{\mathrm{UV+IR}} = 199^{+101}_{-63} \mathrm{{\rm M}_{\odot }} \mathrm{yr}^{-1}$|⁠) galaxy in the Epoch of Reionization. By modelling the kinematics with |$^{\mathrm{3D}}$|BAROLO, we find it has a low-velocity dispersion (⁠|$\bar{\sigma } = 33^{+9}_{-7}$| km s|$^{-1}$|⁠) and a high ratio of ordered-to-random motion (⁠|$V_{\mathrm{rot, ~max}}/\bar{\sigma } = 11 ^{+6}_{-5}$|⁠), indicating that REBELS-25 is a dynamically cold disc. Additionally, we find that the [C ii] distribution is well fit by a near-exponential disc model, with a Sérsic index, n, of |$1.3 \pm 0.2$|⁠, and we see tentative evidence of more complex non-axisymmetric 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 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 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

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 more violent gas accretion (Conselice, Rajgor & Myers 2008; Dekel et al. 2009), as well as more intense active galactic nuclei (AGNs) 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 |$z\gt 1$| galaxies indeed revealed more clumpy and irregular morphologies (e.g. Conselice, Chapman & Windhorst 2003; Elmegreen & Elmegreen 2005; Papovich et al. 2005; Conselice et al. 2008), in comparison to the rotationally supported disc galaxies that make up 50–80 per cent of the |$z\sim$| 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 ionized gas kinematics have found that the turbulence (as measured by the velocity dispersion, |$\sigma$|⁠, of the gas) within galaxies increases from |$z=0$| to |$z \sim 3$|⁠, and that the degree of rotational support (rotational velocity over velocity dispersion, |$V_{\mathrm{rot}}/\sigma$|⁠) decreases over the same redshift range (e.g. Wisnioski et al. 2015, 2019; Di Teodoro, Fraternali & Miller 2016; Johnson et al. 2018; Übler et al. 2019).

While these studies of the |$z\lt 3$| Universe are in line with theoretical expectations, a growing body of work on early (⁠|$z\gt 3$|⁠) massive (⁠|$M_*\gtrsim 10^{9}-10^{10} \mathrm{{\rm M}_{\odot }}$|⁠) galaxies begins to challenge this picture in the more distant Universe. The JWST, with its unprecedented high-angular resolution and sensitivity in the near- and mid-infrared, has now revealed more massive (⁠|$M_*\gtrsim 10^{9}{\!-\!}10^{10} \mathrm{{\rm M}_{\odot }}$|⁠) disc-like galaxies at |$z\gt 6$| (Adams et al. 2022; Akins et al. 2023; Atek et al. 2023; Casey et al. 2023; Labbé et al. 2023; Rodighiero 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 disc galaxies at |$z\gt 3$| (Smit et al. 2018; Sharda et al. 2019; Neeleman et al. 2020, 2021, 2023; Rizzo et al. 2020, 2021; Fraternali et al. 2021; Lelli et al. 2021; Tsukui & Iguchi 2021; de Graaff et al. 2024; Pope et al. 2023; Posses et al. 2023; Roman-Oliveira, Fraternali & Rizzo 2023; Parlanti et al. 2024; Ü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 |$z\lesssim 2$| (Pillepich et al. 2019; Simons et al. 2019).

These 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; Tiley et al. 2021; Gurvich et al. 2022; Kohandel et al. 2023), with gas discs only surviving in |$M_*\gtrsim 10^9 \mathrm{ {\rm M}_{\odot }}$| 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 |$z\gt 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, multiwavelength 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 [C ii] (157.7 |$\mu$|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 (⁠|$T \sim 100$| K) neutral and molecular gas in galaxies (Shibai et al. 1996; Stacey et al. 2010; Pineda et al. 2013; Vallini et al. 2015).

The Reionization Era Bright Emission Line Survey (REBELS; Bouwens et al. 2022) is an ALMA large programme (LP) that has carried out spectral scans for the [C ii] line in ALMA’s Band 6 for |$\sim$| 40 UV-bright star-forming galaxies at |$z\sim 6-8$|⁠. As well as multiple, significant [C ii] 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 programme has highlighted ALMA’s capabilities as a ‘redshift machine’ (see also Smit et al. 2018; Schouws et al. 2022), as well as enabling studies of, for example, dust-obscured star formation (Fudamoto et al. 2021), specific star formation rates (sSFR; Topping et al. 2022), Ly |$\alpha$| emission (Endsley et al. 2022), and dust and ISM properties (Dayal et al. 2022; Ferrara et al. 2022; Sommovigo et al. 2022) at |$z\gt 6$|⁠. Some of these REBELS galaxies show evidence of velocity gradients in the LP data (Schouws et al. in preparation), 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 (⁠|$\gtrsim$|1 arcsec) 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 (⁠|$\sim 1.3$| arcsec |$= 6.7$| kpc) was found to be insufficient to distinguish between a rotating disc and a 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 [C ii] 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 [C ii]-luminous REBELS galaxies (ID 2022.1.01131.S, PI: R. Smit).

In 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 preparation). 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-fitting 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 |$\Lambda$|CDM cosmology with Hubble constant |$H_0=70$| km s|$^{-1}$| Mpc|$^{-1}$|⁠, matter density |$\Omega _m=0.3$| and vacuum energy |$\Omega _{\Lambda }=0.7$|⁠. At the redshift of REBELS-25, this corresponds to a luminosity distance of 72 519 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^{\mathrm{h}} 00^{\mathrm{m}} 32.32^{\mathrm{s}}$|⁠, |$+01^{\circ }44^{\prime }31{_{.}^{\prime\prime}}3$| from rest-frame UV imaging (Stefanon et al. 2019) and a [C ii] spectroscopic redshift of |$z=7.3065\pm 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 |$\sim 1.7$| arcsec and an angular resolution of |$\sim 0.1$| arcsec. In total, 256.67 min 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 |$uv$|-data were therefore used without further modification to the calibration scheme or flagging. This Cycle 8 data set contains a spectral window covering the [C ii] line with a bandwidth of 1875 MHz, and three other spectral windows for detecting the FIR dust continuum surrounding the line (⁠|$\sim 150 \mu$|m), each with a bandwidth of 2000 MHz.

To obtain a [C ii] line-emission cube, we first created a continuum-free 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 [C ii] 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 subparameters initially determined by the recommendations of the ALMA automasking guide¹ 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 subparameters are left at the recommended values (‘lownoisethreshold’ = 1.5, ‘negativethreshold’ = 7, ‘minbeamfrac’  = 0.3). The native channel width of the data is |$\sim$| 5.1 km s|$^{-1}$|⁠. In the final line-cube, we bin the channels by a factor of 3 to a channel width of |$\sim 15.4$| km s|$^{-1}$|⁠, as a compromise between sensitivity and velocity resolution.

All imaging was carried out with a multiscale 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 |$\times$| the beam FWHM, and 3 |$\times$| the beam FWHM. We find that the use of this multiscale CLEAN, in comparison to the default ‘hogbom’ deconvolver, is more successful at masking both compact and extended emission during the cleaning process.

A |$\sim 150\, \mu$|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.

With 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 [C ii] imaging, the natural-weighted beam size is 0.14 arcsec × 0.13 arcsec (710 × 660 pc), and the Briggs-weighted beam size is 0.12 arcsec × 0.10 arcsec (610 × 510 pc). For the dust continuum, the natural-weighted beam size is 0.14 arcsec × 0.13 arcsec (710 × 660 pc), and the Briggs-weighted beam size is 0.11 arcsec × 0.09  arcsec (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 ratio, 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 [C ii] line-cube, the final RMS per 15 km s|$^{-1}$| channel is |$\sim 110\, \mu$|Jy beam|$^{-1}$|⁠. For the dust continuum map, the final RMS is |$\sigma \sim$|5.7 |$\mu$|Jy beam|$^{-1}$|⁠. We find consistent [C ii] and dust continuum flux densities, within the uncertainties, when comparing the high resolution data with the low-resolution REBELS LP data.

Whilst 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 [C ii] 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 [C ii] 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 Fig. 1, we present the [C ii], dust and rest-frame UV morphology of REBELS-25 at sub-kpc resolution. The [C ii] emission reveals a clear extended disc, whereas the rest-frame UV morphology can be separated into clumps that are offset from the [C ii] and dust continuum. Both the [C ii] 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.

Figure 1.

Left: [C ii] moment-0 map from the naturally weighted data cube, with contours showing 2, 3,...14|$\sigma _{\mathrm{RMS}}$| emission (where |$\sigma _{\mathrm{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|$\sigma _{\mathrm{RMS}}$| emission (where |$\sigma _{\mathrm{RMS}}=5.5 ~\mu$|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 [C ii] emission and dust continuum shown by the turquoise and orange contours, respectively.

3.1 [C ii] morphology

To investigate the morphology of the [C ii] 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 [C ii] line emission for the naturally weighted line cube (Fig. 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\pm 28^{\circ }$|⁠, an effective radius, |$r_{\mathrm{e}}$|⁠, of 0.42±0.06 arcsec (⁠|$2.2\pm 0.3$| kpc), and a Sérsic index, n, of |$1.3 \pm 0.2$|⁠, indicating that the [C ii] emission is well fit by a near-exponential disc. The best-fitting ellipticity is 0.13|$\pm$|0.09, which corresponds to an inclination, i, of |$30^{+10}_{-15}$||$^{\circ }$|⁠, assuming the disc is razor thin. However, previous studies have found that high-z 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.

We therefore use the cannubi² package to also fit for the thickness and inclination of the disc. As introduced in Mancera Piña et al. (2020), Fraternali et al. (2021), Mancera Piña et al. (2022a), and Roman-Oliveira et al. (2023), cannubifits the inclination (i), the disc thickness (⁠|$Z_0$|⁠), the morphological centre (⁠|$x_{0,\mathrm{morph}},y_{0,\mathrm{morph}}$|⁠), and morphological position angle (⁠|${\rm PA}_\mathrm{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 [C ii] flux map with resolution-matched 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-fitting 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 autocorrelation times of the parameters (Goodman & Weare 2010). The posterior distributions from the cannubi fitting to REBELS-25 are shown in Fig. B1, and we list the median values and their uncertainties (from the 16th and 84th percentiles) in Table 2.

We find that, overall, the [C ii] morphological parameters are well constrained via the cannubi fitting procedure. One exception is the fitting for |$Z_0$|⁠, from which we can only infer that the best-fitting models have a thickness lower than the |$\sim 0.14\,\rm arcsec$| (710 pc) resolution of our observations. From tests carried out on mock galaxies in Roman-Oliveira 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. To be comparable with other high redshift kinematic studies (e.g. van der Wel et al. 2014; Graaff et al. 2024), we assume a thin disc with an intrinsic axis ratio (⁠|$q_0$|⁠) of 0.2, which corresponds to a disc thickness of |$\sim 460$| pc, although we add if we fix |$Z_0$| to 710 pc in the subsequent analysis we obtain consistent results within the uncertainties. The posterior for the inclination also appears to have a tail extending up to the 15|$^\circ$| 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 in fact lower, this would only increase the derived intrinsic rotational velocities and further strengthen our findings. The |$\rm PA_{\mathrm{morph}}$| posterior also shows a secondary, lower peak at |$\sim 160 ^\circ$|⁠, which is potentially caused by a bright, central, misaligned component which we discuss in Section 3.3.

The morphological centre and |${\rm PA}_{\mathrm{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. We find that this surface brightness profile is best fit with a 1D Sérsic profile with |$n = 1.13 \pm 0.06$|⁠, which is also consistent with the 2D Sérsic fit. The inclination derived from cannubi, |$i=25\pm 6^{\circ }$|⁠, is lower but within error of the |$i=56\pm 29^{\circ }$| value determined from the ratio of the major and minor axes of the low resolution LP data in Hygate et al. (2023), where the [C ii] 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.

In both the cannubi and sersic2d fitting, we see some residual (⁠|$\sim 6\sigma _{\mathrm{RMS}}$| in the cannubi fitting) emission at the centre of the [C ii] map (Fig. 2). This feature is discussed in more detail below, in Section 3.3. Despite these residuals, the reduced chi squared value for the cannubimodel is |$\chi ^2 \sim 0.9$|⁠. For the Sérsic model, |$\chi ^2 \sim 0.18$|⁠, implying that sersic2d is overfitting the data. 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 [C ii] emission returned by cannubi as our fiducial values, and list these in Table 2.

Figure 2.

Left: We show the best-fitting cannubi model, with the observed [C ii] emission from Fig. 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-fitting model, respectively. Right: The residuals from the best-fitting cannubi model, with 2, 3,..., 6|$\sigma _{\mathrm{RMS}}$| contours shown by the solid black lines, and −2, −3|$\sigma _{\mathrm{RMS}}$| contours show by the dashed black lines, where |$\sigma _{\mathrm{RMS}}$| is the local RMS in the [C ii] moment-0 map.

3.2 Dust continuum morphology

In Fig. 1, we also show the dust continuum map of REBELS-25, where we obtain a peak SNR of 10.1. As with the [C ii] morphology, we also fit a 2D Sérsic model to the dust continuum and find |$r_e=0.50$| |$\,\pm 0.08$| arcsec (⁠|$2.5\pm 0.7$| kpc) and |$n=2.2\pm 0.4$|⁠, indicating emission that is more centrally peaked than would be expected for an exponential disc. We find |$\sim 3\sigma$| residuals in clumps to the North and South of the centre from the sersic2d fitting. These clumps are around 0.5 arcsec (⁠|$\sim 2.5$| kpc) from the fitted centre, and each contain around 10 per cent of the flux of the central component.

3.3 Comparing the [C ii], dust, and rest-frame UV morphology

From the lower resolution REBELS LP data, a potential offset of |$0.17 \pm 0.04\,\rm arcsec$| (around |$0.9 \pm 0.2$| kpc) between the [C ii] 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 [C ii] emission are in good agreement, again highlighting the necessity of observations with sufficient angular resolution.

The [C ii] 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 |$PA_{\mathrm{morph}}$| of the [C ii] and dust to differ by |$\sim 60^{\circ }$| at the outer extent of emission.

We also find that the effective radii of the [C ii] and dust emission are comparable, with the |$r_e$| 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\times$| more compact than the [C ii] emission (e.g. Fujimoto et al. 2020; Fudamoto et al. 2022). However, we note that the Sérsic model is a poor fit to the dust emission, and this large |$r_e$| 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.

Notably, both the [C ii] and dust show a central region in the inner |$\sim$| 1 kpc (around 1.4–2 beams wide) that appears misaligned with the |$PA_{\mathrm{morph}}$| of the extended gas disc. We also see this misaligned, central region in the residuals around the fitted centre of the galaxy in both the [C ii] and dust continuum 2D fits, as well as a tail in the posterior distribution for PA|$_{\mathrm{morph}}$| with cannubi. To illustrate this, we fit elliptical isophotes to the [C ii] 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|$\sigma _{\mathrm{RMS}}$|⁠. For the [C ii] 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 leave all parameters free for the dust fitting. The resulting fits are shown in Fig. 3.

Figure 3.

Top left: The elliptical isophote fit to the [C ii] 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 [C ii] (turquoise) and dust (orange) emission.

Within the central 1 kpc, the ellipses fitted to the [C ii] and dust have an average PA of |$\sim 120^{\circ }$|⁠. At roughly the same radii (⁠|$\sim 1$| kpc), the PA for the [C ii] suddenly increases to an average PA of |$\sim 220^{\circ }$| at |$\gtrsim 1$| kpc (consistent with |${\rm PA}_{\mathrm{morph}}$| from cannubi and sersic2d), and for the dust it changes more gradually to |$\sim 190^{\circ }$| (consistent with |$PA_{\mathrm{morph}}$| from sersic2d for the dust). In addition, the ellipticity in the central 1 kpc in both the [C ii] and dust is found to be |$\sim 0.3$|⁠, whereas for the extended [C ii] disc this is found to be |$\sim 0.1$| (near-circular).

Similar 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 |$z=2.467$|⁠, Amvrosiadis et al. 2024 at |$z=3.762$|⁠, Le Conte et al. 2024 at |$1\le z\le 3$|⁠). However, we note that this tentative bar-like central component in REBELS-25 is barely resolved with these observations.

By contrast, the rest-frame UV morphology differs significantly from the [C ii] and dust. As previously identified in Hygate et al. (2023) and Schouws et al. (2022), the HST imaging reveals three clumps that are |$\sim 0.6$| arcsec offset from the centre of both the [C ii] and dust continuum emission. We see in Fig. 1 that these UV clumps lie close to and outside the outer edge of the [C ii] 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 [C ii] 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 |$^{\mathrm{3D}}$|barolo (Di Teodoro & Fraternali 2015, v1.7). |$^{\mathrm{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 channel-by-channel basis, and uses Monte Carlo methods to return the best-fitting 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.

Before the modelling, we extract a 3.5 arcsec |$\times$| 3.5 arcsec subcube centred on the [C ii] emission that covers |$\pm 4\times$| the FWHM of the [C ii] 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\sigma _{\mathrm{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 1.68|$\times 10^{-4}$| Jy beam|$^{-1}$|⁠, which is equivalent to 2|$\sigma _{\mathrm{RMS}}$| of the subcube as calculated by|$^{\mathrm{3D}}$|  barolo, is sufficient to produce a mask that encloses all of the [C ii] emission without introducing too much noise. We set the radial separation of the rings to 0.11 arcsec (⁠|$\sim$|80 per cent of the beam FWHM), resulting in five rings within the mask produced by |$^{\mathrm{3D}}$|barolo. This maximises the number of rings whilst ensuring each ring is nearly comparable to an independent resolution element. We add that when the radial separation is equal to the beam FWHM, we obtain consistent results within the uncertainties. For all fits with |$^{\mathrm{3D}}$|barolo, we use an azimuthal normalisation of the surface density, we fit the kinematics using both sides of the rotation curve, and we minimize the reduced chi squared statistic.

With |$^{\mathrm{3D}}$|barolo, there are nine possible parameters to fit: the rotation velocity (⁠|$V_{\mathrm{rot}}$|⁠), velocity dispersion (⁠|$\sigma$|⁠), radial velocity (⁠|$V_{\mathrm{rad}}$|⁠), systemic velocity (⁠|$V_{\mathrm{sys}}$|⁠), kinematic position angle (⁠|$PA_{\mathrm{kin}}$|⁠), kinematic centre (⁠|$x_{0,\mathrm{kin}},y_{0,\mathrm{kin}}$|⁠), i, and |$Z_0$|⁠. To limit the number of free parameters, we leave only |$V_{\mathrm{rot}}$|⁠, |$\sigma$|⁠, and |$PA_{\mathrm{kin}}$| of each ring free. For the remaining parameters, we fix the centre and inclination to the values returned by cannubi, and we fix |$Z_0$| to 460 pc (corresponding to |$q_0 = 0.2$|⁠). We also fix |$V_{\mathrm{rad}}$| to zero, meaning we assume no radial motion, following preliminary tests described below. For the systemic velocity, (⁠|$V_{\mathrm{sys}}$|⁠), we allow |$^{\mathrm{3D}}$|BAROLO to estimate this from the data cube, which results in |$V_{\mathrm{sys}}=25.76$| km s|$^{-1}$| from the [C ii] line at a redshift |$z=7.3065$|⁠, which is consistent with |$V_{\mathrm{sys}}$| derived from a single component Gaussian fit to the [C ii] spectrum and consistent with the reported uncertainty on the redshift from Table 1.

Table 1.

Table showing the properties of REBELS-25 from previous studies.

Property	Value	Reference
|$z_{\mathrm{[CII]}}$|	7.3065 |$\pm$| 0.0001	1,2
|$M_*$| (⁠|${\rm M}_{\odot }$|⁠)	|$8^{+4}_{-2} \times 10^9$|	2, 3
|$M_{\rm H_2,[CII]}$| (⁠|$\rm {\rm M}_{\odot }$|⁠)	|$5.1^{+5.1}_{-2.6} \times 10^{10}$|	4
SFR|$_{\mathrm{[CII]}}$| (⁠|${\rm M}_{\odot }\,\mathrm{yr}^{-1}$|⁠)	246 |$\pm$| 35	4
SFR|$_{\rm IR}$| (⁠|${\rm M}_{\odot }\mathrm\,\rm {yr}^{-1}$|⁠)	|$185 ^{+101}_{-63}$|	5
|$\rm SFR_{UV}({\rm M}_{\odot }\,\mathrm{yr}^{-1})$|⁠)	|$14 \pm 3$|	3
|$\rm {\it L}_{[CII]}$| (⁠|${\rm L}_{\odot }$|⁠)	1.7 |$\pm$| 0.2 |$\times 10^9$|	4
|$L_{[\rm IR]}$| (⁠|${\rm L}_{\odot }$|⁠)	|$7.1 ^{+3.6}_{-1.5} \times 10^{11}$|	6
|$S_{158\,\mu \rm m}$| (⁠|$\mu$|Jy)	|$260 \pm 22$|	5
FWHM|$_{\rm [CII]}$| (km s|$^{-1}$|⁠)	|$316 \pm 15$|	4
Conversion (kpc arcsec−1)	5.09552

References: [1] Schouws et al. (in preparation), [2] Bouwens et al. (2022), [3] Stefanon et al. (in preparation), [4] Hygate et al. (2023), [5] Inami et al. (2022), [6] Algera et al. (2024).

In Lelli et al. (2023), tests on mock data cubes find that |$^{\mathrm{3D}}$|barolo can reliably constrain intrinsic velocity dispersions that are greater than the instrumental dispersion. For cases where the intrinsic velocity dispersion is smaller than the instrumental dispersion, an exploration of the parameter space with the SPACEPAR function of |$^{\mathrm{3D}}$|barolo finds that the best-fitting velocity dispersion is pushed to the lower boundary of the parameter space. It is therefore reasonable to assume that we cannot constrain velocity dispersions lower than 7 km s|$^{-1}$| (the instrumental dispersion of the data). To be conservative, we have therefore set a lower limit to the allowed intrinsic velocity dispersion of the model (MINVDISP parameter) to 7 km s|$^{-1}$|⁠.

The results from the |$^{\mathrm{3D}}$|barolo fitting are given in Table 2 and shown in Figs 4–6. Whilst we left |$PA_{\mathrm{kin}}$| as a free parameter for each ring, we found that the variation in |$PA_{\mathrm{kin}}$| (of the order of |$\lt 10^{\circ}$|⁠) was negligible in comparison to the uncertainties (⁠|$\sim 30$| for each ring), and we therefore quote only the average |$PA_{\mathrm{kin}}$| across the disc in Table 1.

Table 2.

Table showing the [C ii] morphological and kinematic parameters derived from a 2D Sérsic model, cannubi and |$^{\mathrm{3D}}$|barolo, as described in the text.

Parameter	Value	Model used
n	|$1.3 \pm 0.2$|	2D Sérsic model
|$r_e$| (kpc)	|$2.2 \pm 0.3$|	2D Sérsic model
|$\rm PA_{\mathrm{morph}}$| (⁠|$^{\circ }$|⁠)	|$215 ^{+12}_{-31}$|	cannubi
|$\mathrm{PA}_{\mathrm{kin}}$| (⁠|$^{\circ }$|⁠)	|$218 \pm 6$|	|$^{\mathrm{3D}}$|barolo
i (⁠|$^{\circ }$|⁠)	|$25\pm 6$|	cannubi
|$V_{\mathrm{rot, max}}$| (km s|$^{-1}$|⁠)	|$374^{+86}_{-91}$|	|$^{\mathrm{3D}}$|barolo
|$\bar{\sigma }$| (km s|$^{-1}$|⁠)	|$33 ^{+9}_{-7}$|	|$^{\mathrm{3D}}$|barolo
|$V_{\mathrm{rot, max}}/\bar{\sigma }$|	|$11^{+6}_{-5}$|	|$^{\mathrm{3D}}$|barolo

To further test the robustness of our modelling, we also run |$^{\mathrm{3D}}$|BAROLO with all parameters left free. The resulting inclination (⁠|$i=32\pm 5$|º) is consistent within the uncertainties with the i fitted with CANNUBI, as is the fitted centre (within one beam FWHM of |$x_{0,\mathrm{morph}},y_{0,\mathrm{morph}}$|⁠). The fitted values for |$V_{\mathrm{rad}}$| are consistent with zero, within the uncertainties, for all five rings, and this is also true when a local normalisation is used. The fitted |$Z_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 (⁠|$V_{\mathrm{rot}}$|⁠, |$\sigma$|⁠, and |$PA_{\mathrm{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 |$^{\mathrm{3D}}$|BAROLO in Fig. 4 and the position velocity diagrams (PVDs) in Fig. 5, with a selection of channel maps in Fig. D1. We see a clear velocity gradient in the velocity field maps of Fig. 4 and a near-S-shape curve in the major axis PVD in Fig. 5, which are characteristic features of rotating discs.

Figure 4.

|$^{\mathrm{3D}}$|barolo fitting for REBELS-25. Emission is masked at 2|$\sigma _{\mathrm{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 all plots, the black cross and the ellipse mark the centre and the radial extent of the modelling, respectively. In the first two rows, we also show the position angle of the |$^{\mathrm{3D}}$|barolo model by the dashed line. In the second row, the grey dots give an indication of the separation of each ring along the velocity field (0.11 arcsec). 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.

Figure 5.

Major- (upper) and minor- (lower) axis position–velocity diagrams (PVDs) of REBELS-25, with black contours showing 2|$\sigma _{\mathrm{RMS}}$| and 4|$\sigma _{\mathrm{RMS}}$| emissions, and the dashed grey contours –2|$\sigma _{\mathrm{RMS}}$|⁠, where |$\sigma _{\mathrm{RMS}}=84\,\mu$|Jy^–1 beam and is equal to one standard deviation above the median value as calculated by |$^{\mathrm{3D}}$|barolo. The red contours illustrate the best-fitting model from |$^{\mathrm{3D}}$|barolo. In the major-axis PVD, the yellow markers show the rotation velocities of each ring projected along the line of sight.

In Table 2, we report the final kinematic parameters derived by |$^{\mathrm{3D}}$|BAROLO, including the maximum rotation velocity (⁠|$V_{\mathrm{rot, max}}$|⁠), average velocity dispersion (⁠|$\bar{\sigma }$|⁠) and the |$V_{\mathrm{rot, max}}/\bar{\sigma }$| ratio for REBELS-25. The |$\bar{\sigma }$| gives a measure of the overall disc turbulence, and |$V_{\mathrm{rot, max}}/\bar{\sigma }$| measures the ratio of ordered to turbulent motion. For this analysis, we assume that the |$V_{\mathrm{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.

We also plot the |$V_{\mathrm{rot}}$| and |$\sigma$| values of each ring as a function of radius in Fig. 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, |$V_{\mathrm{circ}}$|⁠. This is a reasonable assumption for REBELS-25 thanks to the high |$V_{\mathrm{rot}}/\sigma$| 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 |$V_{\mathrm{rot, max}}$| value as the maximum circular velocity of the system in subsequent calculations.

Figure 6.

The rotation velocity (⁠|$V_{\mathrm{rot}}$|⁠, dark blue) and velocity dispersion (⁠|$\sigma$|⁠, red) curves for REBELS-25. The size of the band indicates the estimated uncertainty from the |$^{\mathrm{3D}}$|BAROLO fit for |$\sigma$|⁠, and also propagating through the uncertainty on the inclination (⁠|$i=25\pm 6$|°) for |$V_{\mathrm{rot}}$|⁠. The dashed black line indicates the effective radius from sersic2d fitting.

The rotation velocity curve of REBELS-25 is shown in Fig. 6, and we note that the uncertainties in the rotational velocity also include the uncertainty in the inclination fitted by CANNUBI. As we have found a low inclination of |$25 \pm 6{^\circ}$|⁠, this corresponds to large changes in the intrinsic rotational velocity, as indicated by the error bars in Fig. 6. This rotation curve show the velocity rising quickly in the centre and then flattening to the maximum value of |$\sim$|370 km s|$^{-1}$| at around 2 kpc. Similar rotation curves are observed in disc galaxies at lower redshifts (e.g. Noordermeer et al. 2007; Di Teodoro et al. 2016; Lelli, McGaugh & Schombert 2016; Lelli et al. 2023). We also find a potential ‘bump’ feature in the rotational velocity curve within the first |$\sim$| 1 kpc, at a similar radius as the potential gas bar component observed in the [C ii] and dust continuum maps. Whilst the uncertainties are such that this feature is not certain, this bump persists also in the fits with |$^{\mathrm{3D}}$|BAROLO where all the parameters are left free, and could further 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 studies of galaxies at |$z\gt 4$| (e.g. Lelli et al. 2021, 2023; Rizzo et al. 2021; Tripodi et al. 2024; Roman-Oliveira, Rizzo & Fraternali 2024).

A decrease of |$\sigma$| with radius is also observed, as found in some local disc galaxies (e.g. Fraternali et al. 2002; Boomsma et al. 2008; Iorio et al. 2016; Bacchini et al. 2020). 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).

In the velocity dispersion maps in Fig. 4, we see some positive residuals within the radius modelled by |$^{\mathrm{3D}}$|BAROLO. To ensure that we are not underestimating the velocity dispersion, we run tests where we set the intrinsic velocity dispersion of the galaxy to be higher than 80 km s|$^{-1}$|⁠, as detailed in Appendix C. In these tests, we see that the contours of the PVDs do not follow the data well, in comparison to the fiducial model, and that the residuals of the velocity dispersion are now biased negative.

Additionally, some pixels in the observed moment-2 map close to the edge of the galaxy show very low velocity dispersions that are consistent with MINVDISP, which is set to the instrumental dispersion (7 km s|$^{-1}$|⁠). Whilst many of these pixels lie beyond the radius fitted by |$^{\mathrm{3D}}$|BAROLO (indicated by the black ellipse in Fig. 4), some of these pixels are contained within the outer ring. To investigate the effect of these pixels, we test a tighter mask with |$^{\mathrm{3D}}$|BAROLO that excludes all pixels detected in fewer than three channels, which by definition removes most of the pixels with the lowest velocity dispersion values. The details of these tests are discussed in Appendix C. From these tests, we find consistent results within the uncertainties, albeit with a slightly higher velocity dispersion of 38 km s|$^{-1}$| and a |$V_{\mathrm{rot, ~max}}/\bar{\sigma }$| of |$\sim 9$|⁠. However, from an exploration of the parameter space with SPACEPAR, the kinematic parameters are not robustly constrained in the outermost rings.

4.2.1 Non-circular features in the |$^{\mathrm{3D}}$|BAROLO models

Before drawing further conclusions based on the above modelling, we first investigate any features in the velocity field that are not well fit by |$^{\mathrm{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). First, 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. 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 non-axisymmetric structures within galaxies (e.g. Buta 1987; de Naray, Zagursky & McGaugh 2009; Wu et al. 2021). Secondly, in the velocity field map and the residual map of Fig. 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 channel maps at |$-234$| and |$-280$| km s|$^{-1}$| (Fig. D1). 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 |$^{\mathrm{3D}}$|BAROLO fits that incorporate radial motions. As mentioned in Section 4.1, we find |$V_{\mathrm{rad}}$| consistent with zero given the uncertainties, and we also find that these models were not able to reproduce similar high velocity regions.

Additionally, we see some faint emission beyond the rotating disc model in the PVDs. For example, in the upper half of the minor-axis PVD, and also in the upper 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 |$z = 0$| (e.g. Jorsater & Van Moorsel 1995; Fraternali et al. 2001; 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 fig. 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 |$\sim 2\sigma _{\mathrm{RMS}}$|⁠.

Overall, 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 [C ii] component

In Hygate et al. (2023), a potential secondary [C ii] component is identified in the range |$\sim +250$| to |$+650$| km s|$^{-1}$| in the [C II] 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 |$\sim 380$| km s|$^{-1}$|⁠) or a minor merger (modelled as a narrow component with width |$\sim 100$| km s|$^{-1}$|⁠). We have searched for this secondary [C ii] 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 arcsec and attempted to fit a main and secondary [C ii] 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 [C ii] 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 |$\sim 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 disc, as observations 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 [C ii] ALMA-ALPINE data in Jones et al. (2021), and also the PV Split method described in Rizzo et al. (2022).

We 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 |$^{\mathrm{3D}}$|BAROLO. We also make use of the morphological parameters fitted with CANNUBI, and we use |$^{\mathrm{3D}}$|BAROLO’s 2D fit capabilities to an inputted velocity field for the third and fifth criteria.

In Jones et al. (2021), typically at least three of these criteria are met for the ALPINE galaxies classified as rotating discs. For REBELS-25, we find that all five criteria are met:

• We find that a slice along the kinematic major axis reveals a slope in position-velocity space with |$\gt 3\sigma$| significance, implying a smooth velocity gradient.

• For each of the five fitted rings, we find |$V_{\mathrm{rot}}\gt \sigma$| (see Fig. 6).

• We find that the kinematic centre derived by a 2D fit with |$^{\mathrm{3D}}$|BAROLO and the peak in the velocity dispersion map are consistent within the uncertainties.

• The |$\rm PA_{\mathrm{morph}}$| from CANNUBI and the average |$\rm PA_{\mathrm{kin}}$| across all rings are consistent within the uncertainties, as expected for an ideal rotating disc.

• We find that the kinematic centre fitted by |$^{\mathrm{3D}}$|BAROLO and the morphological centre fitted by CANNUBI are consistent within error.

From 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.

In 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, van Albada & Sancisi 2001) in Roman-Oliveira et al. (2023). Using the PVDs returned by |$^{\mathrm{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.

As 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 |$^{\mathrm{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, |$M_{\mathrm{dyn}}$|⁠, enclosed within a radius, R, is given by:

$$\begin{eqnarray} M_{\mathrm{dyn}}(\lt R) =\frac{{k(R)} R v_{\mathrm{circ}}^2(R)}{G}, \end{eqnarray}$$

(1)

where |$k(R)$| is the virial coefficient at radius R. Price et al. (2022) use Sérsic models to derive the virial coefficients for a range of n and intrinsic axis ratios (q). 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 |$k_{\mathrm{tot}}=1.8$|⁠, which is the virial coefficient for |$n=1$| and |$q=0.2$| (typically quoted for thin discs, e.g. Wel et al. 2014) at |$R=r_e$|⁠. We take the average |$V_\mathrm{rot}$| in the outer two rings (⁠|$=359^{+93}_{-95}$| km s|$^{-1}$|⁠), which results in |$M_{\mathrm{dyn, ~tot}} = (1.2^{+1.0}_{-0.6}) \times 10^{11} \mathrm{ {\rm M}_{\odot }}$|⁠. This is larger but within the uncertainties of the value derived in Hygate et al. (2023), which used the low resolution LP data. If we instead take our upper limit on the disc thickness (710 pc), which corresponds to a |$q_0$| of |$\sim 0.32$|⁠, we find a dynamical mass of |$\sim 1.4 \times 10^{11} \mathrm{ {\rm M}_{\odot }}$|⁠, which is within the uncertainties of the quoted |$M_{\mathrm{dyn, tot}}$| above.

Assuming that the contribution of dark matter to the mass budget within |$r_e$| is negligible (as found at high-z in e.g. van Dokkum et al. 2015; Price et al. 2016; Wuyts et al. 2016; Genzel et al. 2017; Gelli et al. 2020, although see de Graaff, Pillepich & Rix 2024 where |$z\gt 6$| galaxies are found to be dark matter-dominated), we can estimate the total gas mass from |$M_{\mathrm{gas}}=M_{\mathrm{dyn}}-M_*$|⁠. Adopting a stellar mass of |$M_* = 8^{+4}_{-2}\times 10^9\, \mathrm{ {\rm M}_{\odot }}$| results in |$M_{\mathrm{gas}}=(1.1^{+1.0}_{-0.7})\times 10^{11}\, \mathrm{{\rm M}_{\odot }}$|⁠, with a gas fraction of |$f_{\mathrm{gas}} = M_{\mathrm{gas, ~tot}}/(M_{\mathrm{gas, ~tot}}+M_*) = 0.93^{+0.03}_{-0.16}$|⁠. However, we note that if we adopt the stellar mass derived in Topping et al. (2022) using a non-parametric star formation history (⁠|$M_*=19^{+5}_{-8} \times 10^9\, \mathrm{ {\rm M}_{\odot }}$|⁠), |$M_{\mathrm{gas,tot}}=(1.0^{+1.0}_{-0.6})\times 10^{11}\, \mathrm{ {\rm M}_{\odot }}$| and |$f_{\mathrm{gas}} = 0.84^{+0.11}_{-0.22}$|⁠.

With these constraints on the gas mass, we can also derive a conversion factor between the [C ii] luminosity and |$M_{\mathrm{gas}}$|⁠, |$\alpha _{\mathrm{[CII]}}$|⁠, assuming that the [C ii] 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 [C ii] spectrum, we find |$L_{\mathrm{[CII]}}=(1.8\pm 0.2)\, \times 10^9 \,{\rm L}_{\odot }$|⁠, which is consistent within the uncertainties with the value derived from the low resolution LP data in Hygate et al. (2023). This results in |$\alpha _{\mathrm{[CII]}}=62^{+68}_{-40} \,\mathrm{{\rm M}_{\odot }/L_{\odot }}$|⁠. This value is a factor of |$\sim 10$| times higher than the median value found for dusty star-forming galaxies at |$z\sim 4.5$| in Rizzo et al. (2021), and a factor of |$\sim$| 2 times higher than the value of 30 |$\mathrm{{\rm M}_{\odot }/L_{\odot }}$| 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 we assume a negligible contribution from dark matter, which may not be appropriate at |$z\gt 6$| (see de Graaff et al. 2024). Studies with a more statistically significant sample would be necessary to investigate |$\alpha _{\mathrm{[CII]}}$| values for star forming galaxies at |$z\gt 6$|⁠.

From the properties derived by |$^{\mathrm{3D}}$|BAROLO, and an estimate for the gas mass, it is also possible to determine the Toomre parameter (Toomre 1964), Q, of REBELS-25. This Q parameter is used to measure the local instability of a disc, with |$Q\lt 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) and the exponential scale length of the disc, we find an average Q of |$0.4^{+1.5}_{-0.3}$| over the radial extent of the [C ii] emission, suggesting it may be locally unstable. However, some studies at |$z\sim 0$| have found that Q does not show any correlation with fragmentation or star formation (e.g. Leroy et al. 2008; Romeo & Falstad 2013; Elmegreen & Hunter 2015; Romeo, Agertz & Renaud 2023), and a recent study also indicates that high-z galaxies are not as locally unstable as expected from the Q 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 (⁠|$V_{\mathrm{rot, max}}/\bar{\sigma }=11^{+6}_{-5}$|⁠) and low overall turbulence (⁠|$\bar{\sigma }=33^{+9}_{-7}$| km s|$^{-1}$|⁠). In Kohandel et al. (2023), galaxy discs are classed according to their |$V/\sigma$| ratios, with |$V/\sigma \gt 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 |$z=7.31$|

In Section 3.1, we found that the [C ii] emission of REBELS-25 is well-fit by a near-exponential disc model, similar to what is found for galaxies in the local Universe and for the dust continuum in some submillimetre galaxies (SMGs) at |$z\sim 1-3$| (e.g. Barro et al. 2016; Hodge et al. 2016; Fujimoto et al. 2018; Rivera et al. 2018; Gullberg et al. 2019). This is in contrast to the clumpy and irregular morphologies observed and anticipated at |$z\gt 6$| (see Section 1). We also find that the gas disc of REBELS-25 is likely thinner than the |$Z_0 \sim 1$| kpc derived from galaxies at |$z\sim 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 |$\Lambda$|CDM galaxy formation simulations (e.g. Brook et al. 2012; Bird et al. 2013; Park et al. 2019; Bird et al. 2021; Renaud et al. 2021).

We 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 [C ii] and dust continuum emission at |$\sim 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 few 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), including gas bars in gas-rich baryon-dominated discs (Bland-Hawthorn et al. 2024). Barred galaxies are common in the nearby Universe (e.g. Knapen, Shlosman & Peletier 2000; Whyte et al. 2002; Laurikainen, Salo & Rautiainen 2004; Buta et al. 2015), but the fraction of barred galaxies has been found to decrease between |$z =$| 0 to 1 (e.g. Sheth et al. 2008). However, some recent observations have identified barred galaxies at |$z=1-4.4$| (Simmons et al. 2014; Hodge et al. 2019; Guo et al. 2023; Huang et al. 2023; Smail et al. 2023; Tsukui 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.

We additionally see some very tentative evidence of arm-like or clump-like extended features at |$\sim$| 3 kpc to the North and South of the peak of emission in the dust continuum map (Fig. 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 programmes (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 ionized gas content may shed more light on whether or not any bar, spiral arm and/or clump structures are present.

Some 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 [C ii] 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 |$z=6.072$| galaxy (Fujimoto et al. 2024).

We also find that REBELS-25’s dynamical mass is a factor of |$\sim 15$| 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 preparation). If we instead adopt the stellar mass derived in Topping et al. (2022)(⁠|$M_*=19^{+5}_{-8} \times 10^9 \mathrm{{\rm M}_{\odot }}$|⁠), |$M_{\mathrm{dyn}}/M_*$| becomes |$\sim 6$|⁠. Graaff et al. (2024) have recently reported |$M_{\mathrm{dyn}}/M_*$| values as high as 30 in high-z galaxies, whereas studies of |$z\sim 2$| galaxies find |$M_{\mathrm{dyn}}/M_* \lesssim 4$| (Maseda et al. 2013; Wuyts et al. 2016). 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 |$z\sim 0$|⁠, gas kinematics of galaxies are typically traced by warm ionized gas tracers, such as H|$\alpha$|⁠. These warm ionized 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.

As discussed in Section 1, the [C ii] 157.7 |$\mu$|m emission line is a useful, strong-line tracer of cold gas kinematics. However, observations of this line become increasingly difficult at |$z \lesssim 3.5$| from ground-based instruments due to the poor transmission through Earth’s atmosphere. For this reason, we compile literature results of both high resolution [C ii] kinematics studies, typically at |$z \gt 3.5$|⁠, and high resolution [CI] and CO kinematics studies, typically at |$0.5 \lesssim z \lesssim 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 |$\gtrsim 3$| independent resolution elements across the semimajor 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.

Rizzo 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 |$z = 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 |$\gtrsim 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 |$M_*$| 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).

Parlanti et al. (2023) collected and analysed all available Band 4, 5, and 6 ALMA archival observations of [C ii] (and [O iii]88|$\mu$|m) emission in galaxies at |$z \gt 3.5$| with angular resolution < 1.5 arcsec 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 |$^{\mathrm{3D}}$|BAROLO for their kinematic analyses, and Neeleman et al. (2020) use QUBEFIT but find that |$^{\mathrm{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 |$^{\mathrm{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 galaxy 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. (2023) are consistent for HZ4, and therefore for HZ4 and the remaining three galaxies (HZ9, J1211, and VR7) selected from Parlanti et al. (2023), we use the kinematic properties derived therein.

To increase the size of our comparison sample, we also consider starburst, quasar host, submillimetre, and gravitationally lensed galaxies that do not show signs of merging activity from Hodge et al. (2012), Rizzo et al. (2020), Fraternali et al. (2021), Rizzo et al. (2021), Amvrosiadis et al. (2023), Lelli et al. (2023), Liu et al. (2023), Neeleman et al. (2023), Roman-Oliveira et al. (2023), Fujimoto et al. (2024), 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 (⁠|$z\gt 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 submillimetre galaxy at |$z=4.6$| with high resolution [C ii] ALMA observations initially analysed by Tadaki et al. (2018), who found a rotating disc with |$V_{\mathrm{rot, max}} \sim 227$| km s|$^{-1}$| and |$\bar{\sigma } \sim 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), Smit et al. (2018), Rybak et al. (2019), Pope et al. (2023), Tamura et al. (2023), 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.

Overall, our selection results in a total of 37 objects with [C ii], [C i], or CO high-resolution kinematic studies from the literature. We plot the stellar mass and SFR for the selected galaxies in Fig. 7 at three different redshift ranges (⁠|$0.5\le z \lt 2.5$|⁠, |$2.5\le z \lt 5$|⁠, |$5\le z \lt 7.5$|⁠). Due to the limited sensitivities and resolution of observations at high redshift, we see that our selection is biased towards massive (⁠|$M_* \gtrsim 10^9 \,\mathrm{ {\rm M}_{\odot }}$|⁠) main-sequence or starburst (SFR|$\gtrsim 10 \mathrm{ {\rm M}_{\odot }}\mathrm{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 Fig. 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.

Figure 7.

Distribution of REBELS-25 and a few comparison samples in the |$M_*$|–SFR plane. We split the sample into three redshift bins and show the main-sequence line from Schreiber et al. (2015) for |$z\sim$| 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|$\sigma$| uncertainty on the main-sequence lines, and the green shaded areas show the region occupied by starbursts according to a threshold of |$4 \times$| the main-sequence lines (e.g. Rodighiero et al. 2011). Galaxies observed with [C ii] 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.

5.3 The cosmic evolution of turbulence within galaxies

In Fig. 8, we plot the average velocity dispersion, |$\bar{\sigma }$|⁠, for REBELS-25 and the literature sources as a function of redshift. We plot also the analytical fit to the large IFU survey KMOS3D at |$z \sim 1-3$| (Übler et al. 2019) for both ionized 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 that turbulence increases with redshift for both warm ionized gas and cold neutral gas tracers. Wisnioski et al. (2015) find that this evolution at |$0\le z \le 3$| 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|$\alpha$| emission of a disc galaxy of mass |$\log (M_*/{\rm M}_{\odot }) \sim 10$| with |$Q=1$| by the shaded grey area in Fig. 8, and we extrapolate this to higher redshifts with the fainter grey region. Additionally, we plot the predictions from IllustrisTNG50 (Nelson et al. 2019) simulations for H |$\alpha$| emission (Pillepich et al. 2019) and from SERRA simulations (Pallottini et al. 2022) for [C ii] emission (Kohandel et al. 2023) for different mass bins. The predictions from the SERRA simulations are specifically for [C ii] emission from mock ALMA observations at different inclinations of 3218 simulated galaxies at |$8\le \log (M_*/{\rm M}_{\odot })\le 10.3$| from |$4\lt z\lt 9$|⁠, with full details of their kinematic analyses given in Kohandel et al. (2023).

Figure 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 [C i] (blue) and [C ii] (orange) observations. The markers are the same as in Fig. 7. Empty markers show where one or more of the kinematic parameters (⁠|$V_{\mathrm{rot}}$| or |$\sigma$|⁠) is given as an upper or lower limit. The grey-shaded region shows the predicted evolution of |$\sigma$| from the semi-analytic model of Wisnioski et al. (2015), which we note is derived for |$0\le z \le 3$| observations, and so we extrapolate this to |$z\gt 3$|⁠. The black and grey lines show fits to observations at |$z\sim 0.6-2.6$| (solid lines) for ionized and molecular gas, respectively, and the extrapolations of these fits to higher redshifts (dashed lines). 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\sigma$| 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.

We find that the average velocity dispersion of REBELS-25 agrees better with the SERRA simulations than the other predictions at high-z. This likely follows from the fact that the Wisnioski et al. (2015) and TNG50 predictions are for warm ionized 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 ionized gas at |$z\sim 1-3$|⁠, although this offset may evolve with redshift or other galaxy properties. Similarly, Kohandel et al. (2023) find that |$\sigma$| measured from H |$\alpha$| emission is |$\sim 2\times$| larger than |$\sigma$| as measured from [C ii] 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 (⁠|$\bar{\sigma } = 33^{+9}_{-7}$| km s|$^{-1}$|⁠) is comparable to |$z\lt 3$| galaxies.³

To investigate why this may be the case, we consider the processes responsible for driving the turbulence, and hence the gas velocity dispersion, within galaxies. First, for REBELS-25, the derived |$f_{\mathrm{gas}}=0.93^{+0.03}_{-0.16}$| 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

$$\begin{eqnarray} \bar{\sigma } = \frac{1}{\sqrt{2}} V_{\mathrm{rot, max}} f_{\mathrm{gas}}(z). \end{eqnarray}$$

(2)

This would predict a |$\bar{\sigma }$| of |$\sim 250$| km s|$^{-1}$| (for a warm ionized gas tracer, so |$\sim 120$| km s|$^{-1}$| for the [C ii] emission according to Kohandel et al. 2023). However, as discussed in Wisnioski et al. (2015), the redshift evolution of the necessary parameters to define |$f_{\mathrm{gas}}(z)$|⁠, such as the depletion time and the sSFR, are uncertain at |$z\gtrsim 3$|⁠, and the extrapolation of this model may therefore not be applicable at the high redshift of REBELS-25.

Secondly, 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. (2023), from which we have selected seven main-sequence 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. This is perhaps due to the lower mass of the galaxies in the Parlanti et al. (2023) 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, we have plotted the predicted velocity dispersion according to both of these Krumholz et al. (2018) models in Fig. 8 (for cool atomic or molecular gas). Here, we see that, as with the galaxies in the Rizzo et al. (2021) and Roman-Oliveira et al. (2024) samples, the models considering only feedback effects better agree with the velocity dispersion obtained for REBELS-25. For REBELS-25, the stellar mass quoted in Table 1 is comparable to those in the Parlanti et al. (2023) 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.

Finally, 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 ionized 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 ionized gas. In high-z 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. 2024; Übler et al. 2024; Fujimoto et al. 2024). We note that upcoming analysis of ALMA Cycle 9 Band 8 high resolution observations of the [O iii]88|$\mu$|m line for REBELS-25 (Rowland et al. in preparation, ID:2022.1.01324.S PI H. Algera) 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, |$V_{\mathrm{rot, max}}/\bar{\sigma }$|⁠. As with Fig. 8, we plot this ratio as a function of redshift for REBELS-25, the literature sources and various predictions in Fig. 9. Here, we see that REBELS-25, and most of the galaxies in our comparison sample, have |$V_{\mathrm{rot, max}}/\bar{\sigma } \gt 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 ionized gas tracers, which are likely to result in higher velocity dispersions and therefore lower |$V_{\mathrm{rot, max}}/\bar{\sigma }$|⁠.

In Fig. 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 (M_*/{\rm M}_{\odot })\gt 9$| sample.⁴ This is also the case for most of the galaxies in our comparison sample at |$z\gt 4$|⁠, although we note that the galaxies from Parlanti et al. (2023) have large uncertainties such that they are consistent with most of the predictions shown here.

Figure 9.

As with Fig. 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.

The 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 Fig. 9 are for kinematics as traced by [C ii] emission, whereas the TNG50 predictions are for H|$\alpha$| 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 other studies, such as 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 (M_*/{\rm M}_{\odot }) \sim 9$|⁠. Other observations at |$z\lt 3$| also indicate a mass dependence on the degree of rotational support, for example Tiley et al. (2021) find that stellar mass, rather than redshift, most strongly correlates with the disc fraction amongst star-forming galaxies at |$z\sim 1.5$|⁠, observing only a modest increase in the prevalence of discs between |$z \sim 1.5$| and |$z \sim 0.04$| at fixed stellar mass. Similar findings are presented in, for example, Simons et al. (2017) and Wisnioski et al. (2019) for |$z\lesssim 2$| observations. 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 Fig. 7, there is a lack of high resolution observations of cold gas kinematics in low-mass (⁠|$M_*\lt 10^9 \mathrm{ {\rm M}_{\odot }}$|⁠) galaxies at high-z to confidently validate these hypotheses.

To investigate this mass dependence further, we therefore also make comparisons to ionized gas kinematic studies at a range of redshifts. At lower redshifts, we compare with results from the KMOS|$^{\mathrm{3D}}$| survey from Wisnioski et al. (2019). At high redshift, we compile results from a recent study of the ionized gas kinematics with JWST NIRSpec Multi-Object Spectroscopy (MOS) observations (Graaff et al. 2024), studies of individual galaxies with NIRSpec IFU observations (including a lensed galaxy at |$z=6.072$| named the Cosmic Grapes, from Fujimoto et al. 2024, ALESS073.1 at |$z=4.76$| from Parlanti et al. 2024, and GN-z11 at |$z=10.60$| from Xu et al. 2024), a study of a |$z=5.3$| galaxy using JWST grism spectroscopy from Nelson et al. (2023), and a study of the |$z=2.259$| galaxy J0901 using SINFONI/VLT H |$\alpha$| observations from Liu et al. (2023). 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 Graaff et al. (2024) 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 Graaff et al. (2024) are amongst the lowest mass galaxies studied kinematically at |$z\gt 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 ([C ii], Lelli et al. 2021) and with JWST (H|$\alpha$|⁠, Parlanti et al. 2024) observations at high resolution, as has the ‘Cosmic Grapes’ |$z=6.072$| galaxy from Fujimoto et al. (2024). In Fig. 10, we plot |$V_{\mathrm{rot, max}}/\bar{\sigma }$| against stellar mass for the same galaxies as in Fig. 9, along with the aforementioned ionized gas observations with magenta markers and lines.

Fig.10 appears to show a tentative positive correlation between |$V_{\mathrm{rot, max}}/\bar{\sigma }$| and |$M_*$|⁠. For the warm ionized 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 |$V_{\mathrm{rot, max}}/\bar{\sigma }$|⁠. 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-z (e.g. Abdurro’uf et al. 2023; Giménez-Arteaga et al. 2023).

Figure 10.

As with Figs 8 and 9, but instead the ratio of ordered-to-random motion is plotted as a function of stellar mass. The magenta markers are from recent studies of ionized gas kinematics at |$z\gt 2$|⁠, and the magenta lines are from H |$\alpha$| observations of 739 galaxies from the KMOS|$^{\mathrm{3D}}$| survey, averaged in three redshift bins, from Wisnioski et al. (2019). The black vertical dashed line indicates a cut-off stellar mass at |$\sim 10^9 \mathrm{ {\rm M}_{\odot }}$| described in Gurvich et al. (2022) as a threshold mass where rotationally supported discs (with |$V_{\mathrm{rot}}/\sigma \gt 1$|⁠) can form.

6 SUMMARY AND CONCLUSIONS

In this paper, we have presented follow-up high-resolution (⁠|$\sim$|710 pc) ALMA [C ii] and dust (⁠|$\sim 150\, \mu$|m) continuum observations of REBELS-25; a massive star-forming galaxy at |$z=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 [C ii] emission in this galaxy is well fit by a near-exponential disc profile, with a Sérsic index of |$1.3\pm 0.2$|⁠. These UV clumps that are offset from the [C ii] 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 that it has a low inclination (⁠|$i= 25\pm 6{^\circ }$|⁠) and potentially a very thin disc (⁠|$Z_0\lt 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 [C ii] and dust continuum emission.

We find the [C ii] kinematics of REBELS-25 are well explained by a rotation-dominated disc using the 3D tilted ring fitting tool, |$^{\mathrm{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-fitting rotating disc model with |$^{\mathrm{3D}}$|BAROLO reveals a low overall turbulence (⁠|$\bar{\sigma }=33^{+9}_{-7}$| km s|$^{-1}$|⁠) and a high ratio of |$V_{\mathrm{rot, max}}/\bar{\sigma } = 11^{+6}_{-5}$|⁠. 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.

Additionally, we find a total dynamical mass of (⁠|$1.2^{+1.0}_{-0.6})\times 10^{11} \mathrm{ {\rm M}_{\odot }}$| for this galaxy. This indicates that REBELS-25 is a highly gas-dominated galaxy, with a gas mass of |$M_{\mathrm{gas, tot}}=(1.1^{+1.0}_{-0.7})\times 10^{11} \mathrm{{\rm M}_{\odot }}$| (if we take |$M_* = 8^{+4}_{-2}\times 10^9 \mathrm{ {\rm M}_{\odot }}$|⁠). However, the total stellar mass of this galaxy is likely very uncertain thanks to dust obscuration, meaning that our estimates of |$M_{\mathrm{gas}}$|⁠, |$f_{\mathrm{gas}}$|⁠, and |$\alpha _{\mathrm{[CII]}}$| are 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.

An increasing number of dynamically cold discs have been identified in the high-z Universe, although these observations are often not sufficiently resolved to confidently distinguish between mergers and discs. With these sub-kpc resolution observations, REBELS-25 is amongst the most distant robustly confirmed cold discs observed to date. This finding of a very distant (⁠|$z=7.31$|⁠), very dynamically cold (⁠|$V_{\mathrm{rot, max}}/\bar{\sigma } = 11^{+6}_{-5}$|⁠) 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 |$z\gt 3$|⁠. By comparing to other cold gas kinematics studies of |$z\gt 0.5$| galaxies with observational data of similar quality, we find that dynamically cold discs seem to be more common in the high-z Universe than predictions based on warm ionized gas tracers (e.g. Pillepich 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 SERRA cosmological simulations (Pallottini et al. 2022) based on mock observations of [C ii] 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 [O iii]88 |$\mu$|m emission (Rowland in preparation, ID:2022.1.01324.S PI H. Algera) and JWST observations of [O iii]5007|$\mathring{\rm A}$| kinematics (GO-6036 PI J. Hodge, and GO-6480 PI S. Schouws) will enable a comparison between cold and warm ionized gas tracers.

REBELS-25 is a massive (⁠|$M_* = 8^{+4}_{-2}\times 10^9 \mathrm{ {\rm M}_{\odot }}$|⁠) 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 (Fig. 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 ionized gas kinematic studies of fainter, lower mass galaxies at high-z now becoming more feasible (e.g. Graaff et al.2024).

Overall, 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-z 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 |$z\sim 6-8$| (Phillips et al. in preparation).

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). MR 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 programme. AF acknowledges support from the ERC Advanced Grant INTERSTELLAR H2020/740120. IDL acknowledges funding support from ERC starting grant 851622 DustOrigin. MA 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. We also thank the anonymous referee for their useful comments and suggestions, which have improved the manuscript.

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.

Footnotes

1

https://casaguides.nrao.edu/index.php/Automasking_Guide

2

https://www.filippofraternali.com/cannubi

3

We add that taking the higher |$\bar{\sigma }$| value obtained from the tighter mask described in Appendix C does not change the conclusions drawn from these findings.

4

We note that using the tighter mask described in Appendix C, which results in a lower |$V_{\mathrm{rot, max}}/\bar{\sigma }$| of |$\sim 9$|⁠, does not change our main conclusions.

APPENDIX A: TESTS USING DIFFERENT [C ii] LINE CUBES

In Table A1, we list the properties of the ALMA [C ii] 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.

Table A1.

Table showing some of the observation details of the different [C II] line cubes tested, and the resulting ratio of ordered to random motion (⁠|$V_{\mathrm{rot}}/\bar{\sigma }$|⁠).

Measurement set	Channel width	Weighting	Robust	Beam	RMS	|$V_{\mathrm{rot, ~max}}/\bar{\sigma }$|
	(km s|$^{-1}$|⁠)			(arcsec|$\times$|arcsec)	(⁠|$\mu$|Jy beam−1 channel−1)
Cycle 8	5	Natural	–	0.14 |$\times$| 0.13	210	10
Cycle 8	10	Natural	–	0.14 |$\times$| 0.13	150	9
Cycle 8	15	Briggs	0.5	0.12 |$\times$| 0.10	140	11
LP |$+$| Cy 8, JvM-corr	20	Natural	–	0.15 |$\times$| 0.14	12|$^{*}$|	8

Note. Note that the residuals for the JvM-corrected cube are scaled by a factor |$\epsilon =0.18$|⁠.

For 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 [C ii] flux from the concatenated dataset that is |$\sim 3.4$| times higher than the [C ii] 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, |$\epsilon$|⁠, 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.

As in Posses et al. (2024), we find that cleaning down to a threshold of |$1\sigma _{\mathrm{RMS}}$| minimizes 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 2 yr 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. We then apply the JvM correction with |$\epsilon = 0.18$|⁠, as determined from the ratio of the CLEAN and dirty beam volumes. Applying this correction results in a consistent total [C II] 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 |$^{\mathrm{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.

We run some initial morphology fitting tests for each data set with Sersic2D and CANNUBI. Notably, we find a Sérsic index of |$\sim 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 |$^{\mathrm{3D}}$|BAROLO. For these kinematic fits, we leave only |$V_{\mathrm{rot}}$|⁠, |$\sigma$|⁠, and PA|$_{\mathrm{kin}}$| as free parameters, as 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 |$V_{\mathrm{rot, ~max}}/\bar{\sigma }$| 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 Fig. 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 |$V_{\mathrm{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.

Figure A1.

The major-axis (top) and minor-axis (bottom) PVDs for the best-fitting |$^{\mathrm{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 |$\sigma _{\mathrm{RMS}}$| as calculated from |$^{\mathrm{3D}}$|BAROLO. Negative contours for the data at −12 and −6 |$\sigma _{\mathrm{RMS}}$| are indicated by the dashed grey contours. Note, however, that this |$\sigma _{\mathrm{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 |$^{\mathrm{3D}}$|BAROLO and their fitted rotation velocities along the line of sight. The asymmetry in the model of the minor-axis PVD is the result of a changing PA across the disc which is preferred for this fit.

APPENDIX B: CANNUBI CORNER PLOTS

In Fig. B1, we show the posterior distributions obtained by applying CANNUBI to our [C II] 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.3410, +1:44:31.153.

Figure B1.

A triangle plot of the posterior distributions for the morphological parameters fitted by CANNUBI.

APPENDIX C: TESTS ON THE GAS VELOCITY DISPERSION

As discussed in the main body of the paper, we see some positive residuals in the moment-2 maps from our best-fitting kinematic model in Fig. 4. We note that these positive residuals do not seem to be axisymmetric, and it is therefore unlikely that they are representative of the overall disc dynamics. One can look at low-redshift galaxies with well-established rotation, where similar positive residuals are also very clear (e.g. the MHONGOOSE survey, de Blok et al. 2024). Despite these residuals, we see that the model PVDs follow the data contours well in Fig. 5. However, if we force the velocity dispersion to be higher, this is no longer the case. To illustrate this, we show the major-axis PVD for a model where we set the intrinsic velocity dispersion of the galaxy to 80 km s|$^{-1}$| in the left-hand side of Fig. C1. Here, we see that this model overestimates the width of the PVD. Additionally, we also show the residual in the moment-2 map (Fig. C2), where we see that the residuals are strongly biased negative.

Figure C1.

Two different major axis position velocity diagrams of REBELS-25, with the data plotted in blue and the models plotted in red, with the contours plotted at the same levels as in Figs 5 and A1. For the model on the left, we have fixed the intrinsic velocity dispersion of the model to 80 km s|$^{-1}$|⁠, but otherwise we follow the same fitting procedure described in Section 4.1. The plot on the right is our fiducial model, where the intrinsic velocity dispersion derived is |$33^{+9}_{-7}$| km s|$^{-1}$|⁠. Here, we see that the contours of our fiducial model better follow the data. In the higher |$\sigma$| model, the widths of the model contours, particularly at high SNR, overestimate the width of the data.

Figure C2.

Dispersion velocity maps (observed: left, model: centre, residuals: right) from |$^{\mathrm{3D}}$|BAROLO fitting for REBELS-25 where we set the intrinsic dispersion velocity to 80 km s|$^{-1}$|⁠.

In addition to these positive residuals in the observed map of the velocity dispersion, we also see some pixels close to the outer edge of the galaxy with |$\sigma$| values that approach the instrumental dispersion, which we have set as the minimum observed velocity dispersion (MINVDISP parameter). These are mostly low SNR pixels that are detected in fewer than three channels in the mask produced by |$^{\mathrm{3D}}$|BAROLO. Many of these pixels lie beyond the extent of the outer ring (black ellipses in Fig. 4). However, some of these pixels still remain within the outer ring, with 26 per cent of pixels in the fifth ring detected in fewer than three channels. We note that the effect of these pixels is likely limited due to the use of an azimuthal averaging, however we investigate if these pixels may be leading to an underestimation of the velocity dispersion in the outer ring by applying an additional step to the mask produced by |$^{\mathrm{3D}}$|BAROLO. To do this, we take the original mask produced by |$^{\mathrm{3D}}$|BAROLO from the SNRCUT and GRWOTHCUT parameters described in Section 4.1, and then we also mask out any pixels that are detected in fewer than three channels. This consequently excludes |$\sim 66~{{\ \rm per\ cent}}$| of pixels with |$\sigma \lt 45$| km s|$^{-1}$|⁠. We then run a fit with |$^{\mathrm{3D}}$|BAROLO again, using this new, tighter mask, and find results that are consistent, within the uncertainties, with our best-fitting fiducial model presented in the main body of the paper (⁠|$\bar{\sigma }=38$| km s|$^{-1}$| and |$V_{\mathrm{rot,~ max}} = 354$| km s|$^{-1}$|⁠). However, when we explore the parameter space of this fit with SPACEPAR (Fig. C3), we see that |$V_{\mathrm{rot}}$| is not well constrained in the outer ring. Therefore, we can conclude that the pixels detected in fewer than three channels are not affecting the derived velocity dispersion, but they are still useful for constraining the rotation of the gas, and so we still include these pixels in our fiducial model.

Figure C3.

The parameter space of the dispersion veloctiy (VDISP) and the rotation velocity (VROT) for each ring fitted with |$^{\mathrm{3D}}$|BAROLO, when we mask out all pixels detected in fewer than three channels, as described in the text. The best-fitting values are indicated with a cross and reported in the top label. The color scale represents the |$\chi ^2$| values, and the white contours correspond to percentage variations (2 per cent, 5 per cent, 10 per cent, and 30 per cent) in |$\chi ^2$| from the best-fitting value. In the fifth and outermost ring at |$\sim$| 0.49 arcsec, we see that the rotation velocity is not well constrained, from which we infer that there are too few pixels to reliably constrain the rotation in this ring. We also see that the velocity dispersion values are pushed to |$\lt 7$| km s|$^{-1}$|⁠.

Finally, we also test the QUBEFIT kinematic fitting tool, as described in Neeleman et al. (2021), and we find consistent results within the uncertainties. In particular, QUBEFIT finds a dispersion velocity of |$40^{+9}_{-8}$| km s|$^{-1}$|⁠.

APPENDIX D: CHANNEL MAPS

In Fig. D1, we show representative channel maps of the [C II] emission of REBELS-25 (top panels) and of the best-fitting kinematic model (bottom panels) obtained with |$^{\mathrm{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|$\sigma _{\mathrm{RMS}}$| and 4|$\sigma _{\mathrm{RMS}}$|⁠, with the negative contours at |$-2\sigma _{\mathrm{RMS}}$|⁠. The green cross shows the kinematic centre fitted by |$^{\mathrm{3D}}$|BAROLO.

Figure D1.

Representative channel maps of REBELS-25, with contours at 2 and 4 |$\sigma _{\mathrm{RMS}}$| for the data in blue (upper panels) and the model in red (lower panels). Negative contours at −2|$\sigma _{\mathrm{RMS}}$| are indicated by the grey dashed contours. As with Fig. 5, |$\sigma _{\mathrm{RMS}}=84\,\mu$|Jy beam|$^{-1}$| and is equal to one standard deviation above the median value as calculated by |$^{\mathrm{3D}}$|BAROLO. The centre of the galaxy as determined from CANNUBI and fixed in the kinematic fit is marked by a green cross.