© ESO 2021

1. Introduction

It is well known that galaxies inhabiting clusters differ from field galaxies, with ellipticals and lenticulars dominating in dense cluster environments and spirals prevailing in the field populations (Hubble & Humason 1931). This observation is the so-called morphology-density relation (Oemler 1974; Dressler 1980). For galaxies of a fixed morphological type, cluster galaxies seem to have older stellar populations (Bower et al. 1990; Rose et al. 1994; Cooper et al. 2010) and suppressed star formation (SF; e.g., Balogh et al. 1997; Gomez et al. 2003; Boselli et al. 2020) relative to field galaxies. Furthermore, galaxies in clusters have a relatively smaller atomic gas reservoir compared to their field counterparts (Haynes et al. 1984; Solanes et al. 2001; Gavazzi et al. 2005). Since almost 40% of galaxies in the universe reside in groups or clusters (e.g., Robotham & Driver 2011), it is clear that the cluster environment plays a crucial role in galaxy evolution in the universe.

Studies of clusters at intermediate redshifts (z <  0.5) show that both the morphologies and the SF rates in cluster galaxies significantly evolve over cosmic time (Poggianti et al. 1999; Dressler et al. 1997; Fasano et al. 2000) with a decreasing fraction of lenticulars and an increasing fraction of spirals towards higher redshifts, while compact clusters at intermediate redshifts have a relatively low fraction of lenticulars with respect to ellipticals and more diffuse clusters have a relatively high fraction of lenticulars (Fasano et al. 2000). In addition, cluster galaxies at intermediate redshifts seem to display a higher SF rate as evidenced by a larger fraction of blue galaxies (Butcher & Oemler 1984), known as the Butcher-Oemler effect, while Abraham et al. (1996) showed that these blue galaxies are preferentially located in the infall region of clusters.

Thus, the dense environments of galaxy clusters somehow expedite the quenching of SF and the morphological transformation of gas-rich, late-type galaxies into gas-poor early-type systems. However, the debate has not been settled yet concerning to what extent the formation (nature) or subsequent evolution of galaxies and the interaction with their surroundings (nurture) give rise to the observed morphology-density relation. Both aspects are likely to play a role but the relative importance of nature and nurture in the formation and evolution of galaxies in dense environments seems to depend on both cosmic time and the environment.

Several gravitational and hydrodynamical mechanisms have been proposed as being responsible for the evolution of galaxies in dense environments such as galaxy groups and clusters (e.g., Gunn et al. 1972; Nulsen 1982; Tinsley & Larson 1979). While gravitational interactions affect both the stellar and the gaseous components of galaxies, hydrodynamical interactions only influence the gaseous component. Gunn et al. (1972) made the first suggestion with their proposed hydrodynamical mechanism of ram-pressure stripping that the interaction between a galaxy’s ISM and the intracluster medium might affect the evolution of a galaxy. Similarly, other mechanisms were proposed such as viscous stripping (Nulsen 1982), starvation (Larson et al. 1980), and thermal evaporation (Cowie & Songaila 1977) considering the effects of viscosity, turbulence and thermal conduction on the depletion of a galaxy’s gas reservoir. Gravitational interactions – such as galaxy-galaxy tidal interactions and mergers (Spitzer & Baade 1951; Tinsley & Larson 1979; Merritt 1983; Springel 2000), galaxy-cluster interactions (Byrd & Valtonen 1990; Valluri 1993), and galaxy harassment (Moore et al. 1996; Jaffé et al. 2016) – mainly occur in denser environments such as galaxy groups and cluster outskirts. In the course of these interactions, the cold gas may also funnel towards the central regions of galaxies and act as fuel for nuclear SF or trigger an active galactic nucleus (AGN; Baldry et al. 2004; Poggianti et al. 2017; Deb et al. 2020). Moreover, cluster galaxies may already have lost their gas and halted their SF activity prior to entering the dense cluster environment, termed as ‘pre-processing’ (Solanes et al. 2001; Tonnesen et al. 2007; Yoon et al. 2017). It is not clear, however, which of these processes are dominant in which cosmic environment, and to what extent they impact the multi-phase ISM and the SF activity of the galaxies.

Since the gas in galaxies is the raw fuel for SF (e.g., Larson 1972; Kennicutt 1998), understanding galaxy evolution requires understanding the prevalence of the various physical processes that are responsible for the accretion and depletion of gas in galaxies in different environments. The extended, fragile, collisional H I gas discs in galaxies serve as excellent diagnostic tracers of environment dependent gravitational and hydrodynamic processes. Observationally, it is well established that the atomic hydrogen discs in galaxies are affected in the dense cluster environments and that cluster galaxies have a smaller H I gas reservoir and often disturbed H I morphologies and kinematics compared to their field counterparts (e.g., Haynes & Giovanelli 1984; Cayatte et al. 1990; Solanes et al. 2001; Gavazzi et al. 2005; Chung et al. 2009).

Until recently, H I studies of the environmental impact on galaxies in clusters have been limited to targeted surveys of galaxies based on their optical properties, typically late-type galaxies, for example, studies of the Virgo cluster (Warmels 1988; Cayatte et al. 1990, 1994; Chung et al. 2009), Abell 262 (Bravo-Alfaro et al. 1997), the Coma cluster (Bravo-Alfaro et al. 2000, 2001), Ursa Major (Verheijen & Sancisi 2001), and Abell 2670 (Poggianti & van Gorkom 2001). In some cases, volume-limited surveys provide a more complete picture such as studies of the Hydra cluster (McMahon 1993), Hercules cluster (Dickey 1997), Abell 2218 (Zwaan et al. 2001), Abell 496, Abell 754, Abell 85 (Van Gorkom et al. 2003), Abell 370 (Lah et al. 2007), Abell 963 and Abell 2192 (Verheijen et al. 2007; Gogate et al. 2020), the ALFALFA Virgo cluster survey (Kent et al. 2008), the Antlia cluster (Hess et al. 2015) and the Fornax cluster (Loni et al. 2021). These studies have provided some insight into the environmental effects on the H I content of cluster galaxies. H I imaging observations show that the H I morphologies in the central cluster galaxies are truncated, asymmetric, or offset from the stellar discs, sometimes trailing the galaxy along its infall trajectory and providing evidence of ram pressure stripping (RPS) as observed in the Coma and Virgo clusters (e.g., Bravo-Alfaro et al. 2001; Chung et al. 2009). For less massive clusters such as the Hydra cluster and Ursa Major, the H I deficiencies in the central galaxies are less dramatic (McMahon 1993; Verheijen & Sancisi 2001; Loni et al. 2021), signifying that those over-densities may be in an earlier stage of formation. Studies have shown (e.g., Solanes et al. 2001; Chung et al. 2009) that there are strong correlations between H I deficiency, proximity to the cluster centre, and the high velocity with which a galaxy moves through the ICM of a cluster.

Wide-area, blind, H I imaging surveys of galaxy clusters are indispensable to reliably distinguish between different environmental mechanisms. By combining the spatially resolved H I morphologies and kinematics with information on the stellar component, SF activity and the global and local environments of a galaxy, it becomes possible to distinguish whether hydrodynamical or gravitational processes dominate the gas depletion in cluster galaxies. With the first pilot observations for massive H I imaging surveys with state-of-the-art telescopes such as MeerKAT (Jonas et al. 2016) and the Australian SKA Pathfinder (ASKAP Johnston et al. 2007, 2008) telescope becoming available, it is timely to broaden our understanding of environmental processes acting on the ISM of galaxies in clusters and their outskirts.

In this paper we present H I and radio continuum imaging data of the galaxy cluster Abell 2626 (hereafter A2626; Abell 1958) obtained with the MeerKAT telescope (Jonas et al. 2016). A2626 is a moderately massive (∼5 × 10¹⁴ M_⊙) galaxy cluster at a redshift of z = 0.055 (Healy et al. 2021). From X-ray observations, A2626 is identified as a cooling-core cluster (Rizza et al. 2000; Wong et al. 2008; McDonald et al. 2018). Extensive X-ray (e.g., Wong et al. 2008; Kadam et al. 2019), optical (e.g. Fasano et al. 2006) and radio continuum (e.g., Rizza et al. 2000; Wong et al. 2008; Gitti 2013; Ignesti et al. 2018, 2020) imaging is available for this cluster, as well as optical spectroscopy (e.g., Cava et al. 2009; Healy et al. 2021) to confirm cluster membership. A2626 is located within the single-dish ALFA Arecibo Legacy Fast survey (ALFALFA; Haynes et al. 2011, 2018) and a few cluster members have been detected in H I by this survey.

