© The Authors 2024

1. Introduction

Low-mass galaxies with stellar masses ≲10⁹ M_⊙ provide a crucial testing ground for the currently favored cosmological model, with dark energy plus cold dark matter (ΛCDM), through detailed comparison of observations and cosmological hydrodynamical simulations. Due to their shallow gravitational potential wells, dwarf galaxies are more sensitive to different baryonic effects, such as stellar feedback (e.g., stellar winds, supernovae, photoionization) and gas turbulence and can hence inform on the current implementation of sub-grid physics in simulations (e.g., Naab & Ostriker 2017). Properly accounting for these effects is important for studying the underlying connection between baryonic matter and dark matter and thereby testing predictions of the ΛCDM model. The significance of the dwarf galaxy population has already been seen in existing discrepancies, such as the large diversity in shapes of rotation curves of observed dwarfs (e.g., de Blok et al. 2008; Oh et al. 2015; Oman et al. 2015) or the too-big-to-fail problem (Ferrero et al. 2012; Papastergis et al. 2015; Papastergis & Shankar 2016) (for a more complete picture, we direct the reader to review articles Bullock & Boylan-Kolchin 2017; Tulin & Yu 2018; and Sales et al. 2022).

Galaxies in relatively isolated environments (field galaxies) make the best targets for tackling the above problems, as their properties are a direct outcome of the initial conditions, dictated by cosmology, and their internal evolution, with minimal environmental impact. Field low-mass galaxies tend to be star forming and gas rich, with gas fractions growing toward lower stellar masses (e.g., Huang et al. 2012a; Catinella et al. 2018). This makes neutral hydrogen (H I) a critical (and often dominant) baryonic component of these systems, and hence H I observations are an extremely useful tool for studying these systems. Furthermore, H I observations obtained with interferometers can provide spatially resolved kinematic information, thereby allowing for the modeling of rotation curves and full insight into the underlying gravitational potential (baryons plus dark matter).

Recent studies on a specific sub-class of low-mass galaxies with low surface brightnesses and large physical sizes, the so-called ultra-diffuse galaxies (UDGs), have pointed toward additional challenges for current models that are struggling to reproduce their observed properties (e.g., Di Cintio et al. 2017; Kong et al. 2022). The UDGs are defined as having a mean effective surface brightness ⟨μ⟩_e, X > 24 mag arcsec⁻² and effective radii R_e, X > 1.5 kpc (van der Burg et al. 2016), where X is an optical photometric band (usually g or r). They were initially identified in large amounts in cluster environments (van Dokkum et al. 2015; Mihos et al. 2015; van der Burg et al. 2016; Venhola et al. 2017; Mancera Piña et al. 2019a; La Marca et al. 2022) but have since been progressively detected also in the field (Leisman et al. 2017; Román & Trujillo 2017; Zaritsky et al. 2023). While the cluster population is gas poor, as expected for cluster environments, field UDGs tend to be gas rich, with gas fractions going up to very high values (⟨M_H I/M_⋆⟩ ≃ 35 for extremely optically faint systems, Leisman et al. 2017; Janowiecki et al. 2019; Poulain et al. 2022). However, it is not yet clear if the two populations have the same origin. The current formation scenarios for the field population (and possibly also the cluster population) include strong stellar outflows (e.g., Di Cintio et al. 2017) and high-spin parameters of host dark matter halos (e.g., Rong et al. 2017), while the cluster population has additional formation mechanisms based on tidal interactions and collisions (e.g., Liao et al. 2019).

Generally, observed dwarf galaxies exhibit a large diversity in their properties when compared to higher mass systems. One such example is the diversity of physical sizes (measured as effective radii R_e) present as a larger scatter in the observed stellar mass–size relation for dwarf galaxies than for higher mass galaxies (e.g., Lange et al. 2016). The increase in the scatter at low masses seems to also be present in the baryonic Tully-Fisher relation (BTFR; see, e.g., Trachternach et al. 2009; Sales et al. 2017; Iorio et al. 2017; McQuinn et al. 2022), which connects the baryonic mass to circular velocity and is a very tight scaling relation for high-mass late-type galaxies (e.g., Lelli et al. 2016b, 2019). It is not yet clear if the observed scatter at low masses is intrinsic or if some of the assumptions (e.g., that the kinematics of H I correctly traces the circular speed of the underlying dark matter halo) do not hold at these scales (Verbeke et al. 2017; Downing & Oman 2023). The H I-rich UDGs further complicate this question, as they seem to be systematically offset from the BTFR toward low circular velocities, as revealed by a spatially resolved H I study of Mancera Piña et al. (2019b, 2020) but as is also supported by unresolved studies (Leisman et al. 2017; Karunakaran et al. 2020; Hu et al. 2023).

High gas fractions of field dwarf galaxies make H I surveys extremely useful for finding these systems. This is especially true for the low surface brightness dwarfs, which are exceptionally hard to find in optical surveys. The Arecibo Legacy Fast ALFA (ALFALFA) H I survey (Giovanelli et al. 2005), a single-dish untargeted H I survey conducted with the Arecibo telescope, has already demonstrated the potential of H I studies in finding galaxies overlooked in optical catalogs (e.g., Cannon et al. 2015; Janowiecki et al. 2015; Leisman et al. 2017). As one might expect, H I-selected galaxies have systematic differences with respect to the optically selected ones, for instance, they tend to be star forming and bluer (e.g., Kovac 2007; Huang et al. 2012a,b; Durbala et al. 2020), as expected from the presence of H I as fuel for star formation.

As mentioned, for a robust determination of kinematic and dynamical properties of galaxies, resolved H I observations have proven crucial, as they allow us to determine how regular a system is and, for those dominated by rotation, extract a rotation curve that can be used for dynamical modeling (e.g., van Albada et al. 1985; de Blok et al. 2008; Read et al. 2016; Mancera Piña et al. 2022a). Previous H I-resolved studies of low-mass galaxies have often been the result of follow-up interferometric observations of preselected samples, selected from either untargeted single-dish H I surveys (e.g., Cannon et al. 2011; Gault et al. 2021) or optical surveys (e.g., Hunter et al. 2012). Being constrained by the preselection of a limited number of candidate sources, these studies focused either on very local targets (e.g., Swaters et al. 2002; Hunter et al. 2012) or on specific types of galaxies, for example, lowest H I masses (Cannon et al. 2011) or UDGs (Gault et al. 2021). Due to these constraints, the intermediate regime between standard dwarf galaxies and extreme UDGs has not been explored in detail. A previous untargeted H I survey conducted with the Westerbork Synthesis Radio Telescope (WSRT; Kovac 2007) had a good spatial resolution (∼30″) for resolving very nearby dwarfs. However, this survey focused on a specific sky region toward the Canes Venatici groups of galaxies and had a limited total sky area (30° × 30°) when compared to single-dish surveys (e.g., ALFALFA footprint of 6900 deg²), thereby potentially providing a less representative sample of the total dwarf galaxy population.

In addition to kinematic modeling, resolved H I data allow for comparison between stellar and H I geometries. While these are generally consistent for most galaxies, UDGs have been shown to often have misaligned H I and stellar morphologies (Gault et al. 2021). Exploration of these misalignments is especially important for the correct determination of rotational velocity, which strongly depends on the disk inclination. Most previous studies have assumed stellar-to-gas disk alignment and have considered only one component for the measurement of disk geometry, with some measurements done on H I total intensity maps (Mancera Piña et al. 2020) and others using optical images (Karunakaran et al. 2020; Hu et al. 2023).

In this work, we draw a sample of low-mass galaxies from an untargeted H I survey undertaken with the phased-array feed, APERture Tile In Focus (Apertif) (van Cappellen et al. 2022; Adams et al. 2022), for the WSRT. With a large sky coverage enabled by Apertif and the good spatial resolution of the WSRT, this allowed us to both find the field population of low-mass galaxies and conduct the kinematic modeling of the galaxies to obtain their rotation curves. It also allowed us to compare the stellar and H I geometries of low-mass galaxies. This enabled us to position our sample in the stellar mass–size relation and the BTFR and work toward linking the standard dwarf population to the UDG population.

Additionally, the Apertif dataset offers a unique opportunity to explore the frequency of dwarf galaxy pairs or multiples using H I observations, as they are easily resolved with Apertif. Dwarf galaxy interactions are thought to play a role in the evolution of dwarf galaxies, for example, by igniting or temporarily suppressing the star-formation process (e.g., Stierwalt et al. 2015; Kado-Fong et al. 2024). However, it is not yet clear how often these interactions take place and what role they have in the evolution of the overall population of dwarf galaxies. Previous studies in this domain have exclusively used optical spectroscopic data (e.g., Sales et al. 2013; Besla et al. 2018). The importance of including H I in such studies can already be seen by the smaller fraction of low-mass galaxies present in spectroscopic surveys when compared to H I surveys (Huang et al. 2012b).

The paper is organized as follows. In Sect. 2, we describe the Apertif H I data, PanSTARRS 1 (PS1) data that we use for the optical counterparts, and the source selection. Section 3 gives an overview of the kinematic analysis done on the H I data and the photometry measurements applied to the PS1 data. In Sect. 4, we present our results. We compare the global properties of our sample with other H I-selected samples in the literature, we compare geometries of the H I and stellar disks in our sample, and we comment on properties of UDGs in our sample. We place our galaxies on the stellar mass–size relation and the BTFR in Sect. 5. In Sect. 6, we discuss the frequency of pairs in our sample and possible biases in our selection procedure; and we further discuss the properties of our H I-selected sample compared to optically selected ones. Finally, in Sect. 7, we state our conclusions.

2. Data

2.1. Apertif H I data

Apertif (van Cappellen et al. 2022) is a phased-array feed receiver system designed for the WSRT. It produced 40 instantaneous compound beams on the sky, thereby increasing the field of view of the telescope up to 8 deg² making it a natural instrument for a wide area survey.

The Apertif imaging survey (Adams et al. 2022) operated between the 1 July 2019 and 28 February 2022. It observed selected regions of the sky above a declination of 30°, simultaneously conducting neutral hydrogen (HI) spectral line, continuum, and polarization imaging surveys. Each individual Apertif observation is 11.5 hours. The H I cubes are produced over the topocentric frequency range 1292.5–1429.3 MHz (Adams et al. 2022) and have spatial resolution of 15″ × 15″/sin δ. The Apertif imaging surveys have two tiers: the Apertif Wide-area Extragalactic Survey (AWES) aimed at covering large area of the sky with ∼2200 deg² of coverage and one or two observations per field; and the Apertif Medium-deep Extragalactic Survey (AMES) aimed to go deeper, targeting specific smaller area of the sky, with up to ten observations per field and the total sky coverage of ∼130 deg² (see Adams et al. 2022).

We use a preliminary H I source list made by running the Source Finding Application (SoFiA) (Serra et al. 2015; Westmeier et al. 2021) on individual H I cubes from single observations taken from August through December 2019 (Hess, priv. comm.). These cubes are separate for individual compound beams and have a typical noise of around 1.6 mJy beam⁻¹ over 36.6 kHz (∼8 km s⁻¹) (Adams et al. 2022). The SoFiA algorithm conducts source finding by using multiple resolution kernels and searches for detections across different kernels based on the S/N. For the production of the source list used in this work, SoFiA was run using three spatial kernels: 15″, 23.4″ and 39″; and three spectral kernels: 7.7 km s⁻¹, 23.2 km s⁻¹ and 54.1 km s⁻¹ (Hess et al., in prep.). The mask of the source was then constructed using all pixels that had S/N > 3.8 in one or more combinations of spatial and spectral smoothing kernels. In the next step, SoFiA corrects for false positive detections by comparing positive and negative detections, with reliability threshold of 0.85 in our case. For sources deemed real, the pipeline proceeds to calculate the H I properties of the detection by producing moment maps and global spectral profiles. Global spectral profiles are constructed using a 2D projection of the initial 3D mask of the source and integrating corresponding values in each channel, both inside and outside the original 3D mask. We note that systemic velocities in the preliminary source list are given in topocentric reference frame. The barycentric correction would be ≲30 km s⁻¹, corresponding to 1.5% at velocities used in this work (> 2000 km s⁻¹). Preliminary distances are calculated from systemic velocities assuming Hubble flow with H₀ = 70 km s⁻¹ Mpc⁻¹.

We use the preliminary source list only for the selection of the sample (described in Sect. 2.3), while the H I cubes used for the kinematic modeling come from a subsequent source finding (currently in progress) that was run on data that have been co-added from all available observations and mosaicked between compound beams within individual fields (full description of the improved processing is in Hess et al., in prep.). Therefore, the data products used in this work have better signal-to-noise ratios than predicted by the preliminary source list.

The properties of the cubes used in this work are given in Table 1. We note that the spectral resolution degrades during the co-adding step due to shifting to a common velocity frame. Noise levels also vary due to differing numbers of observations per field, independent mosaicking of individual fields producing non-uniform noise on the large scale, and instrumental effects (e.g., continuum subtraction) within the Apertif cubes themselves.

2.2. Optical data

