The kinematics of massive high-redshift dusty star-forming galaxies

ABSTRACT

We present a new method for modelling the kinematics of galaxies from interferometric observations by performing the optimization of the kinematic model parameters directly in visibility space instead of the conventional approach of fitting velocity fields produced with the clean algorithm in real-space. We demonstrate our method on Atacama Large Millimeter/submillimeter Array (ALMA) observations of |$^{12}$|CO (2–1), (3–2), or (4–3) emission lines from an initial sample of 30 massive 850 |$\mu$|m-selected dusty star-forming galaxies with far-infrared luminosities |$\gtrsim$||$\, 10^{12}$| L|$_{\odot }$| in the redshift range |$z \sim$| 1.2–4.7. Using the results from our modelling analysis for the 12 of the 20 sources with the highest signal-to-noise emission lines that show disc-like kinematics, we conclude the following: (i) our sample prefers a CO-to-|$H_2$| conversion factor, of |$\alpha _{\rm CO} = 0.74 \pm 0.37$|⁠; (ii) these far-infrared luminous galaxies follow a similar Tully–Fisher relation between the circular velocity, |$V_{\rm circ}$|⁠, and baryonic mass, |$M_{\rm b}$|⁠, as less strongly star-forming samples at high redshift, but extend this relation to much higher masses – showing that these are some of the most massive disc-like galaxies in the Universe; (iii) finally, we demonstrate support for an evolutionary link between massive high-redshift dusty star-forming galaxies and the formation of local early-type galaxies using the both the distributions of the baryonic and kinematic masses of these two populations on the |$M_{\rm b}$| – |$\sigma$| plane and their relative space densities.

galaxies: evolution, galaxies: high-redshift, galaxies: kinematics and dynamics, submillimetre: galaxies

1 INTRODUCTION

The most massive galaxies in the Universe form from the highest density peaks in the primordial matter distribution (White & Rees 1978). Galaxy interactions and mergers are expected to be frequent in such environments and contribute to the growth of the most massive galaxies, as well as potentially imprinting variations in the kinematic structures of galaxies as a function of their mass (and environment): with the most massive galaxies exhibiting pressure-supported spheroidal morphologies (e.g. Ogle et al. 2019). At the present day the majority of these massive systems lack significant on-going star formation; they correspond to the red-and-dead elliptical galaxies that dominate the cores of massive clusters of galaxies (Faber 1973; Dressler 1980).

A number of mechanisms have been suggested to explain why the supply of gas from the intergalactic medium, needed to fuel star formation, is interrupted for many of the most massive galaxies. The accretion of cold gas may cease when the galaxy’s halo becomes massive enough that an accretion shocks develops, interrupting the inflow of cold gas streams necessary to replenish star-forming discs (Dekel & Birnboim 2006). In addition, active galactic nucleus (AGN) feedback from a growing supermassive black hole may heat or expel cool gas, further suppressing star formation (Bower et al. 2006; Hopkins et al. 2006).

Tracing the physical processes driving galaxy evolution is important for understanding the diversity of the galaxy populations we see at the present day. The most massive galaxies may be particularly useful in this regard, as they represent the high-mass limit for individual stellar systems. Connecting massive galaxy populations at different epochs may thus be more straightforward than attempting to track the formation and growth (through mergers and accretion on to larger galaxies) of less-massive systems.

Despite the complex nature of the physical processes that regulate star formation and lead to galaxy morphological transformations, a number of simple scaling relations can be used to investigate the effects of these mechanisms. For example disc galaxies exhibit an empirical correlation between their rotational velocity and their baryonic mass, otherwise known as the Tully–Fisher relation (TFR; Tully & Fisher 1977), while elliptical galaxies exhibit a similar relationship between their baryonic mass and the central stellar velocity dispersion, also known as the Faber–Jackson relation (FJR; Faber & Jackson 1976). Studying these relations for massive galaxy populations across cosmic time can help us understand how these galaxies evolved and potentially link populations at different epochs which are observed at different stages in their evolution.

Among the high-redshift galaxy populations, dusty star-forming galaxies (DSFGs, also referred to as sub-millimetre galaxies, SMGs) originally selected as sources with flux densities |$S \gtrsim 1$| mJy at 850 |$\mu$|m (Smail, Ivison & Blain 1997; Barger et al. 1998; Hughes et al. 1998; Eales et al. 1999), are among the most massive active star-forming systems that have been observed. The redshift distribution of sub-millimetre-selected DSFGs peaks around |$z \sim$| 2–3 (Chapman et al. 2005; Stach et al. 2019; Dudzevičiūtė et al. 2020), where the star formation activity and black hole accretion peak (Madau & Dickinson 2014) and the gas accretion on to galaxies reaches its maximum (Walter et al. 2020). Various studies have now established that these DSFGs have high-stellar masses (⁠|$M_{\star } \sim$| 10|$^{10}$|–10|$^{11}$| M|$_{\odot }$|⁠; Simpson et al. 2014; da Cunha et al. 2015; Smolčić et al. 2015; Miettinen et al. 2017; Dudzevičiūtė et al. 2020), are gas (⁠|$M_{\rm gas} \sim$| 10|$^{10}$|–10|$^{11}$| M|$_{\odot }$|⁠; Bothwell et al. 2013; Birkin et al. 2021) and dust rich (⁠|$M_{\rm dust} \sim 10^{8}$|–|$10^{10} \, \rm M_{\odot }$|⁠; Dudzevičiūtė et al. 2020) and form stars at extreme rates (⁠|$\rm SFR \sim$| 10|$^{2}$|–10|$^{3}$| M|$_{\odot } \rm yr^{-1}$| Swinbank et al. 2014; Miettinen et al. 2017), contributing up to |$\sim$|20 per cent to the total star formation rate density (Swinbank et al. 2014).

Having such high star formation rates, DSFGs can substantially increase their stellar mass on a very short time-scale (⁠|$\sim 100$| Myr; Bothwell et al. 2013; Birkin et al. 2021). Considering that they already have significant mass in stars, that will result in the descendants being very massive systems. It follows that this high-redshift population will evolve to form the most massive galaxies, which are mostly early-type galaxies in our Universe today. Several studies in the literature have proposed an evolutionary link between these two galaxy populations using a variety of arguments: (i) Clustering (e.g. Hickox et al. 2012; Chen et al. 2016; Wilkinson et al. 2017; Amvrosiadis et al. 2018; García-Vergara et al. 2020; Stach et al. 2021), (ii) space densities (e.g. Simpson et al. 2014), and (iii) stellar ages (e.g. Simpson et al. 2014; Dudzevičiūtė et al. 2020; Carnall et al. 2021) among others.

One aspect of the properties of DSFGs that currently we have only limited insight into is the dynamical state of these dusty high-redshift sources. Most of the studies in the literature into the dynamics of high-redshift star-forming galaxies have focused on more ‘typical’¹ less strongly star-forming sources (e.g. Förster Schreiber et al. 2009; Wisnioski et al. 2015, 2019; Price et al. 2016; Tiley et al. 2019). These surveys trace the kinematics of star-forming galaxies from observations of the ionized gas, usually the redshifted H|$\alpha$| emission line. The limitation of using the H|$\alpha$| line as a tracer to study the kinematics is that its susceptible to dust attenuation and the influence of ionized outflows. In addition, these studies typically extend out to |$z \sim$| 2.5, beyond which point, prior to the launch of JWST, we can no longer observe the H|$\alpha$| line from the ground and therefore miss the high-redshift tail of the star-forming galaxies redshift distribution (e.g. Birkin et al. 2023). These factors, combined with the lower surface density of DSFGs compared to typical star-forming galaxies, mean that current studies have not included significant numbers of dust rich systems, which also happen to be the most massive at high redshift (Dudzevičiūtė et al. 2020).

To investigate the kinematics of a representative and complete sample of the more dusty and active DSFG population, we need tracers of the interstellar medium (ISM) that are less influenced by dust obscuration than rest-frame optical emission lines. One such tracer is the emission from carbon monoxide (CO), specifically its low- to mid-J transitions, considered a reliable tracer of the bulk of the molecular gas in these systems. With facilities such as the Atacama Large Millimeter/submillimeter Array (ALMA) or the Northern Extended Millimeter Array (NOEMA), studies of CO in high-redshift dusty galaxies have become feasible in recent years (e.g. Greve et al. 2005; Bothwell et al. 2013; Tacconi et al. 2018; Birkin et al. 2021). In addition to various CO transitions, other far-infrared (IR) emission lines, such as [C ii], [C i], and H|$_2$|O, can be used for dynamical studies. The [C ii] emission line, in particular, has gained attention for being the brightest emission line in the far-IR spectrum. Surveys targetting this line have been undertaken in samples of normal star-forming galaxies up to redshift |$z\sim$| 4–6 (e.g. Le Fèvre et al. 2020), enabling the study of their dynamical properties (e.g. Jones et al. 2021).

While CO and far-IR emission lines are increasingly being used to study the dynamics of DSFGs, many current investigations concentrate on studying a few individual sources, including both lensed and unlensed sources in different environments (field or group/cluster members) and targetting a mixture of emission lines (e.g. Swinbank et al. 2011; Hodge et al. 2012; Dye et al. 2015; Olivares et al. 2016; Chen et al. 2017; Rizzo et al. 2020, 2021, 2023; Fraternali et al. 2021; Hogan et al. 2021, 2021; Lelli et al. 2021, 2023; Tsukui & Iguchi 2021; Xiao et al. 2022; Roman-Oliveira, Fraternali & Rizzo 2023). Due to the complexity of the selection criteria in these studies, it is challenging to generalize their findings to the wider population. However, a trend that is beginning to emerge is that the population of active galaxies, especially those at |$z\gt 4$|⁠, contain a significant fraction of dynamically cold gas discs. To confirm this trend, larger samples with well-defined selection criteria are necessary.

In our work, we study the kinematics of a large sample of unlensed DSFGs, uniformly selected based on their flux at 870 |$\mu$|m, which was originally presented in Birkin et al. (2021). We focus our analysis on sources with CO detections, which is historically considered the best tracer of the molecular gas in these systems. We will use this sample to model the dynamics of these systems and then use these results to study scaling relations between the dynamical properties of the population.

The outline of this paper is as follows: in Section 2, we introduce the sample we will use in this work and discuss its properties. In Section 3, we describe the method we use to model the kinematics for our sources, which is specifically tailored for interferometric observations. In Section 4, we discuss our main scientific results and finally, in Section 5, we give a summary of our findings. Throughout this work, we adopt a spatially flat Lambda cold dark matter (⁠|$\Lambda$|CDM) cosmology with H|$_0=67.8 \pm 0.9$| km s|$^{-1}$| Mpc|$^{-1}$| and |$\Omega _{\rm M}=0.308 \pm 0.012$| (Planck Collaboration XIII 2016).

2 DATA

In this section, we introduce the sample of sources that we use in this work. As already mentioned in the introduction, these sources primarily come from the recent study of Birkin et al. (2021), which presents a large CO survey of massive DSFGs. We focus on the analysis of ALMA observations of sources located in the Extended Chandra Deep Field South (ECDFS) from that study, which are part of the ALESS survey (Hodge et al. 2013), as the observations of sources in other fields are generally at much lower spatial resolution.We complement this sample with three additional sources that were not in the Birkin et al. (2021) sample, but have data available of similar quality in the ALMA archive (ALESS 049.1, 075.1, 122.1).

We begin by discussing our sample selection including defining a signal-to-noise (SNR) cut for our analysis. At the end of this section, we discuss some of the physical properties (i.e. stellar and gas masses) of our sample. We also compare the properties of our sample to other, more ‘typical’, star-forming galaxies at high redshift. We focus the comparison on samples for which large follow-up surveys, targeting the H|$\alpha$| emission line, have been carried out with the aim of modelling their kinematics.

2.1 Sample

The sources used in this work were first discovered in the LABOCA ECDFS Submillimeter Survey (LESS; Weiß et al. 2009) which is a large homogeneous 870 |$\mu$|m survey of the ECDFS conducted with the Large Apex BOlometer Camera (LABOCA) on the APEX telescope. The LESS survey resulted in a sample of 126 sources with flux densities |$S_{870 \rm \mu m} \gt 3.6$| mJy. These sources were subsequently followed-up with ALMA during Cycle 0 in the ALESS survey (Hodge et al. 2013). The sample from the ALESS survey comprises 99 sources from the maps of 69 LESS sources with high quality ALMA observations, with a considerable fraction (⁠|$\gtrsim 35{{\ \rm per\ cent}}$|⁠) of the LABOCA sources are being resolved into multiple sources in the ALMA maps (Karim et al. 2013).