To date, radio continuum studies at various frequencies have focused on the centre of the cluster where a peculiar, extended radio source known as the “kite” is located (e.g., Rizza et al. 2000; Wong et al. 2008; Gitti 2013; Ignesti et al. 2018, 2020). X-ray studies of the cluster (e.g., Rizza et al. 2000; Wong et al. 2008; Kadam et al. 2019), as well as studies of the velocity distribution of galaxies within the cluster based on optical redshifts (Mohr et al. 1996; Mohr & Wegner 1997) have shown evidence of a galaxy group merging with the core of the cluster. A recent, dedicated observing campaign with the multi-object optical spectrograph Hectospec on the 6.5 m MMT facility, provided ∼1900 additional redshifts for galaxies with r <  20.4 mag out to 1.5 R₂₀₀ (Healy et al. 2021). Healy et al. (2021) performed a Dressler-Shectman analysis (Dressler & Shectman 1988) which resulted in the identification of two distinct kinematic substructures of galaxies well within the cluster virial radius (R₂₀₀), as well as four distinct kinematic substructures on the outskirts of A2626, which are thought to be galaxy groups in the process of being accreted by the cluster. Thus, there is evidence of recent and ongoing accretion of galaxies and galaxy groups, making A2626 a particularly interesting laboratory to study multiple, ongoing, environment driven processes of gas depletion, SF quenching and morphological transformations of infalling galaxies. We also note that A2626 is less massive than the well-studied Coma cluster but similar in mass to the nearby Virgo cluster while its dynamical state resembles that of the nearby and ∼10× less massive Fornax cluster (Mohr et al. 1996; Paolillo et al. 2002). Therefore, a comparison of the H I properties of the galaxy populations in these four clusters may reveal whether the dynamical state or the mass of a cluster has an impact on the atomic gas content of its constituent galaxies.

Moreover, A2626 hosts six candidate “jellyfish” galaxies (Poggianti et al. 2016), which are potentially experiencing extreme ram-pressure stripping that is displacing the multi-phase gas into tails trailing the galaxy stellar body along its infall trajectory (e.g., Gullieuszik et al. 2018; Poggianti et al. 2019; Ramatsoku et al. 2019; Deb et al. 2020). The confirmed, striking jellyfish galaxy JW100 has been at the centre of multi-wavelength (radio to X-ray) studies by the GASP (GAs Stripping Phenomena in galaxies, Poggianti et al. 2017) collaboration (e.g., Moretti et al. 2019; Poggianti et al. 2019), aimed at understanding how the stripping process affects the stellar and multi-phase gas components of the galaxy. The H I data presented in this paper will contribute to a better understanding of the jellyfish phenomenon by compelling case studies such as for JW100 (Deb et al. in prep.). The principal goal, however, of the MeerKAT observations of A2626 presented in this work is to provide an H I perspective on the various environment dependent astrophysical processes that govern the evolution of galaxies in a dense cluster and its surroundings, analysing the spatially resolved H I discs of individual galaxies as well as employing the H I stacking technique.

This paper is structured as follows: in Sect. 2, we discuss the MeerKAT observations. We describe the calibration and imaging process, as well as the tests carried out to evaluate the quality of the calibration in Sect. 3. In Sect. 4, we discuss how the source finding was performed on the H I spectral line cube. In the subsequent Sects. 5 and 6, we present the global properties of the detected H I sources. Throughout this paper, we use H₀ = 70 km s⁻¹ and Ω_m = 0.3.

2. MeerKAT observations

MeerKAT observations of A2626 were carried out over the nights of 15–17 July 2019, allocated in response to the first call for MeerKAT-64 open time proposals, which were on a shared risk basis (project code SCI-20190418-JH-01). For these first science observations, the correlator capacity was restricted to provide 4k channels across the entire band pass, resulting in a limited frequency resolution (see below).

The full width half maximum (FWHM) of the MeerKAT primary beam provides a field-of-view (FoV) of 61′ at a frequency of 1346 MHz where H I emission from z ≈ 0.055 is detected. This FWHM is adequate to encompass the R₂₀₀ radius of the cluster, which is R₂₀₀ ≈ 1.59 Mpc, corresponding to 25′ at z = 0.055 (Healy et al. 2021). The sensitivity of MeerKAT, however, enables the detection of H I emission well beyond the FWHM of the primary beam and, therefore, we imaged an area of 2 × 2 deg² centred on the cluster (see Sect. 3.2 for details). Consequently, the H I and radio continuum data presented in this publication may occasionally be corrected for a significant primary beam attenuation.

A2626 attains a maximum elevation of 38 deg above the northern horizon of the MeerKAT site, resulting in a rather elongated synthesised beam and necessitating relatively short tracks. A total integration time of 3 × 5 h was accumulated on-source for A2626. Scans of the flux/bandpass calibrator (J1939−6342) were made at the start and end of each track, while scans of the complex (amplitude and phase) gain calibrator (J2253+1608) were made every 20 min for 2 min. No separate polarisation calibrator was observed. The data were accumulated using the SKARAB correlator in full polarisation and 4k mode, with an integration time of 8 s per visibility record. The total MeerKAT bandwidth is 856 MHz with 4096 channels, which are 207 kHz wide corresponding to 45 km s⁻¹ at z = 0. In this work we only processed two subsets of the total bandwidth as follows: 960–1160 MHz, and 1300–1520 MHz for the continuum imaging (see Sect. 3.1); and 1300–1430 MHz for the H I data (see Sect. 3.2). The velocity range covered by the H I data cube spans from cz = −2011 km s⁻¹ to cz = 27766 km s⁻¹, which corresponds to a volume depth of 420 Mpc.

The large FoV of MeerKAT enabled us to observe the cluster using a single pointing. With 15 h on source, the 3σ H I mass detection limit for a galaxy with a line width of 300 km s⁻¹ is M _H I = 2 × 10⁸ M_⊙, with a 5σ column density sensitivity of N_H I = 2.5 × 10¹⁹ cm⁻² per channel when the data are smoothed to 30″ resolution. The observational parameters can be found in Table 1.

3. Data processing

In order to limit the computing resources required to calibrate and Fourier transform the visibilities into images, we only processed a limited subset of the 856 MHz total bandwidth. For the continuum image (see Sect. 3.1) we processed two frequency ranges separately, 960–1160 MHz and 1300–1520 MHz, in order to avoid prolific Radio Frequency Interference (RFI) from navigation satellites. For the H I data cube (see Sect. 3.2 for details) we only processed the frequency range 1300–1430 MHz. For both the continuum and H I spectral line data, we only included the two orthogonal, horizontal (HH) and vertical (VV) polarisations in our data processing. The visibility data were calibrated on the ilifu computer cluster hosted by the Inter-University Institute for Data Intensive Astronomy (IDIA¹) using the CARAcal pipeline (version 1.0, formerly known as MeerKATHI, Józsa et al. 2020).

We followed the same general strategy as Serra et al. (2019) for the cross-calibration using CARAcal. Given the low altitude of A2626 at the MeerKAT site, the data were flagged for geometric shadowing by nearby dishes. Possible RFI in the calibrator data was flagged using AOFlagger (Offringa et al. 2012), which is bundled with CARAcal. We used the default ‘first pass’ strategy of AOFlagger by only inspecting the Stokes Q amplitudes of the visibilities. The next step was to derive the antenna-based complex bandpass for the primary calibrator (J1939−6342). The bandpass solutions were bootstrapped to the secondary gain calibrator (J2253+1608) to determine its flux scale. The time dependent amplitude and phase corrections, as derived for each scan of the secondary gain calibrator, were applied to the target field by interpolating the gain corrections in time across the 20 min time intervals. It is important to note here, that the gain solutions of the secondary calibrator were determined using only the baselines on which the calibrator was unresolved. This means that all projected baselines longer than 7 kλ were excluded as advised by the NRAO calibrator database. With all cross-calibration solutions applied to the target field, the target field was also flagged for RFI using the same ‘first pass’ AOFlagger strategy.

3.1. L-band continuum imaging

Using the cross-calibrated visibilities, we iteratively imaged and self-calibrated the target field using WSClean (Offringa et al. 2014) and CubiCal (Kenyon et al. 2018). We performed three rounds of self-calibration with solution intervals of two minutes, performing two rounds of phase only, followed by one round of amplitude and phase.

We used multi-frequency synthesis (MFS) in WSClean to image the radio continuum of the target field for three frequency ranges: (i) 960–1160 MHz (LF image), (ii) 1300–1520 MHz (HF image), and (iii) 960–1160 combined with 1300–1520 MHz (full image). Shown in Fig. 1 is the 2 × 2 deg² self-calibrated ‘full’ MFS image constructed with a Briggs robust weighting of r = −1, giving an angular resolution of . The effect of primary beam attenuation can be readily appreciated as well as patches of limited dynamical range around bright sources near the half-power point of the primary beam. The three panels to the right of the main image in Fig. 1 show three extended radio galaxies, including the central “kite” source (middle panel), highlighting the resolution and surface brightness sensitivity of MeerKAT.