For the determination of the stellar properties of our galaxy sample, we use data from the Pan-STARRS 1 (PS1) survey (Chambers et al. 2016), as the only optical survey to date encompassing the whole Apertif coverage. This is a broadband photometric survey made with the 1.8 meter telescope stationed at the Haleakala Observatories in Hawaii. The PS1 survey covers the whole sky region north from δ = −30°, which completely includes the Apertif coverage. In this work, we use g-, r-, and i-band PS1 images with median seeings of 1.31″, 1.19″, and 1.11″. Following the procedure described in Appendix A from (Román et al. 2020), we measured the surface brightness depth of PS1 images at the 3σ noise level over 10″× 10″ area obtaining μ(3σ_10″ ×10″) = (27.43, 27.21, 27.12) mag arcsec⁻² in g-, r-, and i-bands, respectively. In comparison, PS1 images are around 0.5 mag arcsec⁻² deeper than the Sloan Digital Sky Survey (SDSS) (Abolfathi et al. 2018) and about 1 mag arcsec⁻² shallower than The Dark Energy Camera Legacy Survey (DECaLs) (Dey et al. 2019).

2.3. Source selection

As we are interested in gas-rich low surface brightness, low stellar mass galaxies, we base our selection on the preliminary H I source list from Apertif (see Sect. 2.1). This enables us to find galaxies that could otherwise go undetected in optical surveys. In addition, we use the PS1 source catalog as a final check in order to exclude high stellar mass galaxies with low gas mass fractions.

The selection of sources from the preliminary H I source list of Apertif was based on the following constraints:

• (a)

  H I mass: M_H I < 10¹⁰ M_⊙,

• (b)

  width of the global profile at 50% of the peak emission: W₅₀ < 150 km s⁻¹,

• (c)

  average signal-to-noise ratio per channel in the global profile: (S/N)_ch > 3,

• (d)

  systemic velocity: V_sys > 2000 km s⁻¹,

• (e)

  the H I disk resolved by at least 3 Apertif beam elements.

Additionally, we considered only sources for which all observations are fully processed and have been co-added, so the improved data products are available. In the following, we describe the procedure in more detail.

Conditions (a) and (b) are imposed in order to exclude high-mass galaxies. While the latter condition seems to also exclude high inclination galaxies, we point the reader to Sect. 6.2, where we discuss the possible biases of our selection procedure on the inclination.

Condition (c) excludes low S/N detections. As we intend to kinematically model our sample using the total 3D information (see Sect. 3.2), we applied our condition on the average signal-to-noise ratio per channel in the global profile, defined as

$$\begin{array}{c}\hfill {(S/N)}_{\mathrm{ch}}=\frac{F_{H\text{I}}[\mathrm{Jy}\mathrm{Hz}]}{{\mathcal{N}}_{\mathrm{ch}}·\mathrm{\Delta }\nu [\mathrm{Hz}]·\mathrm{rms}[\mathrm{Jy}]},\end{array}$$(1)

where F_H I is the total detected flux, 𝒩_ch the number of channels with detected emission, Δν the channel width, and rms the root mean square noise measured in channels without emission in the global profile of the source (constructed as described in Sect. 2.1). We chose this definition because it gives an estimate of the average level of emission $(F_{H\mathrm{I}}[\mathrm{Jy}\mathrm{Hz}]/({\mathcal{N}}_{\mathrm{ch}}·\mathrm{\Delta }\nu [\mathrm{Hz}]))$ compared to the average level of integrated noise (rms [Jy]) per channel, providing a better estimate for success in producing reliable 3D kinematic models by ensuring sufficient signal in each channel.

Condition (d) aims to minimize the uncertainties introduced by peculiar motions on the estimated distances. The threshold of 2000 km s⁻¹ is chosen so that distance errors are not dominated by typical peculiar velocities of galaxies in the field (one to a few hundred kilometers per second).

Condition (e) aims to ensure that the galaxy is spatially resolved enough for the construction of a reliable 3D kinematic model. Given that H I diameter is not measured in the source finding with SoFiA, we estimated the physical sizes of sources in the preliminary source list using the H I mass – H I diameter relation from (Wang et al. 2016):

$$\begin{array}{c}\hfill log\frac{d_{\mathrm{H}\text{I}}}{\mathrm{kpc}}=(0.506\pm 0.003)log\frac{M_{\mathrm{H}\text{I}}}{M_{\odot }}-(3.293\pm 0.009),\end{array}$$(2)

where d_H I is measured at a surface density of 1 M_⊙ pc⁻². Due to the unknown orientation of the galaxy, we take the geometric mean of the beam axes as an estimate of the beam diameter. We then include sources for which the estimated angular diameter is three or more times larger than the beam diameter. We note that this condition ultimately corresponds to an H I flux cut, but given that we apply it to individual detections (with different beam sizes, see Table 1), it cannot be applied as a universal flux cut on the sample.

After applying the above constraints on the H I content and obtaining a list of candidates, we move on to the optical counterpart. In order to exclude optically bright and high stellar mass sources from our selection, we put a lower limit of −18.5 on the i-band absolute magnitude of candidate sources. This threshold roughly corresponds to stellar masses of ≲10⁹ M_⊙, based on a typical g − r color of 0.2, following the relation from (Herrmann et al. 2016). For this, we search the PS1 source catalog inside the 20″ radius from the H I center position of each source. We excluded the contamination of foreground stars by following the procedure described in Farrow et al. (2014); that is, we put a lower limit threshold of 0.05 mag on the difference between the measured Kron and point spread function (PSF) magnitudes in i-band for each detection. Finally, we transfer the apparent magnitudes in the i-band (m_i) to absolute magnitudes (M_i) using the Hubble flow distance estimate from the H I source list and applied the cut of M_i > −18.5.

After applying the selection criteria described above, we manually checked the PS1 and H I data of the selected candidates. In a few cases, SoFiA detected an H I tail of a larger galaxy as a separate source. In addition, some pairs of galaxies were detected as a single source and the center position of the detection was placed between the two galaxies. These galaxies were in some cases outside of the 20″ radius in which the PS1 catalog was searched for detections. For these cases, we manually searched the PS1 catalog again and excluded the sources in which the corresponding galaxies failed to pass the required selection criteria.

Out of 1231 sources in the source list, 76 detections passed the H I criteria, out of which 47 had fully processed data needed for the production of improved data products (data with co-added observations and mosaicked across compound beams). Out of these, 24 passed the optical criteria.

Properties of H I cubes of the obtained sample are provided in Table 1. Refined physical properties of the selected galaxies (see Sect. 3.4) are listed in Table 2. Note that for galaxies in pairs (see Sect. 4), the cube was separated in two to allow for measurements of properties of individual galaxies. The obtained H I masses range between 8.6 ≲ log (M_H I/M_⊙) ≲ 9.7, while the stellar masses range between 8.0 ≲ log(M_⋆/M_⊙)≲9.7. While the stellar mass range seems to go to higher values than should be permitted by our selection, we note that only three galaxies in the sample have log (M_⋆/M_⊙) > 9.0.

When available, we provide names from the literature obtained by searching the NASA/IPAC Extragalactic Database¹ (NED). For seven galaxies in our sample, there was no recorded entry in NED, and four previously known galaxies have no archival redshift measurement. All of the galaxies previously detected in H I had only single-dish observations and were spatially unresolved. When referring to individual galaxies, we use their literature names when they are shorter and Apertif names otherwise (shortened as AHCJhhmm+ddmm).

3. Analysis

3.1. Summary of H I masks used in this work

In this work, we used several masking techniques for the H I data, depending on the purpose of the task. In this section we briefly describe each mask for easier reference.

As mentioned in Sect. 2.1, SoFiA source finder produces a mask of the source after each detection. This mask is produced using multiple spatial and spectral resolution kernels and extends far enough out to enclose the galaxy emission down the level of noise. Because of this, it allows us to robustly measure the total flux of the galaxy, which we use for the calculation of the H I mass (see Sect. 3.4).

In contrast, when performing kinematic modeling using the whole 3D cube (see Sect. 3.2), we prefer to limit ourselves to regions of high S/N at the full spatial resolution of the data cube in order to minimize the effect of noise on the final kinematic model. In this case, we make use of the default masking scheme from ^3DBarolo (SEARCH option, see ^3DBarolo documentation²) to produce a mask that includes most of the galaxy’s emission, while minimizing the influence of external noise peaks.

And lastly, for the estimation of the geometry of the H I disk using the total intensity map (see Sect. 3.2), we wish to include low level emission to capture the outskirt of the galaxy, but still stay in the mid-high S/N regime in order to limit the influence of the noise on the derived geometry. To produce this mask, we again use ^3DBarolo, but this time we add an option to smooth the cube 1.2 times before selecting regions of mid-high S/N.

3.2. Three-dimensional kinematic modeling

To obtain the rotational velocities of our sample, we adopt the tilted ring model (Rogstad et al. 1974; Begeman 1987b). In this model, the gas is assumed to be in pure circular motion. The galaxy is broken down into concentric rings of different radii, each with its own set of parameters. For simplicity, we group parameters into two categories: geometric (coordinates of the center (x₀, y₀), thickness of the ring (z₀), inclination (i), and position angle (PA)); and kinematic (systemic velocity (V_sys), rotational velocity (V_rot), and velocity dispersion (σ_V)). We make distinction between two position angles (both measured anticlockwise from the north direction): the morphological (PA_H I) describing the geometry of the H I disk as projected on the sky, and kinematic (PA_kin) describing the angle with the steepest gradient in velocity. The tilted ring model was historically widely used for studying galaxy kinematics from 2D velocity field maps (e.g., Begeman 1987a). More recent studies applied the model directly to 3D emission line data cubes (e.g., Józsa et al. 2012; Di Teodoro & Fraternali 2015), which proved essential for accurate determination of underlying kinematics by enabling the correction of the beam smearing effect (e.g., Swaters et al. 2009).

3.2.1. ^3DBarolo setup and assumptions

To derive kinematic properties of our sample, we use the ^3DBarolo (Di Teodoro & Fraternali 2015) software. ^3DBarolo was designed for (and well tested on) poorly resolved data (Di Teodoro & Fraternali 2015; Mancera Piña et al. 2020; Roman-Oliveira et al. 2023), making it a well-suited tool for this study. The software builds 3D model cubes based on the tilted ring model and compares them with the data cube to find the best fitted model. The main strength of this software is the convolution step where the model is convolved with a Gaussian beam before the calculation of residuals between the model and the data cube, thereby minimizing the effect of the beam smearing on the model.

Here we describe a few key parameters when running ^3DBarolo, while the full parameter file is given in Appendix A. As mentioned in Sect. 3.1, we produce the mask using the default masking option by ^3DBarolo. Furthermore, we choose the azimuthal normalization for the flux of the model cube, as our measurements do not have enough spatial resolution to trace gas distribution in much detail. We choose cos²(θ) for the weighing function in the computation of residuals, thereby giving priority to the major axis where rotation is better traced. Finally, we put a limit for the minimum velocity dispersion to half of the spectral resolution to avoid unrealistically low dispersions that can arise in low S/N conditions. We note that in our final kinematic models, the obtained velocity dispersions never reach this floor value.

When running ^3DBarolo, we make a few assumptions about the geometric parameters. We do not allow x₀, y₀, i, and PA_kin to change between different radii, as we cannot track changes in geometric parameters across different rings due to the relatively poor spatial resolution of the galaxies in our sample. The only exception is UGC 12027, which is sampled by eight compound beams along the major axis (for the 1.5 times smoothed cube used for the modeling; see Sect. 3.2.3) and shows signs of a warp. For this reason, we let PA_kin and i change between radii in this case. Additionally, any effect of the disk thickness in the data is dominated by the beam. Physical beam diameters in the sample are all ≥2.3 kpc, while dwarf galaxy disk thickness ranges from a few hundred parsecs to ∼1 kpc (e.g., Banerjee et al. 2011; Das et al. 2020; Bacchini et al. 2020). Even in our physically best resolved case, kinematic models produced with z₀ = 0 pc and z₀ = 500 pc are fully consistent. Therefore, we assume a razor thin rotating H I disk (z₀ = 0 pc).

3.2.2. Kinematic modeling procedure

We constrain the geometric parameters of the H I disk using a publicly available python script CANNUBI³ (briefly described in Mancera Piña et al. 2022b; Roman-Oliveira et al. 2023). CANNUBI uses a Markov Chain Monte Carlo (MCMC) approach where values of parameters are evaluated based on the residuals between the data and the corresponding model produced with ^3DBarolo, made using either total flux maps, or the 3D cubes themselves. We used the total flux maps as the computational time is notably shorter. The mask used in this step is defined in Sect. 3.1. With CANNUBI, we fit x₀, y₀, i, PA_H I, and the extent of the H I disk. We note that, given the angular resolution of our data, our assumption of the razor thin H I disk does not compromise the derivation of the H I inclination as long as true inclination is below ∼80°. In this regime, the residuals of moment 0 maps between models with z₀ = 0 pc and z₀ = 500 pc are lower than the typical rms noise in the Apertif maps, thereby making the two models indistinguishable.

For the final kinematic fit with ^3DBarolo, we use the median values of posteriors from CANNUBI to constrain the geometric parameters. In one case (AGC 239039), we used the optical inclination, as the optical geometry was significantly better aligned with the galaxy kinematics. The morphological position angle is given as an initial guess for the kinematic position angle (PA_kin) (except for three galaxies, UGC 8602, AGC 239039 and AHCJ2218+4059, for which we provide no initial guess on PA_kin), but we let ^3DBarolo fit it as part of a two-stage run. In the two-stage run, ^3DBarolo first fits unconstrained parameters (V_sys and PA_kin in our case) together with V_rot and σ_V in each ring, after which it fits a functional form (a constant value in our case, see Sect. 3.2.1) to the radial distribution of these parameters. In the second run, ^3DBarolo again fits V_rot and σ_V, now keeping all other parameters fixed.

