Robust HI kinematics of gas-rich ultra-diffuse galaxies: hints of a weak-feedback formation scenario

We study the gas kinematics of a sample of six isolated gas-rich low surface brightness galaxies, of the class called ultra-diffuse galaxies (UDGs). These galaxies have recently been shown to be outliers from the baryonic Tully-Fisher relation (BTFR), as they rotate much slower than expected given their baryonic mass, and to have baryon fractions similar to the cosmological mean. By means of a 3D kinematic modelling fitting technique, we show that the HI in our UDGs is distributed in"thin"regularly rotating discs and we determine their rotation velocity and gas velocity dispersion. We revisit the BTFR adding galaxies from other studies. We find a previously unknown trend between the deviation from the BTFR and the disc scale length valid for dwarf galaxies with circular speeds<45 km/s, with our UDGs being at the extreme end. Based on our findings, we suggest that the high baryon fractions of our UDGs may originate due to the fact that they have experienced weak stellar feedback, likely due to their low star formation rate surface densities, and as a result they did not eject significant amounts of gas out of their discs. At the same time, we find indications that our UDGs may have higher-than-average stellar specific angular momentum, which can explain their large optical scale lengths.


INTRODUCTION
In the last five years there have been a significant number of studies aiming to detect and systematically characterize e-mail: pavel@astro.rug.nl a population of low surface brightness (LSB) galaxies with Milky Way-like effective radius, similar to those earlier reported by Sandage & Binggeli (1984) or Impey et al. (1988). Following the work by van Dokkum et al. (2015), who discovered 47 of these so-called ultra-diffuse galaxies (UDGs), different studies have found them in both high-and low-density environments (e.g., van der Burg et al. 2016;Román & Trujillo 2017a,b;Greco et al. 2018;Mancera Piña et al. 2019a;Román et al. 2019, and references therein). Among them, Leisman et al. (2017), hereafter L17, found a population of field galaxies, detected in the ALFALFA catalogue (Giovanelli et al. 2005), that meet the usual optical definition for UDGs ( µ(r, R e ) 24 mag arcsec −2 , R e 1.5 kpc 1 ), but have also large atomic gas reservoirs (≥ 10 8 M ), in contrast to the cluster population. This gas-rich field population is likely to be small in terms of total number. Prole et al. (2019b) estimated that gas-rich UDGs represent about one-fifth of the overall UDG population (cf. Mancera Piña et al. 2018;Lee et al. 2020), and Jones et al. (2018) found that they represent a small correction to the galaxy stellar and H i mass functions at all masses, with a maximum contribution to the H i mass function of 6 per cent at ∼10 9 M . Despite this, their extreme properties make them puzzling and interesting objects to study.
It is well known that resolved 21-cm observations not only reveal interactions between the extended H i galaxy discs and their environments (e.g., Yun et al. 1994;de Blok & Walter 2000;Fraternali et al. 2002;Oosterloo et al. 2007;Di Teodoro & Fraternali 2014), but also allow us to estimate their rotation velocity, angular momentum and matter distribution, key ingredients to understand their formation and evolution (e.g., de Blok 1997;Verheijen 1997;Swaters 1999;Noordermeer 2006;Posti et al. 2018b). Because of these key properties, that may reveal telltale clues about their origins, pursuing studies of UDGs from an H i perspective is potentially very interesting.
From a theoretical perspective, different ideas have been proposed to explain the puzzling nature of UDGs. Di Cintio et al. (2017) presented hydrodynamical simulations where UDGs originate in isolation due to powerful feedback-driven outflows that modify the dark matter density profile allowing the baryons to move to external orbits, increasing the scale length of the galaxies (see also Chan et al. 2018;Cardona-Barrero et al. 2020). On the other hand, Amorisco & Loeb (2016) suggested that the extended sizes of UDGs can be explained if they live in dark matter haloes with high spin parameter (see also Rong et al. 2017;Posti et al. 2018a). While currently those seem to be the most popular ideas, more mechanisms have been proposed in the literature, as we discuss in detail later.
To test these theories, isolated UDGs are very useful. Some of their properties like morphology, circular speed, baryon fraction or angular momentum, can be contrasted with expectations from the above mentioned theories in a relatively straightforward way, since they are not affected by their environments and cannot be explained by interactions with other galaxies (e.g., Venhola et al. 2017;Bennet et al. 2018). Using a combination of H i interferometric data and deep optical images for a sample of six gas-rich UDGs, Mancera Piña et al. (2019b) studied the baryonic masscircular speed plane, finding that these galaxies show a set of intriguing properties: they lie well above the canonical baryonic Tully-Fisher relation (BTFR, McGaugh et al. 2000), in a position compatible with having 'no missing baryons' within their virial radii, and with little room for dark matter inside the extent of their gaseous discs.
In this work we delve into the kinematic properties of the galaxies presented in Mancera Piña et al. (2019b), explaining in detail the methodology used to derive 3D kinematic models. Further, we expand our investigation to other properties of these LSB galaxies, and discuss possible interpretations for our results.
The rest of this manuscript is organised as follows. In Section 2 we describe our sample and give the structural parameters obtained from the optical and H i observations, and in Section 3 we provide details on our methodology and kinematic modelling. In Section 4 we estimate the scale height of the sample and we look briefly into the properties of their interstellar medium (ISM), while in Section 5 we revisit the BTFR, examine the existence of outliers and show that the deviation from the relation at low rotation velocities correlates with the galaxy scale length. A discussion on the implications of our results for proposed UDG formation mechanisms, including the addition of UDGs to the stellar specific angular momentum-mass relation, is given in Section 6. In Section 7 we present our conclusions. Throughout this work magnitudes are in the AB system, and a ΛCDM cosmology with Ω m = 0.3, Ω Λ = 0.7 and H 0 = 70 km s −1 Mpc −1 is adopted.

THE SAMPLE
The sample studied in this work and in Mancera Piña et al. (2019b) consists of six gas-rich UDGs, originally identified by L17, for which dedicated optical and interferometric observations were obtained. The observations and data reduction strategies are explained in detail in L17 and Gault et al. (submitted). We note here that the sample from Gault et al. (submitted) consists of eleven galaxies while ours consists of six. As briefly discussed in Mancera Piña et al. (2019b), we selected the galaxies that were more suitable in terms of data-quality for our kinematic modelling (see below). Figures 1-6 present our data and 3D kinematic modelling, while Table 1 gives the main properties of our galaxy sample.

Optical data
Given the LSB nature of our galaxies, the imaging of widefield public surveys like the Sloan Digital Sky Survey (SDSS, e.g., York et al. 2000) is not deep enough to provide accurate photometric parameters on an individual basis. Because of this, the six galaxies were observed using the One Degree Imager (Harbeck et al. 2014) of the 3.5-m WIYN telescope at the Kitt Peak National Observatory. The g and r bands were used, with a total exposure time of 45 min per filter. The optical image production is described in detail in Gault et al. (submitted). Panel a) in Figures 1-6 shows the r−band optical images of our sample.
As introduced in Mancera Piña et al. (2019b) and shown thoroughly in Gault et al. (submitted), aperture photometry is performed on these images to obtain total magnitudes (and colours) and surface brightness profiles. The central surface brightness and disc scale length (R d ) are obtained 21°20'00"19'30" 00"   . Data and kinematic models for the gas-rich UDG AGC 219533. Panels and symbols as in Figure 1. The H i contours are at 1.1, 2.2 and 4.4×10 20 atoms cm −2 . In this case the data cube is more noisy than in the rest of the sample, and the last contour corresponds to S/N ≈ 5.  from a fit to the observed surface brightness profiles assuming that the light distribution follows an exponential profile, which is a good assumption for these galaxies (see also Román & Trujillo 2017b;Greco et al. 2018;Mancera Piña et al. 2019a).
To derive the stellar masses we employ the mass-tolight-colour relation given by Herrmann et al. (2016): which was specifically calibrated for dwarf irregular galaxies, whose optical morphology is similar to isolated UDGs.  . Data and kinematic models for the gas-rich UDG AGC 749290. Panels and symbols as in Figure 1. The H i contours are at 0.35, 0.7, 1.4 and 2.8×10 20 atoms cm −2 . As shown in the major-axis PV diagram (panel d)), while 3D Barolo models well the main rotating body, there is signal at around velocities of 10 km s −1 and offset 25 arcsec that cannot be reproduced by the best-fit model. The parameters for these galaxy are considered less robust than for the rest of the sample, as we discuss in Section 3.2. In practice, for each UDG we randomly sample Equation 1 using Gaussian distributions on each parameter, to account for the uncertainties in both the relation itself and the photometry. The fiducial value of each parameter is chosen as the mean of its Gaussian distribution while its standard deviation gives the corresponding uncertainty. With this, we obtain a distribution for log(M /L g ) that is then converted to stellar mass using the g−band absolute magnitude distribution. The values for the stellar mass, which we report in Table 1, are the median values of the final distributions for each galaxy, and the uncertainties the difference between these medians and the corresponding 16 th and 84 th percentiles.