A series of campaigns have been conducted in the optical, IR, and millimeter (mm) wavelengths to obtain spectroscopic redshifts for these ALESS sources (e.g. Danielson et al. 2017; Wardlow et al. 2018; Birkin et al. 2021). These redshift identifications come by targetting various CO transitions (⁠|$J_{\rm up} =$| 2–5) and the [C i] (⁠|$\rm ^{3}P_1 - ^{3}P_0$|⁠) line in the mm wavebands (see Wardlow et al. 2018; Birkin et al. 2021) or rest-frame ultraviolet (UV) and optical emission lines (Danielson et al. 2017). From the parent sample of 99 ALESS sources that have been detected in continuum, 30 now have robust CO line detections (Calistro Rivera et al. 2018; Wardlow et al. 2018; Birkin et al. 2021).

The sources that we analyse in this work will be selected from this sample of 30 sources (which were observed as part of the following ALMA proposal IDs: 2016.1.00564.S, 2016.1.00754.S, 2017.1.01163.S, 2017.1.01512.S). As noted earlier we also include ALESS 049.1, 075.1, 122.1, previously discussed in Wardlow et al. (2018) and Calistro Rivera et al. (2018). In Table 1, we list these 30 sources for which we report their redshifts, observed CO transition, the size of the beam² (major |$\times$| minor axis) as well as some other properties which we discuss later in the section.

Table 1.

Properties of the parent sample: (1) ID, (2) spectroscopic redshift, (3) observed CO transition, (4) major |$\times$| minor axis of the synthesized beam, (5) stellar mass, (6) molecular gas mass (derived from the CO luminosity), (7) star formation rate, (8) infrared luminosity, (9) velocity resolution of the cube used in the modelling and producing velocity maps, (10) integrated signal-to-noise ratio of the CO emission line (using the velocity-integrated intensity and its associated error), and (11) classification based on the morpho-kinematical features of the observed 3D cubes (see Section 2.2). The properties of the sources listed here are taken from Birkin et al. (2021), unless the source was not included in that work in which case the properties are taken from various studies in the literature indicated by the footnotes.

ALESS	z	CO	|$I_{\rm CO}$|	Beam	|$\log \big (M_{\star }\big)$|	|$\log \big (M_{\rm gas}\big)$|	|$\rm \log \big (SFR \big)$|	|$\log \big (L_{\rm IR}\big)$|	|$\Delta v$|	SNR	Class
ID		transition	|$\rm Jy \, km \, s^{-1}$|	arcsec	|$\big ( \rm M_{\odot } \big )$|	|$\big ( \rm M_{\odot } \big )$|	|$\big ( \rm M_{\odot } \, yr^{-1} \big )$|	|$\big ( \rm L_{\odot } \big )$|	(km s−1)
003.1	3.375	4–3	|$1.1 \pm 0.1$|	1.77 |$\times$| 1.32	|$11.28^{+0.15}_{-0.25}$|	|$10.96 \pm 0.06$|	|$2.85^{+0.08}_{-0.07}$|	|$12.93^{+0.07}_{-0.06}$|	45	8.7	II
006.1a	2.337	3–2	|$1.7 \pm 0.2$|	0.91 |$\times$| 0.69	|$10.98^{+0.56}_{-0.50}$|	|$11.00 \pm 0.06$|	|$2.32^{+0.07}_{-0.21}$|	|$12.43^{+0.05}_{-0.03}$|	45	17.6	II
007.1	2.692	3–2	|$2.2 \pm 0.3$|	1.79 |$\times$| 1.27	|$11.87^{+0.01}_{-0.02}$|	|$11.22 \pm 0.06$|	|$2.76^{+0.06}_{-0.01}$|	|$12.96^{+0.01}_{-0.01}$|	50	9.9	I
009.1	3.694	4–3	|$1.0 \pm 0.1$|	1.80 |$\times$| 1.63	|$11.86^{+0.18}_{-0.09}$|	|$10.99 \pm 0.05$|	|$2.74^{+0.17}_{-0.15}$|	|$13.01^{+0.09}_{-0.08}$|	96	9.6	II
017.1	1.538	2–1	|$0.6 \pm 0.1$|	1.02 |$\times$| 0.80	|$11.25^{+0.00}_{-0.00}$|	|$10.49 \pm 0.09$|	|$2.14^{+0.00}_{-0.00}$|	|$12.25^{+0.00}_{-0.00}$|	51	9.1	I
022.1	2.263	3–2	|$1.4 \pm 0.1$|	1.80 |$\times$| 1.63	|$11.67^{+0.09}_{-0.07}$|	|$10.89 \pm 0.04$|	|$2.48^{+0.13}_{-0.16}$|	|$12.66^{+0.07}_{-0.06}$|	45	13.5	I
034.1	3.071	3–2	|$0.4 \pm 0.1$|	0.99 |$\times$| 0.81	|$10.66^{+0.07}_{-0.10}$|	|$10.62 \pm 0.08$|	|$2.52^{+0.14}_{-0.18}$|	|$12.66^{+0.16}_{-0.18}$|	55	8.4	II
041.1	2.547	3–2	|$0.85 \pm 0.03$|	1.24 |$\times$| 0.95	|$11.83^{+0.16}_{-0.17}$|	|$10.57 \pm$| 0.08	|$2.44^{+0.15}_{-0.26}$|	|$12.64^{+0.08}_{-0.08}$|	48	29.1	I
049.1b,c	2.945	3–2	|$0.88 \pm 0.03$|	1.24 |$\times$| 0.95	|$10.58^{+0.12}_{-0.22}$|	|$10.89 \pm 0.02$|	|$2.83^{+0.09}_{-0.05}$|	|$12.83^{+0.04}_{-0.07}$|	54	27.9	I
062.2	1.362	2–1	|$1.4 \pm 0.1$|	0.87 |$\times$| 0.69	|$10.68^{+0.02}_{-0.00}$|	|$10.75 \pm 0.07$|	|$2.68^{+0.01}_{-0.01}$|	|$12.52^{+0.01}_{-0.00}$|	48	23.8	I
065.1	4.445	4–3	|$0.9 \pm 0.1$|	0.99 |$\times$| 0.81	|$10.48^{+0.19}_{-0.13}$|	|$11.05 \pm 0.06$|	|$2.48^{+0.10}_{-0.15}$|	|$12.64^{+0.14}_{-0.18}$|	55	10.4	II
066.1a	2.553	3–2	|$1.0 \pm 0.1$|	0.86 |$\times$| 0.69	|$11.42^{+0.33}_{-0.44}$|	|$10.83 \pm 0.04$|	|$2.66^{+0.10}_{-0.13}$|	|$12.78^{+0.05}_{-0.08}$|	48	13.3	I
067.1	2.122	3–2	|$1.2 \pm 0.1$|	2.18 |$\times$| 1.59	|$11.25^{+0.00}_{-0.01}$|	|$10.78 \pm 0.04$|	|$2.19^{+0.01}_{-0.01}$|	|$12.55^{+0.01}_{-0.01}$|	42	9.5	II
071.1	3.709	4–3	|$1.25 \pm 0.07$|	1.15 |$\times$| 0.96	|$12.31^{+0.00}_{-0.03}$|	|$11.20 \pm 0.04$|	|$3.42^{+0.06}_{-0.00}$|	|$13.48^{+0.06}_{-0.00}$|	48	17.6	I
075.1b,c	2.552	3–2	|$0.99 \pm 0.02$|	1.24 |$\times$| 0.95	|$10.48^{+0.00}_{-0.21}$|	|$10.84 \pm 0.02$|	|$2.65^{+0.13}_{-0.14}$|	|$12.58^{+0.17}_{-0.11}$|	48	38.7	I
088.1	1.206	2–1	|$1.4 \pm 0.1$|	0.90 |$\times$| 0.69	|$9.892^{+0.00}_{-0.00}$|	|$10.65 \pm 0.07$|	|$1.83^{+0.00}_{-0.00}$|	|$12.22^{+0.00}_{-0.00}$|	45	15.4	II
098.1	1.374	2–1	|$3.3 \pm 0.2$|	0.86 |$\times$| 0.69	|$11.48^{+0.06}_{-0.05}$|	|$11.13 \pm 0.07$|	|$2.43^{+0.04}_{-0.04}$|	|$12.55^{+0.03}_{-0.01}$|	48	29.7	I
101.1	2.353	3–2	|$2.0 \pm 0.2$|	0.90 |$\times$| 0.69	|$11.20^{+0.18}_{-0.20}$|	|$10.93 \pm 0.04$|	|$2.16^{+0.13}_{-0.18}$|	|$12.30^{+0.09}_{-0.09}$|	46	13.9	II
112.1	2.316	3–2	|$1.8 \pm 0.1$|	0.91 |$\times$| 0.69	|$11.21^{+0.11}_{-0.14}$|	|$11.01 \pm 0.04$|	|$2.44^{+0.07}_{-0.10}$|	|$12.50^{+0.05}_{-0.05}$|	44	16.3	II
122.1b,d	2.024	3–2	|$4.2 \pm 0.8$|	0.45 |$\times$| 0.35	|$10.89^{+0.21}_{-0.21}$|	|$11.30 \pm 0.09$|	|$2.84^{+0.17}_{-0.16}$|	|$12.92^{+0.13}_{-0.25}$|	123	13.9	I
001.1	4.674	5–4	|$1.0 \pm 0.2$|	1.78 |$\times$| 1.33	|$10.93^{+0.12}_{-0.18}$|	|$11.09 \pm 0.09$|	|$2.82^{+0.13}_{-0.16}$|	|$12.94^{+0.17}_{-0.21}$|	90	4.8	III
001.2	4.669	5–4	|$0.8 \pm 0.1$|	1.78 |$\times$| 1.33	|$11.06^{+0.11}_{-0.13}$|	|$11.04 \pm 0.06$|	|$2.56^{+0.17}_{-0.16}$|	|$12.66^{+0.16}_{-0.18}$|	90	7.3	III
005.1	3.303	4–3	|$0.7 \pm 0.1$|	1.80 |$\times$| 1.63	|$11.37^{+0.35}_{-0.01}$|	|$10.78 \pm 0.07$|	|$2.97^{+0.01}_{-0.19}$|	|$12.95^{+0.01}_{-0.13}$|	87	6.6	III
019.2	3.751	4–3	|$0.9 \pm 0.1$|	1.80 |$\times$| 1.63	|$11.43^{+0.17}_{-0.19}$|	|$10.96 \pm 0.07$|	|$2.88^{+0.17}_{-0.07}$|	|$13.02^{+0.17}_{-0.06}$|	96	7.8	III
023.1	3.332	4–3	|$0.7 \pm 0.1$|	1.86 |$\times$| 1.56	|$11.45^{+0.21}_{-0.25}$|	|$10.74 \pm 0.08$|	|$2.75^{+0.11}_{-0.11}$|	|$12.84^{+0.10}_{-0.07}$|	88	6.0	III
031.1	3.712	4–3	|$0.9 \pm 0.1$|	1.80 |$\times$| 1.63	|$11.52^{+0.10}_{-0.14}$|	|$10.93 \pm 0.07$|	|$2.81^{+0.11}_{-0.11}$|	|$12.94^{+0.09}_{-0.09}$|	96	4.8	III
035.1a	2.974	3–2	|$0.9 \pm 0.1$|	2.14 |$\times$| 1.60	|$11.63^{+0.20}_{-0.24}$|	|$10.89 \pm 0.06$|	|$2.64^{+0.12}_{-0.15}$|	|$12.79^{+0.05}_{-0.04}$|	108	7.2	III
061.1	4.405	4–3	|$1.1 \pm 0.1$|	0.98 |$\times$| 0.81	|$10.33^{+0.18}_{-0.01}$|	|$11.14 \pm 0.07$|	|$2.42^{+0.16}_{-0.01}$|	|$12.54^{+0.20}_{-0.01}$|	55	7.9	III
068.1	3.507	4–3	|$0.3 \pm 0.1$|	1.85 |$\times$| 1.56	|$10.97^{+0.26}_{-0.26}$|	|$10.42 \pm$| 0.08	|$2.54^{+0.11}_{-0.11}$|	|$12.65^{+0.09}_{-0.11}$|	92	5.2	III
079.1	3.901	4–3	|$0.4 \pm 0.1$|	1.86 |$\times$| 1.56	|$11.41^{+0.11}_{-0.14}$|	|$10.48 \pm 0.09$|	|$2.32^{+0.12}_{-0.12}$|	|$12.47^{+0.08}_{-0.06}$|	100	3.9	III

a

SED fits are poorly constrained or the SED is potentially contaminated by an AGN.

b

Stellar mass, far-infrared luminosity and SFR from da Cunha et al. (2015).

c

gas mass from Wardlow et al. (2018).

d

gas mass from Calistro Rivera et al. (2018).

2.2 Signal-to-noise ratio selection and classification

We already mentioned that the aim of this work is to model the dynamics of the DSFGs in our sample. In order to perform such an analysis, the data we work with need be of sufficient signal-to-noise ratio (SNR). Here, we refer to the integrated SNR which is defined as the ratio of the velocity integrated flux to its associated error.³ We measured the integrated SNR for each source in our sample, which we report in Table 1.