We propagate errors on the inclination by running ^3DBarolo with the same parameter files as in the original run (with inclination i), but now changing the inclination to i ± σ_i. The error on the final rotational velocity is then taken to be the difference between the original case and the case i + σ_i (i − σ_i) for lower (upper) error estimate when the difference is higher than the statistical error obtained from ^3DBarolo. Otherwise, we adopt the ^3DBarolo error on rotational velocities.

3.2.3. Obtained kinematic models

Out of 24 initial galaxies, 13 showed a velocity gradient that can be interpreted as regular rotation, and had sufficient S/N per channel that enabled us to perform the modeling. In a few cases with low S/N (AHCJ1308+5437, AHCJ1359+3726, UGC 12027, and UGCA 363, Figs. C.13⁴, C.14⁴, C.20⁴, and C.22⁴), we smoothed the cubes 1.5 times the synthesized beam before running CANNUBI to increase the S/N and consequently pick up the faint emission otherwise hidden in the noise. We also used these smoothed cubes for the final kinematic modeling in these cases.

The results of kinematic modeling are listed in Table 3. We report the maximum value of rotational velocity, and the mean value of velocity dispersion across all rings. V_sys is reported in Table 2. For kinematically modeled galaxies, we report V_sys of the best fit model, while for the rest of the sample it is measured as a weighted mean from the global spectral profile. The error corresponds to half of the spectral resolution.

In the description that follows we show the results of kinematic modeling for one galaxy (AHCJ2239+3832), while the rest can be found in Appendix C. Results for AHCJ2239+3832 are shown in the form of channel maps in Fig. 1, where the top panels represent the data and the bottom panels the best fit model. The model is present in all channels where the galaxy emission is detected and represents the data well. Moment maps and position-velocity (PV) slices for AHCJ2239+3832 are shown in Fig. 2. On the leftmost panel of the figure, the H I contours are overlaid on top of the optical image. We note that the noise in the total intensity map is not uniform because the cube is masked with a 3D mask before producing the map, resulting in a different number of channels contributing per pixel. Therefore, the H I contours in the image are the so-called pseudo Xσ contours, obtained by selecting pixels with X-0.25 and X+0.25 values in the S/N map, and taking the average of the flux of the corresponding pixels in the total H I map (see, e.g., Verheijen & Sancisi 2001). The blue ellipse in the image describes the obtained disk geometry from CANNUBI. Contours are spread outside of the ellipse due to the prolongation of the beam in the north-south direction. The velocity field in the middle left panel shows a velocity gradient that suggests the presence of a regularly rotating disk. This is also evident in the PV slices, where data (blue with black contours) shows no major deviation from regular rotation and is well described by the model (red contours).

Fig. 1. Comparison between channel maps of AHCJ2239+3832 and its best-fit model. Top panels and blue contours represent the data, while the lower panels and red contours represent the best fit model. Contours are plotted starting from three times the noise per channel in the cube reported in Table 1, and are spaced by a factor of two in intensity. The black X indicates the center of the galaxy, as determined by CANNUBI (see Sect. 3.2).

Fig. 2. H I kinematics of AHCJ2239+3832. Leftmost: Color-composite image from PS1 g-, r-, and i-bands overlaid with H I contours. Contours are set to start at the level of intensity corresponding to a pseudo 4σ contour (see Sect. 3.2 for explanation) in the total intensity map (white), and grow by a factor of two in intensity toward the redder colors. Lowest contour for AHCJ2239+3832 corresponds to the column density of 7.2 × 1019 cm−2. The overplotted ellipse in blue represents the median geometric parameters obtained from CANNUBI (including the size). Middle left: Velocity field obtained as weighted-mean value from the H I data from 3DBarolo. The black cross represents the best-fit center position and the gray dashed line the kinematic position angle. Contours are given in spacing of 10 km s−1 with the green contour indicating Vsys. Middle right: Position-velocity slice along the major axis with data (blue) and model (red) contours. White points represent the obtained projected rotational velocity of the best fit model. Contours are plotted starting from two times the noise per channel in the cube as reported in Table 1, and grow linearly. Rightmost: Position-velocity slice along the minor axis (perpendicular to the dashed line in the middle left panel) with the same color scheme and contour levels as the middle right panel.

Given the low spatial resolution of most of our sample, it is not straightforward to say whether we are able to trace the flat part of the rotation curve for our galaxies. However, in five cases (UGC 5541, AHCJ1308+5437, AHCJ2207+4008, AHCJ2207+4143, and UGC 12027; Figs. C.12⁴, C.13⁴, C.15⁴, C.16⁴, and C.20⁴), the PV slices along the major axes show a turnover that suggests that the flat part has been reached. While UGC 12027 shows this turnover, its H I inclination shows consistencies with inclinations down to 0°, making it impossible to robustly constrain the rotational velocity. Therefore, for UGC 12027 and the other galaxies in the sample that do not show evidence of a turnover in the PV slices, we consider our derived rotational velocities to be lower limits to the rotational velocity in the flat part of the rotation curve.

From 13 kinematically modeled galaxies, one (UGC 8602; Fig. C.5⁴) seems to be kinematically lopsided, that is, its rotational velocity is higher on one side than the other or one side of the galaxy is not detected far enough out. In either case, this complicates the interpretation of the derived rotational velocity as a tracer of the underlying potential. Another two kinematically modeled galaxies show signatures of warps, with AHCJ2218+4059 (Fig. C.19⁵) showing a warp in position angle, and UGC 12027 (Fig. C.20⁵) showing a warp in both inclination and position angle. Unfortunately, we cannot robustly model the warps due to the low spatial resolution of our sample. For these reasons, we flag these three galaxies in Table 3 and do not include them in the BTFR in Sect. 5.2.

3.2.4. Asymmetric drift correction

As mentioned before, the major contribution to gas kinematics in disk galaxies comes from rotation. However, gas pressure can still have a significant influence in counteracting the pull of gravity. In order to correct for pressure support and determine the circular speed (which allows us to characterize the gravitational potential), we follow the procedure from Iorio et al. (2017). In particular, we made use of their Eq. (9):

$$\begin{array}{c}\hfill V_{\mathrm{A}}^2=-R{\sigma }_{\mathrm{V}}^2\frac{\partial {\sigma }_{\mathrm{V}}{\mathrm{\Sigma }}_{\mathrm{H}\text{I}}cosi}{\partial R},\end{array}$$(3)

where V_A is the asymmetric drift correction (ADC), σ_V the velocity dispersion, Σ_H I the projected surface density of H I, and R the distance from the galaxy center. Given that the rotation curves are sampled by two points for most galaxies in our sample, we fit the term within the derivative with a linear function, and multiply its slope by −R σ_V² in order to obtain the asymmetric drift. Then, the circular velocity (V_circ) is obtained with (e.g., Iorio et al. 2017):

$$\begin{array}{c}\hfill V_{\mathrm{circ}}^2={V_{\mathrm{rot}}}^2+V_{\mathrm{A}}^2.\end{array}$$(4)

The errors on inclination were propagated to the circular velocity the same way as was done for the rotational velocity (see Sect. 3.2).

The maximum circular velocities across the rings and the corresponding ADC terms are given in Table 3. All kinematically modeled galaxies in our sample seem to be highly rotationally supported with little contribution from pressure support. However, we note that given our resolution, we were not able to perform a robust ADC, for example, by fitting a more physically motivated function to the term under the derivative in Eq. (3) (e.g., Bureau & Carignan 2002; Iorio et al. 2017).

3.3. Obtaining optical properties

Our optical analysis consists of two steps: we first conduct isophotal fitting on the i-band image in order to constrain the geometry of the stellar disk, after which we use the obtained stellar geometry to extract surface brightness profiles (in all three bands) following the procedure described in Appendix A of Marasco et al. (2023).

During the isophotal fitting, we mask and slightly smooth the i-band images to lower the influence of small scale overdensities on the obtained geometry of the stellar disk. The size of the smoothing kernel depends on the galaxy and can be found in Table 4. For the fitting we use an Astropy affiliated package photutils (Bradley et al. 2022). The algorithm fits a set of ellipses of increasing semi-major axes using the position of the center, ellipticity, and position angle of the ellipses as free parameters. We run the algorithm two times. The first time we leave all the parameters free to vary, and use the result to constrain the central position by taking the median from all fitted ellipses with semi-major axes larger than half of the full width at half maximum (FWHM) of the PSF. The second time, we fix the central position and fit only for ellipticities and position angles. The final (global) values of ellipticity and position angle were obtained by taking the median value from ellipses, which range in estimated surface brightness between 24 and 27 mag arcsec⁻² (motivated by the classical choice of surface brightness of 25 mag arcsec⁻² as representative of the outer stellar disks) and which have semi-major axes larger than the total FWHM of the PSF in order to exclude the effect of the PSF on the derived geometry. To (at least) partially remove the bias toward the inner regions that were fitted with more ellipses, we calculated the median from semi-equidistant ellipses (along the semi-major axis) satisfying these conditions. The result of isophotal fitting can be found in Fig. 3 for AHCJ2239+3832, and in Appendix C for the rest of the sample.

Fig. 3. Isophotal fitting result for AHCJ2239+3832. Upper panels: The leftmost panel shows the smoothed data overlaid with semi-equidistant ellipses along the semi-major axis length whose parameters are marked with red stars in the lower panels. White regions denote masked pixels in the image. The middle panel shows the model built from all fitted ellipses (blue circles in the lower panels), and the rightmost panel shows the residual of the data image and the model image. Lower panels: Ellipticity and position angle as a function of the semi-major axis length from the second run of the fitting algorithm (see Sect. 3.3) are shown as blue circles, with respect to radius from the fixed center position. Red stars denote parameters of semi-equidistant ellipses plotted in the upper leftmost panel. The black vertical line is located at a distance from the center that corresponds to the FWHM of the PSF after the smoothing, and the gray shaded region corresponds to the region in which the final global geometry was measured (using only the semi-equidistant ellipses, see Sect. 3.3).

To transform from ellipticity (ϵ) to inclination, we used:

$$\begin{array}{c}\hfill {cos}^2(i)=\frac{{(1-ϵ)}^2-q_0^2}{1-q_0^2},\end{array}$$(5)

where q₀ is the intrinsic thickness. We adopt q₀ = 0.3, which is a common value used for dwarf irregular galaxies (e.g., Sánchez-Janssen et al. 2010).

A possible caveat to the accurate determination of the geometry of the stellar disk with isophotal fitting in this mass regime is the influence of clumpy star-forming regions on the obtained geometries of the isophotes. Bright clumpy regions can dominate over the fainter disk component of the galaxy, which we are trying to constrain. Indeed, some galaxies in our sample might be subject to this effect, as they clearly show many bright clumpy regions in their outskirts (e.g., UGC 5541, AHCJ1308+5437, UGC 12005; Figs. C.12⁵, C.13⁵, and C.23⁵). On the other hand, in a few galaxies with a less clumpy morphology (AHCJ2239+3832, UGC 8605, AGC 239039, AHCJ2207+4008; Figs. 3, C.4⁵, C.9⁵, and C.15⁵), there is a clear trend of isophotes becoming more aligned with the H I kinematics with increasing radius. This points toward another caveat in the determination of optical morphology: the depth of optical data. Mancera Piña et al. (2024) has shown that the apparent misalignment between stellar and gas geometry in the UDG AGC 114905 seen in data with surface brightness depths of μ(3σ₁₀″ × 10″)∼28.5 mag arcsec⁻² disappears when analyzing ultra-deep imaging reaching μ(3σ₁₀″ × 10″)∼32 mag arcsec⁻². This demonstrates a need for deeper optical observations in order to trace the underlying stellar disk far enough out for a robust comparison of H I and stellar morphologies.

Using the optical geometric parameters, we conduct surface brightness photometry on full resolution images following the procedure described in Appendix A of Marasco et al. (2023). We input the global geometric parameters and create a set of ellipses of the same geometry, but with increasing semi-major axes. The semi-major axis is sampled so that each ellipse has width equal to the FWHM of the PSF. For each ellipse, we calculate the mean value in image units (counts) and correct for inclination by multiplying the obtained values by cos i. For the estimation of the background and noise, we fit the sky pixel intensity distribution using two-components: a Gaussian function and a Schechter function. The latter accounts for the possible residing contamination from (masked) foreground stars. We obtain the image’s background value and rms noise from the parameters of the Gaussian component. In cases where the fit does not converge, we simply use the mean and the standard deviation of the sky pixel intensity distribution of the masked image. In most cases, we have extracted surface brightness profiles radially until S/N = 1 in a given ring is reached. In some cases⁶ the obtained profiles extended far outside the galaxy due to image artifacts or nearby foreground stars. In these cases we manually truncated the profiles. To convert to magnitudes, we use the calibration from PS1 (Waters et al. 2020):

$$\begin{array}{c}\hfill m[\mathrm{mag}]=-2.5{log}_{10}(\mathrm{\Sigma }\mathrm{counts})+2.5{log}_{10}(T)+25,\end{array}$$(6)