Fig. 1. 2 × 2 deg2L-band multi frequency synthesis continuum image centred on Abell 2626. This image covers a frequency range of 960–1160 MHz and 1300–1520 MHz, and was imaged with a robust weighting of −1. The image has not been primary beam corrected, however the FWHM of the primary beam at 1520 MHz (the high frequency limit of the image) is given by the inner dashed cyan circle and the outer dashed cyan circle represents the primary beam at 960 MHz (the low frequency limit of the image). The black bar in the bottom right corner indicates 2 Mpc at the distance of A2626 (z = 0.055). The synthesised beam () is shown by the small grey ellipse in the bottom left corner. The three stacked panels on the right are zoomed in on two radio galaxies and the central “kite” source. The small filled cyan ellipse in the lower left corner indicates the synthesised beam.

To correct the images for primary beam attenuation, we calculated a model of the primary beam at the centre frequencies using eidos (Asad et al. 2021). The centre frequency used for the LF image is 1060 MHz and for the HF image 1410 MHz. For the full MFS image we used the average of the two primary beams calculated at 1160 MHz and 1320 MHz. We use an averaged primary beam rather than the eidos measured primary beam to approximate the primary beam at the centre frequency of 1240 MHz for the full MFS image because we have effectively masked all data between 1160 MHz and 1300 MHz, and thus do not have data at the centre frequency.

We used PyBDSF (Mohan & Rafferty 2015) to identify continuum sources in the HF image, using the island threshold set to 2.5σ and the pixel threshold to 3.5σ. We cross-matched our catalogue of continuum sources to a subset of unresolved sources in the NRAO VLA Sky Survey (NVSS, Condon et al. 1998), using a search radius of 3″. This cross-match yielded 12 sources within the FWHM of the MeerKAT primary beam. Using the integrated flux density measurements from our PyBDSF catalogue, we compared the fractional difference between the NVSS flux densities and those from our MeerKAT continuum images. We find good agreement between the NVSS and MeerKAT integrated flux densities (see Fig. 2) at all distances from the pointing centre and at all flux density values. The offset and scatter between the two surveys is 0.9% and 6% respectively, this is shown by the grey bar in both panels of Fig. 2. In Fig. 2, the errorbars on each point represent the uncertainty of the integrated flux density in the MeerKAT image.

Fig. 2. Comparison between the integrated fluxes from NVSS and MeerKAT (1300–1520 MHz) for 12 sources unresolved by NVSS. Top: the fractional difference as a function of NVSS integrated flux density. Bottom: fractional difference as a function of radius from the centre of the pointing. The light grey horizontal band in each panel indicates the standard deviation of the differences.

3.2. H I line data reduction

We subtracted the sky-model created during the self-calibration process (see Sect. 3.1) from the calibrated visibilities. Following the subtraction of the sky-model, we also subtracted a 4th order polynomial fitted to the visibilities to remove any residual continuum. The data were Doppler corrected to a barycentric reference frame and then imaged channel by channel using WSClean with a Briggs robust weighting r = −0.5, a pixel scale of , and no further UV-tapering. Our choice of Briggs weighting is motivated by an informed compromise between angular resolution, H I column density sensitivity and the level of sidelobes of the synthesised beam. We imaged the data over a frequency range of 1300–1430 MHz, which covers the cluster at ∼1346 MHz, foreground galaxies (including Milky Way emission), and some cosmic background structures out to z ∼ 0.09.

Before cleaning, we noticed a ripple in the dirty H I cube along the frequency axis (see Fig. 3). Having investigated possible sources thoroughly, we have not been able to identify the origin of this baseline ripple, and it is not something that has been identified in other MeerKAT H I data. In order to remove this ripple, we fitted a 15th order spline to the spectrum of each pixel. This was done over three iterations where the spline was fitted and subtracted from the data, followed by the clipping of outlying pixels including the H I emission from galaxies. We chose the order of the spline to be sufficiently low to avoid fitting H I emission from galaxies. By using a 15th order spline, we remove spectral structures of order 8.7 MHz, or ∼1850 km s⁻¹, which is much larger than the widest H I emission line. We guarded against including H I emission in the baseline fitted by employing an iterative process, clipping any ‘emission’ above and below the baseline subtracted spectrum and excluding those data from a next spline fit. Had this iterative process not worked satisfactorily, there would have been negative ‘bowls’ on either side of the H I emission and this is not what we observe. The initial ripple is evident in the left panel of Fig. 3. The final ‘de-rippled’ cube with a near uniform noise is shown in the right-hand panel of Fig. 3.

Fig. 3. Vertical slice through the dirty data cube, smoothed to an angular resolution of 30″, taken at RA = 23:38:14.65 covering 2 degrees in declination and the full bandwidth. Left: significant spectral ripples are clearly visible after subtracting the continuum model used for self-calibration, hampering faint source detection. Right: illustration of the effectiveness of subtracting an iteratively fitted spline function of order 15 from each spectrum, resulting in uniform noise. H I emission from a galaxy can be observed in the bottom-left quadrant of this position-velocity slice.

From the antenna pattern of the synthesised beam created during the imaging process, we could measure the angular resolution of the H I cube prior to cleaning. The native resolution of the H I cube varies with frequency: to in the north-south direction, and to in the east-west direction.

The deconvolution or cleaning of the H I emission in the ‘de-rippled’ cube is not straight forward. Since we do not know a priori where the H I emission in the cube will occur, the source finding, masking and cleaning are intertwined and iterative processes that will be explained in the following section.

4. H I source-finding and cleaning

The art of source-finding in an H I blind imaging survey is not a uniquely defined problem. The detection criteria depend on instrumental characteristics such as angular and velocity resolution, noise characteristics and uniformity, the presence of side lobes of the synthesised beam and possible residuals after continuum subtraction. The characteristics of the galaxy population and associated H I properties also play a role in setting the optimum detection criteria, and the source-finding method of choice can be tailored according to the science motivation of a particular survey. It should also be noted that balancing the trade off between completeness and reliability of the detected source population is a major challenge that usually requires an empirical approach.

We have used SoFiA (Serra et al. 2015), the Source Finding Application, which is specifically designed for automatic, reliable H I source-finding in interferometric spectral-line data. SoFiA is publicly available at GitHub. We used the SoFiA-2² package (Westmeier et al. 2021) for our analysis and included the corresponding parameter file we used for our source-finding run in Appendix C to allow our results to be reproduced. SoFiA offers several methods to find sources and we have applied the ‘smooth+clip’ algorithm, which smooths the cube spatially and spectrally with multiple, user-defined Gaussian and boxcar kernels respectively. With each smoothing operation the local rms noise is determined and a 3.5-σ clip level is applied. After the smoothing operations, SoFiA evaluates the dimensions of volumes of connected pixels, and estimates the reliability of the detection by statistically comparing the density of positive and negative detections in a three-dimensional parameter space. If no negative detections with similar characteristics are found in the same area of parameter space, then it assigns a high reliability to the detection. We have required a 90% reliability threshold for a detection to be accepted. SoFiA-2 returns three-dimensional binary masks that we subsequently used to define the regions of H I emission that should be cleaned. The final masks were also used to construct data products such as the global H I profiles, total H I maps and H I velocity fields as described in the next section.

The first step was to take our continuum subtracted, ‘dirty’ data cube at its highest angular resolution and to identify and clean the brightest H I sources in order to remove the side lobes of the synthesised beam that could affect the detection of fainter sources. For this purpose we applied SoFiA-2 with seven Gaussian smoothing kernels that resulted in angular resolutions of 15″, 20″, 25″, 30″, 40″, 50″and 60″, and two spectral boxcar kernels of three and five channels wide (∼135 and ∼22 kms⁻¹). The larger angular smoothing kernels exceed the typical H I sizes of the galaxies at the distance of A2626 (∼10″–15″) but enhance the column density sensitivity of the data. This first run of SoFiA-2 resulted in 111 three-dimensional clean masks throughout the entire cube, many of which contained multiple, blended sources due to the relatively large angular smoothing kernels and the sheer number of gas-rich galaxies in the data cube. Nevertheless, we used these large SoFiA-generated masks to clean the H I emission in the ‘dirty’ cube down to 0.3× the rms noise in a channel using the standard Högbom CLEAN algorithm (Högbom 1974) in the GIPSY software package (Groningen Image Processing SYstem, van der Hulst et al. 1992). We restored the clean components with a Gaussian beam of FWHM , which is the FWHM of the synthesized beam at the frequency that corresponds to the redshift of A2626. This resulted in a high resolution data cube in which most H I sources were cleaned.

In the next step, we used this preliminary cleaned data cube to refine the source-finding process and to construct clean masks that are optimised for the individual galaxies. Given the sensitivity, redshift depth and relatively poor spectral resolution of our MeerKAT data, combined with the diversity of the galaxy population in our survey volume, ranging from nearby dwarfs to distant early-type galaxies, and the possibility of detecting tidal H I features, we concluded that with a single, particular set of SoFiA parameters and smoothing kernels we could not identify all H I sources in the cube.