Before we go into the various details of our modelling analysis we want to make some useful clarifications. We attempted to model all 30 sources in our sample, however, we found that the parameters of our model were effectively unconstrained when modelling sources with an integrated SNR |$\lt 8$|⁠. We therefore applied a further selection cut based on SNR, considering only sources with SNR |$\gt 8$|⁠, which results in a sample of 20 sources.

We can now go a step further and divide the sources that satisfy our SNR selection criteria in two classes. We follow a classification approach that previous studies in the literature have used when dealing with 3D integral field spectroscopic data (e.g. Förster Schreiber et al. 2009; Le Fèvre et al. 2020). This classification is based on inspecting both the individual channel maps as well as the velocity maps. The velocity maps of our sources with SNR |$\gt 8$| are shown in Fig. 1. These were produced by fitting a single Gaussian function to the spectrum in each pixel of our 3D cleaned cubes and using the inferred mean of the fitted Gaussian as the value of that pixel in our 2D velocity maps.⁴ The characteristic difference between these two classes is whether the observed emission varies smoothly across the different channels of the cube, resulting in a well defined gradient in the velocity maps (Class I; this is a necessary condition for a source to be considered a rotating disc) or shows a complex behaviour (Class II; which can be considered as an indication that the source is undergoing a minor/major merging event). Two out of these 20 sources, ALESS 003.1/009.1, are placed in Class II due to the lack of sufficient resolution. Finally, all sources that do not satisfy our SNR selection criteria (SNR |$\lt 8$|⁠), for which we lack sufficient SNR to unambiguously determine their kinematic nature, are put in a third class (Class III). These three classes are summarized as follows: Class I – smooth transition between channel maps resulting in a well defined gradient for the velocity map; Class II – complex velocity map or galaxy-pairs or low resolution; Class III – low SNR. The specific class assigned for each source is given in the last column of Table 1.

Figure 1.

Velocity fields for the 20 DSFGs in our sample with CO emission lines that have integrated SNR > 8. The velocity fields are derived from the observed CO emission lines using a pixel-by-pixel spectral fitting method. The kinematic class assigned to each source, from our visual classification, is shown next to the source ID. Sources that fall in Class I display smooth velocity gradients across the source while Class II sources have complex features in the velocity maps. All sources are shown on the same angular scale, where the black bar at the bottom right corner of the figure corresponds to 1 arcsec or |$\sim$| 8.5 kpc at |$z =$| 2. The grey ellipse at the bottom left corner of each velocity map shows the corresponding synthesized beam for that map.

2.3 Properties

The properties of our sample including redshifts, stellar masses (⁠|$M_{\star }$|⁠), gas masses (⁠|$M_{\rm gas}$|⁠), far-IR luminosities (⁠|$L_{\rm IR}$|⁠) and star formation rates (SFRs), were taken from Birkin et al. (2021) or if not available, then from da Cunha et al. (2015), Wardlow et al. (2018), or Calistro Rivera et al. (2018). We briefly summarize here how these quantities were computed, but we refer the reader to Birkin et al. (2021) or the other studies for more details.

Stellar masses, far-IR luminosities and star formation rates were computed using the spectral energy distribution (SED) fitting code magphys (da Cunha et al. 2015; Battisti et al. 2019), keeping the redshifts fixed to the values derived from the observed CO lines. Gas masses were computed from the velocity-integrated intensity, |$I_{\rm CO}$|⁠, of the observed transition lines. This was done as follows: the CO line luminosities were computed using the recipe from Solomon & Vanden Bout (2005). The measured line luminosities of each transition were converted to CO(1–0) luminosities using the excitation corrections that were measured for SMM J2135–0102 (Danielson et al. 2011). Finally, gas masses were calculated from the CO(1–0) luminosities assuming a CO–H|$_2$| conversion factor, |$\alpha _{\rm CO} = 0.74 \pm 0.37$|⁠, which, as we show later in Section 4.2, is the value preferred for our sample and is consistent with previous estimates of this quantity for starburst galaxies (Bolatto, Wolfire & Leroy 2013; Calistro Rivera et al. 2018; Birkin et al. 2021).

We caution that the results from the SED fitting analysis with magphys could potentially be inaccurate for some sources in our sample. Specifically, the presence of a luminous AGN (Wang et al. 2013), which is not accounted for in the SED model, can contaminate the UV to mid-IR part of the SED and lead to a biased estimate of the stellar mass. However, the one source which is confirmed to host an AGN, ALESS 17.1, shows no indication of AGN contamination of their mid-IR SEDs. In addition, we note that there are poor magphys SED fits for two sources: ALESS 66.1 and 75.1, that make their derived properties uncertain [ALESS 66.1 has a nearby foreground quasar (Chen et al. 2020) which could be the cause of the poor SED fit]. We have confirmed that the removal of these sources does not qualitatively change any of our conclusions and so we have retained them in our analysis.

In order to place our sample into perspective, in Fig. 2 we plot our sources on the SFR versus |$M_{\star }$| plane along with other samples of high-redshift star-forming galaxies. We also show the parent sample of ALESS sources as presented in Birkin et al. (2021) and SMGs in a similar redshift range to our sample (⁠|$z\sim$| 1–5) selected from the ALMA survey of the SCUBA-2 CLS UDS field (AS2UDS; Dudzevičiūtė et al. 2020). In addition, we also show SFGs in the redshift range |$z=$| 1.8–2.7 from the KMOS|$^{\rm 3D}$| survey (Wisnioski et al. 2019). This last sample represents more typical population of star-forming galaxies at high redshift, although with a sample selection biased against the most obscured and potentially massive galaxies. Compared to our sample these are less massive and form stars at much lower rates. We chose to use the KMOS|$^{\rm 3D}$| sample for comparison as this is the largest sample with resolved kinematics from observations of the H|$\alpha$| emission line. We will use this sample in later sections for comparing to the kinematical properties of our sample.

Figure 2.

Star formation rate as a function of stellar mass for galaxies in our sample, illustrating that the DSFGs in our study are more massive and more strongly star-forming than ‘typical’ star-forming galaxies at these redshifts. Blue points correspond to sources in our sample, of which those with SNR |$\gt 8$| are shown as open symbols. Black points correspond to typically less active star-forming galaxies from the KMOS|$^{\rm 3D}$| (Wisnioski et al. 2019) in the redshift range |$z=$| 1.8–2.7. The distribution for the SMGs from the AS2UDS survey (Dudzevičiūtė et al. 2020) are shown as the red 2D histogram. The broken power-law model drawn as the black solid line, was optimized using 3D-HST data in all Cosmic Assembly Near-IR Deep Extragalactic Legacy Survey (CANDELS) fields selected in the redshift range |$z=$| 2.0–2.5, where UV+IR SFRs were used (Whitaker et al. 2014). The dashed and dot–dashed curves are |$\times$|4 and |$\times$|10 away (below/above) the mean relation.

As we can see from Fig. 2 our sources have a similar distribution on the |$M_{\star }$|–SFR plane to the much larger AS2UDS sample. In order to determine if our ALESS sample is representative of the DSFGs population, we performed a series of two-sample Kolmogorov–Smirnov (KS) tests for some of the properties listed in Table 1. For these tests we matched the AS2UDS sample based on the flux at 870 |$\mu \mathrm{ m}$| and found probabilities consistent with them being drawn from the same distribution, indicating that our sample is representative of the DSFG population with |$S_{870 \rm \mu \mathrm{ m}} \gtrsim 1$| mJy.

3 METHOD AND ANALYSIS

In this section, we describe the method we use to model the kinematics of galaxies using data that come from interferometric observations. The data we use to constrain the parameters of the models, are called visibilities. Each data point, |$\mathbf {d}_{\rm ij}$|⁠, is a complex number (the |$i,j$| subscript denote the visibility measured by the pair of antennas i and j). Each visibility corresponds to a point in |$uv$|-space (i.e. a sample in Fourier space) which is defined by its |$(u, v)$| coordinates (in units of wavelengths, |$\lambda$|⁠). In the rest of this section, we will describe the dynamical 3D models we use in this work to constrain the kinematics, the method to convert a 3D model cube (x, y, |$\lambda$|⁠) to a set of visibility points (u, v, |$\lambda$|⁠) and the framework we use to perform the parameter inference.

3.1 Data preparation

The raw data were calibrated using the Common Astronomy Software Applications (casa; McMullin et al. 2007), employing the standard pipelines provided with each data set. The un-flagged calibrated visibilities were then exported to a fits format for further analysis with our own python routines.

The modelling analysis, which we will describe later in this section, is performed in the |$uv$|-plane and therefore requires accurate weights associated with each visibility data point. The weights of the visibilities that come us a product of casa’s calibration pipeline are computed in a relative manner and so do not reflect the observed scatter of the data. In order to fit models to the data (i.e. visibilities) we require weights to be computed in an absolute manner. We, therefore, recompute the weights of the calibrated visibilities using casa task STATWT, which re-calculates the weights of the visibilities according to their scatter as a function of time and baseline. In practice what this task does is to select visibilities in a user-defined time interval for each pair of antennas and compute the scatter in the real and imaginary components of the visibilities. The scatter in the two components of the visibilities now corresponds to the error associated with those visibilities that were used to compute it.

3.2 Kinematical model

In order to extract the dynamical properties of our sources we need to fit appropriate models to the data and constrain their parameters. The velocity maps (Fig. 1) of most of our sources indicate the presence of order circular motions, therefore the appropriate choice of a model is that of rotating disc (deviations from ordered circular motions in the data can perhaps result in unphysical model parameters or ‘poor’ fits, which we will discuss later).

We generate 3D parametric kinematical models of a thick rotating-disc using the publicly available software galpak3d⁵ (Bouché et al. 2015). In total, our model is comprised of 10 free parameters: the x, y, and z positions of the source centre in the cube, a scaling factor for the flux, effective radius, inclination, position angle (measured clockwise from East to North), turnover radius, maximum circular velocity (not projected, but asymptotic irrespective of inclination), and velocity dispersion (isotropic and constant over the disc).

In order to generate a model we need to specify the flux, disc-thickness and velocity profiles. For the flux-profile of the disc we use an exponential profile, which is a special case of the more generic Sersic profile (Sérsic 1963) which is given by,

$$\begin{eqnarray} I(r)=I_e \exp \left\lbrace -b_{n} \left[ \left(\frac{r}{r_e}\right)^{1/n} - 1 \right] \right\rbrace \, , \end{eqnarray}$$

(1)

where |$r_e$| is the effective radius of the disc, |$I_e$| is the surface brightness of the disc at |$r_e$|⁠, and |$b_n$| is a constant which given by |$b_n \approx 1.9992n - 0.3271$|⁠, where n is the Sérsic index (the definition of |$b_n$| is such that |$r_e$| is equivalent to the half-light radius, |$R_{1/2}$|⁠). For an exponential profile the Sérsic index is set to |$n=1$|⁠. For the disc thickness profile we use a Gaussian profile,

$$\begin{eqnarray} I(z) \propto \exp \left\lbrace -\frac{z^2}{2 h_z^2} \right\rbrace \, , \end{eqnarray}$$

(2)

which is defined perpendicular to the disc and is added to the disc component. The disc thickness is set to |$h_z = 0.15 R_{1/2}$|⁠, where |$R_{1/2}$| is the half-light radius of the disc.⁶ Finally, for the velocity profile we use an isothermal model, which is given by

$$\begin{eqnarray} V(r) = V_{\rm max} \left[1 - \arctan \left(\frac{r}{r_t} \right) \Big / \left(\frac{r}{r_t} \right) \right] \, , \end{eqnarray}$$

(3)

where |$r_t$| is the turnover radius and |$V_{\rm max}$| is the maximum circular velocity. In what follows, we define the model cube in real-space, given a set of parameters, |$\theta$|⁠, as |$\mathbf {s}(\theta)$|⁠.

3.3 The non-uniform fast Fourier transform

In order to compare our model (i.e. the surface brightness in each channel of the cube) with our data we need to convert our model cubes defined in the real-space to visibilities defined in the Fourier space. An image of the surface brightness of a source, |$I(x, y)$|⁠, in the |$xy$|-space and the visibilities, |$V(u, v)$|⁠, in the |$uv$|-space form a 2D Fourier transform pair,

$$\begin{eqnarray} V(u, v) = \int \int I(x, y) \exp ^{ -2 \pi i (ux + vy) } \mathrm{ d}x \mathrm{ d}y .\, \end{eqnarray}$$

(4)