where T is the exposure time reported in the header. We fit the obtained surface brightness profile with a Sérsic function from which we derive the central surface brightness μ_0, X, Sérsic index n_X, effective radius R_e, X, mean effective surface brightness ⟨μ⟩_e, X, and the apparent magnitude m_X in each band indicated with X (see Appendix B for details). We list m_g, R_e, r and ⟨μ⟩_e, r in Table 4. Magnitudes and mean surface brightnesses are corrected for the Galactic extinction using results from Schlafly & Finkbeiner (2011), which are available through NED. Results of these fits in the i-band are presented in Fig. 4 for AHCJ2239+3832, and in Appendix C for the rest of the sample. We note that in some cases the H I center is offset from the optical center, but given the H I beam diameter size of ≳15″, the two centers are compatible. In two cases (UGC 8605 and AHCJ1308+5437 shown in Figs. C.4⁷ and C.13⁷, respectively), the extracted surface brightness profiles have irregular shapes due to clumpy star-forming regions, which is why we consider their fitted Sérsic parameters to be less reliable. Additionally, the surface brightness profile of AHCJ2207+4008 (shown in Fig. C.15⁷) clearly shows contribution from two components and cannot be well fitted with a single Sérsic function. These galaxies will be flagged in future plots when these parameters (or properties derived from them) are used. The obtained stellar properties for the sample can be found in Table 4.

Fig. 4. Photometry in the i-band of AHCJ2239+3832. Left: Optical i-band image overlaid with measured geometries from H I morphology (blue) and from i-band isophotal fitting (orange). The radius at which the optical geometry is plotted corresponds to the outermost data point in the surface brightness profile. The radius at which the H I geometry is plotted corresponds to the extent of the disk as obtained from CANNUBI. The crosses denote the obtained optical (orange) and H I (blue) centers. Right: Surface brightness profile (orange) and the corresponding best fit with a Sérsic profile (black).

3.3.1. Reliability checks

We found that the Sérsic function gives a good representation of the galaxy light profile, but it assumes the galaxy to extend toward infinity. As galaxies are finite systems, this could result in the overestimation of the total flux of a galaxy. However, the alternative approach of directly measuring the half-light radius (R_e^meas) and total magnitude from the data suffers from the problem of the depth of the data (the comparison of these two approaches in our sample can be found in Appendix B). Trujillo et al. (2001) analyzed this problem and showed that the two approaches would converge to the same value if the surface brightness profile is traced far enough out in the galaxy. To check the reliability of our fitted parameters, we used Eq. (18) from Trujillo et al. (2001) with our fitted R_e^fit values to find the desired outermost radius (in the r-band) for which the two approaches would theoretically give at least 85% compatibility (R_e^meas/R_e^fit ≥ 85%). Galaxies that were traced at least out to this radius are considered to have well constrained Sérsic fit parameters, and consequently, robust photometry. Excluding the galaxies with unreliable Sérsic fit (see Sect. 3.3) out of this analysis, we have 17 galaxies with reliable Sérsic fits and well constrained photometry. The other seven galaxies in the sample are considered to have less reliable derived Sérsic parameters and will be regarded as “optically unreliable”. We note that this only refers to parameters corresponding to the shape of the profile, while the obtained stellar masses are considered to be robustly measured (see Appendix B for the comparison between measured and fitted properties). The optically unreliable galaxies are marked in Table 4 and flagged in plots that use Sérsic parameters or optical properties derived from them.

In some galaxies, the isophotal fitting showed significant changes of geometric parameters with radius. To test how these variations influence the final results, we repeated the extraction and the fitting of surface brightness profiles using the 16th and the 84th percentile of ellipticity and position angle from the same selected region as before (outside the FWHM of the PSF and between surface brightness values of 24 and 27 mag arcsec⁻²). Each time, we changed one of the parameters (either i or PA_op) to a higher/lower value. The resulting difference between R_e of the initial run and R_e of four cases described above (16th or 84th percentile value for either i or PA_op), was shown to be ≲15%, while for the ⟨μ⟩_e the difference was up to 0.4 mag arcsec⁻². These variations were taken into account in the errors by taking the mean difference between the initial run and the four cases, and summing it in quadrature with the error from the fit for all fitted parameters of the Sérsic profile (μ_e, n, and R_e). For apparent magnitudes, the errors were calculated only by taking the mean variation for all different cases. We note that even taking the 16th and the 84th percentile from the whole meaningful range of radii (outside the FWHM of the PSF and surface brightness values lower than 27 mag arcsec⁻²), the variations in parameters were ≲30% in R_e and ≲1.5 mag arcsec⁻² in ⟨μ⟩_e.

3.3.2. The special case of AHCJ2207+4143

Optical images of AHCJ2207+4143 (Fig. C.16⁸) show a strong non-uniform foreground Galactic emission. This emission manifests as the gradient in the background levels across the optical images of the galaxy. In order to correct for this effect, we have produced a 2D image of the background using the Background2D class from the photutils python package. We then subtracted the background from the original images and proceeded with photometry as described in Sect. 3.3.

To characterize the errors of the background subtraction at the position of the galaxy, we additionally extracted surface brightness profiles from the original images (without the subtracted background) and used the difference in the obtained parameters as an additional error estimate. In particular, we updated the errors of Sérsic fit parameters of the profile from the background subtracted image by taking half of the difference between the two cases (subtracted and non-subtracted image) and summing it in quadrature with the statistical error from the fit. For the apparent magnitudes, we have calculated the errors by only taking the half of the difference in the obtained magnitudes between the two cases. This step was performed before the propagation of errors of geometric parameters (see Sect. 3.3.1). The rest of the analysis stays the same as for the rest of the sample.

3.4. Obtaining the distance and masses

3.4.1. Distance

We calculate distances using the Extragalactic Distance Database (EDD) (Tully et al. 2009). The EDD provides two calculators (Kourkchi et al. 2020) based on two different models of local motions: the Numerical Action Methods (NAM) model of Shaya et al. (2017) and the linear density field model of Graziani et al. (2019). As most of our galaxies are outside of the range of the NAM model (≲38 Mpc or ≲3000 km s⁻¹), for consistency, we used the linear model for determination of distances of all our sample. As is reported in Graziani et al. (2019), the linear model has accuracy of 15% or better, depending on the sky region. Unfortunately, the calculators do not provide errors for an individual position, which is why we adopted errors of 15% on all our distances.

Obtained distances are reported in Table 2. In most cases, the Hubble flow distance used for the selection is within a few megaparsecs of the EDD obtained distance, well within the 15% errors. However, in six cases, the difference between the two distances is ∼10 Mpc. This corresponds to ∼20% in absolute error, which we find significant enough to adopt the EDD distances in the rest of the paper. Nonetheless, we note that this difference would not significantly change the outcome of our selection procedure, it would only add two more galaxies to our sample. These galaxies were not initially selected using the Hubble flow distance due to our absolute magnitude cut. We do not include these two galaxies to our sample.

3.4.2. Masses: H I, stellar, and baryonic

We calculated the H I masses using a relation from Kennicutt & Evans (2012):

$$\begin{array}{c}\hfill \frac{M_{\mathrm{H}\text{I}}}{{\mathrm{M}}_{\odot }}=2.343·{10}^5·{\left(\frac{D}{\mathrm{Mpc}}\right)}^2·\frac{F_{\mathrm{H}I}}{{\mathrm{Jy}\mathrm{km}\mathrm{s}}^{-1}},\end{array}$$(7)

where F_H I is the total H I flux and D is the distance to the galaxy. For the calculation of the flux, we used the mask provided by SoFiA, which was produced using multiple spatial and spectral resolution kernels in order to capture all the galaxy emission down to the level of noise. We assume a 15% error on F_H I coming from the calibration of the flux scale, primary beam correction and mosaicking of Apertif compound beams (Kutkin et al. 2022). The obtained H I masses are listed in Table 2.

We obtained stellar masses by applying the mass-to-light color relation for dwarf irregular galaxies from Herrmann et al. (2016):

$$\begin{array}{c}\hfill {log}_{10}{(M/L)}_{g\prime }=-0.601(\pm 0.090)+1.294(\pm 0.401)(g\prime -r\prime ),\end{array}$$(8)

where apparent magnitudes are given in the standard photometric system. We transform the PS1 extinction corrected (see Sect. 3.3) g- and r-band magnitudes by applying the conversion from Eq. (6) in Tonry et al. (2012). For a robust estimation of galaxy color, we extracted additional surface brightness profiles using the largest FWHM of the PSFs between the three bands for each galaxy, truncated them at the same radius, and directly measured apparent magnitudes by integrating the profiles. We prefer this direct approach exclusively for the estimation of color, where only the relative difference in magnitudes is important. In all other situations when we make use of magnitudes, we adopt the value from the Sérsic fit, including the transfer from the g-band magnitude to g-band luminosity in the calculation of the mass. Colors and g-band magnitudes are given in Table 4, while stellar masses are listed in Table 2.

Finally, for the estimation of baryonic mass, we used

$$\begin{array}{c}\hfill M_b=1.33M_{\mathrm{H}\text{I}}+M_{\ast },\end{array}$$(9)

where a factor of 1.33 corrects for the helium contribution to the gas mass of galaxies. Molecular gas mass in low-mass galaxies has been shown to be significantly smaller than H I and stellar masses (around 10% of either) (see e.g., Leroy et al. 2009; Bothwell et al. 2014; Accurso et al. 2017; Ponomareva et al. 2018; Catinella et al. 2018). Therefore, we neglect it in the calculation of our baryonic mass. The obtained masses are listed in Table 2.

4. Properties of the sample

The final sample consists of 24 galaxies whose global properties can be found in Table 2. We discuss and provide specific notes on individual galaxies in Appendix C.

We note that six galaxies are in pairs, with each galaxy having a distinct optical counterpart. We present their results individually throughout this section but note that two pairs, consisting of AHCJ0203+3714-1 with AHCJ0203+3714-2 and CGCG 290-011 with UGC 5480, are clearly interacting as seen from the H I bridges in Fig. C.1⁹. Consequently, their properties might be influenced by the interaction. The third pair (UGC 8605 with UGC 8602) does not show a detectable direct interaction, allowing us to obtain reliable kinematic models for each of them.

One galaxy in the sample (UGC 8503) is positioned at the edge between two separate Apertif fields, and was detected two times (once in every field). Unfortunately, the H I data cubes have low S/N and the two detections show highly inconsistent morphology, which is why we only use these data in order to estimate the H I mass of the source, but do not proceed with the characterization of the H I disk nor the kinematic modeling.

4.1. Global properties of the sample

We put our sample in the context of The Arecibo Legacy Fast ALFA (ALFALFA) H I survey (Giovanelli et al. 2005). We use the completed α.100 (Haynes et al. 2018) catalog for the comparison of H I properties. For consistency, in this comparison we regard our three pairs of galaxies as single sources because they would not be resolved by the Arecibo beam. For comparison of stellar masses, we use the stellar mass of the brighter galaxy in the pair (but the H I mass of the whole pair) in order to be consistent with the ALFALFA-SDSS catalog (Durbala et al. 2020). Additionally, we selected a subsample of ALFALFA galaxies inside the same V_sys range as our sample (i.e., between 2000 < V_sys [km s⁻¹] < 5750) in order to mitigate the bias toward high H I-mass galaxies, which are detectable at greater distances.

The left panel in Fig. 5 shows a histogram of M_H I. Compared to the total α.100 catalog, our sample is probing the regime of lower H I masses, by design with our selection of dwarf galaxies. When considering the ALFALFA subsample of limited volume, our sample peaks at the same H I mass range, but does not show tails toward lower and higher H I masses present in the ALFALFA sample. The absence of the tail toward lower masses is expected as our selection procedure requires galaxies to be resolved, thereby excluding lower H I mass galaxies. The absence of the tail toward larger H I masses is due to the exclusion of high stellar mass galaxies in our selection procedure (by imposing the absolute magnitude cut, see Sect. 2.3) as well as to our (S/N)_ch and W₅₀ cuts, which exclude high mass edge-on galaxies (see Sect. 6.2). Furthermore, by considering only the M_H I range of our sample (8.50 < log (M_H I M_⊙) < 9.75) for ALFALFA, we plot a histogram of W₅₀ shown on the right hand side of Fig. 5. As expected, our sample peaks at lower W₅₀ than the corresponding ALFALFA samples due to the additional W₅₀ cut in our selection. Assuming all samples are randomly oriented, the large tails of the ALFALFA samples point toward the presence of larger total mass galaxies that we have excluded in our selection.

Fig. 5. Comparison of global H I properties between our and the ALFALFA sample. Left: Histogram of MH I normalized to the same maximum bin height with α. 100 sample in gray, the α. 100 volume limited subsample with 2000 < Vsys [km s−1] < 5750 in orange, and our sample in purple. Right: Histogram of log W50 normalized to the same maximum bin height for the H I mass range between 108.5 − 109.75 M⊙. Colors are the same as in the left panel.

The left panel of Fig. 6 shows the M_⋆ − M_H I relation where we compare our sample to the ALFALFA-SDSS catalog (Durbala et al. 2020). They report three different measurements of stellar masses for the sample. We choose the optically based method from Taylor et al. (2011) that (with the translation given in Eq. (3) from Durbala et al. 2020) was shown to be the most consistent with the Spectral Energy Distribution (SED) fitting (see Durbala et al. 2020 for more details). Our sample populates the same area in the graph as the volume limited ALFALFA sample, while the total ALFALFA sample tends to higher H I masses for the same stellar mass due to the bias described before. The right panel of Fig. 6 shows the median M_⋆ for two M_H I bins within 8.8 < log (M_H I/M_⊙) < 9.5. We have nine and eight galaxies in the lower and upper M_H I bin, respectively. As expected from our selection procedure, our sample has slightly smaller median stellar masses than the ALFALFA-SDSS sample for the highest M_H I bin due to our W₅₀ cut as well as the cut on absolute magnitude (a proxy for stellar mass).