H i data
We obtained resolved H i-line observations using the Karl G. Jansky Very Large Array (VLA) and the Westerbork Synthesis Radio Telescope (WSRT). All the galaxies have radio data from the VLA, while AGC 122966 and AGC 334315 have also WSRT data. Details of the data reduction are given in L17 and Gault et al. (submitted). In the case of the two galaxies with VLA and WSRT observations, we use the data with the best quality in terms of spatial resolution and signal-to-noise (S/N), which were the VLA data for AGC 334315 and the WSRT data for AGC 122966. In the rest of this paper we use the parameters derived from these data. For completeness, in Appendix A we present the WSRT data for AGC 334315 2 and the VLA data for AGC 122966, demonstrating the overall good agreement between the different data cubes. We build total H i maps of our sources using the software 3D Barolo (Di Teodoro & Fraternali 2015, see below for more details). These maps are first obtained using a mask that 3D Barolo generates after smoothing the data cubes by a given factor and then selecting those pixels above a chosen threshold in units of the rms of the smoothed cube. Upon inspection of our data, we find sensible values for the smooth-ing factor and the cut threshold around 1.2 and 3.5, respectively. The fluxes of our galaxies are measured from the data cubes using the task flux from GIPSY (van der Hulst et al. 1992;Vogelaar & Terlouw 2001). The measurements of the flux from the VLA and WSRT data cubes are fully consistent with the separate analysis by Gault et al. (submitted) and L17, and in good agreement (within ≈10%) with the values obtained from the ALFALFA single-dish observations, except for AGC 248945, for which we recover ≈30% less flux. Upon inspection of the data cube, we confirm that the emission missing in the VLA profile with respect to ALFALFA is not biased with respect to the velocity extent of the source. Panel b) in Figures 1-6 presents the total H i maps of our galaxies.
We determine the H i mass of our UDGs using the equation with D and F HI the distance and flux of each galaxy, respectively. Distances, taken directly from L17, come from the ALFALFA catalog, which uses a Hubble flow model (Masters 2005). Given the line-of-sight velocities of our sample (see Table 1), and considering that these UDGs live in the field, the possible effects of peculiar velocities are not significant and the Hubble flow distances provide a robust measurement of the 'true' distance, with an uncertainty of ± 5 Mpc.

Interpretation of velocity gradients
As can be seen in the panel c) of Figures 1-6, we observe clear velocity gradients in most of our UDGs; AGC 749290 is the exception and the kinematics of this galaxy is more uncertain, as we discuss below. These gradients are along the morphological H i position angle of the galaxies, and in the following sections we interpret them as produced by the differential rotation of a gaseous disc. Here we briefly discuss other possibilities. One may wonder if the observed velocity gradients could be generated not by rotation but by gas inflow (see Sancisi et al. 2008 for a review) or blown-out gas due to powerful stellar winds (see for instance McQuinn et al. 2019 and references therein). Such winds have been observed in starburst dwarfs traced by Hα emission, where the H i distribution may also be disturbed (e.g., Lelli et al. 2014b;McQuinn et al. 2019). There is, however, clear evidence against these scenarios in the case of our UDGs. First of all, the velocity gradients are aligned with the H i morphological position angle of the galaxies, as happens with normal rotating discs. Further, as we discuss later, our measurements of the gas velocity dispersions point to a rather undisturbed and quiet ISM. Moreover and most importantly, our galaxies have normal-to-low star formation rates (SFR ≈ 0.02-0.4 M yr −1 , see L17), which combined with their extended optical scale length leads to SFR surface densities of a factor about 10 − 20 lower than in typical dwarfs. The fact that our UDGs are gas dominated and there is only one clear velocity gradient implies that if the gradients are due to winds the whole ISM should be in the wind, requiring very high mass loading factors and SFR densities, in contradiction with the information presented above. Based on this discussion we conclude that the possibility of the observed velocity fields being produced by inflows or outflows is very unlikely. In contrast, we show in Section 3 how a rotating disc can reproduce the features observed in our data.

Baryonic mass
The baryonic masses of the galaxies are computed with the equation where the factor 1.33 accounts for the presence of helium. The mass budget of our galaxies is dominated by the gas content, with a mean gas-to-stellar mass ratio (M gas /M ) ≈ 15 (see Mancera Piña et al. 2019b for more details). This ensures that, despite possible systematics when deriving the stellar mass, the estimation of the baryonic mass is robust.
As seen in Eq. 3, we neglect any contribution from molecular gas to the baryonic mass of the galaxies; while the molecular gas mass is indeed often smaller than the stellar and atomic gas ones in dwarfs (e.g. Leroy et al. 2009;Saintonge et al. 2011;Ponomareva et al. 2018, and references therein), it may of course contribute to the total mass budget. This hypothetical baryonic mass gain, however, would place our sources further off the BTFR, only strengthening the results shown in Section 5.

Initial parameters for 3D modelling
Our interferometric observations allow us to estimate rotation velocities for the six galaxies. However, the data have low spatial resolution, with only a couple of resolution elements per galaxy side. Low-spatial resolution observations can be severely affected by beam smearing, which tends to blur the observed velocity fields, and traditional 2D approaches that fit tilted-ring models to beam-smeared velocity fields fail at recovering the correct kinematics, by underestimating the rotation and overestimating the gas velocity dispersion (e.g., Bosma 1978;Swaters 1999;Di Teodoro et al. 2016).
3D Barolo 3 (Di Teodoro & Fraternali 2015) is a software tool which produces 3D models of emission-line observations (e.g., Di Teodoro et al. 2016;Iorio et al. 2017;Bacchini et al. 2019). Instead of fitting the velocity field, it builds 3D tilted-ring realizations of the galaxy that are later compared with the data to find the best-fit model. Thanks to a convolution step before the model is compared with the data, 3D Barolo strongly mitigates the effect of beam smearing, so it is ideal for analysing data like ours. 3D Barolo assumes that the discs are thin; while this is not known a priori, we show in Section 4.1 that the ratio between the radial and vertical extent of our UDGs is large, confirming the validity of our approach.
Due to the small number of resolution elements we prefer to fit only two parameters with 3D Barolo: the rotation velocity and the velocity dispersion. This means that the rest of the parameters need to be determined and fixed, namely the center of the galaxy, its systemic velocity, position angle and inclination.
3D Barolo can robustly estimate the systemic velocity and the centre of the galaxies from the centre of the global H i profile and the total H i map, respectively. We use thus these estimations from 3D Barolo and we keep them fixed while fitting the rings. 3D Barolo can also estimate the position angle and inclination, but for low resolution data these estimates may not be accurate. Therefore, we decide to estimate these two parameters independently and to fix them when fitting the kinematic parameters. The position angle is chosen as the orientation that maximizes the amplitude of the position-velocity (PV) diagram along the major axis. This is done visually using the task kpvslice of the karma package (Gooch 2006). Importantly, we find that in every galaxy the kinematic and morphological (H i) position angle are nearly the same. The exception may seem to be AGC 122966, but as we discuss in Section 3.2 this is an apparent artifact due to the shape of the beam for the WSRT observations of that galaxy.
Estimating the inclination of the galaxies is of crucial importance, as correcting for it can account for a large fraction of the final rotation velocity if the galaxies are seen at low inclinations. Unfortunately, due to the LSB nature of our galaxies, their optical morphologies, often irregular and dominated by patchy regions, provide only an uncertain, if any, constraint on the inclinations (see also Starkenburg et al. 2019 for other limitation of using optical data to determine the inclination of the H i disc.). This, together with the fact that the H i is more extended and massive than the stellar component, motivated us to use the H i maps to estimate the inclinations. We do this by minimizing the residuals between the observed H i map of each galaxy and the H i map of models of the same galaxy but projected at different inclinations between 10 • and 80 • . Such models are produced using the task GALMOD from 3D Barolo, which in turn uses updated routines from the homonym GIPSY task, and takes into account the shape of the beam when generating the models. The centre, surface density and position angle for the models are the same as in the galaxy whose inclination we aim to determine. The inclination of the model that produces the lowest residual when compared with the  . Two of the 40 APOSTLE simulated dwarf galaxies used to test our methods. Top: Total H i maps with their kinematic major axes in magenta. Middle: PV diagrams along the major axis. The black and red contours represent the simulated data and our best-fit model, respectively; grey dashed contours show negative values. The yellow points show the recovered rotation velocities. Bottom: True rotation curves (grey lines) derived directly from the simulations and recovered rotation velocities (red crosses). data (lowest absolute difference between the total H i maps) is chosen as the fiducial inclination.