However, since the samples in the |$uv$|-space are not uniformly distributed this operation reduces to a direct discrete Fourier transform (DFT). The DFT is a memory heavy and time-consuming operation making the analysis of large interferometric data very prohibitive. This problem can be alleviated by the use of a non-uniform fast Fourier transform (NUFFT). The use of NUFFTs for analysing astronomical interferometric data has been discussed in recent works (e.g. Rizzo et al. 2020; Powell et al. 2021) and a plethora of literature exist on the theoretical side (e.g. Beatty, Nishimura & Pauly 2005). In our work, we use the publicly available software pynufft⁷ (Lin 2018) to perform this computation. In order to set-up our NUFFT operator we need to pass it a real-space grid and the |$(u, v)$| coordinates of our Fourier samples.

We determine the number of image-plane pixels and the size of each pixel based on the |$uv$|-coverage and the Nyquist criterion. We require that the size of each pixel is at least half the resolution of our observation which we approximate as the inverse of the maximum |$uv$| distance (in arcsec). The number of image-plane pixels and the size of each pixel are determined by choosing the size of the field-of-view (FoV) and rounding up the number of image-plane pixels to the closest power of 2 so that the pixel scale is at least half the resolution.

Once we have defined the image-plane grid, following the procedure that was described above, we can then initialize our NUFFT operator. The NUFFT operator, |$\mathbf {D}$|⁠, will have dimensions (⁠|$2 N_{\rm vis} \, \times \, N_{\rm p}$|⁠), where |$N_{\rm p}$| is the number of image-plane pixel and |$N_{\rm vis}$| are the number of visibilities that we use in our analysis. To initialize the NUFFT operator, besides the image-plane grid, we also use the |$uv$|-coordinates of our visibilities. Here, we need to note that the |$(u, v)$| position of a visibility recorded at time, |${\bf t}_0$|⁠, shifts radially outwards from the phase-centre as the frequency increases. We therefore define a separate NUFFT operator, |$\mathbf {D}_n$| (where |$n = 1, 2, 3, ..., N_c$|⁠), for each of the |$N_c$| channels of our cube (each of the |$\mathbf {D}_n$| operators on the respective channel image of the cube |$\mathbf {s}_{n}(\theta)$|⁠). Therefore, applying the NUFFT operator to our model cube, |$\mathbf {D} \mathbf {s}(\theta)$|⁠, has an effect equivalent to convolving the model cube with the dirty beam (which can be thought of as the point spread function) in the |$xy$|-plane and the line spread function (LSF) along the z-axis (i.e. the frequency axis).

3.4 Constraining the model parameters

In order to fit our model to the data we use the publicly available software pyautofit⁸ (Nightingale et al. 2021).pyautofit is probabilistic programming language (PPL) that is designed to provide the user with a flexible interface to fit a model to a set of data points by defining the likelihood function of this model given the data.

The optimization of the model parameters, |$\theta$|⁠, is performed in a Bayesian framework. The posterior probability distribution, |$P(\theta | {\bf d})$|⁠, is given by

$$\begin{eqnarray} P(\theta | {\bf d}) \propto P({\bf d} | \theta) P(\theta) \, , \end{eqnarray}$$

(5)

where |$P({\bf d} | \theta)$| is the likelihood and |$P(\theta)$| the prior distribution of our model parameters. It is well-known that the noise dominates over the signal for an individual visibility point and its nature is random (thermal noise), which results in the distributions of both the real and imaginary components of our visibilities being well described by a Gaussian distribution (Thompson, Moran & Swenson 2017). We can therefore assume that the likelihood function is also Gaussian and write it as,

$$\begin{eqnarray} P({\bf d} | \theta) &=& \frac{1}{\sqrt{\det (2 \pi C)}} \\ && \exp \Big \lbrace -\frac{1}{2} \mathop {\sum }_{n = 0}^{N_c} \left(\mathbf {D}_n \mathbf {s}_n(\theta) - \mathbf {d}_n \right)^T \mathbf {C}_n^{-1} \left(\mathbf {D}_n \mathbf {s}_n(\theta) - \mathbf {d}_n \right) \Big \rbrace \, ,\\ \end{eqnarray}$$

(6)

where |${\bf d}_n$| are the observed visibilities, |${\bf C}_n^{-1}$| are the errors of the observed visibilities (see Section 3.1), |${\bf s}_n(\theta)$| are the images for each channel of the model cube given a set of model parameters |$\theta$|⁠, and |${\bf D}_n$| is the NUFFT operator. The index n corresponds to a channel of the cube. The matrix C is a block diagonal matrix where each block correspond to the individual |${\bf C}_n$| matrices.

We use a Gaussian prior for the geometric centre, which is centred on the peak value of the intensity (⁠|$0^{\rm th}$| moment) map and has a standard deviation of 1 arcsec. For the rest of our model parameters we use uniform priors, 0 |$\rm km \, s^{-1}$|<|$V_{\rm max}$|< 600 |$\rm km \, s^{-1}$|⁠; 0 |$\rm km \, s^{-1}$|<|$\sigma$|< 200 |$\rm km \, s^{-1}$|⁠; 0 arcsec <|$r_e$|< 2 arcsec; 0 arcsec <|$r_t$|< 1 arcsec; 0°<|$\theta$|< 90°; and 0°<i< 90°.

We note that for the optimization we actually minimize the log of the likelihood. Having defined our log-likelihood function, pyautofit allows the user to choose between a collection of different non-linear samplers to be used to constrain the parameters of the model by sampling the posterior distribution of these parameters. In our analysis we use a nested sampling Monte Carlo technique of Skilling (2006) implemented in the pymultinest algorithm, which itself is a python wrapper of the multinest (Buchner et al. 2014) algorithm (Feroz & Hobson 2008; Feroz et al. 2009). Finally, in order to speed-up the likelihood evaluation we create manual masks in the velocity axis, ignoring channels that do not contain any emission. This essentially reduces the number of NUFFTs that we have to perform each time we evaluate a likelihood for a given combination of model parameters.

3.5 Best-fit models

As we already discussed in section 2.2, we attempted to model all 30 sources listed in Table 1, but here, we only discuss the results for the 20 sources with SNR |$\gt 8$|⁠. For the remaining 10 sources with SNR |$\lt 8$| the best-fit parameters were largely unconstrained due to the low SNR of our data, effectively not allowing us to characterize their kinematics (rotation/dispersion supported or merger/disrupted). These 10 sources, which are listed at the bottom of Table 1 with Class III, will not be considered as part of our sample for the rest of this work.

In Fig. 3, we show the results from our dynamical modelling analysis for one of the sources in our sample, ALESS 122.1. The different rows correspond to the dirty channel maps of the data, model (ML) and residuals (data-model). This source was previously studied in Calistro Rivera et al. (2018) using the same dynamical modelling software but instead the analysis was carried out in image-plane rather than the |$uv$|-plane. Our inferred values are in good agreement with those reported in Calistro Rivera et al. (2018), but their claimed uncertainties are significantly lower (e.g. |$V_{\rm max} = 564 \pm 8$| and |$\sigma = 129 \pm 1$| km s|$^{-1}$|⁠; Calistro Rivera, private communication). This demonstrates one of the advantages of carrying out the analysis in the |$uv$|-plane where the observational errors (i.e. error on the real and imaginary components of the visibilities) are better defined compared to the uncertainties defined in the image-plane (i.e. the rms noise of cleaned images). Complementary to this figure are the dirty 1st-moment (intensity-weighted velocity; left column) and 2nd-moment (intensity-weighted dispersion; right column) maps, which we show in Fig. 4 for three of the sources in our sample. We reiterate that while our analysis is carried out in the |$uv$|-plane, it is useful to visually inspect both the residual (data versus model) dirty cubes and moment maps to determine if a fit is successful.

Figure 3.

Dirty channel maps of the normalized data, model, and residuals of the CO(3–2) emission line for ALESS 122.1. The velocity of each channel, computed with respect to the redshifted CO emission line (⁠|$z = 2.024$|⁠), is indicated at the top left corner in each channel map of the data panels. The red and blue contours in the data and model panels, respectively, are drawn from |$3\sigma$| to |$6\sigma$| and increase in steps of |$1\sigma$|⁠, where the rms noise is measured considering all channels of the cube. The purple straight line denoted the best-fit major axis. This source was previously analysed in Calistro Rivera et al. (2018), where the modelling was carried out in the real-plane, and found similar values for the model parameters.

Figure 4.

Moment 1 (velocity; left column) and 2 (dispersion; right column) maps for three example sources in our sample (the name of each source is indicated at the top left corner). From left to right in each row we show the observed data, the maximum likelihood model and residual dirty moment maps (the moment maps were computed after masking any emission below |$3\sigma$|⁠). In these example sources (as well as for the other sources regarded as well described by our rotating disc model; see Section 3.5) the model does a good job at reproducing the observed data. The black dotted line in each panel corresponds to the best-fit position angle, |$\theta$|⁠, of the major axis. The black ellipse in the bottom left corner represents the beam.

Among the 20 sources with SNR |$\gt 8$|⁠, our visual classification of their velocity maps places eight of them in Class II (ALESS 003.1, 006.1, 009.1, 034.1, 065.1, 088.1, 101.1, 112.1; see Section 2.2) because they either display complex velocity fields or the resolution is too poor to characterize them. For five of these seven sources – those with sufficient resolution – our model converged to a solution that was rather unphysical, specifically the velocity dispersion converged to values |$\gt 200$| km s|$^{-1}$|⁠. Visually inspecting the residual 3D cubes, using the best-fit model, it is not obvious that the fit was unsuccessful. However, if we instead inspect the 1D spectrum we find that the model is not able to reproduce the data: the wings of the model spectrum are wider compared to the data, effectively trying to model these sources as if they are dispersion dominated. Somewhat similar behaviour of the galpak3d model (i.e. converging to non-physical values, in addition to the velocity dispersion, for the maximum rotational velocity and/or the effective radius as well) was also reported in Hogan et al. (2021) who studied the kinematics of a sample of (ultra) luminous infrared galaxies (U/LIRGs) at |$z\sim$| 2–2.5. For these five sources, the inability of our model to fit the data is not surprising. The model we use does not have the flexibility to account for complex features that deviate from the simple assumption of a regularly rotating disc and can therefore converge to solutions that are not physical, but happen to fit the data better than a model with more physical parameters. To circumvent this we also tried different functional forms for both the flux profile (e.g. Gaussian) and the velocity profile (e.g. Freeman model), but the same behaviour was observed. As the kinematics of these systems are not well-described by a simple rotating disc, we suggest that the complex kinematics in these sources may indicate that they are currently in a merger phase or recently had a merging event which led to their gas dynamics being significantly disturbed from ordered circular motions. Additional support for this hypothesis may be provided by archival Hubble Space Telescope (HST) observations of some of these sources (see Fig. A1), that show clumpy structures.⁹

In addition to these five sources with complex kinematics, our attempt to model a further source, ALESS 67.1, also resulted in large residuals in the velocity map. The kinematics of ALESS 67.1 were previously analysed by Chen et al. (2017) where the authors suggest that this source is consistent with a merger scenario. In addition, the HST image of this source reveals tidal features that considerably strengthen this hypothesis. Therefore, it is not surprising that our thick rotating disc model did not fit the data well. In this work, we also classify ALESS 67.1 as a likely merger and remove it from the sample, leaving 11 galaxies whose kinematics appear to be well described by a rotating disc model. We also note that another source, ALESS 98.1, displays characteristic that might indicate it is a merger (e.g. clumps in the PV diagram shown in Appendix  C). However, as these are not significant to make this claim with certainty we choose to include it in our final sample of discs. We caution that in a recent study by Rizzo et al. (2022) the authors showed that mergers can be missclassified as rotationally supported discs when the resolution of the data is low, which could be the case for some of our source in this work. Nevertheless, as this is the best we can do with the current data at hand for the remainder of this paper we assume that we have correctly classified these 11 sources as rotationally supported discs. We further examine the reliability of our recovered parameters by applying our modelling procedure on realistic simulated observations that mimic our observed data in both resolution and SNR. The results from our simulations are shown in Appendix  B where we find that we can accurately recover the true parameters without any biases.

4 RESULTS AND DISCUSSION

To summarize the modelling in the previous section: from our parent sample of 30 DSFGs, we find that we are unable to constrain the kinematics of the 10 galaxies with integrated CO emission line SNRs of SNR < 8 (Class III). From the remaining 20 DSFGs with CO SNR > 8, we have 11 where our rotating disc model provides a good description of the kinematics, indicating that the molecular gas reservoirs in these systems are likely to be in rotationally supported discs.

4.1 Dynamical parameter relations

In this section, we discuss the dynamical properties of our sample using various scaling relations involving their observable properties. We also compare them with less extreme star-forming samples from the literature to see how they relate to the wider population in terms of their dynamics.

The best-fit maximum posterior (MP) model parameters from the non-linear searches for the 11 galaxies, for which our thick rotating disc model provides a good fit, are presented in Table 2. We used uniform priors for all the model parameters except for the total intensity, I, for which we used a log-uniform prior. The MP solution is computed as the median of the 1D marginalized posterior distributions of each parameter. The upper and lower limits on the MP solution reflect the 1|$\sigma$| error of these parameters and are computed from the 68 per cent percentile of the posterior distributions.