Fig. 6. Comparison of H I and stellar masses between our and the ALFALFA sample. Left: MH I − M⋆ relation. Our sample is shown as purple circles, the total ALFALFA-SDSS sample is given by gray points and contours, and the volume limited subsample with 2000 < Vsys[km s−1] < 5750 is given by orange points and contours. Right: Median values of M⋆ for MH I bins between 8.8 < log (MH I/M⊙) < 9.5. Our sample is represented by purple circles, while the volume limited subsample from ALFALFA-SDSS is represented by orange squares.

4.2. Comparison of optical and H I geometries

We derived geometric parameters (x₀, y₀, i, and PA) for both the stellar (24 galaxies) and gaseous components (17 galaxies) independently. The ellipses describing disk geometries are visually compared in Fig. 4 on the left panel for AHCJ2239+3832, and in Appendix C for the rest of the sample.

The left panel of Fig. 7 shows the comparison between optical and H I inclinations for the galaxies in our sample. For ten galaxies, the optical and H I inclinations are compatible within the errors, while for the other seven galaxies, the discrepancy is significantly larger. Generally, the H I determined inclinations are systematically lower than optically determined ones. A possible explanation for this discrepancy could be the assumption of a finite thickness for the stellar disk when deriving the optical inclination, and a razor thin H I disk when deriving the H I inclination. However, even when assuming a razor thin stellar disk (shown as white circles in the plot), the discrepancy is still present, although mitigated.

Fig. 7. Comparison between H I and optical geometric parameters. Left: Comparison of optical and H I inclinations with purple circles corresponding to optical inclinations obtained by assuming intrinsic thickness of q0 = 0.3, and white circles corresponding to optical inclinations for a razor thin stellar disk. Middle: Difference in optical and H I inclinations with respect to the difference in the optical and H I morphological position angles. Right: Histogram of the cosine values of inclinations obtained from H I in blue and optical in orange. The number in the legend corresponds to the number of galaxies with well constrained inclinations. We note that all six galaxies in the last H I bin to the right show consistency with inclinations down to 0°.

The middle panel of Fig. 7 shows the difference in inclination versus the difference in morphological position angle. For seven galaxies, we observed good consistency between geometric parameters with ΔPA ≲ 25° and Δi ≤ 10°. For cases with large differences in position angle (≳40°), inclinations cannot be directly compared as optical and H I morphologies are misaligned. Therefore, it is not surprising to find large differences in inclinations in these cases. However, even with ΔPA ≲ 25°, we find significant differences in inclinations with four cases having Δi > 10°. We conclude that in these cases it is not straightforward to use one set of parameters for both the gas and stars. These comparisons, however, are subject to some caveats. As mentioned in Sect. 3.3, possible caveats in the isophotal fitting procedure include the influence of clumpy star-forming regions that might dominate over the underlying fainter disk component, as well as the finite depth of the optical images that might not allow us to trace the disk component far enough out for a robust comparison with H I morphology.

On the right panel of Fig. 7, we show a histogram for H I and optically determined inclinations. We point out that all six galaxies in the rightmost H I bin show consistency with inclinations down to 0°, which would potentially flatten the peak at these values. Both distributions seem to lack high inclinations. While H I source finding favors lower inclinations, the majority of this bias is likely a consequence of our selection procedure where we have a cut of W₅₀ < 150 km s⁻¹, and a cut of (S/N)_ch < 3 (see Sect. 6.2 for more details). However, we note that three galaxies (AGC 234932, AHCJ2218+4059, and UGC 12005 in Figs. C.10¹⁰, C.19¹⁰, and C.23¹⁰, respectively) might be edge on. This is seen from the isophotal fitting, which showed an inner region (compared to the one used to infer the global geometry of the system) with ellipticities around 0.7, corresponding to inclinations of 90° for q₀ = 0.3. In comparison, the obtained optical (H I) inclinations are 68° (53°), 63° (55°), and 57° (29°). These galaxies also show a broad H I-emission region in the PV slices along the major axis, which could be a signature of edge-on galaxies, but we note that these cases are not well spatially resolved, which could also produce broad PV profiles. If these galaxies were edge-on, the distribution of inclinations would be flatter than it appears.

4.3. Properties of UDGs in our sample

For the UDG classification, we adopt conditions from Sect. 1 (⟨μ⟩_e, X > 24 mag arcsec⁻², R_e, X > 1.5 kpc) applied to the r-band of PS1 due to better quality of the data compared to the g-band (i-band is not commonly used in the literature for this classification, see Sect. 1). According to this, our sample contains nine UDGs in total, but with four of them being optically unreliable (see Sect. 3.3). The five reliable UDGs are AHCJ2239+3832, UGC 8602, AGC 239039, AGC 239112, and AHCJ1359+3726 (Figs. 4, C.5¹⁰, C.9¹⁰, C.11¹⁰, and C.14¹⁰). They are all marked in Table 4.

Throughout the paper, we compare our UDG sample to the one from Mancera Piña et al. (2019b, 2020) (hereafter M20). The authors of M20 have used similar analysis to ours for both the optical photometry and the H I kinematic modeling. For one of their galaxies (AGC 114905), we use the results from Mancera Piña et al. (2022b), which have been obtained using higher quality H I data. This sample is particularly interesting to use for comparison, as it has been shown to systematically deviate from the BTFR.

In Fig. 8, we show the ⟨μ⟩_e, r − R_e, r relation. We transferred our ⟨μ⟩_e, r to SDSS band (using Eq. (6) in Tonry et al. 2012) in order to compare with the H I selected UDG sample from Leisman et al. (2017) shown as gray circles, and the sample from M20 shown as red hexagons. Leisman et al. (2017) had an additional constraint on the absolute magnitude M_r > −17.6, but all our UDGs (except one optically unreliable) also satisfy that criteria. We note that they do not report the band in which the R_e is measured. Almost all our UDGs populate the same area as UDGs from Leisman et al. (2017), with only one optically unreliable UDG having a significantly higher R_e, r. This outlier is AHCJ2207+4008 whose best fit Sérsic function had a Sérsic index of 2.5, unusually high for previously known UDGs. Compared to the M20 sample, our UDG sample has higher ⟨μ⟩_e, r on average, with optically reliable subsample populating similar area as the three brightest UDGs from M20.

Fig. 8. Mean effective surface brightness versus the effective radius in the r-band. Galaxies with well constrained optical parameters (see Sect. 3.3 for more details) are in color, with UDGs as cyan stars and standard dwarf galaxies as purple diamonds. The rest of the sample with less reliable parameters are shown as white markers, with UDGs as stars and standard dwarfs as diamonds. The H I selected UDG sample from Leisman et al. (2017) is plotted as gray circles, and the sample from M20 as red hexagons. The dash dotted line corresponds to the threshold of surface brightness used for the classification of UDGs.

As rotational velocity of a galaxy strongly depends on its inclination, with lower inclinations introducing larger uncertainties, we compare inclinations of our UDG subsample to inclinations of our standard dwarfs in order to better understand the precision of our derived rotational velocities between the two samples. For H I-derived inclinations, 5 out of 8 UDGs have inclinations lower than 40°, with four of these having inclinations consistent with zero. In comparison, 3 out of 9 standard dwarfs have H I inclination less than 40°, with only two consistent with zero. Taking the optical inclinations, 3 out of 4 galaxies with inclinations lower than 50° are labeled as UDGs. Therefore, UDGs in our sample seem to have lower inclinations than standard dwarfs, possibly leading to higher uncertainties in derived rotational velocities.

5. Scaling relations

5.1. Stellar mass–size relation

We look at the stellar mass–size relation to study the position of our H I-selected galaxies with respect to the optically selected ones. We make our comparison with the sample of late-type galaxies from Fernández Lorenzo et al. (2013). Their sample is part of the Analysis of the interstellar Medium of Isolated GAlaxies (AMIGA) project, specifically selected for high isolation. This ensures low-to-no contribution of environment on the intrinsic galaxy properties, making it a favorable sample for comparison with the H I-selected sample, which naturally contains more isolated systems.

In Fernández Lorenzo et al. (2013), they report effective radii defined in a circular aperture¹¹ (r_e), which we converted to elliptical apertures (R_e, defined along the semi-major axis, i.e., along the disk) using $R_e=r_e\times \sqrt{a/b}$, where a and b are semi-major and semi-minor axes.

The stellar mass–size plot is shown in Fig. 9. To compare where our sample lies with respect to the AMIGA sample, we fit the AMIGA sample with the functional form from Eq. (2) in Fernández Lorenzo et al. (2013). In contrast to their work, we leave all parameters free to fit in order to better follow the trend in their data for our comparison. Our sample is generally more diffuse (tends to larger R_e at the same stellar mass) than the AMIGA sample. This is likely a consequence of basing our selection on H I detections and thereby allowing lower surface brightness (but same stellar mass) galaxies into the selection, as well as of having a cut in absolute magnitude as part of our selection procedure, which could potentially exclude higher surface brightness galaxies of the same stellar mass. This is particularly interesting as the isolated AMIGA sample used for this comparison has been shown to be systematically more diffuse than an optically selected sample without a strict isolation criterion (Shen et al. 2003). Quantitatively, AMIGA galaxies are ∼1.2 times larger, for the same stellar mass.

Fig. 9. Stellar mass–size relation. Our sample is given with the same markers as in Fig. 8 and is plotted with Re values from the i-band. Galaxies from the AMIGA sample with i-band Re values are given in green, and the corresponding fit (see Sect. 5.1) is given in black. The H I-rich UDG sample from M20 is plotted with Re values from the r-band and is denoted as red hexagons. Histogram on the top shows the distribution of stellar masses, with our sample plotted in purple and the AMIGA in green. The histogram on the right shows distributions of effective radii for M⋆< 109 M⊙ in our and the AMIGA samples with the same color scheme.

We also plot UDGs from M20 in our stellar mass–size plot in Fig. 9, in order to compare with our UDG subsample. The M20 sample does not have R_e measurements in the i-band, so we plot their r-band R_e values. On average, our sample has higher stellar masses than M20, but similar R_e as their most massive UDGs. Generally, the M20 sample seems to populate more extreme regime of UDG population having lower surface brightness (see Sect. 4.3), lower stellar masses, and extending to larger effective radii.

5.2. Baryonic Tully-Fisher relation

Here we look at the BTFR connecting baryonic mass and rotational velocity. We want to explore if the UDGs in our sample show offsets from the relation when compared to the rest of the sample, as was indicated by previous works (e.g., Leisman et al. 2017; Mancera Piña et al. 2020). We compare our sample to various samples from the literature. From the SPARC sample (Lelli et al. 2016a), we took 126 galaxies with reliable rotation curves (quality flag Q = 1 or 2) and inclinations i > 25°. From these, we excluded three galaxies that are part of the LITTLE THINGS subsample of 17 galaxies from Iorio et al. (2017) that have more detailed analysis and use a similar approach to ours. Additionally, we add the SHIELD sample (McNichols et al. 2016) of 12 low-mass galaxies. We also compare to two resolved samples of H I-rich UDGs: six UDGs from M20, and 11 edge-on UDGs from He et al. (2019).

The BTFR is shown in Fig. 10. As mentioned in Sect. 3.2, we have used H I geometric parameters (except for the position angle that is derived from the kinematics) for the final kinematic model for all galaxies except AGC 239039. In this case, we used optical parameters due to the high misalignment of H I morphological and kinematic position angles, making the obtained H I (morphological) inclination non-applicable. A more detailed discussion on the inclination impact for our sample in the BTFR can be found at the end of this section.

Fig. 10. Baryonic Tully-Fisher relation comparing the placement of our sample with respect to various samples from the literature. Galaxies in our sample with well constrained rotational velocity (see Appendix C for more details) are denoted as in Fig. 8. We note that the placement on the BTFR is robust for all galaxies (both optically reliable and unreliable) as derived stellar mass is not significantly impacted by the shape of the Sérsic fit (see Sect. 3.3.1 for more details); only the UDG classification of the optically unreliable ones is less certain. For galaxies that show a flattening in the PV slice along the major axis, we provide errors on circular velocities, while others we plot as lower limits. Black dash dotted line is the best fit model to the SPARC sample from Lelli et al. (2016a).

All the standard dwarf galaxies in our sample lie on the relation defined by the SPARC sample, while the UDGs seem to be slightly offset. However, due to the low resolution of our sample, it is not straightforward to say if we are able to trace the flat part of the rotation curve for our galaxies. Lelli et al. (2019) explored the connection of the physical extent at which the rotation is measured and the reached rotational velocity, pointing out that velocities measured at 2R_e (and lower) might not always be a good tracer of the flat part of the rotation curve for low-mass galaxies, systematically underestimating the V_flat. When we apply this threshold to our sample, we find that 7 out of 10 galaxies in the BTFR are traced outside of 2R_e, three of which are UDGs (two of them are optically reliable). As mentioned in Sect. 3.2, we also independently selected four galaxies with potentially flat rotation curve based on the shape of the PVs, all of which also have rotation curves extending beyond 2R_e, and two of which are UDGs. We also note that three UDGs offset from the relation have inclinations consistent with 0°, which could move them to higher rotational velocities in the plot.