Testing on simulated galaxies
We test our method to recover the galaxy initial geometrical parameters using a sample of gas-rich dwarfs from the APOS-TLE cosmological hydrodynamical simulations Sawala et al. 2016). Mock H i data cubes of these galaxies, 'observed' at a resolution and S/N matching our data, are obtained with the software martini 4 (Oman et al. 2019). The simulated galaxies have H i masses of ∼10 8−9 M and rotation velocities around ≈ 20−60 km s −1 . We initially work with four simulated galaxies with similar mass and velocity as our sample, but we project them at different random position angles and inclinations, allowing us in practice to test our methods in 40 different mock data cubes. Figure 7 shows two examples of such simulated galaxies: their H i maps, PV diagrams, and rotation curves.
We treat these mock data in exactly the same way as our UDGs data, using the method described above to derive the position angle and inclination. Figure 8 shows the 4 Version 1.0.2, http://github.com/kyleaoman/martini. true and recovered geometrical parameters. We find that we can consistently recover the position angle of the simulated galaxies and, once this is fixed, the inclination. The mean of the absolute difference between truth and recovered position angles and inclination angles is 8 • and 5 • , respectively. We adopt these values as the uncertainties for these parameters. Note that our method recovers the inclination only for galaxies with inclinations 25 • ; below this all the models look very similar and our method systematically underestimates the inclination by about 10 • . For higher inclinations, in two out of 40 cases we underestimate the inclination by about 15 • . Given the low incidence of this (5 per cent of the times) we do not expect to underestimate the inclination of our real galaxy sample; in any case this would lower the circular velocities of our sample that we report, which would not affect the nature of our results below. The bottom panels of Figure 7 show that our derived rotation velocities well represent the underlying rotation curves after having estimated the position and inclination angles as described above, and then used 3D Barolo.

Running 3D Barolo on the individual systems
After testing our methods we proceed to run 3D Barolo on all our UDGs. As discussed before, we leave the rotation and dispersion as free parameters, and we fix the position angle and inclination; the values of these parameters are given in Table 1. As expected, 3D Barolo is able to estimate the centre and systemic velocity of the sources with good accuracy, as we could verify by visually inspecting the velocity fields and PV diagrams. It is worth mentioning that the systemic velocities in the kinematic fits agree well with the values determined from the ALFALFA global profiles (the mean difference is about 5 km s −1 ). The kinematic centres and systemic velocities used in the models can be found in Table 1. The noise in the data cubes and peak column densities are provided in Gault et al. (submitted). 3D Barolo also applies a correction for instrumental broadening, controlled by the parameter LINEAR, which depends on the spectral smoothing of the data. In our case, we use LINEAR = 0.42 for the VLA data, and LINEAR = 0.85 for the (Hanning-smoothed) WSRT data.
The separation between rings is given by the parameter RADSEP in 3D Barolo, and the value is chosen taking into account the beam orientation and extension for each galaxy. Below we provide information on the value of RADSEP used for each galaxy, as well as some individual comments.
• AGC 114905: The size of the beam is 14.64"×13.31", with a North-West orientation of −5 • . The galaxy position angle is 85 • . The component of the beam projected along the kinematic major axis has a size of approximately the size of the beam minor axis. Given the extension of the galaxy, we use two rings, with RADSEP = 14.5".
• AGC 122966: The size of the beam is 33.16"×18.70", oriented at 15 • . Given the orientation of the galaxy, the component of the beam along its kinematic major axis is ≈ 20.5". Considering the extension of the galaxy (≈ 33") we oversample by a factor ≈ 1.2, using RADSEP = 16.5" and allowing us to have two points. This galaxy gives the impression of having perpendicular morphological and kinematic major axes, but this is an apparent effect due to the elongated shape of the WSRT beam, as it can be seen in Appendix A with the less elongated beam of the VLA data.
• AGC 219533: Beam size of 14.93"×13.62", with orientation at 25.5 • . The projected beam radius along the kinematic major axis of the galaxy is 13.6". We use two independent resolution elements with RADSEP = 14".
• AGC 248945: Beam size of 18.10"×14.11", with a position angle of 57.3 • . The projected beam radius along the kinematic major axis of the galaxy is 14.7". Given the extension of the H i emission we adopt RADSEP = 14.5", ending up with two resolution elements.
• AGC 334315: The galaxy has a beam size of 15.83"×13.94", oriented at −65 • . Along the major axis of the galaxy, the projected size of the beam is ≈ 14". We use two independent resolution elements with RADSEP = 16".
• AGC 749290: Beam size of 21.63"×17.88", oriented at −61 • . The projected radius of the beam along the major axis of the galaxy is ≈ 21.4". Given the extension of the galaxy we oversample by a factor 1.7 to get two resolution elements per galaxy side, with RADSEP = 12", although the resulting two rings are not independent. Because of this, the kinematic parameters of this galaxy are less certain than for the rest of our sample, and we plot the galaxy as an empty symbol when using its kinematic parameters. Nevertheless, we note that the specific values of its circular speed and velocity dispersion are similar to the values of the rest of the sample.

Kinematic models
For all our galaxies the kinematic fits converge and 3D Barolo finds models which are in good agreement with the data. Figures 1-6 show our kinematic models in panels c) to f ): observed and modelled velocity field, and observed and modelled (1 pixel width) PV diagrams.
The PV diagrams and rotation velocities suggest that we are tracing the flat part of the rotation curve, as the two points of the rotation curves are consistent with each other. This may be a possible source of confusion since some PV diagrams, at first-sight, may look like solid-body rotation. However, this is an effect of the beam smearing, and 3D Barolo is able to recover the intrinsic rotation velocities (see for instance figure 7 and 8 in Di Teodoro & Fraternali 2015), although for AGC 749292 this is not possible to establish unambiguously as we oversample the data by a factor 1.7. Moreover, standard rotation curves of simulated dwarf galaxies are expected to reach the flat part well inside our typical values of R out (e.g., Oman et al. 2015). Observed rotation curves do not keep rising after 2-3R d either (e.g., Swaters 1999), which is again inside our values of R out . Yet, higher-resolution and higher-sensitivity observations would be desirable to further confirm this, as well as to trace the inner rising part which we cannot observe at the current resolution.
Taking this into account, and the fact that the inclination is the main driver of uncertainties in the rotation velocity, we estimate the circular speeds and their uncertainties reported in Table 1 as follows: (i) For each galaxy, we run 3D Barolo two more times, but instead of using our fiducial inclination i, we use i−5 and i+5. This means that each ring of a galaxy has three associated velocities, obtained with i, i − 5 and i + 5. 3D Barolo is able to correct for pressure-supported motions, with the so-called asymmetric drift correction (e.g., Iorio et al. 2017), allowing the conversion from rotation velocities to circular speeds. We apply this correction, although it is found to be small, contributing at most ≈ 1 km s −1 .
(ii) For each of the above velocities (at i, i − 5 and i + 5), and for each ring, we generate random Gaussian distributions centred at the value of the velocity, and with standard deviation given by the statistical errors in the fit found by 3D Barolo. A galaxy with two resolution elements has six corresponding Gaussian distributions, three for each ring.
(iii) Finally, we add all these Gaussian distributions in a broader distribution G. For each galaxy, the circular speed (V circ ) corresponds to the 50 th percentile of G, and its lower and upper uncertainties (Table 1) correspond to the difference between that value and the 16 th and 84 th percentiles of G, respectively.

Circular speeds
Our galaxies have circular speeds between 20 and 40 km s −1 . Given their baryonic masses, their velocities are a factor 2 − 4 lower than the expectations from the BTFR (see Section 5 and Mancera Piña et al. 2019b). Our lower-thanaverage circular speeds are consistent with earlier observations by different authors that these kinds of galaxies have narrower global H i profiles than other galaxies with similar masses (L17; Janowiecki et al. 2019).
A question that may arise, given the long dynamical timescales implied by the low rotation velocities of our UDGs, is whether they are in dynamical equilibrium. The average dynamical time for our sample is 2 Gyr. The mean distance from our UDGs to their nearest neighbor, according to the Arecibo General Catalog 5 , is 1 Mpc. If we consider the case where all our galaxies interacted with their nearest neighbor, and we assume that they come from a 10 12 M environment (gas-rich UDGs inhabit low-density large-scale environments, see Janowiecki et al. 2019) with an escape speed of 200 km s −1 , the mean interaction back-time (how long ago did the interaction occur) for them is about 5 Gyr, so the galaxies should have had time to reach a stable configuration, having completed on average more than two full rotations.