Table 2.

The best-fit maximum posterior (MP) parameters, |$\boldsymbol{\theta }$|⁠, of our model and their associated 1|$\sigma$| errors (see the main text). An inclination of |$i = 0$| (deg) means that the source is face-on, while |$i = 90$| (deg) edge-on and the position angle, |$\theta$|⁠, is defined counterclock-wise from North. In the last two columns, we list the circular velocity computed at twice the effective radius (equation 8) and the dynamical mass computed at a radius, |$r=$| 10 kpc (equation 9). We note that we have adopted a minimum 20 per cent uncertainty on the reported |$r_{\rm e}$| in order to reflect assumptions made during fitting.

ID	|$r_{\rm e}$|	|$\theta$|	i	|$V_{\rm max}$|	|$\sigma$|	|$V_{\rm circ}(r=2\, r_{\rm e})$|	|$M_{\rm dyn}(r=10 \, {\rm kpc})$|
	(arcsec)	(deg)	(deg)	(km s|$^{-1}$|⁠)	(km s|$^{-1}$|⁠)	(km s|$^{-1}$|⁠)	(⁠|$\rm 10^{11} \, M_{\odot }$|⁠)
007.1	|$0.61_{\, -0.12}^{\, +0.12}$|	|$153_{\, -8}^{\, +8}$|	|$34_{\, -7}^{\, +6}$|	|$481_{\, -72}^{\, +86}$|	|$42_{\, -18}^{\, +17}$|	399 |$\pm$| 68	3.9 |$\pm$| 1.3
017.1	|$0.41_{\, -0.08}^{\, +0.08}$|	|$32_{\, -15}^{\, +13}$|	|$36_{\, -15}^{\, +14}$|	|$417_{\, -99}^{\, +150}$|	|$73_{\, -16}^{\, +16}$|	255 |$\pm$| 81	2.2 |$\pm$| 1.3
022.1	|$0.44_{\, -0.09}^{\, +0.09}$|	|$57_{\, -7}^{\, +7}$|	|$34_{\, -10}^{\, +8}$|	|$347_{\, -87}^{\, +85}$|	|$15_{\, -12}^{\, +10}$|	334 |$\pm$| 84	2.8 |$\pm$| 1.4
041.1	|$0.29_{\, -0.06}^{\, +0.06}$|	|$75_{\, -2}^{\, +3}$|	|$69_{\, -4}^{\, +4}$|	|$385_{\, -8}^{\, +10}$|	|$50_{\, -7}^{\, +7}$|	400 |$\pm$| 15	4.0 |$\pm$| 0.3
049.1	|$0.16_{\, -0.06}^{\, +0.06}$|	|$115_{\, -20}^{\, +16}$|	|$69_{\, -10}^{\, +15}$|	|$371_{\, -55}^{\, +50}$|	|$78_{\, -12}^{\, +19}$|	346 |$\pm$| 54	5.2 |$\pm$| 1.3
062.2	|$0.28_{\, -0.15}^{\, +0.06}$|	|$166_{\, -6}^{\, +7}$|	|$73_{\, -9}^{\, +17}$|	|$254_{\, -167}^{\, +63}$|	|$99_{\, -66}^{\, +23}$|	334 |$\pm$| 135	4.5 |$\pm$| 3.2
066.1	|$0.43_{\, -0.09}^{\, +0.09}$|	|$109_{\, -46}^{\, +77}$|	|$52_{\, -14}^{\, +15}$|	|$420_{\, -123}^{\, +138}$|	|$92_{\, -18}^{\, +18}$|	344 |$\pm$| 85	3.8 |$\pm$| 1.8
071.1	|$0.48_{\, -0.10}^{\, +0.10}$|	|$138_{\, -4}^{\, +5}$|	|$71_{\, -8}^{\, +7}$|	|$368_{\, -19}^{\, +22}$|	|$67_{\, -24}^{\, +24}$|	381 |$\pm$| 25	3.6 |$\pm$| 0.5
075.1	|$0.32_{\, -0.06}^{\, +0.06}$|	|$115_{\, -2}^{\, +3}$|	|$55_{\, -6}^{\, +4}$|	|$386_{\, -37}^{\, +27}$|	|$108_{\, -10}^{\, +9}$|	389 |$\pm$| 29	4.8 |$\pm$| 0.6
098.1	|$0.22_{\, -0.04}^{\, +0.04}$|	|$100_{\, -3}^{\, +3}$|	|$41_{\, -13}^{\, +7}$|	|$555_{\, -24}^{\, +48}$|	|$26_{\, -17}^{\, +17}$|	519 |$\pm$| 44	7.1 |$\pm$| 1.2
122.1	|$0.62_{\, -0.12}^{\, +0.12}$|	|$86_{\, -6}^{\, +7}$|	|$55_{\, -6}^{\, +8}$|	|$564_{\, -17}^{\, +44}$|	|$157_{\, -18}^{\, +19}$|	533 |$\pm$| 37	6.5 |$\pm$| 0.8

Along with the best-fit model parameters in Table 2, we list some additional quantities, namely the circular velocity computed at a twice the effective radius and the dynamical mass computed within a radius of |$r=$| 10 kpc. The circular velocity, |$V_{\rm circ}$|⁠, is defined as the rotational velocity corrected for asymmetric drift (Burkert et al. 2010). Several studies have suggested that star-forming disc galaxies at high redshift are typically more turbulent (e.g. Wisnioski et al. 2015; Johnson et al. 2016) compared to local analogues. This turbulent motion needs to be accounted for as it contributes to the dynamical support of a system. Assuming a constant velocity dispersion profile and considering that the vertical profile of our thick disc model does not depend on the radius, one has

$$\begin{eqnarray} V_{\rm circ} (r) = \sqrt{ \, V_{\rm rot}^2(r) - \sigma ^2 \frac{\mathrm{ d}ln \Sigma }{\mathrm{ d}lnr}} \, , \end{eqnarray}$$

(7)

which for an exponential profile, |$\Sigma \propto \exp {\left(-r / r_e\right)}$|⁠, reduces to

$$\begin{eqnarray} V_{\rm circ} (r) = \sqrt{ \, V_{\rm rot}^2(r) + 1.68 \sigma ^2 \left(\frac{r}{r_e} \right)} \, , \end{eqnarray}$$

(8)

where |$V_{\rm rot}$| is the rotational speed of the gas (its functional form is defined in Section 3.2), |$\sigma$| is the velocity dispersion, and |$r_e$| is the effective radius of the disc. The dynamical mass, |$M_{\rm dyn}$|⁠, which is a measure of the total mass of the system, is computed as

$$\begin{eqnarray} M_{\rm dyn}(r) = \frac{r V^2_{\rm circ}(r)}{G} , \, \end{eqnarray}$$

(9)

where again we use the circular rotational velocity to account for the effect on turbulent motions.

4.1.1 |$V_{\rm circ} - \sigma$|

We first look at the relation in our sample between the circular velocities, |$V_{\rm circ}$|⁠, and the velocity dispersion, |$\sigma$|⁠, which is shown in Fig. 5 including a comparison sample of star-forming galaxies from KMOS|$^{\rm 3D}$|⁠. The figure demonstrates that our sample is not comparable to‘typical’ star-forming galaxies, but instead represents the most massive sources as judged by their dynamics. The median velocity dispersion of our sample (⁠|$\sim$| 75 km s|$^{-1}$|⁠; blue arrow) is significantly higher than seen in the KMOS|$^{\rm 3D}$| sample (⁠|$\sim$| 40 km s|$^{-1}$|⁠; grey arrow). We argue that this to be a consequence of our selection; the most massive star-forming disc galaxies at any redshift also have the higher velocity dispersions. Indeed, if we apply a cut in circular velocity, |$V_{\rm circ}$|⁠, to the KMOS|$^{\rm 3D}$| sample above 325 km s|$^{-1}$| we find that the median velocity dispersion increases and becomes consistent with our sample (this is important to keep in mind when later we look at the evolution of the velocity dispersion with redshift).

Figure 5.

Circular rotational velocity, |$V_{\rm circ}$|⁠, as a function of the velocity dispersion. The blue points correspond to our sample of massive dusty star-forming galaxies while the grey points correspond to sources from the KMOS|$^{\rm 3D}$| survey which are more representative of the ‘typical’ star-forming galaxies at |$z \sim 2$| (⁠|$V_{\rm circ}$| is computed at twice the effective radius, |$r = 2 \times R_{\rm eff}$|⁠, for both samples). The histograms at the bottom of the figure correspond to the distribution of velocity dispersions for the two populations (arrows indicate the medians). The dotted lines are the |$1:1$| and |$2:1$| relation between the two properties, highlighting that all source shown in this figure are classified as rotationally supported.

We also looked at the evolution with redshift of the ratio of the circular velocity with the velocity dispersion, |$V / \sigma$|⁠. Previous studies have suggested that this ratio decreases with redshift and this evolution is a consequence of the increased velocity dispersion of galaxies at higher redshifts (Price et al. 2016; Turner et al. 2017; Wisnioski et al. 2019; Hogan et al. 2021). However, we do not see a clear trend in our sample (a very mild negative trend is observed when we averaged points in bins of redshift, but this was not statistically significant).

4.1.2 |$\sigma - z$|

We next take a look at the evolution of the intrinsic velocity dispersion, |$\sigma$|⁠, with redshift, z, for our sample. Again, we use sources from the KMOS|$^{\rm 3D}$| survey as the reference sample of typical star-forming galaxies at high redshift (grey points). Previous studies have claimed that the velocity dispersion of star-forming, rotationally supported galaxies increases with redshift (e.g. Übler et al. 2019) and proposed that the increased turbulence is mainly the result of gravitational instabilities (e.g. Übler et al. 2019). In order to quantify the redshift dependence of the velocity dispersion for our sample we fit a linear relation, |$\sigma = a z + b$| and plot this in Fig. 6, where the shaded region corresponds to the 1|$\sigma$| error.

Figure 6.

Velocity dispersion, |$\sigma _0$|⁠, as a function of redshift, z, for the 11 sources with SNR > 8 in our sample that have kinematics that are well described by a rotating disc model. These show a modest increase in |$\sigma$| with redshift, but this is not statistically significant. We also show points from the KMOS|$^{\rm 3D}$| survey (grey; Übler et al. 2019), where the grey line with its corresponding uncertainty (grey shaded regions; 1σ and 2σ) is the best-fit linear relation to these points as derived by Übler et al. (2019). The blue line and its corresponding uncertainty (blue shaded region; 1|$\sigma$|⁠) is the best-fit linear relation to our sample of sources. We also re-calculated averages from the KMOS|$^{\rm 3D}$| points, in the same three redshift bins, using only those sources with |$V_{\rm circ} \gt $| 325 km s|$^{-1}$| to match our sample. These are shown as the square black points and agree with the relation derived from our sample (i.e. consistent with no redshift evolution).

Our fitting suggests no evolution for the velocity dispersion with redshift for our sample (the slope is positive but not statistically significant). Furthermore, we note that our relation is above the one reported with the KMOS|$^{\rm 3D}$| sample, although still consistent within the 1|$\sigma$| errors. If one ignores selection effects, it will seem surprising that our CO-based relation is above the H|$\alpha$|-derived one from Übler et al. (2019). This is because most previous studies claim that velocity dispersions derived from observations of the molecular gas are typically lower than those derived from ionized gas (Levy et al. 2018; Girard et al. 2019, 2021). Here, however, we stress that we are not comparing the same galaxy populations in terms of their intrinsic dynamical properties. As we showed in Fig. 5 our sources are more massive than the typical star-forming galaxy population, and so it is not surprising that their velocity dispersions are also elevated. In order to make a more fair comparison we apply a cut in circular velocity, |$V_{\rm circ} \gt 325$| km s|$^{-1}$| to the KMOS|$^{\rm 3D}$| sample and recalculate the average points in the same bins of redshift as Übler et al. (2019). These are shown as the black points in Fig. 6 and agree with the relation derived from our sample (also suggesting no redshift evolution).

We note that in recent years there have been several high-resolution studies of the dynamics of both SFGs and DSFGs based on observations of various far-infrared emission lines (e.g. Übler et al. 2018; Rizzo et al. 2020, 2021, 2023; Fraternali et al. 2021; Lelli et al. 2021, 2023; Tsukui & Iguchi 2021; Xiao et al. 2022; Roman-Oliveira et al. 2023). Most of these studies report that large fractions of star-forming galaxies have gas kinematics consistent with dynamically cold discs (i.e. low velocity dispersions). However, in our sample, sources with low velocity dispersions (e.g. 007.1, 022.1, 098.1) do not constitute the majority. One potential reason for this difference is the lower angular resolution of the observations modelled in this work. We attempted to address this using simulations in Appendix  B, where we demonstrate that if our sources are indeed well described by a rotating disc model, we can accurately recover their true velocity dispersions without any biases. However, if the intrinsic kinematics of these systems are not disc-like, then this conclusion does not hold. Additionally, the majority of the aforementioned studies involve small samples or a mixed selection of sources (both main-sequence and starburst galaxies), making a direct comparison with the ALESS sample less straightforward. We, therefore, opted not to include sources from the above studies in Fig. 6 for two reasons: first, the majority of these studies focus on individual sources whose selection is different to those from ALESS. Secondly, one of the main aims of this figure is to encourage caution when comparing velocity dispersion estimates for different populations due to the influence sample selection on such comparisons. In our case, the sample we are studying perhaps represents a biased sub-set of the star-forming galaxy populations, specifically the most massive ones.