When compared to UDGs from M20, our UDGs have higher circular velocities for the same baryonic mass. As seen in Sect. 5.1, the M20 sample has on average lower stellar masses than our sample, leading to higher gas fractions for the same baryonic mass. We explored possible correlations between the increase in R_e and M_H I/M_⋆ fraction with the offset from the relation, but no evident trend is seen. However, we note that the M20 sample represents a more extreme UDG population on average, with the main difference between the two samples being the gas fraction. As high gas fractions of UDGs are likely driven by systematic differences in specific angular momentum of these systems (Mancera Piña et al. 2021), they might explain the relative difference in systematic offsets between the two samples, when taken on average. However, more statistically significant samples are needed to confirm this trend.

As shown in Sect. 4.2, our sample shows significant differences in H I and optical inclinations. To explore the impact of this difference on the galaxies’ position in the BTFR, we ran ^3DBarolo again, this time using optical inclinations for galaxies whose optical and H I morphological position angles do not differ by more than 45°. Figure 11 shows the difference of the two cases in the resulting BTFR. While standard dwarf galaxies in our sample are very close to the BTFR from literature in both cases, the UDGs tend to lie off the relation, but their precise locations vary significantly depending on the inclination used. Generally, optical inclinations give systematically lower rotational velocities, as expected from the left panel of Fig. 7.

Fig. 11. Baryonic Tully-Fisher relation where galaxies from our sample are plotted based on the inclination used for kinematic modeling. For both UDGs (stars) and standard dwarfs (diamonds), rotational velocities obtained using H I inclinations are given in blue, while the ones obtained using optical inclinations are given in orange. We note that here we do not make a distinction between the optically reliable vs. unreliable and/or rotationally flagged vs. unflagged galaxies, as long as their PAH I − PAop < 45°.

Taking only the optical or only the H I inclination at the low-mass scales can give a significantly different result in terms of the BTFR. This is especially important as most works up until now have used either the H I (e.g., Iorio et al. 2017; Mancera Piña et al. 2020) or optical (e.g., Karunakaran et al. 2020; McQuinn et al. 2022; Hu et al. 2023) geometry. Our sample demonstrates that at these low masses, the measurements of H I and optical morphologies can give inconsistent results (possibly due to the insufficient depth of optical images, see Sect. 3.3). Therefore, care should be taken when interpreting H I kinematics based on optical morphologies in the low surface brightness regime, and placing galaxies on the BTFR.

6. Discussion

6.1. Incidence of pairs

We found three pairs in our sample of 24 galaxies. This makes the percentage of galaxies in pairs to be 25% (6/24) in our sample. All of our pairs have stellar mass ratios larger than 0.3. In comparison, Sales et al. (2013) studied the frequency of dwarf satellite galaxies as a function of the primary galaxy stellar mass in the Sloan Digital Sky Survey (SDSS).

The fraction of galaxies in pairs for primary masses and absolute magnitude differences of pairs in our sample (up to one magnitude difference) would be ∼5% (see figure 3 from Sales et al. 2013).

Another work by Besla et al. (2018) compared the SDSS catalog with the Illustris-1 cosmological simulation (Vogelsberger et al. 2014; Nelson et al. 2015) and found a similar result (∼4%) for the primary stellar mass range of (2 × 10⁸ − 5 × 10⁹) M_⊙, but could potentially go up to ∼6% for future surveys with better completeness at these masses. As we are finding a significantly larger fraction of dwarf galaxy pairs compared to previous studies on optical spectroscopic data, we explore the potential biases of our selection procedure and the H I source finding technique. If galaxies in pairs were detected as individual objects in the initial source finding catalog from Apertif, there would be no bias in our procedure. However, two pairs in the sample were detected as a single source, possibly introducing the biases. We refer to these cases as close pairs, in contrast to pairs that were not recognized as a single source.

To test if our selection criteria have biased the selection of close pairs in our sample, we apply each criteria again on each galaxy in the two pairs. For consistency, we used the Hubble flow distance and a beam size taken as the geometric mean of the major and minor beam axes (as was used in the selection, see Sect. 2.3) in order to test the M_H I condition, and the condition that the galaxy is resolved by 3 beams. The results are summarized in Table 5. In both pairs, there is one galaxy that passes all criteria. Given that one of the galaxies of the pair always passes the selection criteria, we would have seen the companion galaxy in the data products either way and thus also included it in the sample. We conclude that our selection did not bias us toward finding more pairs than we would obtain if each galaxy was selected as a single source.

On the other hand, the source finding itself could introduce additional biases. We note that the Apertif survey is surface brightness limited for resolved sources, but flux limited for unresolved ones (hence also H I mass limited, see Eq. (7)). As galaxies in our sample are by construction resolved by three or more Apertif beams, they would all be resolved with the largest beam used in the source finding (39″ beam, see Sect. 2.1), thereby putting us in the surface brightness limited regime. Thus there should be no bias in the detection of a galaxy separately or as part of a close pair.

We conclude that our procedure cannot explain the high frequency of dwarf pairs in the sample. However, our sample is small and a more indepth study using a larger sky area from Apertif or other wide-field resolved H I surveys such as WALLABY (Koribalski et al. 2020), is advisable for the robust quantification of the incidence of gas-rich dwarf pairs.

6.2. Impact of our H I-based selection

Basing our selection on the H I properties of galaxies, we are subject to some biases. As seen in Sect. 4.2, we are missing high inclination galaxies in our sample. This is partly a consequence of the H I survey, which is biased toward lower inclinations. The bias comes from the galaxy emission being distributed over a smaller number of spectral channels (increasing the S/N in each individual channel) than for higher inclination galaxies. Additionally, (S/N)_ch threshold in the selection criteria enhances this bias toward lower inclinations for the same reason, as well as the imposed threshold on the W₅₀. We tested the effects of these selection criteria on our sample by exploring how the (S/N)_ch and W₅₀ values obtained from the preliminary Apertif source list would change for a range of inclinations (taking the H I measured inclination as the true one). Only ∼30% of the tested galaxies (17 galaxies, those that have H I inclinations measured) would have been selected for the whole range of inclinations. From the remaining galaxies, ∼30% would be selected only for inclinations ≲40°. We note that in 16/17 cases, the (S/N)_ch criterion was the more stringent one, thereby far more dominant for introducing biases than the imposed W₅₀ cut. This analysis demonstrates that we have a significant bias toward low inclinations in our sample, but a proper quantification is beyond the scope of this work.

Furthermore, as shown in Sect. 5.1, our sample is shifted toward higher effective radii in the stellar mass–size relation compared to the isolated AMIGA sample. This could also be a consequence of the H I-based selection, as it is by construction optically independent and allows us to pick galaxies of lower surface brightness for the same stellar mass. This effect is especially important for galaxies with lower inclinations present in our sample as they appear less bright when projected on the sky and hence have a lower probability of being detected in optical surveys.

We compare the number of UDGs in our sample to expectations based on the UDG H I mass function (HIMF) of Jones et al. (2018). They constructed the HIMF and the velocity width function (made using W₅₀) of UDGs (selected by Leisman et al. 2017) from the α.70 ALFALFA H I untargeted survey catalog (Haynes et al. 2011). They report that the fraction of UDGs to total number of galaxies peaks at log M_H I/M_⊙ ∼ 8.8 with a contribution of 6%, and lowers toward both higher and lower H I masses. In comparison to our optically reliable subsample, we have five UDGs, making ∼30%. However, Jones et al. (2018) do not make a distinction between large stellar mass galaxies with lower gas fractions and dwarf galaxies with lower masses in both stars and gas. Therefore, their comparison for UDG frequency is against a large range of stellar masses and hence total galaxy masses, bringing the final UDG fractions lower. In addition, the contribution of UDGs is seen to steeply rise (more steeply than the total velocity width function) toward the narrower W₅₀ profiles. This goes in accordance with possibly lower rotational velocities as well as low H I inclinations of UDGs compared to standard dwarfs. Given that we have put an upper limit on W₅₀ in our selection, we may have introduced some preference for selecting UDGs in our sample.

7. Conclusions

In this work, we have presented a sample of 24 dwarf galaxies with spatially resolved H I data from the Apertif imaging survey (Adams et al. 2022). Our selection procedure ensured that galaxies are resolved by at least three beams in H I, that the average signal-to-noise ratio per channel in the global profile ((S/N)_ch) is larger than 3, and that their H I and stellar masses are M_H I < 10¹⁰ M_⊙ and M_⋆ ≲ 10⁹ M_⊙, respectively. We measured the geometry of the H I disk for 17 galaxies (out of 24) in the sample and successfully produced kinematic models of 13 of them using ^3DBarolo. We also studied properties of their stellar disks (24 galaxies) by conducting isophotal fitting and extracting surface brightness profiles from g-, r-, and i-bands of the Pan-STARRS 1 photometric survey. We identified five UDGs and four candidate UDGs in the sample. Using the above results, we position our sample with respect to other samples from the literature on the stellar mass–size relation and the BTFR.

In the following, we summarize our main conclusions:

• Our H I-selected sample seems to be more diffuse (has larger R_e, X at given stellar mass) than the optically selected sample of isolated galaxies from the AMIGA project (Fig. 9). This shows that an H I-based selection returns a different population from an optical selection and hence should be taken into account in statistical studies of dwarf galaxies.

• We find apparent misalignments between the derived optical and H I morphologies for nine (out of the 17 measured) galaxies in our sample. However, this comparison is subject to caveats originating from the influence of clumpy star-forming regions on the derived stellar morphology as well as the insufficient depth of PS1 images, which might not trace the underlying stellar disk far enough for a robust comparison with the H I morphology.

• Standard dwarf galaxies in our sample lie on the BTFR determined for higher mass galaxies, but the UDGs in our sample seem to be slightly shifted toward lower rotational velocities (Fig. 10). However, we are limited by a small number of galaxies in our sample, and in most cases, we cannot guarantee that the measured rotational velocities are fully representative of the flat part of the rotation curves.

• Inclination has a large impact on the measured rotational velocity and consequently on the position of a galaxy in the BTFR (Fig. 11). Given the significant misalignments between the measured H I and stellar disk morphologies, care should be taken when deciding which geometric parameters to use for the derivation of the kinematics.

• We find a larger fraction (25%) of dwarf galaxies in pairs than studies based on optical spectroscopic data. We considered the possible biases in our selection procedure and the H I source finding, and find no biases that could account for this. Our findings suggest that H I surveys might be detecting more dwarf galaxy pairs than are found in optical spectroscopic surveys. This can have relevant implications for the formation and evolution of these low-mass galaxies and will be further explored in an upcoming publication.

The above conclusions are based on an early H I source list from Apertif data, which is limited in both sensitivity and coverage. Future work is currently undertaking source finding on more sensitive data and over a larger survey area. This will enable the detection of many more sources that pass our selection criteria and could put stronger constraints on our final conclusions. A more statistical study of the fraction of dwarf pairs and/or UDGs would also benefit from dropping some of our current selection criteria, such as spatially resolving the H I disk or the (S/N)_ch > 3. Furthermore, new upcoming facilities, such as the Square Kilometre Array (SKA), will enable the detection of an even¹² higher number of dwarf galaxies and with superior spatial and spectral resolution. The present study and future follow-ups will clarify the complexity of the present-day dwarf galaxy population, permitting key constraints of theoretical models of galaxy formation and evolution at the low-mass end to be obtained.

Acknowledgments

We thank the anonymous referee for a valuable report that improved this manuscript. We thank Sabrina Stierwalt for pointing us toward previous studies on dwarf galaxy pairs. We also greatly thank Jackson Fuson and John M. Cannon for sharing their asymmetric drift correction code with us. KMH, JG, SSE, and MPR acknowledge financial support from the grant CEX2021-001131-S funded by MCIN/AEI/ 10.13039/501100011033 from the coordination of the participation in SKA-SPAIN funded by the Ministry of Science and Innovation (MCIN); and from grant PID2021-123930OB-C21 funded by MCIN/AEI/ 10.13039/501100011033 by “ERDF A way of making Europe” and by the “European Union”. Authors KMH, BŠ, JG, SSE, and MPR acknowledge the Spanish Prototype of an SRC (SPSRC) service and support funded by the Ministerio de Ciencia, Innovación y Universidades (MICIU), by the Junta de Andalucía, by the European Regional Development Funds (ERDF) and by the European Union NextGenerationEU/PRTR. The SPSRC acknowledges financial support from the Agencia Estatal de Investigación (AEI) through the “Center of Excellence Severo Ochoa” award to the Instituto de Astrofísica de Andaluc ía (IAA-CSIC) (SEV-2017-0709) and from the grant CEX2021-001131-S funded by MICIU/AEI/10.13039/501100011033. KMH and JMvdH acknowledge funding from the ERC under the European Union’s Seventh Framework Programme (FP/2007–2013)/ERC Grant Agreement No. 291531 (‘HIStoryNU’). AAP acknowledges support of the STFC consolidated grant [ST/S000488/1]. PEMP acknowledges the support from the Dutch Research Council (NWO) through the Veni grant VI.Veni.222.364. This work makes use of data from the Apertif system installed at the Westerbork Synthesis Radio Telescope owned by ASTRON. ASTRON, the Netherlands Institute for Radio Astronomy, is an institute of the Dutch Research Council (“De Nederlandse Organisatie voor Wetenschappelijk Onderzoek”, NWO). The Pan-STARRS1 Surveys (PS1) have been made possible through contributions of the Institute for Astronomy, the University of Hawaii, the Pan-STARRS Project Office, the Max-Planck Society and its participating institutes, the Max Planck Institute for Astronomy, Heidelberg and the Max Planck Institute for Extraterrestrial Physics, Garching, The Johns Hopkins University, Durham University, the University of Edinburgh, Queen’s University Belfast, the Harvard-Smithsonian Center for Astrophysics, the Las Cumbres Observatory Global Telescope Network Incorporated, the National Central University of Taiwan, the Space Telescope Science Institute, the National Aeronautics and Space Administration under Grant No. NNX08AR22G issued through the Planetary Science Division of the NASA Science Mission Directorate, the National Science Foundation under Grant No. AST-1238877, the University of Maryland, and Eotvos Lorand University (ELTE). This work made use of Astropy: a community-developed core Python package and an ecosystem of tools and resources for astronomy (Astropy Collaboration 2013, 2018, 2022). This research made use of Photutils, an Astropy package for detection and photometry of astronomical sources (Bradley et al. 2023). This research has made use of the NASA/IPAC Extragalactic Database (NED), which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration.