For example, the highest available angular resolution () is most sensitive to the smallest and unresolved galaxies, albeit with a by-product of yielding more false detections due to a larger number of resolution elements in the cube. However, at that higher resolution, some of the faint extended sources remain undetected due to the limited column-density sensitivity and hence do not exceed the 3.5-σ detection threshold. By smoothing the cube to lower resolutions (50″ or 60″), we improved the column-density sensitivity, thus bringing the faint extended galaxies above the detection threshold as the extended sources better match the larger beam. However, at those large smoothing kernels, we diluted the signal of the smallest galaxies that were barely detected at the highest angular resolution. Figure 4 shows the 5σ column density sensitivity as a function of frequency at the different angular resolutions.

Fig. 4. 5σ column density sensitivity per channel as a function of frequency and the various angular resolutions used in the source finding methodology. Vertical coloured bands indicate redshift ranges corresponding to cosmic over-densities as described in the text. The grey band indicates the frequency of the GPS L3 band, an intermittent RFI source.

Therefore, at first we smoothed the original resolution cube to a circular beam of 15″ using GIPSY to avoid a detection bias in favour of emission that is extended in the north-south direction. Subsequently, we smoothed the 15″ H I cube in SoFiA-2 to different spatial and spectral resolutions and experimented with several detection thresholds and reliability settings in SoFiA-2 in order to explore the trade off between the number of false positives and real sources. Specifically, we smoothed the 15″ H I cube to the resolutions mentioned in Table 2. The normalised spectral noise maps for the native resolution, as well as three circular resolutions (15″, 30″, and 60″) are presented in Fig. 5 and show that the spectral noise is relatively constant across the entire imaged FoV, varying by no more than roughly ±10%.

Fig. 5. Normalised spectral noise maps at various angular resolutions, showing how the noise in the spectrum of each pixel varies across the field of view. The circle indicates the FWHM of the primary beam.

Each of the four SoFiA-2 runs has yielded a list of real and false positive detections, which may or may not be unique for these four different runs, primarily due to stochastic fluctuations. In a blind H I detection experiment, most of the detections are actually very close to the detection threshold due to the slope of the H I mass function. Hence, a small change in the source-finding parameters can result in the removal or addition of detections near the threshold. Thus, by merging the four runs we found some additional genuine detections that would otherwise have been lost due to stochastic fluctuations, thereby increasing the completeness of the catalogue but, at the same time, also accumulating false positives. The final merged list contained 404 unique objects, which may or may not be real H I detections.

We assessed the reliability of each source by searching for a plausible optical counterpart within the SoFiA-2 mask, taking advantage of images from the DECam Legacy Survey (DECaLS Dey et al. 2019), Sloan Digital Sky Survey (SDSS York et al. 2000; Aguado et al. 2019) and Panoramic Survey Telescope and Rapid Response System (Pan-STARRS Chambers et al. 2016). This resulted in 219 sources with a plausible optical counterpart and we considered the other 185 detections as false positives. It is entirely possible, however, that some of these false positives may be actual sources of H I emission associated with tidal debris or remnants of ram-pressure stripping events. Similarly, some of the H I sources may have an optical counterpart as a chance superposition of a background galaxy for which no optical redshift is available. For the 219 sources, confirming optical redshifts are available for 161 sources from the WIde-field Nearby Galaxy-cluster Survey (WINGS Cava et al. 2009), SDSS (Strauss et al. 2002; Aguado et al. 2019) and our MMT/Hectospec survey (Healy et al. 2021).

The next task is to identify the optimum mask for each of the 219 H I sources. If a source is detected in SoFiA in only one of the four runs, then the SoFiA mask generated by the corresponding run is considered as the optimum mask. In case an H I source is detected in multiple runs, the mask from the highest run number is adopted as the optimum mask. Five H I detections turned out to be blended sources, for which we manually de-blended the masks to disentangle the different H I sources in the masks produced by SoFiA. De-blending was done for the following pairs of galaxies: 38, 42, 108, 109, 163, 164, 207, 208, 216, and 217.

The final step is to clean all of the H I sources since for the initial cleaning we have used the masks from a different list of sources (111 sources as mentioned earlier) that does not contain all the 219 H I sources of our final list. For the purpose of cleaning, we have dilated the optimum masks for all 219 sources by 2 pixels (5″) outwards to include emission in the wings of the synthesised beam. We have subsequently used the dilated, optimum masks to define the regions of H I emission in the highest resolution H I cube and cleaned down to 0.3σ. The clean components are restored with a Gaussian beam of similar FWHM as the dirty synthesised beam (), independent of frequency. Finally, we have smoothed the cleaned high resolution cube to a circular beam of 15″ for further analysis.

Figure 6 shows the spatial distribution of the 219 H I detections. The colours of the sources are based to which over-density they belong. Galaxies belonging to A2626 are coloured green, galaxies associated with The Swarm are orange, and the pink points represent galaxies that are associated with Abell 2637. The blue symbols represent galaxies that have the same redshift as the Perseus-Pegasus filament at this location (see Fig. 8 in Healy et al. 2021).

Fig. 6. Sky distribution of the H I detections identified with SoFiA. The inner dashed grey circle represents the R200 of A2626, and the outer light blue dashed circle represents the primary beam FWHM at the frequency of the cluster. The black contours represent the x-ray emission from A2626 taken from ROSAT images. The data points representing the H I detections are coloured by the spectral over-density of which they are a member. The blue indicates likely Perseus-Pegasus filament galaxies (0.035 <  z <  0.0425), green are A2626 galaxies (0.0475 <  z <  0.0615), orange are galaxies in The Swarm (0.0615 <  z <  0.0675), and magenta galaxies represent galaxies part of A2637 (0.0675 <  z <  0.0745); see Sect. 7 for more details. The open squares represent the galaxies also detected by ALFALFA that are used in the H I flux density comparison in Sect. 5.1.

5. H I data products

In this section, we describe the methods we adopted to derive the H I data products and how we determined the H I properties of our sample of 219 H I detected galaxies with optical counterparts.

5.1. H I global profiles

The H I global profiles were constructed by measuring the flux density within the spatially dilated SoFiA masks on a channel-by-channel basis. These are the same masks that were used for the purpose of clean-ing the H I emission. For each galaxy, these 3-dimensional masks vary in shape, size and position from channel to channel due to the rotation of the H I disc. In order to include possible missed emission from the extreme velocity edges of the H I global profiles, for each galaxy we replicated the mask of the first channel to the preceding channel, and the mask of the last channel to the following channel. We used these masks to isolate the data pixels in each channel of the 15″ resolution datacube and applied a channel dependent correction for primary beam attenuation based on Asad et al. (2021). Subsequently, we used the GIPSY task FLUX to measure the spatially integrated flux density of the primary beam corrected data pixels within the masked region of each channel.

Since the size and shape of the masks are different across the frequency channels, the uncertainty in the measured global profile flux at each frequency channel also varies. For calculating the channel dependent uncertainties in the global profile flux, we have taken an empirical approach by projecting the H I mask for a particular channel of a galaxy at 24 different line-free, nearby locations in that same channel. After applying the primary beam correction, the flux inside each offset mask was measured and we calculated the variance of the 24 measurements as the uncertainty on the flux from the galaxy for that particular frequency channel.

We measured the observed widths of the global H I profiles at the 20% and 50% levels of the peak flux. For Gaussian and asymmetric profiles we adopted the single, highest flux of the global profile. For clearly double-horned profiles, however, we adopted separate peak flux levels for the high and low frequency sides of the profile. The frequencies of the two data points above and below the 20% and 50% flux levels at both sides of the profile were linearly interpolated to determine the frequencies , , and on the high and low frequency sides of the global H I profile. We note that these frequencies were determined ‘outside-in’, which implies that we may have overestimated the width of the H I profile in case the flux of the profile edges is not declining monotonically. The uncertainties on these frequencies were determined by considering the uncertainties in the channel-dependent flux levels.

The observed central frequency of the global profile is calculated according to:

(1)

This central frequency is indicated with a black downward arrow in the global profiles as presented in the atlas pages (see Figs. B.1 or B.2 for example). From this observed frequency the galaxy’s systemic recession velocity is calculated as

(2)

where ν₀ is the rest frequency of the H I emission line. The frequency widths of the global profiles are measured according to:

(3)

These frequency widths are indicated by grey horizontal double arrows in the global profiles as presented in the atlas pages (see Figs. B.1 or B.2 for example). From this, the rest frame observed velocity widths at the 20% and 50% levels were calculated according to:

(4)

The calculated values of V_sys, and are tabulated in Appendix A along with their uncertainties based on the uncertainties in . It should be noted that the reported line widths are measured as is and not corrected for the poor instrumental velocity resolution.