4.2 Dynamical constraints on |$\alpha _{\rm CO}$|

In addition to the dynamical scaling relations discussed above, our kinematic estimates of the masses of these galaxies also allow us to constrain a critical parameter used to estimate the gas masses of these systems: |$\alpha _{\rm CO}$|⁠.

The bulk of the molecular gas in the ISM of galaxies is in the form of molecular hydrogen, H|$_2$|⁠. Unfortunately, due to H|$_2$| being symmetric molecule, the conditions needed to observe its transitions in emission are extreme (⁠|$T \gt 500$| K or a strong UV radiation field), compared to typical ISM conditions (⁠|$T \sim 20$| K). In order to overcome this limitation most studies typically use CO, which is the second-most abundant molecule, as an indirect tracer of the total molecular gas in the ISM. The problem with using the CO molecule to estimate molecular gas masses is that one has to assume a conversion factor, otherwise known as the CO-to-|$H_2$| conversion factor, |$\alpha _{\rm CO}$| (Bolatto et al. 2013).

There have been several attempts in previous studies to estimate the value of |$\alpha _{\rm CO}$|⁠. The general picture that emerged from these studies is that there is a dichotomy in the measured value of |$\alpha _{\rm CO}$| between ‘normal’ star-forming galaxies (e.g. Daddi et al. 2010) and ‘starburst’ galaxies (e.g. Downes & Solomon 1998; Hodge et al. 2012; Bothwell et al. 2013). In addition, the value of |$\alpha _{\rm CO}$| has been shown to strongly depend on the gas phase metallicity (e.g. Bolatto et al. 2013; Combes 2018), among other factors (Tacconi, Genzel & Sternberg 2020).

Here, we attempt to estimate |$\alpha _{\rm CO}$| following the same approach as previous studies in the literature (e.g. Calistro Rivera et al. 2018). Specifically, we use the dynamical mass within a fixed radius, |$M_{\rm dyn}$|⁠, and equate it to the sum of the masses of the different galactic components (dark matter and baryons),

$$\begin{eqnarray} M_{\rm dyn} = M_{\rm gas} + M_{\star } + M_{\rm DM}. \end{eqnarray}$$

(10)

Expressing the molecular gas mass as, |$M_{\rm gas} = \alpha _{\rm CO} \, L_{CO(1-0)}$|⁠, and the dark matter mass as, |$M_{\rm DM} = f_{\rm dm} * M_{\rm dyn}$|⁠, the mass equation above, becomes,

$$\begin{eqnarray} \alpha _{\rm CO} = (M_{\rm dyn} - M_{\star } - f_{\rm dm} * M_{\rm dyn}) / L_{\rm CO(1-0)} \, , \end{eqnarray}$$

(11)

where |$L_{\rm CO(1-0)}$| is the CO (1–0) line luminosity and |$f_{\rm dm}$| is the dark matter fraction.

Before we calculate the value of |$\alpha _{\rm CO}$| we need to make some assumptions about the various terms on the right-hand side of equation (11). First, we chose a fixed radius, |$r =$| 10 kpc, within which we calculate the dynamical mass. This value is close to twice the median effective radius, |$r_e$|⁠, for the sample and thus should contain the bulk of the baryonic material in these systems. Secondly, we assume a fixed dark matter fraction, |$f_{\rm dm} = 0.25$|⁠, which is a reasonable value at the radius we choose to estimate the dynamical mass (e.g. Genzel et al. 2017; Sharma, Salucci & van de Ven 2021). In theory, instead of adopting a fixed value for the dark matter fraction, one could directly constrain this value by incorporating the dark matter contribution into the dynamical modelling analysis, as is done in some recent studies (e.g. Genzel et al. 2020; Bouché et al. 2022). However, our observations lack the necessary resolution and SNR to be able to independently constrain the contribution from both the baryonic and dark matter components. Finally, we use the stellar masses obtained via SED-fitting (da Cunha et al. 2015) which are listed in Table 1. We note, however, that stellar mass estimates can be uncertain, especially in starburst systems for which star formation histories (among other ingredients of the SED models) are poorly constrained. One approach to account for our lack of knowledge in the assumptions we make during the SED-fitting is to substitute the stellar mass in our mass equation with the term |$L_{\rm H} \left(\rm M_{\star } / L_{\rm H}\right)$|⁠, where |$L_{\rm H}$| is the rest-frame H-band luminosity corrected for dust obscuration. If we were to make this substitution, we could allow the ratio, |$M_{\star } / L_{\rm H}$|⁠, to be a free parameter and constrain it simultaneously with the |$\alpha _{\rm CO}$| as was done in a recent study by Calistro Rivera et al. (2018). Exploring this possibility is outside the scope of this paper, however, we note that this ratio is degenerate with the |$\alpha _{\rm CO}$| and will result in larger uncertainties for each individual source.

We compute the value of |$\alpha _{\rm CO}$| and its associated error (via bootstrapping) for each of the sources in our sample. We show these as posterior distributions in the upper panel of Fig. 7. We then combined the individual distribution to compute the joint constraint on |$\alpha _{\rm CO}$| for our sample, shown in the same figure. Our approach yields an estimate of |$\alpha _{\rm CO} = 0.74 \pm 0.37$|⁠, for a dark matter fraction of, |$f_{\rm dm} = 0.25$|⁠, which is consistent with estimates from previous studies (e.g. Tacconi et al. 2008; Hodge et al. 2012; Bolatto et al. 2013; Bothwell et al. 2013; Calistro Rivera et al. 2018; Birkin et al. 2021). We note, however, that in a recent study by Dunne et al. (2022) the authors report a value of |$\alpha _{\rm CO} = 4.0 \pm 0.1$|⁠, which is in tension with the limits we place here. Given the differences in sample selection between our work and theirs, where the later included both lensed and unlensed sources, understanding the origin of this difference is not straightforward and outside the scope of this work.

Figure 7.

The CO–H|$_2$| conversion factor, |$\alpha _{\rm CO}$|⁠, as a function of the dark-matter fraction, |$f_{\rm DM}$| (within a radius of |$r \lt 10$| kpc). The zoom in plot at the top of the figure shows the individual posteriors distributions (grey), for a dark matter fraction of the |$\alpha _{\rm CO}$| value for each source in our sample as well as the combined distribution (blue) which is computed by multiplying together all the individual distributions.

We also explored how the value of |$\alpha _{\rm CO}$| changes when varying the adopted dark matter fraction, which we show in the bottom panel of Fig. 7. This indicates that we can place a firm upper limit on |$\alpha _{\rm CO}$| of |$\alpha _{\rm CO} \lesssim 2$| (assuming no dark matter within 10 kpc radius). Similarly adopting a dark matter fraction, |$f_{\rm dm} \gt 0.5$|⁠, results in the value of |$\alpha _{\rm CO}$| becoming negative.¹⁰

4.3 Tully–Fisher relation at high redshift

Next we use the results we obtained above from our dynamical modelling analysis of the 11 DSFGs with disc-like kinematics, specifically the rotational velocity corrected for inclination, the velocity dispersion and the global constraint on |$\alpha _{\rm CO}$| and hence the gas masses, to study the stellar and baryonic mass Tully–Fisher relations (sTFR/bTFR; McGaugh et al. 2000) for this population.

The TFR holds important information about the interplay between the build-up of galaxies and the dark matter haloes in which they grow. It connects an observable of the baryonic mass, where in this case we consider both the stellar mass alone and the stellar + gas mass, with a proxy for the total potential of the halo, the circular velocity. Following Übler et al. (2017) we use the circular velocity as it is directly connected to the total potential of the halo for sources with elevated velocity dispersions, as are the sources in this study. In the left and right panels of Fig. 8, we show the sTFR and bTFR, respectively, and compare these to the less active galaxies from the KMOS|$^{\rm 3D}$| survey.

Figure 8.

The stellar and baryonic TFR, sTFR (left) and bTFR (right), respectively for the 11 disc-like sources with SNR > 8 in our sample (blue points) as well as for |$z \sim 2.3$| sources from the KMOS|$^{\rm 3D}$| (black points). The best-fit power-law models to our sample are shown as the blue solid lines. For the sTFR and bTFR, compared to the local relations from Reyes et al. (2011) and Lelli et al. (2016), we find offsets in the normalization of |$\Delta = -0.28 \pm 0.24$| and |$\Delta = -0.28 \pm 0.18$|⁠, respectively, indicating modest evolution of the masses at fixed circular velocity. These offsets are consistent with those seen by Übler et al. (2017) for less active galaxies.

In order to quantify the relation between the observed mass (stellar or baryonic) and the circular velocity, we fit our data points using a linear regression model of the form:

$$\begin{eqnarray} M = V_{\rm circ}^a 10^b \, , \end{eqnarray}$$

(12)

where a and b are constant parameters, which we would like to constrain, and M corresponds to either the stellar or the baryonic mass, where the latter is the sum of the stellar and gas mass. For the fitting, we use an orthogonal distance regression method taking into account the error in both the x and y coordinates (we use the scipy.odr package). During the fitting process we actually fix the value of the slope, a, to values found in local studies and only optimize the offset, b. This is a common practice in high-redshift studies of the TFR (e.g. Cresci et al. 2009; Price et al. 2016; Tiley et al. 2016; Übler et al. 2017) simply because the dynamic range in mass is typically very limited (which is particularly evident for our sample of sources, mainly probing the high-mass end of the TFRs) and is therefore not easy to robustly constrain. The results are then quoted in terms of a zero-point offset compared to the |$z\sim$| 0 value, |$\Delta b$|⁠. For the sTFR and bTFR we adopt the local slopes of |$\alpha = 3.60$| (Reyes et al. 2011) and |$\alpha = 3.75$| (Lelli, McGaugh & Schombert 2016), respectively. These are the same values that were used in Übler et al. (2017) with which we want to eventually compare our results.

The best-fit offsets, b, we obtain from our fitting are |$b = 2.08 \pm 0.24$| and |$b = 1.90 \pm 0.18$|⁠, translating to zero-point offsets relative to the |$z\sim$| 0 values of |$\Delta b = -0.28 \pm 0.24$| and |$\Delta b = -0.28 \pm 0.18$| for the sTFR and bTFR, respectively. The power-law models are shown in each of the panels in Fig. 8, with the shaded area corresponding to the 1|$\sigma$| errors. The values we obtain are in good agreement with those found for less-active galaxies in Übler et al. (2017) for both the sTFR (⁠|$\Delta b = -0.42 \pm 0.05$|⁠) and bTFR (⁠|$\Delta b = -0.27 \pm 0.05$|⁠).¹¹

We caution here that our sample of sources spans a wide range in redshift compared to the sample in Übler et al. (2017). This could potentially introduce a bias in the zero-point offsets that we infer and will definitely contribute to the scatter of the relation. When splitting our sample in redshift, below and above |$z \sim$| 2.5 (roughly the median for our sample), and repeating the power-law fit we find no significant differences in the inferred zero-point offsets; however, the errors are significant due to the small number of data points in each redshift bin.

The interpretation of our results is rather straightforward. Our population of massive dusty star-forming galaxies with disc-like kinematics appears to represent the high-mass end of the TFR traced by the more typical star-forming galaxies at high redshift. While the DSFGs have higher star formation rates, higher gas and stellar masses, than the more typical galaxies surveyed by KMOS|$^{\rm 3D}$|⁠, they also reside in more massive dark matter haloes, which ultimately translates to baryon fractions (in the case of the bTFR for example) that are similar to these more typical star-forming galaxies. Hence, the key conclusion from this comparison is the massive nature of these disc-like DSFGs, that appear to just extend the scaling relations found in less massive galaxies.

4.4 The descendants of massive DSFGs

In this final sub-section, we use our kinematic modelling of the DSFG galaxies to test if there is an evolutionary link between early-type galaxies in the local Universe and massive dusty star-forming galaxies at high redshift by comparing the dynamical properties of these two populations.

In order to investigate the evolutionary link between between these populations we compare their distribution in terms of their baryonic mass, |$M_{\rm b}$|⁠, and central velocity dispersions (i.e. velocity dispersion, |$\sigma _e$|⁠, within the half-light radius, |$R_e$|⁠), frequently referred to as the |$M_{\rm b}$|–|$\sigma$| relation. The connection between these two galaxy populations, using their distributions on this plane, was previously studied in a statistical manner by Birkin et al. (2021). Here, instead we will use the results from the dynamical modelling analysis which allow us to infer both the rotational speed of the gas as well as the velocity dispersion of individual galaxies.