References

Appendix A: ^3DBarolo parameter file example

Here we give an example of the parameter file for kinematic modeling used in this work. For additional information, please visit the ^3DBarolo documentation.¹³

FITSFILE...: data cube name

3DFIT True: perform 3D fitting

NRADII...: number of rings in the model

RADSEP...: separation between rings in arcseconds

VSYS...: systemic velocity in kilometers per second

XPOS...: X axis position of the center of the galaxy in pixels

YPOS...: Y axis position of the center of the galaxy in pixels

INC...: inclination

PA...: kinematic position angle

Z0 0: thickness of the disk in arcseconds

DISTANCE...: distance to the galaxy in megaparsecs

FREE VROT VDISP VSYS PA: parameters for fitting

NORM AZIM: azimuthal normalization of the model cube

MASK SEARCH: method for producing a 3D mask

THRESHVELOCITY 1: keyword for masking, which ensures not taking into account emissions with two or more channels apart

FLAGROBUSTSTATS True: use median and MADFM for calculating cube statistics

SNRCUT 3.3: initial signal-to-noise ratio cut for producing a mask

FLAGGROWTH True: enable enlarging the initial mask

GROWTHCUT 2.5: enlarge the mask until the S/N reaches 2.5

BWEIGHT 0: not penalizing models extending outside the mask

WFUNC 2: choosing cos²θ as a weighing function

LINEAR ...: 0.5 times spectral resolution, in units of channels

NOISE... : rms noise per channel in the cube

TWOSTAGE True: enable two-stage fitting

POLYN 0: fit functional form to selected parameters across the rings, in this case VSYS and PA are fitted with a constant value

FLAGERRORS True: calculate and report errors on rotational velocity and velocity dispersion

MINVDISP...: minimal value of velocity dispersion, in our case half of the spectral resolution of the cube

SIDE B: fit both receding and approaching sides of the galaxy

Appendix B: Truncated versus interpolated photometry

In this work, we have obtained optical parameters (R_e, X, m_X and ⟨μ⟩_e, X) from the fits of the extracted surface brightness profiles for which we used the Sérsic function defined as

$$\begin{array}{c}\hfill {\mu }_X(r)={\mu }_{e,X}+\frac{2.5b}{ln(10)}[{\left(\frac{r}{R_{e,X}}\right)}^{1/n_X}-1],b=2n_X-\frac{1}{3}+\frac{4}{405n_X},\end{array}$$(B.1)

where μ_e, X is surface brightness at R_e, X, and n_X the Sérsic index. The m_X and ⟨μ⟩_e, X are then obtained using equations (see e.g., Graham & Driver 2005)

$$\begin{array}{c}\hfill {⟨\mu ⟩}_{e,X}={\mu }_{e,X}-2.5log\left(\frac{n_X{e}^b}{b^{2n_X}}\mathrm{\Gamma }(2n_X)\right),\end{array}$$(B.2)

$$\begin{array}{c}\hfill m_X={⟨\mu ⟩}_{e,X}-2.5log(2\pi {(R_{e,X})}^2),\end{array}$$(B.3)

where Γ(x) is the gamma function. R_e, as defined by the Sérsic function, is theoretically the same as half-light radius. However, the Sérsic function extends to infinity in r, while galaxies do not. Because of this, Eq. B.3 could be overestimating m_X and thus also overestimating R_e, X, or rather R_e^fit > R_e^meas, where R_e^meas is the half-light radius.

On the other hand, we can directly measure m_X and R_e, X from the profile. To do this, we sum the intensity in all rings that were used in the fitting of the profiles to get the total flux and convert it to m_X using Eq. 6. We provide these measured values in Table B.1. The comparison of values obtained by the fit and those from direct measurement can be found in Fig. B.1. Similarly, we find the half-light radii by taking the radius containing half of the total flux. The comparison of obtained measured and fitted R_e, X can be found in Fig. B.2. We computed ⟨μ⟩_e, X by summing the flux inside the measured R_e, converting it to a magnitude and using Eq. B.3 in reverse. Finally, we compare the obtained stellar masses in Fig. B.3, computed using either the fitted magnitude, or the measured magnitude in the g− band (the color is always computed in the same way). The obtained masses are consistent within the errors for all galaxies, except UGC 8605 and AHCJ2207+4008 for which they are compatible within 2σ.

Fig. B.1. Comparison of apparent magnitudes in different bands for values obtained from the fit of the Sérsic profiles and the directly measured ones.

Fig. B.2. Difference of effective radii obtained from the fit of Sérsic profiles and the directly measured ones for the r−band, given in percentages. The points are color coded using the Sérsic index, with optically reliable galaxies represented as circles and optically unreliable as squares.

Fig. B.3. Comparison of stellar masses calculated using either the magnitude obtained from the fitted Sérsic function or the directly measured magnitude from the extracted surface brightness profile (corresponding to $M_{\star }^{\mathrm{meas}}$).

Directly measured magnitudes are systematically fainter than ones obtained by the fit. This is not surprising as m_X^fit are obtained by integration of the Sérsic function to infinity, as discussed before. Similarly, measured R_e, X are systematically smaller than fitted values with differences being more pronounced for Sersic indices above 1.

Appendix C: Notes on individual galaxies

In this section we provide detailed comments for each galaxy considering both gas and stellar properties.

AHCJ2239+3832. The galaxy is shown in Figs. 2, 3, and 4. The H I disk shows an apparent decrease in inclination with radius, as seen from the total intensity map in Fig. 2. Therefore, the output from CANNUBI showed the galaxy being consistent with a wide range of inclinations down to the inclination of 0°. The kinematics are regular with no apparent signs of non-circular motions. The stellar counterpart is irregular with a slight excess of emission toward the west compared to the east side of the galaxy in the outer parts. Interestingly, fitted isophotes seem to align more with the H I geometry as they go further out of the galaxy, both by the shifting of the PA_op and by slightly decreasing ellipticity.

AHCJ0203+3714-1 and AHCJ0203+3714-2. These two galaxies are part of a pair shown on the uppermost panel of Fig. C.1¹⁴ with AHCJ0203+3714-1 being eastward of the two and AHCJ0203+3714-2 westward. The two galaxies are connected by a gas bridge as seen from the 3D cube, and consequently on the total intensity and velocity maps. The internal kinematics of both galaxies seem to be synchronized with their common motion as seen by a continuous velocity gradient in the velocity map. The stellar counterpart of AHCJ0203+3714-1, shown in better detail in Fig. C.2¹⁴, is partly obscured by a nearby bright star, which might have slightly influenced the isophotal fitting of the galaxy and the extracted surface brightness profile. On the other hand, the stellar counterpart of AHCJ0203+3714-2, shown in Fig. C.3¹⁴, shows variation in geometric parameters with radius. However, the change is small in our region of interest and the resulting geometry seems to well represent the orientation of the disk when compared to the orientation of H I (Fig. C.1¹⁴).

UGC 8605 and UGC 8602. These two galaxies also form a pair as is shown in the middle panel of Fig. C.1¹⁴ with UGC 8605 being the eastward and UGC 8602 the westward galaxy. These galaxies do not show a detectable gas bridge and are rotating in opposite directions as seen from the velocity map. For this pair, we were able to produce independent kinematic models. The obtained H I geometric parameters of UGC 8605, seen in the first two panels of Fig. C.4¹⁴, seems to well describe the H I morphology. However, CANNUBI was not able to fully constrain the inclination of the disk showing consistencies down to the inclination of 0°. The kinematics of the galaxy seems regular and the kinematic model well reproduces the data. The stellar counterpart shows complex morphology with the three apparent stellar arms going out of the inner disk toward the southeast, northwest, and northeast. This features influenced the isophotal fitting as seen from the large change in the position angle and the increase of ellipticity with radius. As in the case of AHCJ2239+3832, the fitted isophotes seem to align more with the H I geometry as they extend further to the ourskirts of the galaxy. The extracted surface brightness profile could not be robustly fitted with a Sérsic function due its irregular shape. In the case of UGC 8602, shown in Fig. C.5¹⁴, the H I geometric parameters could also not be well constrained, showing consistency with inclinations down to 0° and thereby making the morphological position angle less constrained. The galaxy also seems to be kinematically lopsided as seen from its asymmetric PV diagram along the major axis. This means that the derived rotational velocity might not be a reliable tracer of the total underlying potential. Its stellar counterpart has irregular morphology seen by the large change in ellipticity and position angle with radius.

CGCG 290-011 and UGC 5480. These galaxies make the third pair in our sample and are shown in the lowermost panel on Fig. C.1¹⁴ with CGCG 290-011 being the eastward and UGC 5480 the westward galaxy. The two galaxies are not directly connected in the 3D cube, due to a separation along the spectral axis. UGC 5480 is the larger of the two galaxies and shows a clear velocity gradient in the velocity map. The stellar counterpart of CGCG 290-011, shown in Fig. C.6¹⁴, is close to a nearby star, which might have slightly influenced the isophotal fitting and the extracted surface brightness profile. The stellar counterpart of UGC 5480, shown in Fig. C.7¹⁴, shows a drop in ellipticity in the outskirts, but an almost constant ellipticity in the rest of the galaxy.

UGC 8503. The galaxy is shown in Fig. C.8¹⁴. This galaxy is positioned at the edges of two Apertif observed fields. As the fields are not being mosaicked together at the present moment, the galaxy was detected two times in the new source finding and was assigned two names. One of them reported in Table 2, and the second being AHC133045.2+324548. Due to being at the edge of the primary beam, the H I data quality is bad and unreliable for this galaxy. Additionally, in the i−band of PS1, there is a drop in background counts eastward from the galaxy that was masked for the optical analysis. The stellar counterpart shows a slight twist in position angle, but an almost constant ellipticity with radius. There is an indication of a faint emission beyond the fitted ellipses, which would probably result in a decrease in ellipticity at larger radii.

AGC 239039. The galaxy is shown in Fig. C.9¹⁴. The kinematics shows a twist along the minor axis (as seen from its velocity field and from the velocity gradient in the PV slice along the minor axis) and an asymmetric PV slice along the major axis. The twist along the minor axis could be a sign of non-negligible radial motions. The asymmetry in the PV slice is caused by higher velocity dispersion in the approaching side of the galaxy indicating higher disturbance, but with no effect on the underlying rotation. There is also a large offset between the morphological and kinematic position angles of the H I disk. In this case, CANNUBI was unable to constrain the inclination, showing consistency with inclinations down to 0°. Stellar morphology shows a twist in position angle as a function of radius.

The inner region is almost perpendicular to the northeast-southwest orientation of the outskirts. However, the kinematic position angle from H I is also oriented in the northeast-southwest direction indicating the orientation of the disk, which is why we adopt optical geometric parameters for the kinematic model of this galaxy.

AGC 234932. As shown in Fig. C.10¹⁴, the galaxy shows signs of disturbance in the disk seen from the asymmetry in PV slices in the form of a seemingly separate peak in emission in the receding part. Interestingly, the stellar disk displays a peculiar curvature in the same part of the disk. These features indicate that the galaxy might be a merger in a late stage. Even with these features, the galaxy seems morphologically regular and geometry was successfully constrained in both the optical and H I. In the inner region of the optical counterpart, the ellipticity goes up to 0.7, which would be equal to the inclination of 90° for q₀ = 0.3. There also seems to be a broad H I emission region in the PV slice along the major axis, which could be another signature of edge-on galaxies, but we note that the galaxy is not spatially well enough resolved to distinguish this from the effect of the beam.

AGC 239112. The galaxy is shown in Fig. C.11¹⁴. The H I total intensity map shows a very irregular morphology so the geometry of the H I disk could not be constrained. The twist along the minor axis seen in the velocity field could be a sign of non-negligible radial motions. The stellar counterpart shows a twist in the position angle with radius, but this time toward higher misalignment with the kinematics of the H I.

UGC 5541. The galaxy is shown in Fig. C.12¹⁴. The H I geometry was well constrained and galaxy shows a well ordered rotation. However, the galaxy has an H I tail extending toward the south indicating a possible interaction. Going 250 kpc eastward, there is the NGC 3182 galaxy with stellar mass of 10^10.3 M_⊙ (Pak et al. 2023), which could potentially be interacting with UGC 5541. The optical counterpart shows high disturbance with many bright star-forming regions in the outskirts, another possible signature of an interaction. Some of these bright regions were masked for the isophotal fitting, but the obtained position angle was still slightly influenced by the region in the northeast direction. The surface brightness profiles could not be robustly fitted with a Sérsic function for this galaxy due to its irregular shape.