We integrated the global H I profiles to obtain the integrated H I flux densities S_V = ∫S_νdV in units of Jy km s⁻¹ and tabulated these values in Appendix A. From these integrated flux density values we infer the total H I gas masses in units of M_⊙ according to:

(5)

where ∫S_νdV is the integrated H I flux density in Jy km s⁻¹ and D_L is the luminosity distance to the galaxy in units of Mpc. The luminosity distance for each galaxy is determined using the galaxy redshift. The inferred H I masses are presented in Fig. 11b where the M_H I uncertainties are determined from the flux uncertainties in each channel of the global profile.

H I flux and profile comparison. In the same way that we have compared the integrated flux densities of a number of sources in our continuum image (see Sect. 3.1) to previously published flux densities of the same sources, we compare a sample of the integrated H I flux densities to previously published values. In this case, we have an overlap of 12 sources with the Arecibo Legacy Fast ALFA survey (ALFALFA; Haynes et al. 2018). The 12 sources have no known nearby galaxies that could possibly confuse the H I spectrum, and are all classified as ‘code 1’ detections in the ALFALFA catalogue (S/N >  6.5). All 12 sources are located outside the FWHM of the MeerKAT primary beam. The ALFALFA and MeerKAT H I spectra for the 12 galaxies are shown in Fig. 7, and the measured line widths and integrated flux densities are presented in Table 3.

Fig. 7. Global profiles for the 12 galaxies identified in the ALFALFA survey that are also identified in the MeerKAT H I data. The MeerKAT spectra are shown in orange, the ALFALFA spectra in light grey (original spectral resolution) and dark grey (MeerKAT spectral resolution). The red and blue horizontal dot-dash lines indicate where the W20 and W50 respectively are measured for the MeerKAT spectra.

Figure 8 shows the comparison between ALFALFA and MeerKAT for the properties measured from the H I global profiles (see Sect. 5.1). The H I redshift and integrated flux density measurements for the ALFALFA galaxies are taken directly from the α.100 catalogue (Haynes et al. 2018) as listed Table 3. Using the same procedure as discussed in Sect. 5.1, we re-measured the W₂₀ and W₅₀ line widths from the ALFALFA spectra (R = 11 km s⁻¹) after they were re-sampled to the same spectral resolution as the MeerKAT spectra (R = 45 km s⁻¹). The re-measured line widths are presented in Table 3, and are used in the comparison of line widths between the two surveys.

Fig. 8. Comparison between measured properties from the global profiles of galaxies in both the MeerKAT and ALFALFA surveys. (a) the difference between the recessional velocities calculated from the H I redshifts. (b) the log ratio of the line widths, open circles represent the W50 and filled circles represent the W20. (c) log ratio of the integrated flux densities.

Panel a of Fig. 8 shows the comparison between the measured H I redshifts, shown here as the recession velocity cz_H I. The offset between the MeerKAT measurements and ALFALFA is 4 km s⁻¹ with a scatter of 11 km s⁻¹. Given that both the offset and the scatter are much less than the velocity width of a single channel (45 km s⁻¹), we conclude that the measured H I redshifts for the two surveys are equivalent.

In the lower two panels of Fig. 8, we present the log ratios of the measured profile line widths (b) and the integrated flux densities (c). The W₂₀ measurements are plotted as filled circles, and the W₅₀ measurements as open circles. The offset for the W₂₀ and W₅₀ measurements are 1.8 per cent (0.008 dex) and 4.7 per cent (0.02 dex) respectively, while the scatters are 9.6 per cent (0.04 dex) and 12 per cent (0.05 dex) respectively. For the smallest galaxies in the overlap (∫Sdv^MK <  1 mJy km s⁻¹), ALFALFA measures slightly higher integrated flux density than MeerKAT. The offset and scatter between the two surveys for galaxies with ∫Sdv^MK >  1 mJy km s⁻¹ is 1.6 per cent (0.007 dex) and 12 per cent (0.05 dex).

5.2. Total H I maps

The total H I maps display the integrated H I column densities of the H I gas along the line of sight. These maps are constructed by applying the same dilated, optimum 3-dimensional masks used to derive the global H I profiles, in order to isolate the primary beam corrected data pixels in the clean-ed data cube at 15″ resolution that contain the H I signal. The pixels outside the masks are set to zero and the pixel values inside the mask are summed along the frequency axis to obtain the 2-dimensional total H I maps. The resulting map is then multiplied by the rest frame channel width in km s⁻¹ to set the units of the pixel values to Jy km s⁻¹ beam⁻¹. The advantage of this method is that we avoid adding noise from pixels without H I signal and thereby obtain a higher signal-to-noise for certain pixels in the H I maps. A complication of this method is that the noise becomes non-uniform across the H I map. Consequently, the 3σ column density level in the H I maps is not well defined and will be different for each galaxy in our sample. Next, we explain how the pseudo-3σ contour level for the total H I map is estimated.

To define the approximate 3σ column density level in a total H I map, we constructed a signal-to-noise map for each galaxy. We took an empirical approach similar to how we estimated the uncertainties in the global H I profiles. We shifted the 3-dimensional mask to 24 nearby, line-free locations in the 15″ data cube and followed the procedure above to construct 24 2-dimensional maps that contain only noise. Subsequently, for corresponding pixels in these emission-free maps we calculated the variance of the 24 values to create a representative noise map for each galaxy. The total H I map was divided by this noise map to obtain a signal-to-noise map. From this signal-to-noise map we selected the pixels with and calculated the average value of the corresponding pixels in the total H I map. We adopt this average value as the representative 3σ H I column density level, being the lowest contour in the total H I map of each galaxy that is reliable. We note, however, that this pseudo-3σ contour may run through areas in the total H I map where the signal-to-noise value may deviate significantly from 3. The signal-to-noise maps and the pseudo-3σ contour are included in the atlas pages presented below.

Subsequently, the pixel values in the total H I map are converted from Jy km s⁻¹ beam⁻¹ to H I column densities in units of cm⁻² using:

(6)

where N_HI is the H I column density in cm⁻², dV is the rest frame velocity width in km s⁻¹ over which the emission line is integrated at each pixel in the map and T_b is the brightness temperature in K, which is calculated according to:

(7)

where S_ν is the flux density in units of mJy/beam, ν₀ and ν are the rest frequency and the observed frequency of the H I emission respectively, and Θ_x and Θ_y are the major and minor axes of the Gaussian synthesized beam in arcseconds.

5.3. H I velocity fields

A galaxy’s H I velocity field displays the rest frame line-of-sight component of the velocity vector along which the H I gas moves in the gravitational potential of the galaxy. Usually the systemic heliocentric recession velocity, cz, is added to this projected velocity, indicated by the iso-velocity contour along the kinematic minor axis of the galaxy’s rotating H I gas disc.

Various methods can be adopted to construct a 2-dimensional H I velocity field, ranging from calculating an intensity weighted first moment of the H I spectrum at each pixel position, fitting (multiple) Gaussians, a one-sided Gaussian or a Gauss-Hermite function in an attempt to counteract the effects of beam smearing. The method of choice typically depends on the resolution and signal-to-noise of the data. We decided to fit a single Gaussian to the H I spectrum at each pixel position in the total H I map and used initial estimates for the centroid and dispersion of the Gaussian fit from an intensity weighted first and second moment calculation of the H I spectra censored by the 3-dimensional mask. The peak flux density of the censored spectrum was adopted as the initial estimate for the amplitude of the Gaussian. With these initial estimates, we fitted single Gaussians to the full, uncensored spectra but forced the baseline to zero. Uncertainties on the Gauss parameters were calculated on the basis of the noise in a spectrum.

Gaussian fits are accepted only when satisfying the following conditions: (1) the amplitude is larger than 3 times the rms noise in that profile, (2) the uncertainty in the velocity centroid of the Gaussian is less 1/3 of the width of a velocity channel, (3) the velocity dispersion is within the range 15–100 km s⁻¹ considering cases of edge-on galaxies suffering from significant beam-smearing. We have to keep in mind that the velocity resolution is comparatively poor (45 km s⁻¹) and the H I emission for most galaxies is limited to only a few velocity channels (3–5). Consequently, for many galaxies, the velocity fields are patchy with holes due to limited signal-to-noise and, therefore, a lack of acceptable fits.

5.4. H I position-velocity diagrams

Position-velocity diagrams of inclined rotating H I gas discs are generally constructed by extracting a two-dimensional slice from the data cube, centred on the dynamical centre of a galaxy and in the direction of the kinematic major axis. Such position-velocity diagrams demonstrate the projected shape of a galaxy’s rotation curve and may indicate the presence of any kinematic asymmetry of the rotating gas disc. Usually, the kinematic major axes and dynamical centres are calculated from the two-dimensional velocity fields of spatially well resolved galaxies. Since most galaxies in our sample are not well resolved and the H I velocity fields are not of the highest quality due to the poor velocity resolution of our data, we have adapted a different strategy for determining the central position and position angle of the position-velocity diagrams.