Velocity dispersion
The narrowness of the (beam-smeared) PV diagrams of our galaxies, shown in Figures 1-6, suggests a rather low gas velocity dispersion for most of them. This is indeed confirmed by the best-fit models of 3D Barolo. The mean velocity dispersion for AGC 114905, AGC 219533, AGC 248945 and AGC 749290, with a channel width of ∆ v ≈ 4 km s −1 , is σ = 3 ± 2 km s −1 , which is below ∆ v . However, based on tests using artificial data cubes, we find that, for data like ours, 3D Barolo cannot recover the exact velocity dispersion if this lies below ∆ v , but it tends to find σ ≈ ∆ v . Therefore, for these data cubes we assume an upper limit of σ 4 km s −1 . For AGC 334315, with the same ∆ v ≈ 4 km s −1 , we find σ = 7±2 km s −1 , and we adopt this value. The WSRT data for AGC 122966, which is Hanning smoothed and of lower spatial resolution (∆ v ≈ 6 km s −1 ), has σ = 7 ± 2 km s −1 .
The observed gas velocity dispersions are lower than observed in typical spiral and dwarf galaxies. The upper limit in the velocity dispersion of the VLA cubes is indeed similar to the velocity dispersion of the 'cold' neutral medium of Leo T (Adams & Oosterloo 2018). For comparison, Iorio et al. (2017) in their reanalysis of the kinematics of dwarf galaxies from LITTLE THINGS (Hunter & Elmegreen 2006) found σ ∼ 9 km s −1 , similar to the σ ∼ 10 km s −1 of both the more massive spirals of Tamburro et al. (2009) and Bacchini et al. (2019) and the regularly-rotating starburst dwarfs from Lelli et al. (2014b). In the next section we explore the repercussions of these results.

Thickness of the gas disc
Given the gravitational potential and gas surface density of a galaxy, the value of its velocity dispersion can be used to estimate its gas disc scale height h (see for instance § 4.6.2 in Cimatti, Fraternali & Nipoti 2019). Since our galaxies are dominated by the gas component rather than the stellar and dark matter components, at least up to the outermost measured point (see figure 3 in Mancera Piña et al. 2019b), we can consider the simple case of a self-gravitating disc with constant circular speed in hydrostatic equilibrium (e.g., van der Kruit 1988; Marasco & Fraternali 2011). This exercise sources detected in the ALFALFA survey plus all the galaxies with optical spectroscopic detections within the ALFALFA footprint. It is compiled and maintained by Martha Haynes and Riccardo Giovanelli. only provides an indicative value for h, but it is still instructive as this measurement has not been yet carried out for UDGs. The scale height of such discs is given by the equation with σ the gas velocity dispersion, G the gravitational constant and Σ gas the gas surface density. Assuming a mean velocity dispersion constant with radius and the mean surface density of the disc 6 , we obtain a mean (median) disc scale height of h = 260 (150) pc. Note that these values may in reality be smaller, as i) we are adopting an upper limit in the velocity dispersion for most galaxies, and ii) we completely neglect the potential provided by the stars and dark matter, which, even if small, would contribute to flatten the disc. We can conclude that our galaxies do not appear to have H i discs significantly thicker than other disc galaxies. For reference, the H i discs studied in Bacchini et al. (2019) have mean values for h between 130 − 540 pc, depending on the galaxy, and the dwarfs from Banerjee et al. (2011) have h ≈ 500 pc. Note that the differences in the assumed shape of the vertical profile are not very big: for instance, the correction for using sech 2 instead of a Gaussian function (as in Bacchini et al. 2019) is less than 10% of the value of h.
It is also worth highlighting that given the H i radius (R HI ) of our sample (Gault et al. submitted), the values of h indicate that they have relatively 'thin' discs: the extension of their discs (using for reference R HI ) is on average 50 times larger than the size of the scale height. This result also confirms that our approach of using 3D Barolo (where galaxies are modelled as thin discs) is perfectly adequate.

Turbulence in the ISM
According to the Field (1965) criterion for thermal instabilities, the ISM should only exist in stable conditions in two well-defined phases. These two phases correspond to the cold (CNM) and warm neutral media (WNM), with temperatures of ∼ 70−100 K and ∼ 6000−8000 K, respectively, although in realistic conditions gas in the interfaces of both media exists at intermediate temperatures (e.g., Heiles & Troland 2003). These temperatures imply a thermal speed of 0.75−1 km s −1 for the CNM and 7 − 8 km s −1 for the WNM.
In this context, our UDGs are an intriguing case because the observed intrinsic velocity dispersions are lower than the expected thermal speed of the WNM. Assuming that indeed the galaxies lack of a significant amount of WNM, the velocity dispersion can be then attributed entirely to the thermal broadening of the CNM plus turbulence in the disc 7 . By further assuming that turbulence is driven entirely by supernova explosions, we can compute the supernova efficiency in transferring kinetic energy to the ISM (see for instance § 8.7.4 in Cimatti, Fraternali & Nipoti 2019 or § VI in Mac Low & Klessen 2004). We find that efficiencies between 2 and 5 per cent are enough to reproduce the observed low gas velocity dispersions. While these values are limited by all our uncertainties and are valid only within about one order of magnitude, they indicate that the supernova efficiency in our UDGs is likely similar to the expectations for disc galaxies from different theoretical papers like Thornton et al. (1988);Fierlinger et al. (2016)