For comparison we use samples of local massive early-type galaxies from two surveys, ATLAS|$^{\rm 3D}$| (Cappellari et al. 2011) and MASSIVE (Ma et al. 2014). These two samples together span a wide range in stellar mass |$M_{\star } \sim 10^{10}$|–|$10^{12}$| M|$_{\odot }$|⁠.The ATLAS|$^{\rm 3D}$| survey (Cappellari et al. 2011) is a volume-limited (⁠|$D \lt $| 42 Mpc) study of 260 early-type galaxies selected to have absolute K-band magnitude of |$M_K \lt -$|21.5. Stellar central velocity dispersion measurements for the ATLAS|$^{\rm 3D}$| sources are taken from Cappellari et al. (2013). The MASSIVE survey (Ma et al. 2014) is a volume-limited (⁠|$D \lt $| 108 Mpc) survey of 116 early-type galaxies where galaxies were selected to have absolute K-band magnitude of |$M_K \lt -$|25.3 mag. Stellar central velocity dispersion for these are taken from Veale et al. (2018). For the samples of early-type galaxies, stellar masses are computed from their K-band magnitude which provides a fairly robust approximation of their stellar mass (the K band is less sensitive to dust absorption compared to optical wavelengths and has relatively small mass-to-light ratio variations with star formation history).

In order to place our sample of massive DSFGs at high redshift on the |$M-\sigma$| plane, along with the samples of early-type galaxies, we first need to make three assumptions. The first assumption that we make is that by |$z \sim$| 0 all of the available gas that our sources have will be depleted in order to form stars, without any loss of mass. Therefore, their total baryonic mass at |$z \sim 0$| will be the sum of two components, their stellar and gas masses at the redshift of detection. We note, however, that some of their gas might be expelled due to their strong star formation activity or if they undergo a Quasi Stellar Object (QSO) phase, although it may also subsequently be re-accreted or they may accrete stars and gas though future mergers.

The second assumption is that the energy of a system will be conserved when transitioning from the active star-forming phase where the system is dominated by rotation, to a dispersion-dominated phase at |$z = 0$|⁠. Under this assumption, we can use the virial theorem (Epinat et al. 2009) and equate their virial masses at these two redshift epochs. This allow us to estimate the velocity dispersion of the system in this dispersion-dominated phase from,

$$\begin{eqnarray} \beta \frac{\sigma ^2_{\rm e, \, z_0} R_{\rm e, \, z_0}}{G} = \frac{V^2_{\rm circ, max, \, z} R_{\rm e, \, z}}{G}, \end{eqnarray}$$

(13)

where |$V_{\rm circ, max, \, z}$| is the maximum circular velocity of the system at the observed redshift (these are the values we inferred from our dynamical modelling analysis; see Table 2), |$\beta$| is a constant which for a spherical system takes a value of |$\beta \sim 5$| (Cappellari et al. 2013), |$R_{e, \, z}$| and |$R_{\rm e, \, z_0}$| are the radii before and after the system transitioned to a dispersion-dominated phase.

The final assumption is that the radius of the system after it transitions to a dispersion-dominated phase will be roughly equal the radius it had when it was still in the rotation-dominated phase. We make this assumption to further simplify equation (13) but we note that these systems potentially do evolve in size (e.g. through minor mergers, which do not significantly increase their stellar mass). However, since the velocity dispersion in early-type galaxies is measured only in the central region, our assumption can be justified.

Using the assumptions described above, we can place our sample of DSFGs on the same M–|$\sigma$| plot along with the samples of local early-type galaxies in Fig. 9. The estimated baryonic masses of the DSFGs place them at the highest masses seen for early-type (or any) galaxies at |$z\sim$| 0, |$M_{\rm b} \gtrsim$| 10|$^{11}$| M|$_\odot$| (e.g. Ma et al. 2014), and the dynamical masses we infer for the DSFG support these high masses. Indeed, based on the gas dynamics, the DSFGs scatter around the trend line from Cappellari et al. (2013) that describes the relationship between the total stellar mass and the velocity dispersion of stars in local early-type galaxies (see Fig. 9). The median values of our sample, |$M_{\rm b}=(3.6\pm 2.9)\times 10^{11}$| M|$_\odot$| and |$\sigma =$| 186|$\pm$|59 km s|$^{-1}$| are within |$\sim$|1|$\sigma$| of the |$z\sim$| 0 trend. We therefore conclude that there is good agreement on the M–|$\sigma$| plane between the predicted properties of the DSFGs (subject to the assumptions listed above) and those found for the most massive early-type galaxies found in the Universe at low redshifts. This provides further support for the existence of an evolutionary link between the formation of the most massive early-type galaxies in the local Universe and massive gas-rich and disc-like DSFGs at high redshift.

Figure 9.

The M–|$\sigma$| relation for DSFGs (blue unfilled points) and early-type galaxies from two surveys: ATLAS|$^{\rm 3D}$| (black points) and MASSIVE (purple points). For early-type galaxies, |$\sigma$| correspond to the central velocity dispersion and |$M_b$| is the stellar mass (the gas mass in these systems is negligible compared to their stellar mass, where the latter is computed from their K-band luminosity). For DSFGs on the other hand, the quantity |$\sigma$| is computed from equation (13) and |$M_b$| is the sum of the gas and stellar masses of the galaxies (this assumes that by |$z \sim$| 0 these galaxies will have converted all of their gas to stars). For comparison, the shaded region corresponds to median trend found by Birkin et al. (2021) using the parent CO survey sample from which the data analysed in this work were selected (see Section 2). The larger blue filled point represents the median of all the blue unfilled points where the errors were calculated via boostrap. The solid line shows the trend followed by the local early-type galaxies and is computed using equation (5) from Cappellari et al. (2013). The parameters of that equation were optimized from fitting directly to the early-type galaxies plotted and as can be seen from the figure it adequately describes points at the high-mass end of the relation.

Finally, we note that a crude estimate of the space density of the sources in our sample yields a rough volume density of |$\sim$| 1 |$\times$| |$10^{-5}$| Mpc|$^{-3}$| (based on a sub-set of |$\sim$|110 sources from the parent ALESS survey over 0.35 deg|$^2$|⁠, which matches the median |$S_{\rm 870\mu m}=$| 4.5 mJy for our sample, and adopting a full width at half-maximum (FWHM) of the redshift distribution for our sources spanning |$z=$| 2.0–3.9, giving a survey volume of 8.3 |$\times$| 10|$^6$| Mpc|$^3$|⁠). Taking the median lifetime of |$\sim$| 400 Myr (twice the gas depletion time-scale) for SMGs from (Birkin et al. 2021), we need to apply a duty cycle correction of |$\sim$| 4 |$\times$| to account for the duration of the SMG phase over the 1.7 Gyr corresponding to |$z=$| 2.0–3.9, indicating a descendant volume density of |$\sim$| 4–5 |$\times$| 10|$^{-5}$| Mpc|$^{-3}$|⁠. This falls between the volume density of galaxies in the MASSIVE survey of |$\sim$| 2 |$\times$| 10|$^{-5}$| Mpc|$^{-3}$|⁠, which are typically more massive than our sources, and |$\sim$| 10 |$\times$| 10|$^{-5}$| Mpc|$^{-3}$| for the sub-set of the ATLAS|$^{\rm 3D}$| galaxies more massive than |$M_\ast \ge$| 10|$^{11}$| M|$_\odot$|⁠. This suggests that a substantial fraction of local early-type galaxies more massive than |$M_\ast \sim$| 10|$^{11}$| M|$_\odot$| could be formed through the massive DSFG population in our sample (see McAlpine et al. 2019; Dudzevičiūtė et al. 2020).

5 SUMMARY AND CONCLUSIONS

We have presented a method to model the dynamics of sources in 3D, from interferometric observations of emission lines, that builds upon existing methods (i.e. galpak3d; Bouché et al. 2015) with the added extension that the analysis is performed directly in the |$uv$|-plane. Performing the modelling in the |$uv$|-plane results in more realistic estimates of errors on the inferred best-fit parameters of the model we are fitting to the data (e.g. ALESS 122.1; see Section 3.5), because the errors on the visibilities are well defined (i.e. Gaussian errors). Our method does not require any image pre-processing, such as the CLEAN task from casa (McMullin et al. 2007), which can depend on user assumptions (e.g. masking, threshold, etc.), a process that is not reversible (see Tsukui et al. 2023, for an alternative approach to dealing with correlated noise). Our main science results are summarized below:

We find that 11 of the 20 sources that satisfy our SNR selection criteria, SNR > 8 in their integrated CO emission, can be reasonably well-fit by a simple thick rotating disc model. We classify all of these sources as rotationally supported discs, suggesting that |$\sim$| 55 per cent of massive DSFGs fall in this category. This fraction is consistent with the behaviour seen at high masses in previous studies of more typical star-forming galaxies (e.g. Wisnioski et al. 2015). We note, however, that at this stage, this fraction should only be considered an upper limit because the modest resolution of our data can lead, in some cases, to a misclassification of mergers as discs (see Rizzo et al. 2022). Among the remaining eight sources in our sample, six are classified as potential mergers, with the majority of them displaying complex features in their velocity maps, and for two sources we are not able to assign a classification due to resolution limitations.

We used the best-fit parameters of our kinematic model to investigate the dynamical state of our sources and compare them with samples in the literature. First, we looked at the relation between the circular velocity, |$V_{\rm circ}$|⁠, and the velocity dispersion, |$\sigma$|⁠. Comparing our sample with more typical star-forming samples from the literature (e.g. KMOS|$^{\rm 3D}$|⁠) we find that the DSFGs occupy the high velocity part of the distribution on this plane, with velocity dispersions that partly overlap but extend to higher dispersions than typical star-forming galaxies. Although, we find that our sample has an elevated median velocity dispersions compared to the ‘typical’ star-forming galaxy population, their similarly higher circular velocities means that the ratio of these two quantities, |$V/\sigma$|⁠, has values consistent with the less actively star-forming population. We also looked at the variation in the typical velocity dispersion, |$\sigma$|⁠, with redshift and found little evidence for evolution. We highlighted that, although our inferred |$\sigma$|–z relation lies above that derived from H|$\alpha$| gas kinematics, this may be a consequence of the higher masses of our sample, although, the low-angular resolution of our observations may also play a part. Nevertheless, the relation we infer is not statistically different from that reported by Übler et al. (2019) and so we can not draw any conclusions about the difference in velocity dispersion between molecular and ionized gas kinematics.

Combining the dynamical and physical properties of our sample we make the following conclusions:

• We were able to constrain the median |$\alpha _{\rm CO}$| conversion factor, a quantity that is used to convert CO(1–0) line luminosities to gas masses, to have a value of |$\alpha _{\rm CO} = 0.74 \pm 0.37$|⁠. This measurement is consistent with previous estimates based on samples of similarly far-IR luminous galaxies (e.g. Tacconi et al. 2008; Hodge et al. 2012; Bolatto et al. 2013; Bothwell et al. 2013; Calistro Rivera et al. 2018; Birkin et al. 2021).

• We studied the stellar and baryonic TFRs, using our sample of massive DSFGs at high redshift, and compared it to a sample of more ‘typical’ star-forming galaxies at |$z \sim$| 2.3 from the KMOS|$^{\rm 3D}$| survey. Our DSFGs occupy the high-mass end of these relations, but have normalizations that are consistent with those measured for the KMOS|$^{\rm 3D}$| sources. This shows that these DSFGs represent some of the most massive disc-like galaxies that have existed in the Universe.

• Finally, we have shown that under three reasonable assumptions (see Section 4.4) the most massive DSFGs at high redshift (⁠|$z \sim$| 1.2–4.7) have remarkably similar distributions on the baryonic mass versus velocity dispersion plane to the most massive (⁠|$M_{\rm b} \gtrsim 10^{11}$| M|$_{\odot }$|⁠) early-type galaxies in the local Universe. The apparent agreement between the distributions of these two populations, and their similar space densities, adds further evidence to support the hypothesis that there is an evolutionary link between these two galaxy populations (e.g. Dudzevičiūtė et al. 2020; Birkin et al. 2021).

ACKNOWLEDGEMENTS