PAs were initially determined by visual inspection of the optical image while the receding side with respect to the systemic recession velocity of the galaxy was determined from the H I kinematics, either from the two-dimensional velocity field or from the three-dimensional cube if the former was poorly defined. For nearly face-on galaxies with a poorly defined optical PA, the PA of the receding side was decided solely from the H I velocity field. In case of ambiguity of both the optical image and the H I velocity field, the PA was based on the elongation of the H I map. In few cases, the PA was chosen to cut through interesting H I morphological features that were not aligned with the optical or kinematic major axes of the galaxy. Instead of using the dynamical centres derived from the velocity fields of the rotating discs, we have adopted the optical centres of the galaxies from the SDSS source name. The position-velocity diagrams presented in the atlas pages (see Sect. 6 and Appendix B) were extracted from the 15″ resolution cube.

6. The H I catalogue and H I Atlas

In the previous section, we described the methodology used to produce the H I data products. In this section, we discuss how the H I atlas pages for the 219 H I detections with optical counterparts are presented.

6.1. The H I catalogue

The observed H I parameters of all 219 galaxies detected in our MeerKAT survey are listed in Table A.1. These properties are derived from the global H I profiles as presented in the atlas pages. Information that can be derived from the two-dimensional H I column density maps and velocity fields, such as the size and asymmetry of the H I gas disc and its kinematic position angle, will be presented in a forthcoming paper. The column entries of Table A.1 are as follows:

Column (1): The assigned H I identifier based on the ascending Right Ascension (J2000.0) of their optical counterparts as provided by SDSS (DR14).

Column (2): The SDSS identifier for the optical counterpart of the H I detected galaxies, based on their Right Ascension and Declination (J2000.0)

Columns (3), (4): Heliocentric, systemic recession velocity (cz) based on the H I redshift measured as the midpoint of the 20% line width of the H I global profile.

Columns (5)–(8): The observed rest frame line width of the H I global profiles and their uncertainties in km s⁻¹, computed at the 20% and 50% level of the H I peak flux respectively, measured at a velocity resolution of 45 (1 + z)  km s⁻¹. See Sect. 5.1 for details.

Column (9) and (10): The measured, primary beam corrected, integrated H I flux densities and their uncertainties in Jy km s⁻¹.

Column (11) and (12): The radio continuum flux density for the H I detected sources and their uncertainties based on the primary beam corrected continuum maps with robust weighting, r = −0.5.

6.2. The H I atlas

The H I data products as discussed in Sect. 5 are shown as atlas pages for all 219 H I detected galaxies with fully reduced data. A schematic of the atlas pages is shown in Fig. 9, each of the following subsections correspond to one of the panels on the atlas page.

Fig. 9. Layout of the H I atlas for the 219 galaxies with fully reduced data. See Appendix B.

6.2.1. Optical image

A DECaLS r-band greyscale image is presented for each galaxy. The SDSS identifier of the galaxy associated with the H I detection is shown above the image. The physical size of the image in terms of kpc at the H I redshift of the galaxy is denoted by a scalebar in the bottom right of the panel. The orange open circle indicates the optical centre from the SDSS.

6.2.2. H I column density contours on optical image

Contours of the H I column density distribution are overlaid on the DECaLS r-band image represented in greyscale in the background. For each galaxy, the H I map is constructed based on the optimum H I mask, which is the SoFiA mask at the corresponding resolution of the H I detection or the lowest resolution mask in case the H I source was detected at multiple resolutions. The H I column density contours are drawn at the levels with n = 0,1,2,3,… while , the average 3σ column density for each galaxy, is printed in the top-left corner of the signal-to-noise panel on the atlas page (see below). The FWHM of the beam size for each H I source is indicated by the hatched circle in the bottom-left. The physical size of the map in units of kpc is presented as a scalebar in the bottom-right corner of the panel.

6.2.3. Total H I map

The total H I map is presented in greyscale with the H I column density contours superimposed. All positive pixels inside the total H I map have an assigned greyscale value. The H I column density contour levels are the same as the neighbouring aforementioned panel. Greyscale pixels outside the outermost contour typically have a signal-to-noise below 3. The purple line indicates where we have extracted the position velocity diagram, starting at the receding side of the galaxy (see Sect. 6.2.7). The orange open circle signifies the optical centre from SDSS.

6.2.4. Signal-to-noise map

The H I signal-to-noise map is represented as a colour map with the different colours representing different levels of signal-to-noise as indicated by the vertical scale bar on the right of the panel. The average 3σ H I column density for each galaxy is calculated and printed in the top-left corner of the panel. The H I contours are repeated in white, allowing to evaluate the signal-to-noise along the lowest contour and thereby its reliability.

6.2.5. Radio continuum map on optical

The cleaned 21 cm radio continuum map is plotted as contours superimposed on the DECaLS r-band image in greyscale in the background. The continuum contours are drawn at the levels of 2ⁿ× the rms noise, where n = 0,1,2,3,.. . Dashed contours are drawn at −2× the rms noise. The rms noise levels are reported in the upper-left corner of the continuum map after correcting for primary beam attenuation and, therefore, may vary significantly for different sources. The FWHM of the beam size for each radio continuum source is indicated by the hatched circle in the bottom-left.

6.2.6. H I velocity field

The greyscale image indicates where radial velocities could be measured. Lighter greyscale and blue iso-velocity contours signify the approaching side of the galaxy while darker greyscale and red iso-velocity contours denote the receding side of the H I disc. The white contour indicates the systemic velocity derived from the H I global profile. The contours are drawn at intervals of ±n × 25  km s⁻¹, both for the approaching and receding sides. The orange open circle in the middle indicates the optical centre from the SDSS.

6.2.7. Position-velocity (PV) diagram

The position-velocity diagram is made along the major axis of the optical image, starting at the receding side of the galaxy. The coordinates of the optical centre are adopted from the SDSS. The strategies we adopted for determining the central positions and the PAs are mentioned in Sect. 5.4.

The PA of each galaxy is printed in the top-left corner of the panel. The directions of the start and end points of the slices are indicated respectively in the bottom left and right corners of each panel.

For all the galaxies, the PV diagrams were extracted from the cleaned cube at 15″ resolution and at the 45 km s⁻¹ velocity resolution. The contour levels are (−2, 2, 3, 4.5, 6, 9, 12, 15, 20, 25) × rms noise. The orange contour outlines the dilated, optimum mask as automatically generated by SoFiA and used for each galaxy as described in Sect. 4. The vertical dashed line corresponds to the optical centre of the galaxy. The horizontal dashed line indicates the central frequency as derived from the global profile.

6.2.8. H I global profile

The H I spectrum or global profile of each galaxy is constructed by applying the dilated, optimum H I mask for each galaxy to the 15″ resolution datacube as described in Sect. 5.1. The method of deriving the errorbars as shown on the global profile is also described in Sect. 5.1. The vertical black downward arrow in the middle shows the central frequency, corresponding to the systemic velocity, which is the midpoint of the 20% velocity width (W₂₀). The H I redshift is calculated from this central frequency and reported in the top right of the panel. For the galaxies with measured optical redshifts, an orange upward arrow is indicating their optical redshifts. The grey, horizontal double-arrows indicate the measured line widths at the 20% and 50% level of the peak intensity as tabulated in the Appendix A.

7. Sample properties of the H I detections

In this section we briefly discuss the distributions of the measured parameters, the inferred H I masses as a function of redshift and the redshift distribution of galaxies within the surveyed volume.

7.1. Distributions of detection related parameters

The detectability of a source depends on its distance (redshift) and H I mass, resulting in an integrated H I line flux density (S_int = ∫Sdv), as well as its H I line width (W₂₀), peak flux density (F_peak), angular size and distance from the centre of the primary beam, which attenuates the flux density from sources away from the pointing centre. In this subsection we explore the distributions of several of these observables (e.g., Kovač et al. 2009) to gain insight in the overall properties of our H I detections.

Figure 10 shows the W₂₀ line width of each galaxy as a function of the log of the peak flux density of the global profile (panel a) and the log of the integrated flux density of the galaxy (panel b). The distributions show that there are no galaxies with large W₂₀ with small integrated flux densities and low peak flux densities in the observed volume. This is because galaxies with a small integrated or peak flux densities that is spread over a large velocity width would be very difficult to detect and might lie below our detection threshold. By smoothing in velocity and hence increasing the sensitivity, we could detect some galaxies with low flux densities and large line widths, still our list of H I sources is not complete. The detection thresholds based on the three velocity smoothing kernels used for the source-finding are represented by the black lines in Fig. 10b. The grey histograms atop panels a and b, and to the left of panel b in Fig. 10 show the one-dimensional distributions of log F_peak, log S_int, and W₂₀ respectively.