UNDERSTANDING THE DEVIATION FROM THE BTFR
As discussed in Mancera Piña et al. (2019b), our UDG sample lies off the canonical BTFR, with circular speeds 2 − 4 times lower than galaxies with similar masses or, equivalently, with 10 − 100 times more baryonic mass than galaxies with similar circular speeds. This result holds after taking into account all the different possible systematics while deriving the circular speeds and baryonic masses.
Mancera Piña et al. (2019b) postulate that it may not be surprising that no other galaxies have been found to lie on a similar position off the BTFR as interferometric observations are usually targeted based on optical detections and the UDGs are a faint optical population. Some galaxies in the literature (e.g., Geha et al. 2006;Kirby et al. 2012;Oman et al. 2016) also appear to be outliers from the BTFR, although concerns regarding their kinematic parameters have been raised (see discussion in Oman et al. 2016). In this section we study in more detail the existence of outliers from the BTFR, using more galaxies with resolved 21-cm observations from the literature than in Mancera Piña et al. (2019b). We plot our UDGs (stars) and the different comparison samples in Figure 9, together with the best-fit line to the SPARC galaxies from Lelli et al. (2016b), extrapolated towards the low-circular speed regime, and the expected relation for galaxies with a baryon fraction equal to the cosmological mean (see Mancera Piña et al. 2019b). The reader interested in the main results of this comparison, without delving into the details, may wish to go ahead to Figures 9, 10 and Section 5.
We start by considering all the galaxies from the SPARC sample (Lelli et al. 2016a). From the 175 galaxies listed in the data base 8 , 135 have an available measurement of their asymptotically flat rotation velocity. Five out of these 135 galaxies are included in the LITTLE THINGS sample of Iorio et al. (2017), and since their analysis is more detailed and similar to ours we do not use the SPARC values for these galaxies. From the remaining 130 galaxies we select those with inclinations i ≥ 30 • and good quality flag on their rotation curve (Q = 1, 2, see Lelli et al. 2016a for details), ending up with 120 galaxies, shown in Figure 9 as cyan circles.
We consider also the LITTLE THINGS galaxies from Iorio et al. (2017), shown in Figure 9 as blue pentagons, and the SHIELD galaxies ), plotted as green octagons. Additionally, we include UGC 2162 (red hexagon), a UDG with resolved GMRT data presented in Sengupta et al. (2019), and a sample of nearly edge-on 'H i-bearing ultra-diffuse sources' (HUDs, see L17) with AL-FALFA data from He et al. (2019), shown as magenta diamonds.
While we restrict our comparison to samples with resolved H i data and the sample of UDGs from He et al. (2019), some other studies based on unresolved H i data are also worth briefly mentioning in the context of the BTFR. For example, the sample from Geha et al. (2006) shows a number of low-mass dwarf galaxies that increase the scatter of the BTFR at low velocities. In particular, for V rot 40 km s −1 , most of their dwarfs rotate too slowly for their baryonic mass (see their figure 7). Also, Guo et al. (2019) have recently used unresolved data from ALFALFA and faint SDSS imaging to suggest that 19 dwarfs they observe show similar properties as those discussed in Mancera Piña et al. (2019b): the galaxies seem to rotate too slowly for their masses and to have less dark matter than expected. The result from Guo et al. (2019) by itself, however, may be subject of concern as they used unresolved data to estimate the rotation velocities, and lack information on the inclination of the H i disc as they derive an inclination from shallow SDSS data that may not inform us on the actual orientation of the disc (see for instance Starkenburg  The outcome of including all the different samples can be seen in Figure 9. Appendix B provides comments on the most interesting individual galaxies from each of the samples discussed above. In general, Figure 9 suggests that it is likely that our UDGs are not the only outliers from the canonical BTFR at low circular speeds, although they may be the most extreme cases. In this context, we examine the deviation from the SPARC fit as a function of central surface brightness and disc scale length; in Appendix B we provide the references from which the structural parameters of the galaxies in Figure 9 are obtained. We realize that those galaxies above the SPARC fit usually have lower surface brightness than galaxies in the relation, as expected for a constant M/L (see discussion in Zwaan et al. 1995 andde Blok 1998). However, this is not true for all the galaxies, and the analysis may be significantly influenced by the different strategies employed to derive the surface brightness in the literature (e.g., if values are corrected or not for inclination, dust reddening and Galactic extinction, and if different filters were used). Instead, measuring the radius of galaxies is more straightforward, as it has been shown to be less dependent on the different optical and infrared bands used to derive it (see for instance figure  . Circular speed vs. baryonic mass plane for different galaxy samples. Gas-rich UDGs are shown as orange stars, except for AGC 749290, whose circular speed is less robust than for the rest of the sample, and it is shown as a white star. The black dotted like shows the fit to the SPARC galaxies from Lelli et al. (2016b), extrapolated towards low circular speeds. The grey solid line is the expectation for galaxies that have a baryon fraction equal to the cosmological mean. When including more galaxies from the literature in this plane, the apparent gap between the canonical BTFR and our UDGs is populated. See the text for details.  2008), with a larger fraction of LSB galaxies, reported that the scale lengths of their sample do correlate with the residuals of the BTFR, with smaller galaxies deviating towards higher velocities at fixed baryonic mass (note however that they looked at the BTFR using V max , instead of V flat like in the other two mentioned studies). Apart from the existence of these discrepancies, those works include only a few galaxies with circular speeds similar to those of our UDGs (V circ ≈ 20 − 40 km s −1 ), so it is interesting to re-consider the possible existence of correlations within the same range in velocity as our sample, using the compilation of galaxies that we have shown in this section.
Given the values of V circ for our sample, we consider galaxies with 15 km s −1 < V circ < 45 km s −1 , and we use them in Figure 10 to build the R d − ∆ M bar plane. Here R d is the stellar scale length and ∆ M bar the vertical distance of the galaxies from the BTFR, defined as the logarithmic difference between the observed baryonic mass and the value expected from the extrapolated SPARC BTFR, ∆ M bar ≡ log(M bar,obs /M bar,BTFR ). A clear trend is found: at these low circular speeds, larger (more diffuse) galaxies lie systematically above the SPARC BTFR while the more compact ones lie below; the correlation has a slope around 1.5.
Spearman tests tell us that the correlation is significant at a 99.9 per cent confidence level (p−value ≈ 10 −8 ) when all the samples are considered. This holds even if we exclude our UDGs (p−value ≈ 0.0003). The correlation is less robust but still significant at the 95 per cent level when considering exclusively the SPARC and LITTLE THINGS galaxies (p−value ≈ 0.02).
This trend is potentially of great importance because it provides evidence supporting the idea that the deviation from the BTFR at low circular speeds is driven by physical processes related to the optical size of the galaxies (which is independent of the kinematics), and that it is not only an effect produced by observational biases.
One may wonder whether it is possible to interpret the trend as a spurious relation due to a severe underestimation of the circular speed of the galaxies: if the galaxies that deviate from the SPARC BTFR have wrong measurements but actually have V circ ≈ 80 − 100 km s −1 , then they would be expected to have larger scale lengths, giving rise to the trend observed in Figure 10.
We find this unlikely for several reasons. First, as discussed in Mancera Piña et al. (2019b), a significant underestimation of the circular speeds of our sample is extremely unlikely. Further, since galaxies from several independent samples analyzed with independent techniques all seem to follow the trend in Figure 10, so the circular speeds of all the other galaxies would need to be underestimated in precisely the same way, which seems very unlikely. Finally, let us consider, ad absurdum, the following scenario. If we assume that all galaxies that deviate from the SPARC BTFR have wrong measurements, but they actually lie on it with V circ ∼ 80 km s −1 , then those galaxies should have higher surface brightness than dwarfs with V circ ∼ 20 − 40 km s −1 . So, if the trend in Figure 10 is spurious, we should also find that Symbols are as in Figure 9 and the dashed line represents no offset from the SPARC BTFR. A correlation between both parameters is observed, with larger galaxies falling systematically above the BTFR. Some samples have no reported uncertainty in R d , so we do not plot any horizontal error-bar for consistency.
galaxies which (apparently, due to wrong measurements) deviate from the BTFR have higher surface brightness than the dwarfs (V circ ∼ 20 − 40 km s −1 ) in the BTFR, which is not observed. Based on this we are led to believe that the correlation in Figure 10 is real, and it provides an extra parameter to explain the deviation from the canonical BTFR and its larger scatter at the low-velocity regime. The vertical offset from the BTFR can also be seen as a progression in the baryon fraction of the galaxies: at fixed V circ , the more baryonic mass a galaxy has, the higher its baryon fraction is. This, coupled with our results above, implies that at low circular speeds R d plays a role affecting the baryon fraction of galaxies (see Section 6.5 for more details).
Based on our discussion, we outline a possible interpretation of our results regarding the phenomenology of the BTFR. Perhaps, the SPARC BTFR holds at low circular speeds (V circ 50 km s −1 ), but the distribution of galaxies in the V circ − M bar plane may be more complex than a welldefined and tight relation as the one established at larger circular velocities, with dwarf galaxies showing baryon fractions above the one implied by the canonical BTFR but still below the cosmological limit. From our analysis, it looks very possible that the disc scale length is an important parameter regulating the deviation from the canonical BTFR. A more extreme scenario would be one where the canonical BTFR breaks down at low circular speeds, being replaced by a 2D distribution. Given the selection biases and small statistics of the samples being analysed, we cannot discern among these options, and this should be addressed with more complete and representative samples.

DISCUSSION: THE ORIGIN OF GAS-RICH UDGS
Using the kinematic information derived in the previous sections, here we discuss how our results compare with predictions from some of the main theories that have been proposed to explain the origin and properties of UDGs.
6.1 Brief comparison with NIHAO simulations: formation via feedback-driven outflows?
Di Cintio et al. (2017) studied simulated dwarf galaxies from the Numerical Investigation of a Hundred Astrophysical Objects (NIHAO) simulations (Wang et al. 2015), and found a subset of them with properties similar to observed UDGs in isolation. They found that intermittent feedback episodes associated with bursty star formation histories modify the dark and luminous matter distribution, allowing dwarf galaxies to expand, as their baryons move to more external orbits (see also e.g., Navarro, Eke & Frenk 1996;Read & Gilmore 2005;Pontzen & Governato 2012or Read et al. 2016 for further considerations). Because of this, some of their simulated dwarf galaxies become larger, entering in the classification of UDGs. Chan et al. (2018) reported similar results with the FIRE simulations (e.g., Hopkins et al. 2018). We have observational evidence suggesting that our galaxies have low velocity dispersions and thus a low turbulence in the ISM. In principle, this seems at odds with models that require stellar feedback strong enough to modify the matter distribution. A detailed comparison between our observations and this kind of simulations it is beyond the scope of this paper. Yet, it is interesting to make some brief comments on some apparent similarities and discrepancies between the simulated NIHAO UDGs and our sample. By inspecting the optical scale lengths, we see that our largest galaxies have no counterparts among the NIHAO UDGs (their largest simulated UDG has R d ≈ 2 kpc). In general, the mean values differ by a factor 2.5 (ours being larger), but strong selection effects are at play so this should be studied with a complete sample. The gas mass of our galaxies and NIHAO UDGs largely overlap, but our distribution has a sharp selection cut around M HI < 10 8.5 M .
The UDGs formation mechanism proposed by Di Cintio et al. (2017) can also be contrasted with the observational results of Mancera Piña et al. (2019b), in particular the baryon fraction of the galaxies with respect to the cosmological mean and the inner amount of dark matter. Di Cintio et al. (2017) mention that their simulated UDGs show a correlation between their optical size and baryon fraction, with their largest UDG having a baryon fraction up to 50 per cent of the cosmological value, with a mean of 20 per cent for the whole sample. Our UDGs have ≈ 100 per cent of the cosmological value. Nevertheless, as mentioned before, our galaxies have also larger scale lengths than the NIHAO UDGs, so one may wonder whether their higher baryon fraction is just a consequence of this. Extending our sample to include UDGs with smaller R e may shed light on the connection between them and the simulated NIHAO UDGs. The inner dark matter content is a major discrepancy between our observations and the UDGs that the NIHAO simulation produces: our galaxies show very low dark matter fractions within their discs (measured within ∼ 2 R e on average), while Jiang 2019b found that the NIHAO UDGs are centrally dark-matter dominated (measuring the dark-matter content within 1 R e ). Related to this, Di Cintio et al. (2017) reported that their UDGs have dark matter concentration parameters typical of galaxies with similar halo masses. This does not seem to be the case in our sample: preliminary attempts of rotation curve decomposition of our UDGs show that if they inhabit 'normal' NFW (Navarro, Frenk & White 1996) dark matter haloes (i.e. with a halo mass typical of galaxies with their stellar mass), their concentration parameters need to be extremely low (see also Sengupta et al. 2019), far off expectations of canonical concentration-halo mass relations (e.g., Ludlow et al. 2014). This should be investigated further with data at higher spatial resolution, but it opens the exciting possibility of providing clues on the nature of dark matter itself (e.g. Yang et al. 2020). Producing artificial data cubes of the NIHAO UDGs to explore their H i kinematic parameters (like their position with respect to the BTFR), as well as obtaining SFR histories for our sample would also allow an interesting and conclusive comparison, although the latter has been proved to be challenging even for closer UDGs (e.g., Ruiz-Lara et al. 2018;Martín-Navarro et al. 2019). Stellar kinematics seems to be a promising tool as well (Cardona-Barrero et al. 2020).

High angular momentum
Angular momentum is a fundamental quantity to understand the origin of high surface brightness and LSB galaxies (e.g., Dalcanton et al. 1997b;Di Cintio et al. 2019), and it is usually studied via the so-called spin or λ-parameter for dark matter haloes (e.g., Mo et al. 1988;Dutton & van den Bosch 2012;Rodríguez-Puebla et al. 2016;Posti et al. 2018a) or with the specific angular momentum-mass relation for the stellar component (e.g., Fall 1983;Romanowsky & Fall 2012;Fall & Romanowsky 2018;Posti et al. 2018b).
One of the main ideas to explain the large scale lengths and faint luminosities of UDGs is that they are dwarfs living in high spin (high-λ) dark matter haloes. This is supported by some semi-analytical models and hydrodynamical simulations, where the size of a galaxy is set by its λ, that seem to reproduce different observational properties of the (cluster ) UDG population like abundance, colours and sizes (e.g., Amorisco & Loeb 2016;Rong et al. 2017;Liao et al. 2019). Some other simulations, however, do not find anything atypical in the angular momentum content of UDGs (e.g., Di Cintio et al. 2017;Tremmel et al. 2019).
In this Section we investigate the angular momentum content of our sample, looking separately at the specific angular momentum of gas and stars.

H i specific angular momentum
Based on H i observations, L17 and  suggested that gas-rich UDGs could indeed have higher λ-parameter than other galaxies of similar mass. However, these results are derived from the relation given by Hernandez et al. (2007) to estimate λ from observations, which is highly assumption-dependent, as discussed in detail in Dutton & van den Bosch (2012). In particular, our galaxies do not follow the same BTFR nor seem to have the same disc mass fraction as the galaxies used by Hernandez et al. (2007) to calibrate their relation. Therefore, we decided not to estimate the λ-parameter in that way, and we emphasize that the calibration of Hernandez et al. (2007) should be used with caution, as also mentioned in L17,  and Sengupta et al. (2019).
Unfortunately, we cannot robustly estimate the angular momentum of the gas component of our galaxies as we lack the resolution needed to determine the shape of the surface density profile (e.g., Obreschkow & Glazebrook 2014; Kurapati et al. 2018), so we can only make qualitative statements (see also Sengupta et al. 2019). In this context, the fact that our gas-rich UDGs lie on the H i mass-size relation (Gault et al. submitted) is useful, as it tells us that their H i discs have normal sizes for their H i mass. Additionally, we have shown that the gas rotates at velocities much lower than other galaxies with the same H i mass. Together, these results suggest that our UDGs have low-to-normal gas specific angular momenta compared with galaxies of similar H i mass.

Stellar specific angular momentum
As mentioned before, the stellar specific angular momentum-mass relation (sometimes called 'Fall' relation, see Fall 1983;Romanowsky & Fall 2012) is often used as a more direct way to study the angular momentum of galaxies. To compute this relation, high-resolution (stellar) rotation curves and stellar surface density profiles are needed. However, it is common (see discussion in Romanowsky & Fall 2012;Rizzo et al. 2018) to adopt the approximation where j is the stellar specific angular momentum, R d the optical disc scale length and V rot, the stellar rotation velocity. This approximation has been proved to work very well, and it is valid for galaxies with exponential light profiles and flat rotation curves. Thus, to use this simplified version to compute j for our sample, we have to assume flat rotation curves, which seems at least tentatively supported by our data, and exponential profiles, that describe reasonably well the stellar profile of our galaxies (see Gault et al. submitted).
As the rotation velocity of the stars is needed, the next step is assuming that their rotation can be inferred from the circular velocity of the galaxies, by means of the stellar asymmetric drift correction (V AD, ), via the equation To compute V AD, , we follow the approach described in Posti et al. (2018b), using the equation: with σ 0,z the vertical stellar velocity dispersion 9 . For simplicity we use only the outermost point of the rotation curve, so log(j / kpc km s 1 ) Disc galaxies (Posti + 18b) HI rich UDGs Mean HI rich UDGs Figure 11. Stellar specific angular momentum-mass relation. Orange stars show our UDGs (AGC 749290 is in white as in Figure 9) and the red square their mean position. Blue circles show the sample analyzed by Posti et al. (2018b), while the black dashed line and the pink band are their best-fit relation and its 1σ scatter, respectively.
effectively R = R out . As discussed by Posti et al. (2018b), the value of σ 0,z depends on the central surface brightness (Martinsson et al. 2013), and for our sample it is about 5 km s −1 .
The results of this exercise are shown in Figure 11, where we plot our galaxies in the M vs. j plane. We compare our UDGs with the galaxies studied in Posti et al. (2018b), showing also the best-fit relation (dashed line) and scatter (pink band) that they obtain. We stress here that the assumptions that we have made in our analysis are the same as in Posti et al. (2018b), making the comparison in Figure 11 as fair as our data allow.
Three of our galaxies (AGC 114905, AGC 248945 and AGC 749290) lie within the 1σ scatter of the relation. AGC 219533 and AGC 334315 have a j about 3-4 times larger than the best-fit line, although the observational scatter is relatively large at those values of M * . The outlier with the highest j is AGC 122966, which has both the largest optical disc scale length and highest rotation velocity of our sample, resulting in a j about 9 times larger than expected. Note, however, that the Fall relation is not well-constrained at M < 10 8 M . A caveat to bear in mind regarding AGC 122966 is that it has the lowest surface brightness of our sample (see Table 1 in Mancera Piña et al. 2019b), so its scale length is relatively more uncertain than for the other UDGs. The mean (median) ratio between the measured and expected j of our galaxies is 3.3 (2.5). These numbers are of course dependent on the several assumptions we have made, and need to be confirmed with better data and more accurate calculations (for instance by obtaining high-resolution stellar rotation curves to formally compute j instead of using the approximation of Eq. 5 with Eq. 7). Yet, our simple analysis indicates that, as a population, our gas-rich UDGs may have a j a factor ≈ 3 higher than the expectations for dwarf irregular galaxies with similar M , as shown graphically with in Figure 11, where the mean value for our sample is indicated in red. This larger j may help explaining why UDGs have more extended optical disc scale length/effective radius than other dwarf irregulars.
The H i component of our galaxies is both more massive and more extended, so one may speculate that its specific angular momentum is likely to be more representative of the spin of the dark matter halo. If this is the case, our UDGs would be galaxies that inhabit dark haloes with normal-tolow λ but with higher-than-average j , meaning that they would be galaxies with a higher-than-average 'retained' fraction of angular momentum ( j / j halo ), as suggested by Posti et al. (2018a).

'Failed' Milky Way galaxies
Another mechanism proposed to explain the nature of UDGs is that they could be 'failed' Milky Way-like galaxies, with massive dark matter haloes that for different reasons (e.g., strong supernova feedback or gas stripping) failed at converting their gas into stars (van Dokkum et al. 2016). This idea is mainly motivated by the high velocity dispersions of globular cluster observed around a few UDGs (e.g., van Dokkum et al. 2016;Toloba et al. 2018) and by the large effective radius of UDGs (but see Chamba, Trujillo & Knapen 2020). However, several other studies, both observational and theoretical, show that most UDGs should reside in dwarf galaxysized dark matter haloes (e.g., Beasley & Trujillo 2016;Rong et al. 2017;Amorisco 2018;Ruiz-Lara et al. 2018;Kovacs et al. 2019;Prole et al. 2019a;Tremmel et al. 2019).
A dark matter NFW halo of virial mass 10 12 M , following a standard concentration-mass relation, is expected to have a maximum circular speed around 170 km s −1 , with a value of about 120 km s −1 at our typical R out , much larger than the velocities observed in our sample. Thus, our data can safely exclude the 'failed' Milky Way scenario as an origin for our gas-rich UDGs.
Actually, the fact that our galaxies have a baryon fraction close to the cosmological mean (see Figure 9) suggests that they may reside in ∼ 10 10 M dark mater haloes 10 , although this should be confirmed with a detailed mass decomposition and will be the subject of future work.

Ancient tidal dwarfs
Tidal dwarf galaxies (TDGs) are self-gravitating systems formed from the collapsed tidal debris of interacting galaxies. Because of this, one expects to find them inhabiting the chaotic environments near their progenitors, even after some Gyr (e.g., Duc et al. 2014;Lelli et al. 2015). In fact, TDGs in the nearby universe are usually found within 15 R e,p of their progenitors (R e,p being the effective radius of the progenitor, Kaviraj et al. 2012). As they form from pre-enriched material they are expected to show high metallicities (e.g., Duc & Mirabel 1998, but see also Hunter et al. 2000), and due to their weak gravitational potential they should be free of dark matter (see Hunter et al. 2000;Braine et al. 2001;Lelli et al. 2015 and references therein).
From an observational point of view, TDGs present some properties similar to those in our sample of UDGs.
For instance, it has been argued that they are comparable in terms of effective radius and surface brightness (Duc et al. 2014). Perhaps more intriguingly, they also share dynamical properties: they lie off the BTFR and show dark matter fractions close to zero within their H i discs (Lelli et al. 2015). Given all this, it is pertinent to ask whether or not our gasrich UDGs may be TDGs. In this Section we explain why we consider this scenario unlikely.
The strongest evidence against a tidal origin is the environment of our UDGs. They cannot be young tidal dwarfs given their totally different environments: they are isolated with the mean distances to their nearest, second-nearest and third-nearest neighbors being 1 Mpc, 1.4 Mpc and 1.7 Mpc, respectively. This is completely different from the expected progenitor-TDG separation of 15 R e,p , which would typically be around 100 kpc. It would be required that all our galaxies formed as TDGs some Gyr ago, and the separation between all of them and the other interacting galaxy increased up to at least 1 Mpc, which seems unlikely (see also the text in Section 3.3.1).
One may also argue that we do not find remaining signs of tidal interactions around our galaxies, but this could be just because detecting such interactions is hard (e.g., Holwerda et al. 2011;Müller et al. 2019). Perhaps more importantly, from cosmological simulations one would expect TDGs to have smaller sizes than typical dwarfs (see the discussion in Haslbauer et al. 2019), while UDGs are exactly the opposite, a population of galaxies with much larger optical radii than other dwarfs; Duc et al. (2014), however, argue that some old TDGs are larger than normal dwarfs.
Currently, our interpretation is that our galaxies live in ∼ 10 10 M dark mater haloes and the dark-matter deficiency is restricted to the extent of the H i disc. However, our current data cannot unambiguously distinguish between this scenario and one where our UDGs have very little, if any, dark matter in their whole extension, as expected in TDGs. The mass decomposition in our galaxies would then conclusively tell us if in fact there is room to accommodate (low-concentration) 10 10 M dark mater haloes or if their lack of dark matter is analogous to that in TDGs.

Weak feedback producing little mass losses
Different cosmological hydrodynamical simulations predict large mass loading factors in dwarf galaxies (e.g., Vogelsberger et al. 2013;Ford et al. 2014), although other theoretical studies show that big mass losses do not necessarily take place (see for instance Strickland & Stevens 2000;Dalcanton 2007;Emerick et al. 2018;Romano et al. 2019). Supporting this latter scenario, there is recent observational evidence (McQuinn et al. 2019, see also Lelli et al. 2014a) that the mass loading factors in dwarf galaxies are indeed relatively small, as often the outflows do not reach the virial radii of the galaxies and are kept inside their haloes, available for the regular baryon cycle. These results suggest that some dwarfs may have baryon fractions larger than expected, qualitatively in agreement with our gas-rich UDGs with 'no missing baryons'.
In this work, along with Mancera Piña et al. (2019b), we suggest that a scenario where feedback processes in our UDGs have been relatively weak and inefficient in ejecting gas out of their virial radii could explain their quiescent ISM (inferred from the velocity dispersion) and high baryon fractions (as derived from the BTFR). In Figure 10, we found a significant trend for low-mass galaxies (15 km s −1 < V circ < 45 km s −1 ) where they deviate more from the canonical BTFR, towards higher masses (and thus have a larger baryon fraction), if they have large disc scale lengths. As most of these galaxies have normal or low SFRs with respect to their stellar masses (e.g., Teich et al. 2016; L17, but more robust measurements would be desirable), their large disc scale lengths then imply that they have lowerthan-average (about an order of magnitude) SFR surface densities, and we speculate that this affects their capability to drive outflows powerful enough to eject baryons out of their virial radius, and thus allowing the galaxies to retain a higher-than-average baryon fraction. Recently, using high-resolution hydrodynamical simulations, Romano et al. (2019) have shown, for an ultra faint dwarf galaxy, that this indeed may be the case: the gas removal of a galaxy depends on the distribution of supernovae explosions and thus in the distribution of star formation. Those authors find that the more evenly spread or 'diluted' the distribution of OB associations is, the higher is the gas fraction that the galaxy keeps, in agreement with our interpretation.

CONCLUSIONS
In this work we present the 3D kinematic models of six dwarf gas-rich ultra-diffuse galaxies (UDGs). By analysing VLA and WSRT 21-cm observations with the software 3D Barolo, we derive reliable measurements of the circular speed and gas velocity dispersion of our sample galaxies. Our models have been used by Mancera Piña et al. (2019b) to show that the galaxies lie significantly above the baryonic Tully-Fisher relation.
Our main findings can be summarised as follows: • We have shown that the kinematic models are robust (Figures 1-8) and our galaxies have circular speeds of 20 − 40 km s −1 .
• Our UDGs exhibit low gas velocity dispersions, lower than observed in most dwarf irregular galaxies before. Their H i layers have the typical thicknesses observed in other dwarfs and disc galaxies, and gas turbulence appears to be fed by supernovae with efficiencies of just a few percent.
• We have reviewed the canonical baryonic Tully-Fisher relation, showing that these gas-rich UDGs are likely not the only outliers, although they may represent the most extreme cases ( Figure 9).
• At circular speeds below ≈ 45 km s −1 the vertical deviation from the canonical baryonic Tully-Fisher relation correlates with the disc scale length of the galaxies (Figure 10).
• The low velocity dispersions observed in our sample seem at odds with models where UDGs originate from feedback-driven outflows. Our galaxies tend to have larger scale lengths than the simulated NIHAO UDGs and to have higher baryon fractions, but they share some other structural properties. The most important discrepancy is that, unlike NIHAO simulated UDGs, our galaxies have little dark matter in the inner regions (within ∼ 2R e ).
• We find indications that the gas specific angular momentum of our sample is similar or slightly lower than that in other galaxies of similar H i mass. However, the specific angular momentum of the stellar component may be (a factor ≈ 3) higher-than-average at given M ( Figure 11). This can help to explain the large optical scale length of UDGs.
• The measured low circular speeds rule out the possible origin as 'failed' Milky Way galaxies for our UDGs.
• Fully testing the idea that all our six galaxies are old tidal dwarf galaxies is not possible with our observations, but their isolation seems to go against this possibility.
• To explain the high baryon fractions and low turbulence in the discs, Mancera Piña et al. (2019b) have suggested that these galaxies experienced weak feedback, allowing them to retain all of their baryons. Here we have shown that this idea is consistent with our findings: the larger the optical disc scale length of dwarf galaxies is, the more they depart from the canonical baryonic Tully-Fisher relation towards higher baryonic masses. Their extended optical sizes coupled with normal star formation rates result in very low star formation rate surface densities, impacting their capability to lose mass via outflows.

ACKNOWLEDGEMENTS
We thank an anonymous referee for the comments provided on our manuscript. We would like to thank Arianna di Cintio, Vladimir Avila-Reese, Anastasia Ponomareva, Nushkia Chamba and Francesco Calura for fruitful discussions regarding our data and results. We also acknowledge the work done by Michael Battipaglia, Elizabeth McAllan and Hannah J. Pagel on acquisition and reduction of data used in this work. P.E.M.P. thanks Giuliano Iorio  This work uses observations from the VLA and the WSRT. The VLA is operated by the National Radio Astronomy Observatory, a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. The WSRT is operated by the ASTRON (Netherlands Institute for Radio Astronomy) with support from NWO. Our work is also based in part on observations at Kitt Peak National Observatory, National Optical Astronomy Observatory, which is operated by the Association of Universities for Research in Astronomy (AURA) under a cooperative agreement with the National Science Foundation.
We have made extensive use of SIMBAD and ADS services, as well as of the Python packages NumPy (Oliphant 2006), Matplotlib (Hunter 2007), SciPy (Jones et al. 2001) and Astropy (Astropy Collaboration et al. 2013), and the tool TOPCAT (Taylor 2005), for which we are thankful.
APPENDIX A: WSRT VS. VLA DATA As mentioned in the main text, two galaxies in our sample, AGC 122966 and AGC 334315, have both VLA and WSRT data. Given the spatial resolution and S/N of the data, we adopt as our 'fiducial' data the WSRT observations for AGC 122966 and the VLA observations for AGC 334315.
Here we delve into the reasons in which we based these choices, and we show the complementary data: the VLA data for AGC 122966 as well as the WSRT data for AGC 334315. We also make the comparison between the kinematic and geometrical parameters of the different data.
The VLA data cube of AGC 122966 has a lower S/N than the WSRT one. Moreover, the emission is less extended, and we have to oversample by a factor 2 in order to fit the galaxy with 3D Barolo. Figure A1 shows the stellar image, H i map, and observed and modelled velocity fields and PV diagrams. From inspecting the major-axis PV, we can see that the emission connecting approaching and receding sides is not detected. More importantly, it looks likely that there is emission missing at the end of the approaching and residing sides too, which would significantly affect the recovered value of the rotation velocity. In fact, while for the WSRT cube 3D Barolo finds a projected velocity of 26 km s −1 , it finds 16 km s −1 for the VLA cube. Given all this, we decided to use the WSRT data for this galaxy (Figure 2).
Despite the problems deriving the kinematics, we can still use the VLA data to independently estimate the inclination of the galaxy as described in the main text. Our method finds an inclination of 44 • ± 5 • using the WSRT H i map and of 40 • ± 5 • with the VLA data, meaning that the inclination estimates from the different data cubes are consistent with each other. The VLA data of AGC 122966 is not only useful to confirm our inclination measurements, but also because the morphology of the galaxy can be better appreciated without the elongated beam of the WSRT observations ( Figure 2). Finally, the models for the different cubes have the same systemic velocity and physical center.
On the other hand, for AGC 334315 the VLA data ( Figure 5), have better quality (although a factor 2 less extended) than the WSRT data, shown in Figure A2: the spatial and spectral resolution are better, and the beam is more circular. With this, one can more clearly appreciate the H i structure as well as the intrinsic shape and velocity field of the galaxy. Apart from this, the parameters found for the two different data cubes are perfectly compatible with each other: we find an inclination of 45 • ± 5 • for the VLA map and of 52 • ± 5 • with the WSRT one, in agreement within the reported uncertainties. The final circular speed is also the same within the uncertainties: 25 km s −1 at 8.5 kpc for the VLA data, and 26 km s −1 at 18 kpc for the WSRT observations. The models for the different cubes share also the same systemic velocity and physical center.

APPENDIX B: COMPARISON SAMPLES
B1 The BTFR Figure 9 shows the circular speed-baryonic mass plane, with galaxies scattering around the BTFR. The main result from the figure is that the region between our extreme UDGs and galaxies following the extrapolation of the SPARC BTFR at low masses, is populated by other dwarf galaxies. When discussing specific galaxies which are particularly interesting given their position with respect to the SPARC BTFR, we provide their coordinates on the V circ − M bar plane.

B1.1 SPARC
Within this sample (cyan circles in Figure 9), two galaxies deviate significantly from the SPARC BTFR: UGC 7125 (65 km s −1 , 10 9.86 M ) and UGC 9992 (34 km s −1 , 10 8.77 M ). UGC 7125 has an inclination of 90 • and Q = 1, so its position seems robust, even though, as Lelli et al. (2016b) mention, the galaxy has a relatively large uncertainty on its mass due to a high distance-correction for Virgocentric infall.
UGC 9992 has an inclination of 30 • ± 10 • and Q = 2. Its PV diagram (Swaters 1999) is regular and traces the flat part of the rotation curve. To bring this galaxy back to the BTFR an inclination of ∼ 17 • would be needed. For this hypothetical inclination to be correct, an error of 13 • in the measured inclination is needed, a bit larger than the quoted uncertainty of 10 • but within 1.5σ. It is worth mentioning that two other galaxies with i = 30 • ± 10 • (F571-V1 and UGC 7261) lie much closer to the BTFR, as well as, for instance, the galaxies UGC 11914, UGC 6930 and UGC 10310, with i = 31 • ± 5 • , 32 • ± 5 • and 34 • ± 6 • , respectively.

B1.2 LITTLE THINGS
Most LITTLE THINGS galaxies lie around the extrapolation of the SPARC BTFR, even when some of their rotation curves have not clearly reached their flat part (see . DDO 50 (also known as UGC 4305 or Ho II) is an outlier, in a position very close to our UDGs (39 km s −1 , 10 8.95 M ). Its rotation curve has clearly reached the flat part but the inclination of the galaxy is relatively low, 30 • . Different values have been proposed, ranging from 30 • to 50 • (see Oman et al. 2016;Iorio et al. 2017 and references therein), so the exact value of its circular speed remains somewhat uncertain. However, an inclination of ∼ 18 • is needed for the galaxy to lie directly on the canonical BTFR, which is outside the wide range proposed in the literature.
IC 1613 (also known as DDO 8, not shown in Figure 9), is another well known candidate to be an outlier from the BTFR with a baryonic mass of 10 7.9 M and rotation velocity of about 20 km s −1 (Oh et al. 2015;Oman et al. 2016). It is part of the LITTLE THINGS galaxies studied in Oh et al. (2015), but was excluded in the work of Iorio et al. (2017). While the galaxy is potentially interesting, there may be issues with its inclination and equilibrium state, as mentioned in Read et al. (2017), Because of this, as in Iorio et al. (2017), we do not consider this galaxy in Figure 9.
In the sample from Iorio et al. (2017), DDO 101 (59 km s −1 , 10 7.94 M ) stands out but for being an outlier that rotates too fast for its baryonic mass. However, Read et al. (2017) have demonstrated that this discrepancy is largely mitigated if the galaxy is farther away than the distance used in Iorio et al. (2017), so it is possible that the uncertainties in distance (and thus in baryonic mass) have been underestimated for this galaxy.

B1.3 SHIELD
Due to the very low spatial resolution data the analysis of this sample was very challenging, as discussed in McNichols et al. (2016). Five out of twelve galaxies in the sample have rotation velocities estimated from fitting tilted-ring models to the observed low-spatial resolution velocity fields, and at least three of the resulting rotation curves have no indication of a flattening. The remaining velocities were derived from the visual inspection of different PV slices (see McNichols et al. 2016 for details). In addition to this, the asymmetric drift correction was not applied to the galaxies, and the inclination comes from optical data rather than from the H i itself. Despite these differences with other analyses, the twelve galaxies lie around the extrapolation of the canonical BTFR.

B1.4 Edge-on HUDs
He et al. (2019) selected a sample of nearly edge-on HUDs that were not originally selected in the catalogue of L17 due to selection effects (given their high inclination the galaxies have apparent higher surface brightness). As the optical morphology indicates that the galaxies are edge-on, corrections for inclination are negligible. Their optical images, shown in He et al. (2019), suggest that their stellar structure is not exactly the same as in our sample, as theirs look more regular and thinner, but it is hard to unequivocally judge this as our sample is not edge-on. However, the stellar component of our sample is both more extended (by a factor about 2.5) and of lower surface brightness.
Most of these galaxies (magenta diamonds in Figure 9) show velocities nearly compatible with the SPARC BTFR. The clear outlier is AGC 202262, that shows the narrowest velocity width, with a rotation of about 30 km s −1 . An inclination of ∼ 40 • , totally incompatible with the optical image (but see the caveats on the optical-gas misalignment reported in Starkenburg et al. 2019 or Gault et al. submitted), would be required to put the galaxy on the canonical BTFR.

B1.5 UGC 2162
UGC 2162 is an UDG which has been studied using 3D Barolo in resolved GMRT observations (Sengupta et al. 2019). Two caveats exist regarding this galaxy. First, Sengupta et al. (2019) mention that it has large uncertainties in its baryonic mass (probably because of the uncertainty associated with its distance; the galaxy is much closer than those in our sample), and second, its rotation curve does not show signs of flattening (the emission is less extended than in our data), although it is presumably close to reaching the flat part (see figure 2 in Oman et al. 2015).
While these caveats should be taken into account, the position of UGC 2162 11 in Figure 9 (red hexagon) seems in line with our results. Note also that UGC 2162 is less 'extreme' than our UDGs: it is less massive, smaller and slightly brighter. Assuming that the amplitude of the rotation curve does not increase significantly beyond the outermost measured radius, the galaxy would need an inclination of ≈ 39 • to be on the SPARC BTFR, far from the measured inclination of 55 • in Sengupta et al. (2019).
11 UGC 2162 has no reported uncertainties in rotation velocity nor H i mass by Sengupta et al. (2019), so we assume an uncertainty of 7 km s −1 for the velocity (the spectral resolution of the data, which is also larger than the difference in the velocity obtained from the best-fit model of 3D Barolo and the global H i profile) and a typical value of 20 per cent in the H i mass.

B2 Disc scale length and central surface brightness values
In Figure 10 we study the deviation from the BTFR as a function of the disc scale length; the central surface brightness of the galaxies was also inspected. We have obtained these structural parameters from the following sources in the literature. For the SPARC galaxies, scale lengths are taken directly from Lelli et al. (2016a) and surface brightnesses from Verheijen (1997);de Blok (1997) and Swaters (1999). For LITTLE THINGS, scale lengths come from Read et al. (2016) and surface brightnesses from Hunter & Elmegreen (2006). Regarding SHIELD, size and surface brightness come from Teich et al. (2016) and Haurberg et al. (2015), respectively. He et al. (2019) provides both size and surface brightness for their sample, and the values used for UGC 2162 come from Trujillo et al. (2017).
This paper has been typeset from a T E X/L A T E X file prepared by the author.