AA was supported by ERC Advanced Investigator grant, DMIDAS [GA 786910], to C.S. Frenk. AA, JEB, IRS, AMS, and JN acknowledge support from STFC (ST/T000244/1 and ST/X001075/1). This paper makes use of the DiRAC Data Centric system at Durham University, operated by the Institute for Computational Cosmology on behalf of the STFC DiRAC HPC Facility (www.dirac.ac.uk). This equipment was funded by BIS National E-infrastructure capital grant ST/K00042X/1, STFC capital grants ST/H008519/1 and ST/K00087X/1, STFC DiRAC Operations grant ST/K003267/1, and Durham University. DiRAC is part of the National E-Infrastructure. This paper makes use of the following ALMA data: ADS/JAO.ALMA # 2016.1.00564.S, 2016.1.00754.S, 2017.1.01163.S, 2017.1.01471.S, and 2017.1.01512.S. ALMA is a partnership of ESO (representing its member states), NSF (USA), and NINS (Japan), together with NRC (Canada), MOST and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO, and NAOJ. This work was performed using the Cambridge Service for Data Driven Discovery (CSD3), part of which is operated by the University of Cambridge Research Computing on behalf of the STFC DiRAC HPC Facility (www.dirac.ac.uk). The DiRAC component of CSD3 was funded by BEIS capital funding via STFC capital grants ST/P002307/1 and ST/R002452/1 and STFC operations grant ST/R00689X/1. DiRAC is part of the National e-Infrastructure. This work used the DiRAC@Durham facility managed by the Institute for Computational Cosmology on behalf of the STFC DiRAC HPC Facility (www.dirac.ac.uk). The equipment was funded by BEIS capital funding via STFC capital grants ST/P002293/1 and ST/R002371/1, Durham University and STFC operations grant ST/R000832/1. DiRAC is part of the National e-Infrastructure.

Software:numpy (van der Walt, Colbert & Varoquaux 2011), scipy (Jones, Oliphant & Peterson 2001), and astropy (Astropy Collaboration 2013).

DATA AVAILABILITY

The ALMA data that were used in this study are available from the ALMA Science Archive at https://almascience.eso.org/asax/.

Footnotes

1

We use the term ‘typical’ to characterize star-forming galaxies that have lower dust masses and star formation rates compared to the average population of DSFGs selected at sub-millimetre wavelengths.

2

Some sources were observed in multiple projects (e.g. ALESS 041.1, 075.1). In these cases, we report the synthesized beam of the best available data set (i.e. highest resolution and/or signal-to-noise ratio), which are the ones we use in our modelling analysis.

3

The error on the velocity integrated flux is computed as the quadratic sum of the errors in all channels of the cube that were used to measure the flux (i.e. the width of the spectrum). The error on the flux in each channel is computed as the standard deviation of fluxes estimated in |$N =$| 100 regions, that have the same size as the aperture that was used to measure the flux, which do not contain any of the emission (otherwise referred to as the ‘random aperture method’; Tsukui et al. 2023).

4

In order to determine if a pixel is included in our velocity maps, we compute two individual |$\chi ^2$| values: one for the Gaussian model, |$\chi ^2_{\rm Gauss}$|⁠, and another for a flat line going through zero, |$\chi ^2_{0}$| (i.e. no signal). If |$\sqrt{\chi ^2_{0} - \chi ^2_{\rm Gauss}} \gt 5$| for a pixel, then it is included in our final maps. The |$\sigma$| value that is used to calculate the |$\chi ^2$| is computed as the standard deviation of pixel fluxes excluding the frequency range of the emission line.

5

Available at http://galpak.irap.omp.eu/.

6

This is hard-coded in the source code and can not be used as a free parameter in our model. We note that we have not explored how this parameter affects our inference for the rest of the model parameters.

7

Available at https://github.com/jyhmiinlin/pynufft.

8

Available at https://github.com/rhayes777/PyAutoFit.

9

We caution that dusty galaxies showing clumpy structure in the rest-frame UV/optical does not necessarily mean that they are in the process of merging or recently had a merging event, but this morphology may instead reflect structured dust extinction within these systems.

10

We note that there is a significant scatter among the constraints from the individual sources with some of them centred on negative values. Negative |$\alpha _{\rm CO}$| values occur because the stellar masses, which are estimated via SED-fitting, are in some cases higher than the estimated dynamical masses, which can occur due to uncertainties in these measurements (e.g. Wuyts et al. 2016).

11

We note that the scatter in the sTFR appears larger than the individual errors, suggesting there may be a second parameter at work. However, we have searched for correlations of the offsets from the relation with other properties derived from the SED fitting (e.g. age, |$A_v$|⁠) but found no significant correlations.

APPENDIX A: HST MORPHOLOGIES

In Section 3.5, we discussed the results from our dynamical modelling analysis for the 20 sources with SNR|$\gt 8$|⁠. We argued there that when our non-linear search convergences to unphysical values, |$\sigma \gt 200$| km s|$^{-1}$|⁠, this indicates that these sources can not be described by our simple rotating disc model. For three of the sources we have archival HST observations in the |$H_{160}$| band (ID 12866; Chen et al. 2015; Hodge et al. 2016, 2019), which we show in Fig. A1. These sources appear to be breaking into multiple components in the HST images (5–6 individual clumps from fitting the images with GALFIT; Chen et al. 2015), a feature that can be used to strengthen our interpretation that these sources are mergers.

We note however that previous studies have shown that in some star-forming galaxies at |$z\sim$| 1–3 even though their stellar morphologies appear clumpy, their kinematics are consistent with ordered circular motions (e.g. Genzel et al. 2011). We therefore stress that we are not claiming that clumpy looking sources in the rest-frame UV are all necessarily mergers, but when both their kinematics and their stellar morphologies appear irregular this is strong evidence in favour of a merger/interactions scenario.

Figure A1.

HST images (greyscale) of three of the Class II sources in our sample (ALESS 088.1, 101.1, 112.1). The contours show the distribution of the CO emission, while the HST imaging suggests potentially complex morphologies for these sources.

APPENDIX B: EXPLORING THE EFFECTS OF RESOLUTION AND SNR

In this section, we explore the effects of resolution and SNR on the inferred best-fit parameters from our modelling analysis using both simulations and observations.

B1 Simulations

We start by outlining the steps for creating simulated observations which are used to evaluate the reliability/accuracy of the best-fit model parameters. Specifically, we want to examine how well we can recover the true parameters of our model for data with different SNR and resolution. In order to create simulated visibilities for this exercise we use the observed uv-coverage of ALESS 122.1, which achieves a resolution of |$\sim 0.5$| arcsec, the highest out of all data sets used in this study. This gives us the flexibility to taper the simulated data to produce lower resolution simulations.

To create our simulations we randomly draw values for each model parameter from a wide enough range to include most of the best-fitting values we infer for the sources in our sample (e.g. 200 <|$V_{\rm max}$|< 600, 50 <|$\sigma$|< 200). After creating the model cube and Fourier transforming it to compute the model visibilities, we need to add noise. To ensure the noise properties of our simulated data are as realistic as possible we use the observed visibilities for ALESS 122.1 in a spectral window where no emission line is observed. We treat these as the noise, which is then added to our simulated visibilities.

Next, we need to ensure that we create simulations of a given SNR. For each combination of model parameters we first create a simulation with arbitrarily high intensity, and therefore SNR. We then estimate the SNR of this data set following the same approach as we used on the real observations (see Section 2.2). Having determined the SNR of this reference data set we can then scale the intensity of the model sources to produce data sets with the desired SNR. In Fig. B1, we show position velocity (PV) diagrams along the major/minor axis for one set of simulations (i.e. for different resolution and SNR) with the same model parameters.

Figure B1.

Position velocity (PV) diagrams along the major/minor axis, shown in the left and right panels, respectively, for one simulated data set. The top and bottom rows show the PVs for the 0.5/1.0 arcsec resolution simulations, respectively. The different columns correspond to simulations with different SNR, indicated at the bottom left corner. The contours are drawn at 3σ, 6σ, and 9|$\sigma$|⁠.

Finally, in Fig. B2, we show fitting results to our simulations with different resolutions (0.5 and 1.0 arcsec) and SNR (10, 15, 20). We focus on two parameters, the maximum velocity and the velocity dispersion, which are the ones that are mostly used in the main body of this work. In general, we find that we always recover the true input parameters within 3|$\sigma$|⁠. The errors on the output parameters are larger for data set with poorer resolution and lower SNR. In addition, we colour-coded the points in these graphs according to their true inclinations to highlight that sources with lower inclination typically have larger errors.

Figure B2.

Output versus input values for the two main parameters of interest, the maximum rotation velocity (⁠|$V_{\rm max}$|⁠) and the velocity dispersion (⁠|$\sigma$|⁠), colour-coded by the true inclination. The different columns corresponds to data sets with different SNR, indicated at the top left corner of each panel. The top and bottom figures correspond to data sets with different resolution, 0.5 and 1.0 arcsec, respectively.

B2 Observations: ALESS 073.1

In this section, we explore how resolution affects the inferred parameters of our model using observed data instead of simulated data. For this exercise, we utilize data for ALESS 073.1, originally presented in Lelli et al. (2021), achieving a resolution of approximately 0.1 arcsec (2017.1.01471.S).

The dynamical modelling of this source strongly suggests it is a rotating disc galaxy (Lelli et al. 2021). Initially, we modelled this source following the procedure outlined in Section 3 and used the observed data (i.e. visibilities) at the native resolution. Subsequently, we repeated the fitting process, this time tapering the data to a resolution of about 1 arcsec, similar to the average resolution of data in our main ALESS sample. In Fig. B3, the corner plot displays our fitting results for both the native (blue) and tapered (red) resolution data. As evident from this figure, the best-fit parameters inferred from these two data sets are consistent at 3|$\sigma$|⁠. Notably, constraints tighten when using higher resolution data.

Figure B3.

Corner plot for ALESS 073.1 using the data presented in Lelli et al. (2021). The blue and red contours, drawn at 2|$\sigma$| and 3|$\sigma$|⁠, show the posterior distribution using the data at the native resolution (0.1 arcsec) and after tapering to a resolution of |$\sim 1$| arcsec, respectively. The best-fitting values for all the inferred parameters from these two data sets are consistent and also in agreement with the values reported in Lelli et al. (2021).

Lastly, the best-fitting values we infer for our model parameters align with those reported in Lelli et al. (2021). Due to differences in the modelling software used in that work, an exact direct comparison is not possible. Nevertheless, we make a rough comparison. For instance, the inclination reported in Lelli et al. (2021) (25°|$\pm$| 3°) closely matches our finding 29°|$\pm$| 2°. The velocity dispersion of the tilted rings in Lelli et al. (2021) varies from |$\sim 50-60$| km s|$^{-1}$| in the inner part to |$\sim 10$| km s|$^{-1}$| in the outer part of the disc, while we find a value of |$41 \pm 3$| km s|$^{-1}$| using a uniform velocity dispersion across the disc.

APPENDIX C: FITTING DIAGNOSTICS

In this section, we show various fitting diagnostic figures for all Class I sources which are considered to be well described by a disc kinematic model. These include 0th moment maps (Fig. C1), channels maps (Fig. C2), PV diagrams along the major/minor axis (Fig. C3), and rotation curves extracted from the velocity maps (Fig. C4). All these figures (except for the rotation curves) were produced from dirty images. However, as discussed in the main text, the analysis is carried out directly in the uv-plane.

Figure C1.

Normalized dirty images of the 0th moment (intensity) of the observed, model, and residual maps for each source in our sample. The red and blue contours in the first two panels show contours from 3|$\sigma$| to 6|$\sigma$| levels in increments of 1|$\sigma$|⁠. The normalized residuals contours levels are drawn at |$\pm 2 \sigma$|⁠, |$\pm 3 \sigma$|⁠.

Figure C2.

Dirty channel maps for the data, model, and residuals, normalized by the rms noise of the cube. Contours for the data and model are drawn from |$3\sigma$| to |$6\sigma$| in increments of |$1\sigma$|⁠, while contours in the residual panels are drawn at |$\pm 2\sigma$| and |$\pm 3\sigma$|⁠.

Figure C3.

Position velocity (PV) diagrams along the major/minor axis for the 11 Class I sources in our sample which are well described by a disc kinematic model. The two left panels in each figure show the PV diagrams for the data and model. The red and yellow contours correspond to the data and model, respectively, and are drawn at 3σ, 6σ, and 9|$\sigma$|⁠. The right panel in each figure shows the residuals, where the colourbar goes from −3|$\sigma$| to 3|$\sigma$| (007.1, 075.1, and 098.1 have |$\gt 3\sigma$| residuals).

Figure C4.

Rotation curves (RC) extracted from velocity maps for the 12 Class I sources in our sample which are well described by a disc kinematic model. The RC is sampled in steps of |$\rm FWHM_{\rm maj} / (2 \sqrt{2 \ln 2})$|⁠.

The rotation curves are extracted from our observed velocity maps shown in Fig. 1. We note that these can be subject to artefacts, at this resolution, that can arise when performing the imaging with casa (we tested this by generating model visibility data sets from a smooth rotating disc model and repeated the rotation curve extraction in the same manner as the data). Nevertheless, one thing that we take away from this figure is that not all rotation curves flatten at large radii which can bias some of our estimates of the maximum rotational velocities.