Fig. 10. Bivariate distributions of measured properties from the global H I profiles of the H I detected sources. Peak flux density (a) and integrated flux density (b) as a function of W20 line width, derived from the global H I profiles. The detection thresholds based on the three velocity smoothing kernels used for the source-finding are represented by the black lines in Fig. 10b. The colours of the data points represent the over-density to which the galaxy belongs: A2626 sources are green, The Swarm sources are orange, A2637 sources are pink, and the sources coloured blue belong to the Perseus-Pegasus filament. Grey data points are galaxies in the volume that do not belong to a spectral over-density. The histograms along the top and right side of the panels show the 1D distribution of the observable along that axis.

The distribution of integrated flux density with redshift (z), presented in Fig. 11a, shows a segregation of detections in the various over-densities mentioned earlier: A2626, The Swarm and A2637. We observe a lack of galaxies with small S_int (<  − 1.5logS_int Jy⁻¹ km s⁻¹) at low z (z <  0.03). This is likely due to two reasons, partly due to the smaller surveyed cosmic volume at low z, and partly since the low mass dwarf galaxies which we have the sensitivity to detect at this low redshift have velocity widths of the order of the channel width, thus escaping detection. In Fig. 11a, lines of a constant H I mass in M_⊙ with a logarithmic step of 1 dex are shown as dotted lines. The H I mass is calculated using Eq. (5) where the distance is the luminosity distance calculated assuming H₀ = 70 km s⁻¹ Mpc⁻¹ and Ω_m = 0.3. We observe that most of the galaxies have H I masses in the range of 10⁹ − 10¹⁰ M_⊙. At low redshifts, the cosmic volume we survey is much smaller and large galaxies are atypical. Furthermore, the large velocity channel width makes it very difficult to detect galaxies with narrow line widths (e.g., the low-mass dwarf population).

Fig. 11. (a) Integrated flux versus redshift. The dotted lines show lines of constant H I mass. (b) Measured H I mass as a function of redshift. The coloured round data points in both panels indicate the 219 H I detections; the colours correspond to the spectral over-densities highlighted in the histogram in panel b, and the black data points correspond to the grey histograms that include sources not associated with any over-density. The light grey open histogram represents the optical redshift catalogue from Healy et al. (2021). The dashed grey lines indicate the 5σ H I mass limits at the three velocity resolutions used in the source finding: dashed for 45 km s−1, dot-dashed for 135 km s−1, and dotted for 225 km s−1.

7.2. H I mass versus redshift

As mentioned in Sect. 2, the expected 3σ H I mass limit for a galaxy with an H I line width of 300 km s⁻¹ was M_H I = 2 × 10⁸ M_⊙ at the distance of the cluster and the centre of the field-of-view. Figure 11b shows the distribution of the H I masses of the detected galaxies as a function of redshift. The dashed and dotted grey lines indicate the redshift dependent mass limits at the centre of the primary beam, based on the source finding criteria.

The redshift distribution of the H I detected galaxies is shown by the filled histogram at the bottom of Fig. 11b, while the unfilled histogram represents the optical redshift catalogue covering the field from Healy et al. (2021). The optical catalogue is 70% complete across the surveyed area for galaxies brighter than r = 19.1 mag. Most of the H I detected galaxies fall in the interval 0.0475 <  z <  0.0745. The green over-density corresponds to the redshift range of A2626 (0.0475 <  z <  0.0615), and the orange and magenta to respectively the two background over-densities The Swarm (0.0615 <  z <  0.0675), an association of galaxies and galaxy groups (Healy et al. 2021) and the galaxy cluster A2637 (0.0675 <  z <  0.0745, Abell et al. 1989; Healy et al. 2021). The blue histogram identifies a distinct over-density in the foreground, likely associated with the Perseus-Pegasus filament (0.035 <  z <  0.0425; Batuski & Burns 1985; Healy et al. 2021).

8. Conclusions and future work

We have conducted a blind 21 cm H I imaging survey with a single pointing of MeerKAT in 4k mode targeting the Abell 2626 galaxy cluster at cz = 16576 km s⁻¹, over an effective area of ∼4 deg² with an angular resolution of , a velocity resolution of 45 km s⁻¹ and average noise of 80 μ Jy beam⁻¹ channel⁻¹ over the radial velocity range of 0–25 000 km s⁻¹. The main summary points of this work are:

1. In line with our expectations, our MeerKAT data have enabled us to make 219 direct detections of galaxies with H I mass M_H I > 10⁸ M_⊙. Only 12 of the galaxies were previously detected in H I in the ALFALFA (Haynes et al. 2011, 2018) single-dish survey. We have presented the H I properties of the 219 H I detected galaxies as a data catalogue and as an atlas with total H I map, H I velocity field, signal-to-noise H I map, position-velocity diagram along optical major axis, H I global profile or spectrum, and an overlay with optical DECaLS r-band image for each galaxy.

2. We identified the H I sources in the datacube using SoFiA-2 with the smooth+clip finder algorithm with four different combinations of spatial and spectral smoothing and then identified the unique sources combining those four catalogues. The source finding was non-trivial, since the number of detections depend on determining a trade-off between reliability and completeness. Because of our MMT/Hectospec redshift survey and existing ancillary data from DECaLS, WINGS, SDSS, and PanSTARRS, we could finally settle at 219 H I sources with reliable optical counterparts.

One of the principal objectives of observing A2626 was to image the striking jellyfish galaxy JW100 (Moretti et al. 2019; Poggianti et al. 2019) and a few other jellyfish candidate galaxies (JW98, JW99, JW101, JW102, JW103, Poggianti et al. 2016) in H I. Further work on the jellyfish galaxies in A2626 is presented in Deb et al. (in prep.). There are three main over-densities in the line-of-sight of the observation: A2626 itself, The Swarm (Healy et al. 2021) and A2637. Identification of clusters and sub-structures in this pointing area has been done and presented in Healy et al. (2021). Detailed analysis of H I morphology and kinematics of the direct H I detections will be presented in Deb et al. (in prep.) and analysis of H I stacking of non-detections in Healy et al. (in prep.).

Motivated by the unexpected high number of direct H I detections due to MeerKAT’s unprecedented sensitivity, our H I study of A2626 will be expanded by four additional MeerKAT pointings in 32k-mode covering the ENE-SSW sector of A2626. The observations have been allocated time and will be used to investigate the H I properties of galaxies before they have entered the cluster environment. These higher spectral resolution observations will complement the large collection of H I detected galaxies presented in this paper, will allow detailed investigations of the kinematics of detected galaxies, and also enable the detection of low mass dwarf galaxies that escaped detection in the 4k data.

Supplementary Material

Full Figs. B.1 and B.2 (atlas). Access here

Acknowledgments

We thank the anonymous referee for their useful comments that have improved this work. We thank J Delhaize for comments on this manuscript. JH acknowledges research funding from the South African Radio Astronomy Observatory. This paper makes use of the MeerKAT data (Project ID: SCI-20190418-JH-01). The MeerKAT telescope is operated by the South African Radio Astronomy Observatory, which is a facility of the National Research Foundation, an agency of the Department of Science and Innovation. MV acknowledges support by the Netherlands Foundation for Scientific Research (NWO) through VICI grant 016.130.338. MR’s research is supported by the SARAO HCD programme via the “New Scientific Frontiers with Precision Radio Interferometry” research group grant. We acknowledge the use of the ilifu cloud computing facility – www.ilifu.ac.za, a partnership between the University of Cape Town, the University of the Western Cape, the University of Stellenbosch, Sol Plaatje University, the Cape Peninsula University of Technology and the South African Radio Astronomy Observatory. The ilifu facility is supported by contributions from the Inter-University Institute for Data Intensive Astronomy (IDIA – a partnership between the University of Cape Town, the University of Pretoria and the University of the Western Cape), the Computational Biology division at UCT and the Data Intensive Research Initiative of South Africa (DIRISA). This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement no. 679627; project name FORNAX). The data published here have been reduced using the CARACal pipeline, partially supported by ERC Starting grant number 679627 “FORNAX”, MAECI Grant Number ZA18GR02, DST-NRF Grant Number 113121 as part of the ISARP Joint Research Scheme, and BMBF project 05A17PC2 for D-MeerKAT. Information about CARACal can be obtained online under the URL: https://caracal.readthedocs.io. This project has received funding from the European Research Council (ERC) under the Horizon 2020 research and innovation programme (grant agreement N. 833824). We acknowledge financial contribution from the contract ASI-INAF n.2017-14-H.0, from the grant PRIN MIUR 2017 n.20173ML3WW_001 (PI Cimatti) and from the INAF main-stream funding programme (PI Vulcani). We acknowledge support from the Italian Ministry of Foreign Affairs and International Cooperation (MAECI Grant Number ZA18GR02) and the South African Department of Science and Innovation’s National Research Foundation (DSI-NRF Grant Number 113121) as part of the ISARP RADIOSKY2020 Joint Research Scheme. This research made use of Astropy (http://www.astropy.org), a community-developed core Python package for Astronomy (Astropy Collaboration 2013, 2018).