AHCJ1308+5437. The galaxy is shown in Fig. C.13¹⁴. Due to low S/N, the cube was smoothed 1.5 times for the estimation of H I geometry and the kinematic modeling. The galaxy shows regular rotation. The PS1 data in the i-band for this galaxy has a region with very low background counts in the northeast corner that was masked for the optical analysis. The stellar counterpart shows a large star-forming region almost completely outside of the inner disk that was partially masked for the isophotal fitting but unmasked for the extraction of the surface brightness profile. The obtained geometry of the inner disk was well constrained, but the derived surface brightness profile cannot be well fitted with the Sérsic function.

AHCJ1359+3726. The galaxy is shown in Fig. C.14¹⁴. Again, due to low S/N, the cube was smoothed 1.5 times for the estimation of H I geometry. The H I kinematics are very irregular as seen in the PV slices. Therefore, we could not produce reliable kinematic model for this galaxy. Optical counterpart shows significant variations of geometric parameters with radius. Specifically, the ellipticity grows from the central part and falls off at the outer part. The obtained optical geometry is in accordance with the kinematic position angle from H I, but not with H I geometry.

AHCJ2207+4008. The galaxy is shown in Fig. C.15¹⁴. The H I geometry and kinematics of the galaxy are well constrained. Its stellar counterpart shows a variation in position angle and ellipticity with radius. Interestingly, the evolution of position angle with radius is moving toward the position angle seen in H I.

AHCJ2207+4143. The galaxy is shown in Fig. C.16¹⁴. Its kinematics shows a very well ordered disk, but a significant non-orthogonality of major axis with its systemic velocity contour seen from the velocity field map. This also causes the velocity gradient in the PV along the minor axis and is a likely consequence of radial motions of gas. Unfortunately, the galaxy is positioned close to Galactic cirrus, making the optical analysis less reliable. Otherwise, the galaxy seems to be regular.

AHCJ2232+3938. The galaxy is shown in Fig. C.17¹⁴. The H I geometry was well constrained and the H I kinematics shows a well ordered disk. The optical counterpart shows little variation of geometric parameters with radius and is well constrained.

AHCJ2216+4024. The galaxy is shown in Fig. C.18¹⁴. The H I geometry was well constrained and the kinematics seems regular. The optical counterpart shows some variation in position angle with radius, but it is not unexpected due to low inclination.

AHCJ2218+4059. The galaxy is shown in Fig. C.19¹⁴. The geometry of the H I disk was well constrained. However, the galaxy shows a warp in the position angle seen as a twist in the velocity field, which is also apparent in the H I morphology. Its optical counterpart also has a well constrained morphology and little variation of geometric parameters with radius. We note that in the very inner part, the ellipticity is around 0.7, corresponding to the inclination of 90° for q₀ = 0.3. As before, there also seems to be a broad H I emission region in the PV slice along the major axis, which could be another signature of edge-on galaxies. Again, we note that the galaxies in this sample are not spatially well enough resolved to distinguish this broadening from the effect of the beam.

UGC 12027. The galaxy is shown in Fig. C.20¹⁴. Due to low S/N, the H I cube of this galaxy was smoothed 1.5 times for the estimation of H I geometry and the kinematic modeling. The galaxy shows a drop in rotational velocity at larger radii, which is also apparent in the velocity map. This might be a sign of a warp in inclination and possibly the position angle. The change in inclination (or ellipticity) also seems to be present in its optical counterpart, with higher inclination in the inner part and a drop in inclination in the outskirts.

UGC 8363. The galaxy is shown in Fig. C.21¹⁴. The H I geometry was well constrained and the kinematics shows a well ordered rotating disk. The optical counterpart shows a change in ellipticity from a higher constant value in the inner part to a slightly lower constant value in the outer parts. The change is minimal in our region of interest for the estimation of geometric parameters and is in accordance with the H I geometry.

UGCA 363. The galaxy is shown in Fig. C.22¹⁴. Due to low S/N, the H I cube was smoothed 1.5 times for the estimation of H I geometry. However, kinematic modeling was not possible for this galaxy due to its high disturbance and the apparent complete misalignment between the morphological and kinematic position angle in H I as well as the misalignment between the optical position angle and the kinematic one. The stellar counterpart shows signs of three stellar arms toward the northeast, south and southwest directions. The obtained position angle is in accordance with the orientation of the inner disk.

UGC 12005. The galaxy is shown in Fig. C.23¹⁴. The geometry of the H I disk was well constrained with CANNUBI. Interestingly, the obtained morphological position angle is almost perpendicular to the kinematic one (70° difference). Because of this and the fact that the obtained optical morphology is also misaligned with the kinematic position angle, the kinematic modeling was not possible for this galaxy. Additionally, the PV along the major axis shows the galaxy is kinematically lopsided meaning its rotational velocity might not be a good tracer of its underlying potential. The stellar counterpart looks disturbed with a high number of star-forming regions in the outskirts. These regions cause the drop of ellipticity in the outermost part of the disk. In the inner part, the galaxy has ellipticity of around 0.7, corresponding to inclination of 90° for q₀ = 0.3. As before, there is an indication of a broad H I emission region in the PV slice along the major axis, which could be another signature of edge-on galaxies, but is indistinguishable from the effect of the beam due to having too few resolution elements.

AHCJ2249+3948. The galaxy is shown in Fig. C.24¹⁴. The H I geometry of the disk could not be constrained for this galaxy due to its peculiar morphology with a significant extension of H I toward the south. Additionally, the velocity field map shows kinematic disturbances along the disk that could be a sign of interaction with another galaxy. However, there is no recorded galaxy with a redshift measurement in NED within 600 kpc from AHCJ2249+3948 that could confirm this assumption. The stellar counterpart shows two stellar arms in the northeast and southwest direction causing the variation of position angle with radius. The H I kinematics and optical morphology in the outer parts are completely misaligned.

All Tables

All Figures

Fig. 1. Comparison between channel maps of AHCJ2239+3832 and its best-fit model. Top panels and blue contours represent the data, while the lower panels and red contours represent the best fit model. Contours are plotted starting from three times the noise per channel in the cube reported in Table 1, and are spaced by a factor of two in intensity. The black X indicates the center of the galaxy, as determined by CANNUBI (see Sect. 3.2).

In the text

Fig. 2. H I kinematics of AHCJ2239+3832. Leftmost: Color-composite image from PS1 g-, r-, and i-bands overlaid with H I contours. Contours are set to start at the level of intensity corresponding to a pseudo 4σ contour (see Sect. 3.2 for explanation) in the total intensity map (white), and grow by a factor of two in intensity toward the redder colors. Lowest contour for AHCJ2239+3832 corresponds to the column density of 7.2 × 1019 cm−2. The overplotted ellipse in blue represents the median geometric parameters obtained from CANNUBI (including the size). Middle left: Velocity field obtained as weighted-mean value from the H I data from 3DBarolo. The black cross represents the best-fit center position and the gray dashed line the kinematic position angle. Contours are given in spacing of 10 km s−1 with the green contour indicating Vsys. Middle right: Position-velocity slice along the major axis with data (blue) and model (red) contours. White points represent the obtained projected rotational velocity of the best fit model. Contours are plotted starting from two times the noise per channel in the cube as reported in Table 1, and grow linearly. Rightmost: Position-velocity slice along the minor axis (perpendicular to the dashed line in the middle left panel) with the same color scheme and contour levels as the middle right panel.

In the text

Fig. 3. Isophotal fitting result for AHCJ2239+3832. Upper panels: The leftmost panel shows the smoothed data overlaid with semi-equidistant ellipses along the semi-major axis length whose parameters are marked with red stars in the lower panels. White regions denote masked pixels in the image. The middle panel shows the model built from all fitted ellipses (blue circles in the lower panels), and the rightmost panel shows the residual of the data image and the model image. Lower panels: Ellipticity and position angle as a function of the semi-major axis length from the second run of the fitting algorithm (see Sect. 3.3) are shown as blue circles, with respect to radius from the fixed center position. Red stars denote parameters of semi-equidistant ellipses plotted in the upper leftmost panel. The black vertical line is located at a distance from the center that corresponds to the FWHM of the PSF after the smoothing, and the gray shaded region corresponds to the region in which the final global geometry was measured (using only the semi-equidistant ellipses, see Sect. 3.3).

In the text

Fig. 4. Photometry in the i-band of AHCJ2239+3832. Left: Optical i-band image overlaid with measured geometries from H I morphology (blue) and from i-band isophotal fitting (orange). The radius at which the optical geometry is plotted corresponds to the outermost data point in the surface brightness profile. The radius at which the H I geometry is plotted corresponds to the extent of the disk as obtained from CANNUBI. The crosses denote the obtained optical (orange) and H I (blue) centers. Right: Surface brightness profile (orange) and the corresponding best fit with a Sérsic profile (black).

In the text

Fig. 5. Comparison of global H I properties between our and the ALFALFA sample. Left: Histogram of MH I normalized to the same maximum bin height with α. 100 sample in gray, the α. 100 volume limited subsample with 2000 < Vsys [km s−1] < 5750 in orange, and our sample in purple. Right: Histogram of log W50 normalized to the same maximum bin height for the H I mass range between 108.5 − 109.75 M⊙. Colors are the same as in the left panel.

In the text

Fig. 6. Comparison of H I and stellar masses between our and the ALFALFA sample. Left: MH I − M⋆ relation. Our sample is shown as purple circles, the total ALFALFA-SDSS sample is given by gray points and contours, and the volume limited subsample with 2000 < Vsys[km s−1] < 5750 is given by orange points and contours. Right: Median values of M⋆ for MH I bins between 8.8 < log (MH I/M⊙) < 9.5. Our sample is represented by purple circles, while the volume limited subsample from ALFALFA-SDSS is represented by orange squares.

In the text

Fig. 7. Comparison between H I and optical geometric parameters. Left: Comparison of optical and H I inclinations with purple circles corresponding to optical inclinations obtained by assuming intrinsic thickness of q0 = 0.3, and white circles corresponding to optical inclinations for a razor thin stellar disk. Middle: Difference in optical and H I inclinations with respect to the difference in the optical and H I morphological position angles. Right: Histogram of the cosine values of inclinations obtained from H I in blue and optical in orange. The number in the legend corresponds to the number of galaxies with well constrained inclinations. We note that all six galaxies in the last H I bin to the right show consistency with inclinations down to 0°.

In the text

Fig. 8. Mean effective surface brightness versus the effective radius in the r-band. Galaxies with well constrained optical parameters (see Sect. 3.3 for more details) are in color, with UDGs as cyan stars and standard dwarf galaxies as purple diamonds. The rest of the sample with less reliable parameters are shown as white markers, with UDGs as stars and standard dwarfs as diamonds. The H I selected UDG sample from Leisman et al. (2017) is plotted as gray circles, and the sample from M20 as red hexagons. The dash dotted line corresponds to the threshold of surface brightness used for the classification of UDGs.

In the text

Fig. 9. Stellar mass–size relation. Our sample is given with the same markers as in Fig. 8 and is plotted with Re values from the i-band. Galaxies from the AMIGA sample with i-band Re values are given in green, and the corresponding fit (see Sect. 5.1) is given in black. The H I-rich UDG sample from M20 is plotted with Re values from the r-band and is denoted as red hexagons. Histogram on the top shows the distribution of stellar masses, with our sample plotted in purple and the AMIGA in green. The histogram on the right shows distributions of effective radii for M⋆< 109 M⊙ in our and the AMIGA samples with the same color scheme.

In the text

Fig. 10. Baryonic Tully-Fisher relation comparing the placement of our sample with respect to various samples from the literature. Galaxies in our sample with well constrained rotational velocity (see Appendix C for more details) are denoted as in Fig. 8. We note that the placement on the BTFR is robust for all galaxies (both optically reliable and unreliable) as derived stellar mass is not significantly impacted by the shape of the Sérsic fit (see Sect. 3.3.1 for more details); only the UDG classification of the optically unreliable ones is less certain. For galaxies that show a flattening in the PV slice along the major axis, we provide errors on circular velocities, while others we plot as lower limits. Black dash dotted line is the best fit model to the SPARC sample from Lelli et al. (2016a).

In the text

Fig. 11. Baryonic Tully-Fisher relation where galaxies from our sample are plotted based on the inclination used for kinematic modeling. For both UDGs (stars) and standard dwarfs (diamonds), rotational velocities obtained using H I inclinations are given in blue, while the ones obtained using optical inclinations are given in orange. We note that here we do not make a distinction between the optically reliable vs. unreliable and/or rotationally flagged vs. unflagged galaxies, as long as their PAH I − PAop < 45°.

In the text

	Fig. B.1. Comparison of apparent magnitudes in different bands for values obtained from the fit of the Sérsic profiles and the directly measured ones.
In the text

	Fig. B.2. Difference of effective radii obtained from the fit of Sérsic profiles and the directly measured ones for the r−band, given in percentages. The points are color coded using the Sérsic index, with optically reliable galaxies represented as circles and optically unreliable as squares.
In the text

	Fig. B.3. Comparison of stellar masses calculated using either the magnitude obtained from the fitted Sérsic function or the directly measured magnitude from the extracted surface brightness profile (corresponding to $M_{\star }^{\mathrm{meas}}$).
In the text