References

Appendix A: Catalogue of H I detections

Appendix B: Examples of H I atlas pages

Fig. B.1. Atlas page for H I ID = 45. See section 6 for details. Clockwise from the top left panel: DECaLS r-band image, total H I map on optical image, total H I map on greyscale, signal-to-noise map, radio continuum map on optical, H I velocity field, position-velocity (PV) diagram, H I global profile. The full collection of 219 atlas pages is available online.

Fig. B.2. Atlas page for H I ID = 128. See section 6 for details. The full collection of 219 atlas pages is available online.

Appendix C: SoFiA parameter file used for the H I source finding

### SoFiA 2.1.1 (default\_parameters.par)      ###
### Source Finding Application                               ###
### Copyright (C) 2020 Tobias Westmeier                      ###

# Global settings
pipeline.verbose = false
pipeline.pedantic = true
pipeline.threads = 0

# Input
input.data = A2626.norip.cl.bs15.fits
input.region =
input.gain =
input.noise =
input.weights =
input.mask =
input.invert = false

# Flagging
flag.region =
flag.auto = false
flag.threshold = 5.0
flag.radiusSpatial = 0
flag.log = false

# Noise scaling
scaleNoise.enable = true
scaleNoise.mode = local
scaleNoise.statistic = mad
scaleNoise.fluxRange = negative
scaleNoise.windowXY = 225
scaleNoise.windowZ = 1
scaleNoise.gridXY = 75
scaleNoise.gridZ = 1
scaleNoise.interpolate = true
scaleNoise.scfind = false

# S+C finder
scfind.enable = true
scfind.kernelsXY = 0
scfind.kernelsZ = 0, 3
scfind.threshold = 3.5
scfind.replacement = 2.0
scfind.statistic = mad
scfind.fluxRange = negative

# Threshold finder
threshold.enable = false
threshold.threshold = 5.0
threshold.mode = relative
threshold.statistic = mad
threshold.fluxRange = negative

# Linker
linker.radiusXY = 3
linker.radiusZ = 1
linker.minSizeXY = 3
linker.minSizeZ = 1
linker.maxSizeXY = 0
linker.maxSizeZ = 0
linker.keepNegative = false

# Reliability
reliability.enable = true
reliability.threshold = 0.90
reliability.scaleKernel = 0.15
reliability.fmin = 5.0
reliability.plot = true

# Mask dilation
dilation.enable = false
dilation.iterations = 1
dilation.threshold = -1

# Parameterisation
parameter.enable = true
parameter.wcs = true
parameter.physical = false
parameter.prefix = SoFiA
parameter.offset = false

# Output
output.directory = /data/users/deb/Cube_bs15
output.filename = A2626_bs15
output.writeCatASCII = true
output.writeCatXML = false
output.writeCatSQL = false
output.writeNoise = false
output.writeFiltered = false
output.writeMask = true
output.writeMask2d = false
output.writeRawMask = false
output.writeMoments = true
output.writeCubelets = true
output.marginCubelets = 5
output.overwrite = true

All Tables

All Figures

Fig. 1. 2 × 2 deg2L-band multi frequency synthesis continuum image centred on Abell 2626. This image covers a frequency range of 960–1160 MHz and 1300–1520 MHz, and was imaged with a robust weighting of −1. The image has not been primary beam corrected, however the FWHM of the primary beam at 1520 MHz (the high frequency limit of the image) is given by the inner dashed cyan circle and the outer dashed cyan circle represents the primary beam at 960 MHz (the low frequency limit of the image). The black bar in the bottom right corner indicates 2 Mpc at the distance of A2626 (z = 0.055). The synthesised beam () is shown by the small grey ellipse in the bottom left corner. The three stacked panels on the right are zoomed in on two radio galaxies and the central “kite” source. The small filled cyan ellipse in the lower left corner indicates the synthesised beam.

In the text

	Fig. 2. Comparison between the integrated fluxes from NVSS and MeerKAT (1300–1520 MHz) for 12 sources unresolved by NVSS. Top: the fractional difference as a function of NVSS integrated flux density. Bottom: fractional difference as a function of radius from the centre of the pointing. The light grey horizontal band in each panel indicates the standard deviation of the differences.
In the text

Fig. 3. Vertical slice through the dirty data cube, smoothed to an angular resolution of 30″, taken at RA = 23:38:14.65 covering 2 degrees in declination and the full bandwidth. Left: significant spectral ripples are clearly visible after subtracting the continuum model used for self-calibration, hampering faint source detection. Right: illustration of the effectiveness of subtracting an iteratively fitted spline function of order 15 from each spectrum, resulting in uniform noise. H I emission from a galaxy can be observed in the bottom-left quadrant of this position-velocity slice.

In the text

	Fig. 4. 5σ column density sensitivity per channel as a function of frequency and the various angular resolutions used in the source finding methodology. Vertical coloured bands indicate redshift ranges corresponding to cosmic over-densities as described in the text. The grey band indicates the frequency of the GPS L3 band, an intermittent RFI source.
In the text

	Fig. 5. Normalised spectral noise maps at various angular resolutions, showing how the noise in the spectrum of each pixel varies across the field of view. The circle indicates the FWHM of the primary beam.
In the text

Fig. 6. Sky distribution of the H I detections identified with SoFiA. The inner dashed grey circle represents the R200 of A2626, and the outer light blue dashed circle represents the primary beam FWHM at the frequency of the cluster. The black contours represent the x-ray emission from A2626 taken from ROSAT images. The data points representing the H I detections are coloured by the spectral over-density of which they are a member. The blue indicates likely Perseus-Pegasus filament galaxies (0.035 <  z <  0.0425), green are A2626 galaxies (0.0475 <  z <  0.0615), orange are galaxies in The Swarm (0.0615 <  z <  0.0675), and magenta galaxies represent galaxies part of A2637 (0.0675 <  z <  0.0745); see Sect. 7 for more details. The open squares represent the galaxies also detected by ALFALFA that are used in the H I flux density comparison in Sect. 5.1.

In the text

Fig. 7. Global profiles for the 12 galaxies identified in the ALFALFA survey that are also identified in the MeerKAT H I data. The MeerKAT spectra are shown in orange, the ALFALFA spectra in light grey (original spectral resolution) and dark grey (MeerKAT spectral resolution). The red and blue horizontal dot-dash lines indicate where the W20 and W50 respectively are measured for the MeerKAT spectra.

In the text

	Fig. 8. Comparison between measured properties from the global profiles of galaxies in both the MeerKAT and ALFALFA surveys. (a) the difference between the recessional velocities calculated from the H I redshifts. (b) the log ratio of the line widths, open circles represent the W50 and filled circles represent the W20. (c) log ratio of the integrated flux densities.
In the text

	Fig. 9. Layout of the H I atlas for the 219 galaxies with fully reduced data. See Appendix B.
In the text

Fig. 10. Bivariate distributions of measured properties from the global H I profiles of the H I detected sources. Peak flux density (a) and integrated flux density (b) as a function of W20 line width, derived from the global H I profiles. The detection thresholds based on the three velocity smoothing kernels used for the source-finding are represented by the black lines in Fig. 10b. The colours of the data points represent the over-density to which the galaxy belongs: A2626 sources are green, The Swarm sources are orange, A2637 sources are pink, and the sources coloured blue belong to the Perseus-Pegasus filament. Grey data points are galaxies in the volume that do not belong to a spectral over-density. The histograms along the top and right side of the panels show the 1D distribution of the observable along that axis.

In the text

Fig. 11. (a) Integrated flux versus redshift. The dotted lines show lines of constant H I mass. (b) Measured H I mass as a function of redshift. The coloured round data points in both panels indicate the 219 H I detections; the colours correspond to the spectral over-densities highlighted in the histogram in panel b, and the black data points correspond to the grey histograms that include sources not associated with any over-density. The light grey open histogram represents the optical redshift catalogue from Healy et al. (2021). The dashed grey lines indicate the 5σ H I mass limits at the three velocity resolutions used in the source finding: dashed for 45 km s−1, dot-dashed for 135 km s−1, and dotted for 225 km s−1.

In the text

	Fig. B.1. Atlas page for H I ID = 45. See section 6 for details. Clockwise from the top left panel: DECaLS r-band image, total H I map on optical image, total H I map on greyscale, signal-to-noise map, radio continuum map on optical, H I velocity field, position-velocity (PV) diagram, H I global profile. The full collection of 219 atlas pages is available online.
In the text

	Fig. B.2. Atlas page for H I ID = 128. See section 6 for details. The full collection of 219 atlas pages is available online.
In the text