© ESO 2019

1. Introduction

Hickson compact groups (HCGs, Hickson 1982) of galaxies are systems characterised by a high local density while being located in low-density environments when viewed at larger scales. This high density combined with low-velocity dispersions (Hickson et al. 1992) in many cases leads them to exhibit multiple physical processes associated with galaxy–galaxy interactions: tidal tails and bridges visible optically, in atomic gas (H I), or both (e.g. Verdes-Montenegro et al. 1997, 2005a; Sulentic et al. 2001; Serra et al. 2013; Konstantopoulos et al. 2013); intragroup diffuse X-ray emission (Belsole et al. 2003; Desjardins et al. 2013; O’Sullivan et al. 2014a); shock excitation from starburst winds or galaxy–tidal debris collisions (Rich et al. 2010; Vogt et al. 2013; Cluver et al. 2013); anomalous star formation (SF) activity, molecular gas content, and morphological transformations (Tzanavaris et al. 2010; Plauchu-Frayn et al. 2012; Alatalo et al. 2015; Eigenthaler et al. 2015; Zucker et al. 2016; Lisenfeld et al. 2017), among others.

Single-dish H I studies of HCGs (Williams & Rood 1987; Huchtmeier 1997) revealed that most are deficient in H I. Verdes-Montenegro et al. (2001) expanded on this discovery by performing a comprehensive study of the total H I content of 72 HCGs observed with single-dish telescopes, together with an analysis of the spatial distribution and kinematics of the H I gas within a subset of 16 HCGs observed with the Very Large Array (VLA). As a result of the analysis the authors proposed an evolutionary sequence in which compact group galaxies become increasingly H I deficient as the group evolves. In phase 1 of the sequence the H I is relatively unperturbed and found mostly in the discs of the galaxies, with the remaining gas found in incipient tidal tails. In phase 2, 30–60% of the total H I mass forms tidal features. In phase 3a, most if not all of the H I has been stripped from the discs of the galaxies and is either found in tails or is not detected. The least common phase, 3b, involves groups where the H I gas seems to form a large cloud with a single velocity gradient that contains all the galaxies. However, of the four groups proposed to be in this phase, HCGs 18 and 54 are now thought to be false groups (Verdes-Montenegro et al. 2001, 2002) and the H I distribution in HCG 26 probably does not fulfil the necessary criteria (Damas-Segovia et al., in prep.), leaving only HCG 49 and raising the question of whether phase 3b is a genuine phase of CG evolution. A slightly different evolutionary sequence was proposed by Konstantopoulos et al. (2010), where the evolution follows a similar sequence, but all groups are split into two categories: (a) those where the gas is mostly consumed by SF in the galactic discs before major interactions can strip it, leading to late-time dry mergers, and (b) those where the gas is removed from the galaxies through tidal stripping early on in the evolution of the group, leading to a hot, diffuse intra-group medium (IGrM) at late times.

Borthakur et al. (2010, 2015) compared single-dish H I spectra of HCGs obtained with the Green Bank Telescope (GBT) and VLA H I maps to demonstrate that some HCGs have a diffuse H I component that was not detected by the VLA and can extend to up to 1000 km s⁻¹ in velocity width. The fraction of H I in this component seemed to be greater for groups with larger H I deficiencies, and thus makes up some, but not all, of the “missing” H I reported by Verdes-Montenegro et al. (2001). The connection between the H I content and distribution, SF activity, and X-ray emission has been the subject of numerous studies (e.g. Ponman et al. 1996; Rasmussen et al. 2008; Bitsakis et al. 2011; Martinez-Badenes et al. 2012; Desjardins et al. 2013; O’Sullivan et al. 2014b; O’Sullivan et al. 2014a), however, how the observed H I depletion occurs and more generally how the groups might evolve from phase 2 to 3, remains far from understood.

Assuming the proposed evolutionary scenario is correct, detailed studies of phase 2 groups are of special relevance for addressing the unknowns above, because in these groups the processes that drive the transformation to phase 3 HCGs should be at work. HCG 16 is a prototypical example of this intermediate phase of evolution. Its H I gas is in the process of leaving the discs of the galaxies and filling the intragroup medium with significant amounts of H I in tidal tails, but the group has not yet become globally H I deficient. HCG 16 also hosts an array of other ongoing processes that will likely shape its future evolution: active galactic nuclei (AGNs), a new member, starburst events, and the accompanying winds and shocks. Many of these have been studied in detail in an extensive set of papers focusing on the group or a small sample of groups including HCG 16 (Ribeiro et al. 1996; Mendes de Oliveira et al. 1998; Rich et al. 2010; Vogt et al. 2013; Konstantopoulos et al. 2013; O’Sullivan et al. 2014b; O’Sullivan et al. 2014a), however, to date there has been no study specifically targetting the extraordinary H I component of the group, which is the focus of the current work.

The aim of this paper is to shed light on how the final stages of evolution in HCGs are reached by performing a census of the ongoing physical processes in HCG 16, identifying those that could be influencing the fate of the H I in the group and its evolution towards a phase-3 morphology.

In the following section we give a brief overview of HCG 16 and in Sect. 3 we describe the observations and standard data reduction. Section 4 covers the separation of H I into galaxies and tidal features and in Sect. 5 we present our results for the group as a whole and the individual galaxies. In Sect. 6 we discuss their interpretation and attempt to construct a coherent picture of the evolution of the group. Throughout this paper we assume a distance of 55.2 ± 3.3 Mpc for HCG 16 and all its constituent galaxies.

2. Overview of HCG 16

HCG 16 was first identified in the Atlas of Peculiar Galaxies (Arp 1966), Arp 318, and later classified as a compact group by Hickson (1982). Since then it has been referenced in approximately 100 journal articles and is thus an extremely well-studied group for which there is a large amount of multi-wavelength data. However, this work represents the first focused investigation of its H I component.

The core group contains four disc galaxies, each with stellar mass of the order 10¹⁰–10¹¹ M_⊙, that all fall within a projected separation of just 7′ (120 kpc). There is a fifth, similarly sized member of the group (NGC 848) to the southeast that was identified as being associated with the core group through optical spectroscopy (de Carvalho et al. 1997), and was later shown to also be connected in H I (Verdes-Montenegro et al. 2001). de Carvalho et al. (1997) identified two further dwarf galaxy members of HCG 16, PGC 8210 to the southwest and 2MASS J02083670−0956140 to the northwest. The latter is not considered in this work as it falls outside of the primary beam of the H I observations. The basic optical properties of the others are summarised in Table 1.

In Sect. 5 we discuss the details of each galaxy individually, but here we provide a brief overview of their properties for readers unfamiliar with this group. Figure 1 shows an optical grz image (from the Dark Energy Camera Legacy Survey) of the group. From northwest to southeast the galaxies are HCG 16b, a, c, d, and NGC 848. PGC 8210 is to the southwest of the core group. The four galaxies in the core group form two interacting pairs, HCG 16a and b, and HCG 16c and d. In the first pair, both galaxies host an AGN, but have limited SF activity, while the second pair does not host AGNs and both galaxies are currently undergoing nuclear starburst events. NGC 848 is physically connected to the core group by an enormous H I tail, while PGC 8210 appears quite separate and shows no evidence in H I for a past interaction with the core group.

Fig. 1. DECaLS grz colour image of HCG 16 with the member galaxies labelled.

3. Observations and data reduction

3.1. H I data

HCG 16 was mapped with the VLA in C and D configurations in 1999 and 1989 respectively. These data were reduced using AIPS (Astronomical Image Processing System) and presented in Verdes-Montenegro et al. (2001) and Borthakur et al. (2010). In this work we have re-reduced the raw data using CASA (Common Astronomy Software Applications, McMullin et al. 2007)¹ and re-imaged them using multi-scale CLEAN in the CASA task tclean. The H I line emission was imaged over the velocity range 3246 km s⁻¹ – 4557 km s⁻¹ with a resolution of 21 km s⁻¹. The dataset was imaged twice to generate two cubes using Brigg’s robust weighting parameters of 2 and 0. While almost all of the following analysis relies on the robust = 2 cube, the robust = 0 is useful to see some parts of the highest-column-density gas with finer angular resolution. The multi-scale CLEAN angular scales used in these two cubes are 0, 8, 16, 24, and 40 pixels in the robust = 2 cube, and 0, 4, 8, 16, 24, and 40 pixels in the robust = 0 cube, where each pixel was 4″ across in both cases. The resulting beam sizes of these two cubes were 37.2″ × 30.3″ and 19.4″ × 14.8″ respectively. At the assumed distance of 55.2 Mpc, 20″ corresponds to a projected distance of 5.4 kpc. The robust = 2 and robust = 0 cubes have rms noises of 0.36 mJy beam⁻¹ and 0.40 mJy beam⁻¹ respectively, which correspond to 3σ H I column density sensitivities of 2.2 × 10¹⁹ cm⁻² and 9.6 × 10¹⁹ cm⁻² at a velocity resolution of 21 km s⁻¹. An interactive, 3D figure displaying the robust = 2 cube created using the X3D pathway introduced in Vogt et al. (2016) is available online².

To create a source mask within which the zeroth moment and total integrated flux could be calculated we made use of the SoFiA package (Source Finding in Astronomy, Serra et al. 2014, 2015), using smoothing kernels over spatial scales approximately equal to one and two times the (robust = 2) beam size, over one and three channels, and clipping at 3.5σ_rms (shown with contours in Fig. 2, right panel). A reliability threshold of 100% was set to remove spurious noise spikes, the sources were merged and the final mask dilated using SoFiA’s mask optimisation tools in order to include all flux associated with the group. An equivalent procedure with a threshold of 5σ_rms and without dilating the mask was used to produce the first-moment map of the group H I emission (Fig. 3). In addition we made a more traditional source mask based on a 3σ_rms clipping in each channel (of the robust = 2 cube) using CASA. A comparison of the moments generated by these two masks can be seen in Fig. 2.

Fig. 2. Left: moment-zero map (primary beam corrected) of the H I emission of HCG 16 calculated using a 3σrms mask in each channel (this mask is intended only for visual purposes, the whole analysis uses the SoFiA-generated mask shown in the right panel). The black ellipse in the lower right corner indicates the beam size (37.2″ × 30.3″ or 10.0 kpc × 8.1 kpc) and the small red cross shows the centre of the GBT pointing from Borthakur et al. (2010). The extended features that we separated from the galaxies are indicated by arrows, with the exception of the emission joining galaxies HCG 16c and d. The filled black circle and the white plus sign indicate the pointing centres of the VLA D and C array data respectively. Right: moment -zero contours (uncorrected for primary beam) overlaid on a DECaLS r-band image. In this case the map was generated using the SoFiA-generated mask described in Sect. 3.1, which includes more extended emission. The galaxies in the main band of H I emission going from top-right (NW) to bottom-left (SE) are: HCG 16b, a, c, d, NGC 848, and the single galaxy to the lower right of the group is PGC 8210. Contour levels: −2.45, 2.45, 9.80, 24.4, 49.0, 73.5, and 98.0 × 1019 cm−2, where 2.2 × 1019 cm−2 corresponds to the 3σ sensitivity in one channel. In order of increasing flux the contours are coloured: black (dashed), black, blue, purple, red, orange, and yellow.

The spatial and spectral smoothing performed by SoFiA results in a more extended mask (even though the threshold is 3.5σ_rms rather than 3σ_rms), which should include more of the low-column-density emission. The standard 3σ_rms mask is only used for visual representation as it more clearly separates the higher-column-density features (precisely because it excludes the fainter emission). In all of the following sections and analysis we use the robust = 2 cube and the SoFiA source mask, unless explicitly stated otherwise.

3.2. Optical images

Throughout this paper we compare H I features with optical images from the Dark Energy Camera Legacy Survey (DECaLS³). The three DECaLS bands (g, r, and z) have surface brightness limits in the field of interest of 28.5, 28.7, and 28.0 mag arcsec⁻² (3σ in 10″ × 10″ boxes), respectively. This is the deepest image that covers the entire field of which we are aware; therefore, we focus on this image to look for faint optical features which may be associated with extended H I features. However, Hubble Space Telescope and Spitzer images of the group have been published by Konstantopoulos et al. (2013).

We used the DECaLS images as published, which were processed with the automated Dark Energy Camera Community Pipeline. This processing slightly over-subtracts the sky in the vicinity of large galaxies, negatively impacting the sensitivity for large-scale faint features near large galaxies (like those in HCG 16). As this work is focused on the H I component of the group, we note this issue, but do not reprocess the images.

4. Separation of H I features

The H I content of HCG 16 is enormously complicated, with multiple blended galaxies and tidal features. Therefore, in order to study the properties of each galaxy and tidal feature, the H I emission must be separated into these individual objects wherever possible. By doing this we can assess whether each galaxy has a high, normal, or low H I content (irrespective of the global group H I content), if they have well-ordered rotation and morphology or have been disrupted by interactions, and where the gas in tidal features likely originated. The answers to these questions will in turn help to build a consistent picture of how the group has evolved to date and how it might proceed in the future.

To perform this separation we used the SlicerAstro package (Punzo et al. 2016, 2017). The SoFiA mask was imported into SlicerAstro and divided into sub-regions corresponding to galaxies or specific tidal features. By importing the same mask we ensure that all of the flux included in the integrated measurement is assigned to a particular galaxy or feature. This is especially important for obtaining the H I masses of individual galaxies. Although it is more straightforward to separate the higher-column-density regions into distinct sources, this would lead to the H I flux being underestimated as low-column-density emission would be systematically missing, while the scaling relations we use to define the deficiency were calculated from single-dish observations that include all the H I flux of target (isolated) galaxies. The final separation was made through an iterative comparison of the channel maps, optical images, 3D visualisations, and moment maps of the individual galaxies. This process was unavoidably subjective to some degree, particularly in the region of HCG 16c and HCG 16d where emission smoothly changes from gas associated with the galaxies to multiple high-column-density tidal features. However, despite the resulting large uncertainties in the galaxy parameters this is still an instructive tool for assessing the likely history of the group, as discussed in the following sections. More specific details of this process are described in Sect. 5.3, along with the results for individual galaxies once they have been separated from surrounding features.

5. Results

In this section we present the results of our H I analysis, first for the group as a whole and then for each galaxy.

5.1. Global H I morphology

Figure 2 (left) shows the zeroth-moment map of the H I emission created using a standard 3σ_rms clipping in each channel. This map excludes much of the lowest-column-density emission, making many features easier to distinguish by eye. The galaxies HCG 16a, b, c, d, NGC 848, and PGC 8210 are all detected along with tidal features across the whole group that appear to connect HCG 16a and b to HCG 16c, HCG 16c to d, and HCG 16d to NGC 848. The most striking tidal feature is the southeast tail, which stretches over a projected distance of ∼160 kpc, connecting NGC 848 to the core group. There are no indications of an interaction between PGC 8210 and the rest of the group.

Figure 2 (right) shows the zeroth moment again, but this time made using the SoFiA-generated mask (Sect. 3.1) and displayed as contours overlaid on the DECaLS grz image. This mask includes more of the low-column-density emission than the previous mask and thus provides a more complete measurement of the H I content, but makes most individual features more difficult to identify by eye. This masking is used throughout the following analysis as we want to include low-column-density emission because a large fraction of the group’s H I is in tidal features.

In Fig. 3 we show the first-moment map of the entire group (using a SoFiA mask with a 5σ_rms threshold). Due to the high spatial density of the group, the galaxies and tidal features are superposed and confused in this image, complicating its interpretation. Without separating out individual objects and features, clear signs of ordered rotation are only visible in NGC 848 and PGC 8210, although in the latter cases it is barely larger than the beam. The velocity field of each galaxy is presented in Sect. 5.3 after separating them from the rest of the emission in the group.

Fig. 3. First-moment map of the H I emission calculated using the SoFiA-generated mask with a threshold of 5σrms and without dilation. The rainbow colour scale indicates the recessional (optical) velocity in km s−1. The black ellipse in the lower right indicates the beam size, the grey lines are isovelocity contours separated by 40 km s−1, and the small black star symbols indicate the locations of the optical centres of the galaxies in the group.

5.2. H I flux and mass

Integrating the entire H I emission shown in Fig. 2 (right) gives the total H I flux of HCG 16 as 47.2 Jy km s⁻¹ and its mass as log M_HI/M_⊙ = 10.53 ± 0.05 (using the standard formula, M_HI/M_⊙ = 235 600 × [D/Mpc]² × [S_HI/Jy km s⁻¹]). This is considerably higher than the values of Borthakur et al. (2010), 21.5 Jy km s⁻¹ and 10.19 ± 0.05 × 10¹⁰ M_⊙, based on the single-dish spectrum taken with GBT. However, as pointed out in that work, this difference arises because the GBT HPBW is 9.1′ and the pointing was centred on the group core (red x in the left panel of Fig. 2). This means that the majority of the flux from the tail extending to the SE, NGC 848, and PGC 8210 is missing from the GBT spectrum. The moment-zero map was weighted with a 2D Gaussian window (FWHM of 9.1′) centred on the GBT pointing centre to estimate the fraction of the total emission that the GBT observations would have been able to detect. This gives an H I mass of log M_HI/M_⊙ = 10.24 ± 0.05 M_⊙ and a flux of 24.2 Jy km s⁻¹, which is about 10% higher than the GBT measurement (Fig. 4). This slight difference could arise from calibration or baseline uncertainties, or simply because a Gaussian is not a completely accurate representation of the beam response.

Fig. 4. Integrated VLA spectrum of the entire group (thick black line) compared to the GBT spectrum of the core group (Borthakur et al. 2010, high spectral resolution green line), the VLA spectrum within the GBT beam area (thick black dashed line), and the spectrum extracted from the HIPASS cube (orange line). Vertical dashed lines indicate the optical redshifts of the member galaxies.

We compare our total flux measurement (47.2 Jy km s⁻¹) to that in the H I Parkes All Sky Survey (HIPASS, Barnes et al. 2001; Meyer et al. 2004). As HCG 16 is extended over many arcmin it is a marginally resolved source even for the Parkes telescope. Therefore, we cannot use the HIPASS catalogue values, which assume it to be a point source. Using the HIPASS cube in this region of the sky we perform a source extraction with SoFiA, applying a threshold limit of 4.5σ over smoothing kernels of 3 and 6 spatial pixels and 3 and 7 velocity pixels (each pixel is 4′ × 4′ × 13 km s⁻¹). The resulting spectrum is shown in Fig. 4. The integrated flux in HIPASS is 43.8 Jy km s⁻¹, which is within 10% of our measurement with the VLA. Given the difficulty in absolute calibration (e.g. van Zee et al. 1997), this is about the level of agreement that is to be expected. However, there are some additional discrepancies with this comparison which we discuss further in Appendix B.

5.3. Individual galaxies

Table 2 shows the basic H I properties of the individual galaxies and tidal features based on the separation performed in Sect. 4. We have chosen to use the flux-weighted mean velocity and flux-weighted velocity dispersion to quantify the centre and width of each emission profile, rather than more common measures such as for example W₅₀, because some of the profiles, particularly those of the strongly disturbed galaxies or tidal features, do not follow typical profile shapes of H I sources (Figs. 5 and 6). The measurements of the B-band luminosity, L_B, of each galaxy (Table 1) were used as inputs to the H I scaling relation of Jones et al. (2018) in order to estimate their expected H I masses if they were in isolation, and thus their H I deficiencies. We note that although it has been fairly common in past works to consider the morphology of a galaxy when calculating H I deficiency, here we choose to ignore it. There are two main reasons for this. Firstly, Jones et al. (2018) suggest that morphology should be ignored unless the sample has a large fraction of ellipticals, because their piece-wise scaling relations (split by morphology) are quite uncertain due to the small number of galaxies that are not Sb-Sc in the Analysis of the interstellar Medium of Isolated GAlaxies (AMIGA; Verdes-Montenegro et al. 2005b) reference sample. Therefore, using these piece-wise relations in the case of HCG 16 would trade a small bias for a large uncertainty that is difficult to accurately quantify. Secondly, galaxies are expected to undergo morphological transformation in compact groups as they are stripped of their gas and perturbed by interactions. Therefore, it is perhaps misguided to use their present-day morphologies in these scaling relations.

Fig. 5. Spectral profiles of each of the six galaxies detected in H I near the core group calculated from the separation of tidal features and galaxies performed in Sect. 3.1. The black vertical dashed lines show the optical redshifts of each galaxy as given in Table 1. The dot-dashed and dotted vertical lines show the central and extreme velocities, respectively, from the rotation curve measurements of Rubin et al. (1991) in red and Mendes de Oliveira et al. (1998) in blue. These measurements are only available for the four core galaxies; also Rubin et al. (1991) do not specify Vmax values for HCG 16b or c due to the peculiar shape of their rotation curves. The left panels have vertical scales that go to 11 mJy, whereas the scales on the right are a factor of approximately six higher (60 mJy).

Fig. 6. H I spectral profiles of each of the eight separate tidal features in the group calculated from the separation of tidal features and galaxies performed in Sect. 3.1. Each row has a different vertical scale.

Our analysis indicates that only galaxies HCG 16a and b are missing neutral gas compared to their expected H I content if they were in isolation, and that all the remaining galaxies are consistent with having a normal quantity of H I. The group as a whole (galaxies and extended features) has an H I mass of (3.39 ± 0.29)×10¹⁰ M_⊙ (accounting only for distance uncertainty) and a global H I deficiency of −0.12 ± 0.09⁴. Thus the group as a whole does not appear to be missing H I gas (it is marginally H I-rich), but the galaxy pair HCG 16a and b has likely lost the vast majority of its original H I content. The probable fate of this lost gas is discussed in the following sections.

Here we note that had we used the Haynes & Giovanelli (1984) scaling relation instead of the updated Jones et al. (2018) relation, then we would have found the global H I deficiency to be 0.05, which is approximately 2σ higher (more deficient). The differences between these relations are discussed in detail in Jones et al. (2018), but in this case the most important point is that the uncertainties in the measurements of L_B for the reference sample used to calibrate the relation results in the Haynes & Giovanelli (1984)M_HI–L_B scaling relation overestimating the H I deficiency of galaxies.

The uncertainty estimates for the H I deficiencies in Table 2 are dominated by the scatter in the M_HI–L_B scaling relation and also have small contributions due to the uncertainty in L_B and the group distance, 55.2 ± 3.3 Mpc, which was estimated using the Mould et al. (2000) local flow model as described in Jones et al. (2018), which has corrections for flow towards the Virgo cluster, the Great Attractor, and the Shapley supercluster. The uncertainties on the masses of each component are challenging to estimate because these values strongly depend on where subjective boundaries are drawn between tidal gas and gas in a galactic disc. As H I deficiency is the quantity of interest for the galaxies and the scatter in the scaling relation (0.2 dex) results in an uncertainty of 40−60%, which is expected to dominate the error budget, we do not attempt to estimate the H I mass uncertainty due to the subjective measurements. The same issue applies to the tidal features, but is more severe as they are mostly made up of lower-column-density gas, and therefore their individual mass measurements should be treated with caution.

5.3.1. HCG 16a

HCG 16a is an Sab galaxy that has been classified as an active star forming galaxy in infrared (IR; Zucker et al. 2016) and hosts an AGN (Seyfert 2) that has been confirmed with both optical line ratios (Véron-Cetty & Véron 2010; Martínez et al. 2010) and X-rays (Turner et al. 2001; Oda et al. 2018). The galaxy also hosts a ring of SF (star formation) that is clearly visible in the IR and UV (Tzanavaris et al. 2010; Bitsakis et al. 2014). In terms of its H I, HCG 16a is the second most deficient galaxy in the group, apparently missing approximately 5.5 × 10⁹ M_⊙ of H I.

The discs of HCG 16a and b overlap in the plane of the sky which complicates the separation of their H I emission. However, in velocity the gas that is co-spatial with HCG 16a forms a clear, continuous, and steep gradient, as would be expected from an inclined disc galaxy. This is most readily seen in the X3D interactive figure which accompanies this paper. The emission associated with HCG 16b forms a small, offset clump at approximately the same velocity as the lower-velocity “horn” of the HCG 16a profile (Fig. 5). The two galaxies are connected in H I by a faint bridge of emission, visible in the channel maps (Figs. D.1 and D.2) between velocities 3858 and 4006 km s⁻¹. Given the size of the beam and the fact that the optical discs of HCG 16a and b overlap, it is not possible to reliably assign the emission to either galaxy. Therefore, we simply split it approximately half way between the two sources. The uncertainty in where this bridge should be split introduces a negligible error to the H I properties of HCG 16a as the entirety of the emission that we attribute to HCG 16b would only increase the H I mass of HCG 16a by 0.13 dex.

When viewed projected on the plane of the sky there is an H I tail which connects HCG 16a to HCG 16c (the NW tail, see Fig. 2, left panel). There is also an accompanying stellar tail (identified by Rubin et al. 1991) that extends in the same direction away from HCG 16a (clearly visible in the Figs. 1 and 7). However, when the H I data cube is studied in 3D it is clear that the H I tail does not form a kinematic connection with HCG 16a and the apparent connection is the result of a projection effect. In Fig. 7 it can be seen that as the NW tail is traced away from HCG 16c its velocity decreases from ∼3800 km s⁻¹ to ∼3500 km s⁻¹, whereas the H I profile of HCG 16a covers the approximate range 3850−4300 km s⁻¹ (right panel of Figs. 5 and 8). This tail then ends in a dense clump, which may have formed a tidal dwarf galaxy (TDG), discussed further in Sect. 5.4.

Fig. 7. Top: segmented position–velocity diagram showing the spatially and kinematically continuous H I emission spanning HCG 16. The vertical dashed lines show the locations of the nodes making up the segmented slice through the data cube. Galaxies are labelled in red and tidal features in blue. We note that the noise in this plot increases significantly near the left edge as this region is approaching the primary beam edge. Bottom, left panel: position–velocity slice as a segmented red line on top of the moment-zero map (using the standard 3σ threshold mask as in the left panel of Fig. 2) and right panel: same line but overlaid on the DECaLS r-band image.

Figure 8 shows the zeroth and first moments of the H I emission. From these two maps it can be seen that the high-column-density gas in the immediate vicinity of the optical disc of HCG 16a is relatively undisturbed, despite the large amount of missing H I. The moment-zero map has a regular oval shape and is centrally peaked (almost coincident with the optical centre). While the velocity field appears to be relatively regular in the centre, the line of nodes forms an “S” shape in the outer regions, likely indicating the presence of a warp in the disc (e.g. Bosma 1978). Given the stellar tails in the optical and the clear ongoing interaction with HCG 16b, such an asymmetry is not unexpected.

Fig. 8. Left: moment-zero map of HCG 16a. The white and black stars indicate the optical centres of HCG 16a and b respectively, and the black ellipse shows the beam. Right: first-moment contours for HCG 16a overlaid on the DECaLS grz image of HCG 16a and b. The green lines are separated by 20 km s−1.

Mendes de Oliveira et al. (1998) also found that the rotation in the inner region of HCG 16a is very regular, rising quickly and becoming flat well within one of our beam widths. The velocity extent of our H I spectrum (Fig. 5) approximately agrees with that implied by their Hα rotation curve. However, Rubin et al. (1991) found the conflicting result that the Hα rotation has an anomalous structure, declining on one side of the galaxy at large radii. As our H I map extends to much larger radii we would expect to see a continuation of this in our data, but we do not. Therefore, the H I velocity field appears more consistent with the rotation curve of Mendes de Oliveira et al. (1998). Nevertheless, the major axis position angle given in that work is almost exactly aligned N-S, whereas in H I it is clearly rotated (counter-clockwise). We expect this is due to a combination of beam smearing and disturbances in the outer regions of the disc.

5.3.2. HCG 16b

HCG 16b is an Sa galaxy and the most H I -deficient galaxy in the group, probably having lost about 4.6 × 10⁹ M_⊙ of its H I⁵. The galaxy has a central source identified as an AGN with X-rays and line ratios (Turner et al. 2001; Martínez et al. 2010; Oda et al. 2018), it has also been classified as a LINER (Véron-Cetty & Véron 2010). This is the only galaxy in the group that was classified as quiescent in terms of its WISE IR colours (Zucker et al. 2016).

As described above, HCG 16b has an H I bridge connecting to the much more H I -massive HCG 16a. In the absence of more information, the faint bridge is simply split approximately half-way between the centres of emission of each object in each channel. One may argue that this bridge should be classified as a separate tidal feature. However, given the fact that the outer regions of the optical discs of the two galaxies blend together and taking into account the resolution of the H I cube, this is not a practical suggestion. While it is difficult to quantify the resulting uncertainty in the H I mass of HCG 16b due to the simplistic separation, it is likely that the procedure assigns too much emission to HCG 16b rather than HCG 16a, because the latter is both more optically luminous and has a higher H I mass, and therefore more of the bridge is probably gravitationally bound to it. Therefore, the key result that HCG 16b has lost around 90% or more of its expected H I content would be unchanged.

Figure 9 shows that the little remaining H I in HCG 16b is very off-centre compared to its optical disc. The velocity field shows a small but clear velocity gradient aligned with the major axis of the optical disc, suggesting that at least some of the remaining H I is likely still rotating with the optical disc. However, this gradient is disrupted on the eastern side of the galaxy where there is a gas extension to the NE that connects to emission around the NW tail and that around HCG 16a.

Fig. 9. Left: moment-zero map of HCG 16b. The black stars indicate the optical centres of HCG 16a and b, and the black ellipse shows the beam. Right: first-moment contours for HCG 16b overlaid on the DECaLS grz image of HCG 16a and b. The green lines are separated by 20 km s−1. The signature of a small, consistent velocity gradient across the optical disc is visible, but the H I distribution is strongly off centre and disturbed.

For HCG 16b there is dramatic disagreement between the rotation curves of Rubin et al. (1991) and Mendes de Oliveira et al. (1998), with them even disagreeing on the direction of rotation. The direction that we see in H I agrees with that of Rubin et al. (1991), with redshifted emission occurring to the NE of the galaxy centre and blueshifted to the SW. If we take the Rubin et al. (1991) systemic velocity then it appears that we are mostly detecting gas on the redshifted side of the disc, with little contribution from the blueshifted side. This would also imply that the entire velocity gradient seen in Mendes de Oliveira et al. (1998) is on the redshifted side of the galaxy.

5.3.3. HCG 16c

HCG 16c is classified as an S0-a galaxy and a luminous infrared galaxy (LIRG) and is the most actively star-forming galaxy in the group. It has been classified as both a pure starburst based on X-ray observations (Turner et al. 2001) and a composite object with both AGN and SF activity using optical line ratios (Martínez et al. 2010). The overall morphology of the galaxy is highly reminiscent of M82, with the nuclear starburst driving a bipolar galactic wind that was studied in detail by Vogt et al. (2013).

HCG 16c is apparently involved in a number of ongoing interactions: the NW tail extends between HCG 16c towards HCG 16a, the NE tail extends from HCG 16c away from the rest of the group, and there is an H I bridge connecting HCG 16c and d, all of which can be seen in Fig. 2, the X3D plot, and in the channel maps (Figs. D.1 and D.2). All three of these features contain more than 10⁹ M_⊙ of H I, that is, they are comparable in H I mass to the galaxy itself, yet only considering the H I emission in the region we identified as the H I disc of HCG 16c the galaxy is not found to be H I deficient. This strongly suggests that the majority of the neutral gas surrounding HCG 16c likely originated elsewhere (discussed further in Sect. 6).

We trace the NW tail from the central velocity and position of HCG 16c extending towards HCG 16a in the plane of the sky, but away from it in velocity space, ending when separated from HCG 16a by over 300 km s⁻¹. The NE tail appears to originate from the receding edge of HCG 16c and extends out to the NE, but then curves back towards HCG 16a, overlapping both spatially and in velocity with the NW tail. In this region there is a great deal of extended emission of uncertain origin and it is impossible to reliably separate the two features anywhere except at their bases where the column density is high. Therefore, we assign emission at different velocities that is mostly co-spatial (on the plane of the sky) with the bulk of each feature. This means the fluxes of these two tails should be treated with great caution, however the measurement of the overall flux of the extended emission is not impacted by this issue. Finally, there is the cd bridge, which is predominantly made up of the high-column-density H I emission that forms a bridge between galaxies HCG 16c and d. However, all the remaining emission in the vicinity of galaxies HCG 16c and d that had not been assigned to any of the galaxies or other features listed above was also assigned to the cd bridge and therefore it is more poorly defined than the other extended features we describe.

Despite the evidence for interactions listed above, the centre of the H I disc of HCG 16c still has a consistent velocity gradient across it, indicating that it is rotating and has not been completely disrupted (Fig. 10). However, the H I distribution is clearly extended (asymmetrically) in the direction of HCG 16d, and the outer regions of the velocity field also trail off in that direction, indicating that there are major disturbances in the outer parts of the disc.

Fig. 10. Left: moment-zero map of HCG 16c. The white and black stars indicate the optical centres of HCG 16c and d respectively, and the black ellipse shows the beam. Right: first-moment contours for HCG 16c overlaid on the DECaLS grz image of HCG 16c and d. The green lines are separated by 20 km s−1.

The Rubin et al. (1991) and Mendes de Oliveira et al. (1998) rotation curves are again quite different for this source, but this is because the spectral slit that Rubin et al. (1991) used was aligned with the R-band photometric major axis, which is rotated about 40° with respect to the kinematic axis identified by Mendes de Oliveira et al. (1998). The kinematic major axis of the H I data (Fig. 10, right) appears approximately consistent with this position angle (120°) obtained from the Hα line. The Mendes de Oliveira et al. (1998) velocity field is more or less regular, but the range of velocities seen is considerably larger than we see in H I. This may indicate that the rotation curve declines beyond the central region that they measure, or that their rotation curve was contaminated by Hα associated with the central outflow, as suggested by Vogt et al. (2013).

Figure 7 shows the peculiar kinematics of the NW tail. Rather than connecting to the outer edge of the disc of HCG 16c, as is typical for tidal tails (in fact this is visible on the opposite side of the same plot, around NGC 848), the NW tail appears to intersect HCG 16c at its central velocity. This behaviour leads us to consider three competing hypotheses which we will discuss in the following section: (1) the gas is being accreted directly to the central regions of the galaxy and is providing fuel for the nuclear starburst, (2) the gas is being ejected from the central regions of the galaxy by the nuclear starburst, or (3) the entire feature is superposed with HCG 16c and the coinciding velocities do not correspond to a spatial connection.

5.3.4. HCG 16d

Similarly to HCG 16c, HCG 16d is a LIRG, classified as S0-a, and is morphologically similar to M 82. It contains a central source classified as both a LINER and as a LINER/Seyfert 2 double nucleus (Ribeiro et al. 1996; de Carvalho & Coziol 1999) with optical line ratios, or as an AGN in X-ray observations with XMM (Turner et al. 2001). However, HST observations show that the suggested double nucleus is instead a group of star clusters (Konstantopoulos et al. 2013) and Chandra X-ray observations indicate that the hard X-rays might originate solely from X-ray binaries and not from an AGN (O’Sullivan et al. 2014b; Oda et al. 2018), while IFU observations suggest that the LINER line ratios in the optical may arise from shock excitation driven by the ongoing starburst (Rich et al. 2010).

Upon first inspection, the integrated H I emission (Fig. 11) appears almost like a face-on galaxy with a central H I hole; however this is quite misleading and does not agree with the optical image, which shows a highly inclined disc. The apparent central hole is instead H I absorption in front of a central continuum source. The H I depression is of the same shape as the beam in both the robust = 2 and robust = 0 cubes (i.e. at two different resolutions) and the spectral profile at that position switches from a flux density of approximately 1.5 mJy to −1.5 mJy at the central velocity of the galaxy (Fig. 12).

Fig. 11. Left: moment-zero map of HCG 16d. The white and black stars indicate the optical centres of HCG 16d and c, respectively, and the black ellipse shows the beam. A central depression is clearly visible around the white star marking the centre of HCG 16d. Right: first-moment contours for HCG 16d overlaid on the DECaLS grz image of HCG 16d and c. The central region of the galaxy is removed as the H I absorption feature would contaminate the first-moment map there. The green lines are separated by 20 km s−1. While the velocity field is highly irregular, there is a slight gradient roughly aligned with the minor axis of the optical disc.

Fig. 12. H I spectral profile of the centre of HCG 16d. The dashed vertical line shows the systemic velocity as determined from optical lines. The profile was extracted from the robust = 0 cube (higher resolution) by placing a double Gaussian window at the position of the optical centre of the galaxy with the shape and orientation of the synthesized beam.

This absorption feature is redshifted relative to the central velocity (from optical spectra) by about 100 km s⁻¹. A regular rotating disc would form a symmetric absorption feature about the central velocity (Morganti & Oosterloo 2018), so the fact that this feature is redshifted suggests it could be gas falling towards the centre and fuelling the starburst event. However, a 100 km s⁻¹ velocity shift is comfortably within what might be expected from orbiting gas in a galaxy of this size. Given the resolution of the data it is not possible to distinguish between these two possibilities, or indeed that the absorption may be due to an intervening clump of stripped gas in the group.

The absorption also means that the integrated flux for HCG 16d will be underestimated. To estimate how much the integrated flux is reduced we made a crude linear interpolation across the absorption feature between 3879 and 4006 km s⁻¹. This gives the missing flux as 0.1 Jy km s⁻¹, which is barely more than 1% of the total H I emission flux of HCG 16d. Indeed, the absorption feature is not apparent in its integrated spectrum (Fig. 5). Therefore, we ignored this feature for our analysis of the H I deficiency.

When considering the kinematics of the gas in HCG 16d there is an apparent contradiction between H I and Hα. In H I there is only a weak velocity gradient in approximately the north–south direction, almost perpendicular to the optical major axis (Fig. 11). In contrast, the existing Hα rotation curve of the inner regions of the galaxy (Rubin et al. 1991; Mendes de Oliveira et al. 1998; Rich et al. 2010) shows significant rotation (V_rot ≈ 100 km s⁻¹) and has its major axis almost aligned with the optical disc, that is, its axis is approximately perpendicular to the slight gradient in the H I emission. The full extent of the Hα velocity field in Mendes de Oliveira et al. (1998) is about 15″, which is less than the beam size of the H I data, and therefore it is possible that the gas could undergo a major kinematic warp (or other distortion) beyond this region, meaning the H I and Hα results are not necessarily in conflict. However, as the Hα kinematic major axis is aligned with the optical disc of the galaxy (which extends to much larger radii) this seems unlikely.

The H I data thus demonstrates two key points. The first being that there is a continuum source at the centre of HCG 16d, although the H I data offer no information on its nature,. Secondly, the H I gas appears to be completely kinematically disconnected from the optical disc of the galaxy. These points raise the question of whether the H I is truly associated with the galaxy or is simply superposed on it, which we discuss in the following section.

There are also a number of tidal features in the vicinity of HCG 16d, most notably the SE tail that extends for about 10′ (160 kpc) towards NGC 848. Although this feature connects HCG 16d and NGC 848 in both velocity and the plane of the sky (Fig. 7) there is no evidence for an accompanying optical tail emanating from HCG 16d. Most of the emission from the tail is contained within only four channels (a range of 84 km s⁻¹), although there are some clumps of emission which we have assigned to this feature near HCG 16d that extend further in velocity space such that the entire profile of the SE tail covers ten channels (Fig. 6). The small velocity range along the length of the tail likely indicates that it is probably almost aligned with the plane of the sky. HCG 16d is also accompanied by two small, dense clumps on the east and south, and at approximately the same velocity as HCG 16d. Each is several times 10⁸ M_⊙ in H I mass. These clumps are discussed further in Sect. 5.4. There is also a high-H I -column-density bridge connecting HCG 16d and c. The complexity of this region of the group makes separation of these features highly subjective and thus their H I properties should be treated with caution.

5.3.5. NGC 848

NGC 848 is a barred spiral galaxy (SBab) that is separated from the rest of the group by about 10′, or around 160 kpc, and is another galaxy in the group undergoing a starburst (Ribeiro et al. 1996; de Carvalho & Coziol 1999). Although this galaxy was not included in the original Hickson catalogue, once spectra were obtained for the galaxies it was immediately noticed that NGC 848 was at the same redshift as the core group (Table 1) and therefore likely associated (Ribeiro et al. 1996). However, it was not until there were H I interferometric observations that the physical connection was discovered in the form of a 160 kpc tidal tail (Verdes-Montenegro et al. 2001).

Despite its connection to this enormous tidal feature (visible in Fig. 2) the galaxy itself has an entirely normal global H I content (Table 2). This alone indicates that the gas in the tail probably did not (for the most part) originate from NGC 848, but from somewhere else in the group, especially as the total H I mass in the tail is approximately the same as the H I mass of NGC 848. There is another tail connected to NGC 848, which appears to originate from its receding side and extends to the southeast; we refer to this as the NGC 848S tail (Fig. 2). At least part of this feature is included in the SoFiA mask, however it may extend much further, looping back around towards the core group. This is at low signal-to-noise ratio however, and (if real) the emission in adjacent channels also shifts in position, such that smoothing spatially and in velocity does little to improve its signal-to-noise ratio. Some parts of this faint feature are visible in the left panel of Fig. 2 (as well as the X3D plot) slightly to the north of NGC 848 and the SE tail. We include this extended loop in Table 2 and Fig. 6, but it is not included in the total integrated emission of the group as it does not have sufficient significance to be included in the SoFiA mask.

The H I distribution in the galaxy itself (Fig. 13) appears mostly regular, with a centroid coincident with the optical centre and a mostly uniform velocity field. However, the kinematic major and minor axes are not quite perpendicular, which is probably due to the presence of the bar (e.g. Bosma 1981). The H I gas is also slightly extended towards the NW and SE, that is, in the direction of the tidal tail, while in the optical image the spiral arms appear very loosely wrapped (Fig. 13), probably indicating that the outer disc is beginning to be unbound.

Fig. 13. Left: moment-zero map of NGC 848. The white star indicates the optical centre, and the black ellipse shows the beam. Right: first-moment contours for NGC 848 overlaid on the DECaLS grz image. The green lines are separated by 20 km s−1.

5.3.6. PGC 8210

The final galaxy detected within the primary beam of the VLA observations is PGC 8210. This is another spiral galaxy but it is considerably smaller and less massive than the core members of the group. Its B-band luminosity is almost 1 dex lower than the galaxies in the core group and its H I mass is considerably less than all but HCG 16b, which is extremely H I deficient. Because the velocity of PGC 8210 is coincident with that of the group, it is assumed to be part of the same structure, or at least about to join it. However, as shown in Fig. 14 the H I distribution appears undisturbed (with the caveat that it is hardly larger than the beam), as does its optical disc, and its total H I content is entirely normal. Therefore, it seems very unlikely that this galaxy has had any meaningful interaction with the group to date.

Fig. 14. Left: moment-zero map of PGC 8210. The white star indicates the optical centre, and the black ellipse shows the beam. Right: first-moment contours for PGC 8210 overlaid on the DECaLS grz image. The green lines are separated by 20 km s−1.

5.4. Tidal dwarf galaxy candidates

We have identified three dense clumps of H I emission in the group that are not associated with any of the main galaxies. The first two are listed in Table 2, the E and S clumps in the vicinity of HCG 16d, and the third is a dense clump at the northwestern end of the NW tail, which is clearly seen as an excess at low velocities in the spectrum of the NW tail (Fig. 6). These could be candidate low surface brightness (LSB) dwarf galaxies or TDGs, both of which can occur in galaxy groups and can be rich in H I (Lee-Waddell et al. 2012, 2016; Román & Trujillo 2017; Spekkens & Karunakaran 2018), or transient clumps of stripped gas. Their moment-zero and first-moment maps are shown in Fig. 15 (NW, S and E clumps, left to right).

Fig. 15. Top: moment-zero maps of the NW, E, and S clumps (left to right) overlaid on the DECaLS grz image. Contour levels: 0.98, 1.47, 1.96, 2.45 × 1020 cm−2. Bottom: first-moment maps of the same features with the beam shown as a black ellipse in the corner.

All three H I clumps have masses well above 10⁸ M_⊙, the minimum mass thought to be needed to form a long-lived TDG (Bournaud & Duc 2006). Their masses are 3.9, 5.2, and 9.0 × 10⁸ M_⊙, respectively (going left to right in Fig. 15). While the S clump is the most massive (in H I) it is also the most spatially extended, has the lowest peak column density, shows little evidence for rotation in its first-moment map, and has no apparent optical counterpart. Therefore, despite its large mass it seems unlikely that the S clump is a long-lived TDG or a LSB galaxy, and instead is probably a transient feature associated with the SE tail.

The E clump is in between the other two clumps in terms of H I mass and column density, however its first-moment map shows a clear velocity gradient across it indicating that it may be rotating. Taking a very approximate H I diameter of 1.5′ and a rotation velocity of 25 km s⁻¹, the dynamical mass (assuming equilibrium) would be 1.7 × 10⁹ M_⊙, which is about three times its measured H I mass. This simple estimate implies that some dark matter component may be required to explain the velocity gradient if it is due to rotation. However, there are also a number of reasons why this estimate might be inaccurate. For example, the assumption of dynamical equilibrium is unlikely to be correct and even if it were, the source is only three beams across, which means both its velocity field and spatial extent are likely to be heavily affected by beam smearing.

The NW clump has the smallest H I mass of the three clumps, but the highest column density, peaking at over 2.45 × 10²⁰ cm⁻². The first-moment map does not show clear evidence for rotation, however it is barely two beams across and any potential signature of rotation may be confused with gas in the NW tail. Tidal dwarf galaxies are expected to be found at the tip of tidal tails, as this clump is (e.g. Bournaud et al. 2004), and unlike for the other two clumps there is a clear but very LSB-optical counterpart visible in the DECaLS image at RA = 2h 9m 26s, Dec = −10° 7′ 9″ (top left panel of Fig. 15 and right panel of Fig. 8) and it is also just visible in the Galaxy Evolution Explorer (GALEX, Martin et al. 2005) survey. The blue colour of this optical counterpart suggests that it is the result of in situ SF, not old stars which might have been stripped along with the gas. This strengthens the case for this candidate TDG, but there remains the possibility that this is a gas-rich LSB dwarf which is having its H I stripped, rather than a new galaxy that has formed as a result of the tidal interactions in the group. While with the current information it is not possible to rule out either of these hypotheses, we favour the TDG interpretation because this dwarf appears at the end of an enormous tidal feature and it seems implausible that such a small galaxy stripped this gas from the other (much more massive) group members or that it could have originally been sufficiently gas-rich to be the origin of the observed tidal gas.

6. Discussion

6.1. Gas consumption timescales

One way that a galaxy can become gas deficient is by converting its gas into stars without replenishment. Ideally the detailed star formation history (SFH) of a galaxy would be compared to its present day gas content, but in the absence of a SFH the gas consumption time given by the current star formation rate (SFR) can be used instead. If the galaxy is in the middle of converting its gas reservoir into stars then this should be apparent.

To estimate the gas consumption timescales of the cold gas in the galaxies, we took stellar masses and SFRs from Lenkić et al. (2016) and molecular gas mass measurements from Martinez-Badenes et al. (2012). Lenkić et al. (2016) use Spitzer IRAC (Infrared Array Camera) 3.6 and 4.5 μm photometry to estimate stellar masses based on the prescription of (Eskew et al. 2012). To estimate SFRs these authors use the UV (2030 Å, from Swift) and IR (24 μm) luminosities as proxies for the unobscured and obscured (re-radiated) emission due to SF, and combine the two to estimate the total SFR, which in theory allows them to avoid correcting for internal extinction. Martinez-Badenes et al. (2012) use CO data from Boselli et al. (1996), Leon et al. (1998), and Verdes-Montenegro et al. (1998) to estimate molecular gas masses and use a standard constant value of the CO-to-H₂ conversion factor, X = N_H2/I_CO = 2 × 10²⁰ cm⁻² K⁻¹ km⁻¹ s. As all the CO observations are single pointings, Martinez-Badenes et al. (2012) also extrapolate for emission beyond the primary beam by assuming a pure exponential disc with a scale length of 0.2r₂₅, where r₂₅ is the major optical 25 mag arcsec⁻² isophotal radius.

Table 3 shows the gas consumption timescales of the four core galaxies of HCG 16, estimated by combining our measurements of the H I masses with measurements of the H₂ masses (Martinez-Badenes et al. 2012) and the SFRs (Lenkić et al. 2016). The total hydrogen mass is multiplied by 1.3 to account for all other elements, giving T_con = 1.3(M_HI + M_H2)/SFR. The galaxies are roughly split into two categories in terms of their gas consumption time: those that at their current SFR will exhaust their existing gas reservoirs in about a gigayear (HCG 16c and d) and those that will take several gigayears to do the same (HCG 16a and b).

HCG 16c and d are both starbursting LIRGs, so it is unsurprising that their gas consumption timescales are short. It is tempting to think that HCG 16a and b probably looked much the same approximately a gigayear in the past and now have slowed SFRs and have become gas deficient. However, as the group is not globally H I deficient it seems unlikely that their atomic gas was converted to H₂ and consumed in SF, unless the group was gas-rich to begin with. However, even if we allow for that possibility this scenario does not seem to agree with their gas and stellar content. To become depleted in H I but not in H₂ via SF would have required them to have undergone a SF episode, which would consume much of the molecular gas, and then for the remaining H I to condense into H₂. The stellar populations (see Sect. 6.3) do not support this scenario and tidal stripping is a more natural explanation as it preferentially removes the most loosely bound gas, which is typically H I, not H₂. However, this would imply that the encounter(s) responsible for stripping the H I gas did not trigger a major SF event in these galaxies, despite the presence of a considerable amount of molecular gas.

6.2. Hot gas in the IGrM

Belsole et al. (2003) and O’Sullivan et al. (2014a) measured the hot diffuse gas component of the IGrM of HCG 16 with the XMM-Newton and Chandra satellites respectively. The fact that this hot diffuse medium is detected at all is already unusual for a spiral-rich HCG, but in addition it is also quite massive. Belsole et al. (2003) estimate the total hot gas component of the IGrM to be 4.5 × 10¹⁰ M_⊙ and O’Sullivan et al. (2014a) estimate 5.0 − 9.0 × 10¹⁰ M_⊙ (after adjusting to a distance of 55.2 Mpc). As thoroughly discussed in O’Sullivan et al. (2014a) the origin of such a large amount of hot gas is difficult to fully explain. The group is not virialised, so it is unlikely that the gas has a cosmic origin and has simply fallen into the group halo and been heated. The group is also not globally deficient in H I, meaning that stripped gas cannot be the main source either. O’Sullivan et al. (2014a) conclude that the most probable origin of the majority of the hot gas is the galactic winds of HCG 16c and d.

The hot gas in the vicinity of the group core is of course co-spatial with a large amount of H I in tidal structures, demonstrating that the IGrM is multi-phase. This hot gas is unlikely to act as a source of additional cool gas due to its long cooling timescale (7−10 Gyr, O’Sullivan et al. 2014a), but it could negatively impact the lifetime of the H I, which we discuss in Sect. 6.6.

6.3. Stellar populations, star formation rates, and outflows

O’Sullivan et al. (2014b) used the STARLIGHT code and spectra from SDSS to model stellar populations in HCG 16b and c, the only two of the galaxies with spectra in SDSS. These latter authors concluded that HCG 16b is entirely dominated by an old stellar population with a characteristic age of ∼10 Gyr. They also find some evidence of a very minor SF event occurring at some point in the past few hundred million years. This event was likely triggered by the ongoing interactions with HCG 16a, but it represents a negligible fraction of the total stellar population. For HCG 16c the results were heavily dependent on the choice of stellar population models, but the general finding was that a significant minority of the stellar mass of HCG 16c was formed in a starburst event during the past few hundred million years, with the rest of the population being made up of old stars (5−10 Gyr).

As mentioned in the previous section, Lenkić et al. (2016) estimated the SFRs of all the core galaxies. Although there is no stellar population estimate for HCG 16a, its central region appears similar in colour to HCG 16b, suggesting it is made up of an evolved stellar population. However, it is surrounded by a ring of SF that is bright in GALEX UV bands and greatly elevates the estimated SFR. In the case of HCG 16d the estimated SFR is very high (∼17 M_⊙ yr⁻¹) as is expected for a LIRG.

Using the Wide Field Spectrograph (WiFeS) on the ANU 2.3 m telescope Rich et al. (2010) studied the biconical outflow from HCG 16d. This outflow, driven by the ongoing nuclear starburst, contains ionised and neutral gas, traced by Hα and Na D lines. They also find A-type stellar absorption features throughout the stellar disc, and suggest that starbursting galaxies like HCG 16d might be progenitors for E+A galaxies. This A-type stellar population in the stellar disc indicates that the galaxy has undergone a global star formation event less than a gigayear ago, in addition to the current nuclear starburst (although they may represent different phases of the same sustained event). Despite this evidence of significant recent star formation, HCG 16d still has sizeable reservoirs of both molecular and atomic gas (Table 3). However, it is not clear whether the H I component of this gas is truly associated with the galaxy, or is a chance superposition.

One reason to believe that the H I gas might not be associated with the galaxy is because of its peculiar velocity structure, which is very disorderly and any potential gradient appears to be almost perpendicular to the stellar disc (Fig. 11). Here we draw a comparison with the H I distribution around M 82, one of the best-studied starbursting galaxies. The H I distribution around M 82 appears somewhat similar to that around HCG 16d, with the velocity gradient in H I aligned with the outflow rather than with the stellar disc (e.g. Martini et al. 2018). In the case of M 82, the H I distribution is interpreted as gas that is entrained in the hot wind, although the details of the exact mechanism are uncertain. This lends support to the interpretation that this H I gas observed in the vicinity of HCG 16d really is associated with it and that its anomalous velocity structure is a result of the current wind. However, it is also possible, and even likely, that both proposed scenarios are somewhat true. There is a great deal of extended H I in the IGrM so it is very plausible that we have mistakenly attributed some of this to HCG 16d.

HCG 16c is another LIRG in the group, and like HCG 16d has a high SFR. Vogt et al. (2013) studied the M 82-like wind that also exists in this galaxy, also with the WiFeS instrument. In striking similarity to HCG 16d, the wind also appears to be a nuclear-starburst-driven phenomenon and the galactic disc shows signs of an A-type stellar population throughout. This indicates that the recent past has been very similar for these two galaxies, with each experiencing a global SF event within the last gigayear and both currently undergoing a nuclear starburst that is powering a galactic wind. The simplest explanation for this synchronised evolution is that it is driven by their tidal interaction with each other. However, there is another plausible explanation, that the passage of NGC 848 through the group triggered these events in both galaxies at approximately the same time (we consider this timescale in the following section). Vogt et al. (2013) argue that although NGC 848 could have triggered the event responsible for the A-type population, the timescales are not compatible for the ongoing starbursts; however it is possible that the events were initially global and have since been funnelled to the nuclear regions.

Despite the apparent similarity in their recent SFHs and the presence of winds, the H I properties of HCG 16c and d are quite disparate. While HCG 16d shows no signs of rotation and a possible velocity gradient along the minor axis, HCG 16c has a mostly regular H I velocity field (Fig. 10) except in its outskirts. Vogt et al. (2013) discussed the different nature of the galactic winds in the two galaxies. The wind emanating from HCG 16c is still (mostly) confined to two bubbles above and below the disc within the H I surrounding envelope, indicating that it is young (only a few million years old). On the other hand the HCG 16d wind is biconical and apparently free streaming. These authors also argue that the primary driving mechanisms of the winds differ, with one being shock-excited (HGC 16d) and the other photoionised (HCG 16c). Given that these two galaxies are in the same environment and appear to have had similar recent interactions, the differences in these winds are probably due to pre-existing differences in the host galaxies or simply the different phases we are currently observing them in. We refer the reader to Rich et al. (2010), Vogt et al. (2013), and references therein, for further discussion on this topic.

6.4. Tail age estimates

The SE tail, which links the core group to NGC 848, was most likely formed by NGC 848 passing very close to the core group, unbinding (or attracting already loosely bound) H I gas and stretching it out to form the ∼160 kpc-long tail. Across most of its length the SE tail is visible in just four velocity channels (3942–4006 km s⁻¹), which suggests the feature is approximately aligned with the plane of the sky and that we can consider its projected length as almost equal to its physical length. Given the large separation between NGC 848 and the core group it is also reasonable to assume that when NGC 848 passed through the group it was travelling at approximately the escape velocity. Summing the stellar masses on the four core galaxies gives the total stellar mass of the core group as 4.42 × 10¹⁰ M_⊙ (Table 3). Combining this with the stellar mass–halo mass relation from Matthee et al. (2017, their Eq. (3) and Table 2) calculated from the Evolution and Assembly of GaLaxies and their Environments (EAGLE) simulations, we estimate the dark matter halo mass of the core group as 2.8 × 10¹² M_⊙. This corresponds to an escape velocity of ∼400 km s⁻¹ at the present separation. Assuming this velocity, NGC 848 would have passed by the group approximately 400 Myr ago. It should be noted that this value is relatively uncertain as the simplistic argument above hides many complexities. However, it is still useful as an order-of-magnitude guide.

The only disturbance to NGC 848 that is visible in the optical image is that its spiral arms appear quite loose and are extended to the northwest and the southeast. However, over the majority of the extent of the H I tail there is no detectable optical counterpart. Konstantopoulos et al. (2010) estimate that optical tidal features in groups will be dispersed within about 500 Myr. Given our age estimate of the H I feature, this may explain why no optical counterpart is seen. However, it is also possible that the tail is formed of H I gas that was already loosely bound and did not host any significant stellar population. In this case there may never have been an optical counterpart. We discard the possibility of in situ SF within the SE tail because at no point along its length does the column density rise to 10²¹ cm⁻² (the peak value is 1.3 × 10²⁰ cm⁻²).

Following the arguments in Sect. 3.5 of Borthakur et al. (2015) we can also make an estimate of how long the H I content of the SE tail can persist into the future. The peak column density along the spine of the SE tail is about 10²⁰ cm⁻² and it is about 1′ (16 kpc) in width. If we assume that the tail is a cylinder of gas then the average density is about 0.016 cm⁻³. Borthakur et al. (2015) find that H I becomes susceptible to ionisation from background radiation below densities of about 10⁻³ cm⁻³ at column densities of about 10¹⁹ cm⁻² (their Fig. 6). If we assume that the H I clouds making up the SE tail are expanding with a fiducial velocity of 20 km s⁻¹ then it will reach this threshold column density after about 1 Gyr and the threshold density after about 1.7 Gyr. We therefore expect this feature to survive for at least another gigayear.

Konstantopoulos et al. (2013) use the colours (from SDSS images) of the optical tail extending eastward of HCG 16a to estimate an age of between 100 Myr and 1 Gyr. While in projection this tail appears to be associated with the NW tail (Fig. 2), as discussed previously this H I feature does not connect kinematically to HCG 16a. Instead, if this optical tail has an H I counterpart it is probably the much smaller H I feature visible on the eastern side of HCG 16a at 4027 km s⁻¹ (see the channel maps Figs. D.1 and D.2). Given the loose constraint on the age of this tail and the density of the core group it is difficult to say what interaction is responsible for this tail. It could simply be the ongoing interaction between HCG 16a and b, or an interaction with HCG 16c and d, or perhaps even due to the recent passage of NGC 848.

6.5. The NW tail: accreting gas, tidal tail, or outflow?

The NW tail intersects HCG 16c at its centre on the plane of the sky and at its central velocity in the spectral direction (Fig. 7). As HCG 16c is currently undergoing a SB event it is worthwhile considering if this could be a sign of low-angular-momentum, cool gas accreting onto the centre of HCG 16c and fuelling its starburst. However, before asserting such an exceptional hypothesis we should first examine and attempt to eliminate other more mundane scenarios. As mentioned previously, two other competing hypotheses are that this feature could be the result of an outflow or the chance superposition of a tidal tail.

First, we consider the outflow hypothesis. We have already argued that H I gas in HCG 16d is being disrupted by its galactic wind and that there is also a wind emanating from HCG 16c. However, the NW tail is a well-collimated feature, albeit with a pronounced curve, whereas the gas in the vicinity of HCG 16d is disordered and even the H I velocity field of M 82 (Martini et al. 2018) does not show collimated features like this. Without invoking a mechanism to collimate the ouflowing neutral gas over many tens of kiloparsecs it would not be possible to form such a feature with a galactic wind. We therefore discard the possibility of the NW tail being an outflow.

Next, we discuss the possibility of a chance superposition. Given that the feature intersects the centre of HCG 16c in both velocity and position, a chance superposition seems, at first, unlikely. However, we know already that NGC 848 must have passed extremely close to HCG 16d (and therefore HCG 16c) in order to form the SE tail. Also the SE tail and the NW tail appear as though they may be part of one continuous feature which passes through the core group (Fig. 7, top panel). This feature may trace the path of NGC 848 when it passed through the core group, with the NW tail being the leading tail formed from the gas surrounding HCG 16c and d as NGC 848 approached the group, and the SE tail being the trailing tail which became very extended as NGC 848 exited the opposite side of the group. In this scenario a chance superposition of the centre of HCG 16c and the NW tail is not nearly as contrived as it might otherwise be.

In summary, while it remains possible that the NW tail is accreting onto the centre of HCG 16c, given the other information about the likely recent past of the group it seems that a chance superposition is the most likely explanation, that is, that the NW tail is a tidal tail with the same redshift as the centre of HCG 16c, but not at the same line-of-sight distance (owing to peculiar motions).

6.6. Survival of H I in the IGrM

We observe a considerable amount of H I gas in HCG 16 that is not associated with any individual galaxy. Here we discuss the stability of this gas considering the ongoing processes within the group. As noted by Verdes-Montenegro et al. (2001), HCGs seem to have two final-stage morphologies⁶; those with almost no remaining H I, and those with H I found only in extended features, not associated with individual galaxies. Which of these will be the fate of HCG 16?

Borthakur et al. (2015) discussed the distribution and fate of diffuse H I gas in four compact groups. These authors formulated an analytic approximation for the minimum distance an H I cloud of a given column density can be from a starburst event forming stars at a given rate (their Eq. (3)). The lowest-column-density contours in Fig. 2 are 2.45 and 9.80 × 10¹⁹ cm⁻², and these enclose much of the area around the galaxies in the core group. Using the fiducial column density of 5 × 10¹⁹ cm⁻² and the expressions from Borthakur et al. (2015) we estimate that such H I clouds should not be stable within ∼100 kpc of either HCG 16c or d, yet this would rule out most of the core region of the group, where we already know there is H I gas.

Borthakur et al. (2015) found a similar apparent contradiction in HCG 31, but reasoned that the diffuse H I probably followed a similar distribution to the higher-column-density H I and was thus likely shielded from ionising photons. In HCG 16 it may be that the apparently low-column-density features are really made up of dense clumps, smaller than the resolution of our images (∼30″ or about 8 kpc), which have their emission smeared out by the beam. Furthermore, the energy output from the starbursts in HCG 16c and d will be highly non-isotropic, which may mean that a small faction of the core group is very hostile to H I clouds, but that the majority is not.

Using deep Chandra observations O’Sullivan et al. (2014a) estimated the temperature of the hot, diffuse IGrM in HCG 16 as 0.3 keV (3.5 × 10⁶ K) and its number density as around 1 × 10⁻³ cm⁻³. Following Borthakur et al. (2010) we estimate that the critical radius of H I clouds to prevent evaporation due to conductive heating is about 2 kpc. Given the spatial resolution of our VLA data we cannot verify this directly, other than to say that the persistence of H I in the IGrM of HCG 16 implies that the H I clouds are larger than this limit. We also note that lack of correlation between the H I properties and hot IGrM properties found in other HCGs (Rasmussen et al. 2008) disfavours conductive heating as being a key H I removal mechanism in HCGs.

The final source of energy which could have the potential to remove the H I gas from the IGrM is the starburst-driven galactic winds emanating from HCG 16c and d. Rich et al. (2010) and Vogt et al. (2013) describe these winds that are driving out neutral gas from the discs of these two galaxies. Despite the dramatic nature of these winds, Borthakur et al. (2013) argue that the energy release rate of an M 82-like wind is similar to the cooling rate of the surrounding medium, so while a starburst can maintain a galactic-scale wind while it is active, it is unlikely to have a major long-term effect on the temperature of the IGrM.

To summarise, we are not aware of any source that is likely to destroy a considerable fraction of the H I currently observed in the IGrM on short timescales (< 1 Gyr). Therefore, we expect the dominant effect to be gradual evaporation due to background UV radiation.

6.7. Correspondence of H I and optical extended features

The DECaLS images were masked using the Noisechisel software (Akhlaghi & Ichikawa 2015). We then rebinned the images to a pixel size of 1″ followed by smoothing with a kernel of 2 pixels. This enhances the diffuse emission of the data, allowing the lowest-surface-brightness features of the data, not detectable at the original pixel scale of ≈0.27″, to be identified. By comparing images in the g, r, and z bands we identified four new LSB features (Fig. 16) which appear to be associated with HCG 16 (there are other LSB features which were attributed to background galaxies or clusters). The over-subtraction and fluctuation in the reference background level, especially close to extended galaxies, makes the surface brightness of the features unreliable, but we estimate the features to be between 27 and 28 mag arcsec⁻² in r-band.

Fig. 16. Top: enhanced DECaLS gr image of the core group with LSB features highlighted. The black regions around the galaxies are due to over-subtraction of the sky that occurs in the DECaLS pipeline. The grey regions are the masked high-surface-brightness emission. Bottom: as above but for the area surrounding NGC 848. This image also includes H I contours from the combination of the three channels where the NGC 848S tail is most pronounced. Levels: −1, 1, 2, 3, 4, 5 × 0.05 Jy km s−1 per beam, ×4.9 × 1019 cm−2, or ×0.39 M⊙ pc−2.

In the core group we were unable to confidently identify any H I structures which correspond to these LSB features. As stellar tails are generally thought to be shorter lived than H I tails this likely indicates that these LSB features have a gas-poor origin. For example, it is possible that the most western feature originates from either HCG 16a or b, which are both H I deficient. An alternative source for these features could also be disrupted dwarf members of the group that have been accreted by the larger members, leaving only faint stellar tails.

Near NGC 848 there is a very linear LSB feature, which was almost obscured by a star (Fig. 16, bottom). This feature is co-spatial with the H I emission of the NGC 848S tail and is most likely its optical counterpart. It has a g − r colour of 0.4 ± 0.1. The contours in the lower panel of Fig. 16 are from the combination of three channels (4027−4069 km s⁻¹) where the NGC 848S tail is most prominent. The peak-column-density contour within the tail is 1.5 × 10²⁰ cm⁻², and therefore in situ star formation is unlikely and we are probably detecting the old stellar population in this tail. As optical tidal tails are expected to be observed for only ∼500 Myr after formation, a barely visible stellar component is consistent with our age estimate of 400 Myr for the SE tail and for the interaction of NGC 848 with the rest of the group. If this tidal feature formed at approximately the same time then today it would be expected to be barely visible in the optical.

Even in these enhanced images there is still no indication of a stellar counterpart for the SE tail. The deepest of the three bands is r with a 3σ surface brightness sensitivity of 28.7 mag arcsec⁻² in a 10″ × 10″ box. This null detection supports the hypothesis that the SE tail was formed from loosely bound H I gas that had no associated stellar counterpart.

6.8. Morphological transformation

Hickson compact groups are known to host an excess of early-type galaxies, particularly S0s, and have a corresponding shortage of late-type galaxies relative to the field (Hickson et al. 1988). This shortage of spirals and excess of S0s has led various authors to suggest that as spiral in CGs are stripped of gas they may evolve into S0s (e.g. Verdes-Montenegro et al. 2001; Sulentic et al. 2001). In addition, CGs are found to have a pronounced gap in IR colours (Johnson et al. 2007), referred to as the “canyon” or the infra-red transition zone (IRTZ), where galaxies appear to be in the late stages of shifting their morphology after passing the optical green valley (Alatalo et al. 2014) and beginning to lose their molecular gas (Lisenfeld et al. 2017). Furthermore, in HCGs the IR canyon is dominated by lenticulars and early-types, again suggesting that spirals may be transitioning to later types in CGs.

HCG 16 contains two S0a galaxies, HCG 16c and d, which are both undergoing starbursts, as would be expected in the model proposed by Bekki & Couch (2011), in which spirals transform into S0s through repeated close interactions and the associated SF events. This raises the question of whether or not HCG 16c and d are currently undergoing morphological transformation from spirals to lenticulars.

In the DECaLS image (Fig. 1) faint stellar shells are visible around both HCG 16c and d. Shells in massive galaxies have been widely regarded as the result of interactions with low-mass companions on approximately radial orbits (e.g. Quinn 1984; Dupraz & Combes 1986; Hernquist & Quinn 1988, 1989) and have generally been found around early-type rather than late-type galaxies (e.g. Malin & Carter 1983; Atkinson et al. 2013). However, more recently Pop et al. (2018) performed a study of the occurrence of shells based on the Illustris cosmological simulations, finding that they can be caused by more equal-mass interactions as well. In addition, they find that although for 1:10 mass ratio interactions (the minimum found to produce shells) the orbit must be almost purely radial, the more equal the mass ratio of the progenitors the wider the range of impact parameters that can produce shells. Therefore, the simplest and most likely explanation for both HCG16c and d displaying shells is that one leads to the formation of the shell of the other and vise versa. However, this does not rule out the possibility that they both independently formed shells in mergers with now unseen companions prior to their present interaction with each other.

The presence of shells in HCG 16c and d also suggests that they are likely to have earlier morphologies in the near future, as most shells are found around early-type galaxies. However, the formation mechanism of S0s in groups and the field is still a topic of debate, and evidence exists that supports either a secular (Guérou et al. 2016; Rizzo et al. 2018) or violent (Laurikainen et al. 2010; Querejeta et al. 2015; Tapia et al. 2017; Eliche-Moral et al. 2018) evolutionary pathway from spirals to lenticulars, or perhaps both (Fraser-McKelvie et al. 2018). It has even been suggested that spirals are not the progenitors at all and that lenticulars represent an old population distinct from both ellipitcal and spiral galaxies (Gao et al. 2018).

It should also be noted that HCG 16c and d are not in the IRTZ/canyon because they are actively forming stars, which IRTZ/canyon galaxies are not, almost by definition. These latter are however already classified as S0a, indicating the lenticular-like appearance arose before any potential future transition. The only galaxy in the group that has apparently traversed the canyon already is HCG 16b (Zucker et al. 2016), which has clear spiral features. However, just because crossing the canyon was not accompanied by a morphological change for this galaxy does not mean it cannot be for others, but it does indicate that caution is needed when trying to interpret small samples.

Given this lack of consensus regarding the origin of S0s and what relation they may have to the IRTZ, it is not possible at present to say whether HCG 16c and d are likely undergoing a transformation from spiral to lenticular morphology, or whether they were lenticular already. However, this is a possibility that merits further consideration both in HCG 16 and other CGs.

7. Summary and outlook

In this final section we first review our findings of the present state of the group as evidenced by the H I data and the rich multi-wavelength data set described in the literature. To conclude we then propose a scenario which explains these findings and discuss the probable end state of the group.

7.1. Summary

Overall, HCG 16 is not deficient in H I, it has the expected quantity of H I gas given the B-band luminosity of its members. However, while the total amount may be equivalent to that found in isolated galaxies, the gas is unevenly distributed throughout the group and its members. The northwestern pair, HCG 16a and b, have both lost the vast majority of their expected H I content, while the other members all have normal H I masses. The remaining gas is spread out through the IGrM in a tangle of tidal tails, bridges, and dense clumps, some of which may be TDGs.

Despite being spread out across the entire group, the most plausible origin for the majority of the extended H I gas is the HCG 16a and b pair. This pair is undergoing a strong interaction resulting in multiple optical and H I tidal features between and around them. However, despite this interaction (and the presence of significant molecular gas reservoirs) they do not appear to have highly elevated SFRs (although HCG 16a does have a ring of SF in its outer disc) and their optical colours and stellar population models do not point to bursts of SF in the recent past. Given that these galaxies are both H I deficient and have not recently converted large amounts of gas into stars, we conclude that the H I gas must have been tidally stripped by interactions without triggering starburst events.

The pair at the centre of the group, HCG 16c and d, is also undergoing a strong interaction, but with quite the opposite outcome. Both galaxies are currently starbursting and have galactic winds powered by these events. They are also embedded in a tangled web of extended H I features which form a high-column-density bridge between the two galaxies, tidal tails to the NE and NW, dense knots, and an enormous tail to the SE. The H I kinematics of HCG 16c are mostly regular in its inner region, although the outer disc is somewhat disturbed. The NW tail connects in position and velocity to the centre of HCG 16c, presenting the possibility of cool gas with low angular momentum accreting directly to its centre and fuelling the starburst. However, given the interactions we traced through the group we favour an interpretation of this as a chance superposition in velocity space that does not correspond to the two objects being truly co-spatial in three dimensions. In HCG 16d, the H I kinematics appear to be completely disrupted; there is only a faint indication of a gradient in its velocity field and this is approximately perpendicular to the major axis of the disc. We interpret this as the starburst event disrupting the H I and it becoming entrained by the galactic wind. The rapid SFRs of these two galaxies mean that they will exhaust their gas reservoirs within about 1 Gyr, leaving them gas deficient without external replenishment.

The final large galaxy in the group, NGC 848, is separated from the other galaxies in optical emission by approximately 160 kpc, but this distance is traversed by the enormous SE tail, which forms an H I connection between NGC 848 and the core group. Other than the unwinding edges of the disc of NGC 848, this tidal tail has no apparent optical counterpart along its entire length, indicating that it is either too old for any accompanying stellar component to visibly persist to the present day, or that it was formed from loosely bound gas that did not have an associated stellar component to begin with. A simplistic estimate of the age of the tail, based on the assumption that NGC 848 is travelling at approximately the escape velocity of the group, gives 400 Myr, which is on the upper end of how long a counterpart is expected to survive. Given the approximate nature of this estimate, it is hardly conclusive either way. However, due to the abundance of neutral gas in the core group that is not associated with any galaxy, we favour the latter interpretation.

Therefore, the dominant processes modifying the H I content of HCG 16 are tidal stripping and star formation. Tidal interactions removed the majority of the H I content of HCG 16a and b, and spread it out across the entire group, while SF in HCG 16c and d is rapidly consuming molecular gas (which will presumably be replenished from the available H I reservoirs) and disrupting the H I disc of HCG 16d. If this interpretation is correct then it would contradict the Konstantopoulos et al. (2010) modified evolutionary sequence for HCGs because the HCG 16a and b pair would have been a gas-rich, strong interaction which did not result in gas being consumed by SF, whereas HCG 16c and d are another gas-rich pair which clearly has resulted in elevated SF, demonstrating not only that these two different results are possible, but that they are even possible in the same group. Therefore, the group cannot be classified as either a case where gas is mostly consumed by SF before major interactions occur, or as a case where tidal interactions remove the gas before it can be consumed by SF (the two distinct pathways in that scheme).

After reviewing the potential mechanisms for ionising the H I currently in the IGrM we find that the hot component of the IGrM is not energetic enough to evaporate large (> 2 kpc) H I structures, that the ongoing SF in the group does not appear to be strongly effecting the existing H I in the IGrM, and that while the galactic wind of HCG 16d may currently be ejecting gas and energy into the IGrM, the effects of this are unlikely to persist once the starburst event has ended. Thus, evaporation by the UV background will likely be the principal mechanism for removing H I from the IGrM on a long timescale (> 1 Gyr).

7.2. Global picture and future outlook

With all of the above results in mind we attempted to construct a coherent picture of the past evolution of the group that fits with all the available evidence: Strong tidal interactions involving HCG 16a and b likely unbound much of their H I gas without triggering a major SF event. This unbound gas was then dragged through the centre of the group by the passage of NGC 848 about 0.5 Gyr ago. This close passage also started SF episodes in HCG 16c and d, generating the E+A-like spectra their discs have today. This passage of NGC 848 is traced by the SE tail and the NW tail, which together form a continuous structure spanning the entire group, and the latter of which may have formed a TDG at its tip. At present the SFRs of HCG 16c and d are highly elevated, now driven by their interaction with each other, while the SFRs in HCG 16a and b remain more modest, aside from a ring of SF activity in HCG 16a.

Over the following gigayear HCG 16c and d will likely convert, consume, or expel much of their gas supply through SF leaving themselves H I deficient like HCG 16a and b, though by different means. Meanwhile, the H I in all the galaxies will continue to be stripped by tidal interactions. It is unclear whether NGC 848 is travelling fast enough to escape the group, or whether it will fall back to be the sole H I-normal large galaxy in the group. The extended H I features in the group are expected to persist for several gigayears as they are gradually evaporated by the UV background. This will result in HCG 16 resembling a phase 3a group where there is little or no H I remaining in the galaxies (with the possible exception of NGC 848), but extended H I features are still visible in the IGrM.

Acknowledgments

MGJ is supported by a Juan de la Cierva formación fellowship. We also acknowledge support from the grants AYA2015-65973-C3-1-R and RTI2018-096228-B-C31 (MINECO/FEDER, UE). This work has been supported by the Spanish Science Ministry “Centro de Excelencia Severo Ochoa” program under grant SEV-2017-0709. MGJ wishes to thank B. Koribalski, K. Lee-Waddell, and S. Cazzoli for helpful discussions. We also thank the referee for his thorough comments which helped to improve this paper. This project used archival data from the VLA. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. This project used data obtained with the Dark Energy Camera (DECam), which was constructed by the Dark Energy Survey (DES) collaboration (full acknowledgement at legacysurvey.org/acknowledgment). 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. We also acknowledge the use of the HyperLeda database (Makarov et al. 2014). This research made use of APLpy, an open-source plotting package for Python (Robitaille & Bressert 2012; Robitaille 2019), astropy (Astropy Collaboration 2013, 2018), Aladin (Bonnarel et al. 2000), mayavi (Ramachandran & Varoquaux 2011), and SAOImageDS9 (Joye & Mandel 2003).

References

Appendix A: Full workflow and reproducibility

Aside from the scientific aims of this work we also aimed to follow the scientific method as fully as possible by making the reduction and analysis completely reproducible. Our intention is to provide other astronomers with access to our data and methodology such that they are able to understand what was done to the data (at any desired level), to verify our analysis, or to reuse any of our outputs. We were not able to fully meet these goals for every aspect of our analysis, but have done so for the majority. The exhaustive findings of this study will be detailed in another publication, but in this section we provide a brief description of our methodology and experiences with regard to reproducibility. Our full workflow is available online⁷, while the final plots can be recreated and modified⁸ using the default version of the final data product (as in this paper) without needing to execute the whole workflow.

A.1. Need and motivation

Reproducibility is a key pillar of the scientific method, but in the current age of digital research it is frequently not possible to reproduce or reuse previous results. Although we could point to a number of studies indicating that science (of all sub-disciplines) is currently in a reproducibility crisis, indicating the necessity of methodological improvements, it is perhaps more informative to highlight this point with some concrete examples from our own project.

– Figure 4 shows that the spectral profile and total flux measurements do not quite agree with those from HIPASS. Although we have identified the baseline subtraction around this source as a potential problem, without a substantial investment of time and an expert from within the HIPASS team it would not be possible to fully understand this difference. While the HIPASS reduced data are publicly available and the general methodology has been published, it is difficult to know exactly what has happened to a specific part of the data. This meant that our investigation of this issue was only possible from one side, which was ultimately not completely successful.

– The VLA data used in this paper have been used in previous publications (Verdes-Montenegro et al. 2001; Borthakur et al. 2010). We had access to one of the previous data cubes (reduced using AIPS) that one of the authors still had on their hard drive. Initially we were making comparisons between that cube and the one from this work to verify our results and investigate the effect of using multi-scale CLEAN. However, these comparisons were eventually abandoned because even with the file history in the fits file, it was difficult to know exactly what had been done previously. Furthermore, if we disagreed with a step from the previous analysis it would have been a considerable amount of additional work (if even possible) to reproduce the previous reduction and modify that step.

Our goal is to take a step forward in this work to show how such barriers to reproducibility can be minimised in future. In recent years this difficulty with reproducibility in modern science has received considerable attention (including by members of the research team at the Instituto de Astrofísica de Andalucía). One such recent work, Wilkinson et al. (2016), describes a set of generic principles regarding how scientific data should be stored and documented in order to make them as reusable as possible, in particular, for automated processing. These are the FAIR principles (Findable, Accessible, Interoperable, and Reusable)⁹. Ourselves, Wilkinson et al. (2016), and others argue that these principles should not only be applied to data, but also to algorithms and tools in order to enable full transparency and reproducibility in practice. We have used these principles to help guide our approach, which we discuss below.

A.2. Discussion of our approach

Before embarking on creating a workflow framework it is necessary to decide at what level a project should be reproducible. This ranges from providing just the raw and reduced data, to allowing every value in the final paper to be re-generated in a completely automated way, such that if any error were to be identified in the reduction or analysis (or that of previous works on which we depend) then all dependent values throughout the paper would be automatically changed once it was corrected. Unfortunately the latter is not yet fully achievable (discussed further below), but we also wish to go beyond providing only the data by allowing other astronomers to fully understand our methodology and to modify or reuse it. Therefore, we have elected to make every digital object created in this work a reproducible product of the workflow. In this case this corresponds to every figure and table which displays results of this work. This includes Table 2 as the properties it contains are all non-trivial results of this work, however, Tables 1 and 3 are not included as they are either values taken from the literature or the result of trivial calculations with literature values. Unfortunately, these aims have not been fully realised and there are a few results which are not included in our workflow, which we discuss below.

To fully describe a workflow it is not sufficient to provide only the data (raw or reduced) and the reduction scripts as different versions of the same software might not produce identical results, even given the same input, and installing unfamiliar software can often be a barrier. Therefore, it is also necessary to regularise the software environment. To do this we made extensive use of docker¹⁰ containers and Conda¹¹ environments. Each piece of software used, in this case CASA and SoFiA, was placed inside a docker container so that exactly the same software environment can be used on any machine. A Conda environment was created to regularise the Python environment which was used for generating the final plots. All the steps of the reduction, analysis, and plotting were chained together in a single workflow using the CGAT-core¹² workflow management software. This allows the final workflow to be executed using a single script which downloads the raw data and scripts, builds the Conda environments, pulls the docker containers, runs the data reduction and analysis steps, and finally generates the plots displayed in the paper. Furthermore, any one of these steps can be executed individually (provided the prerequisite steps have already been executed) and the reduction and analysis scripts themselves are stored separately from the software environment so that they can be modified according to the wishes of any end user and the workflow can be re-run. For example, if a particular user wanted the data to be CLEANed to a lower threshold than was used, this would require only the threshold value in the imaging script to be changed and the workflow executed again¹³.

As mentioned above, our original goal has not been fully achieved and some figures are missing from the workflow. Figures 7 and 16 are not included and the velocity labels in the first moments of the galaxies (Figs. 8–11, 13, and 14) had to be added manually, so will not be updated if the iso-velocity contours change. The X3D interactive plot was also not included in the workflow as at present this still requires manual alterations to circumvent bugs in the software.

Figure 7 was not included as this was the only figure for which the tools in SlicerAstro were required (unless we were to completely reinvent the “wheel”). While SlicerAstro was a key tool for separating emission from features and galaxies, this process ultimately results in a set of masks which define the regions assumed to contain each source. As these masks are created in a manual way, the result can only be preserved rather than the process being made completely reproducible. Hence, the final masks were imported into CASA and used to generate the moments, which meant that creating a docker container for SlicerAstro was not necessary, with the exception of this figure which uses its segmented PV slice feature. Although it is possible to script processes in SlicerAstro, the parent Slicer software is not designed in a way that makes this straightforward for general users; for example it is not possible to record steps performed in the GUI as a script. Therefore we decided to omit this plot from the final workflow due to the time investment it would take to include it. Figure 16 was not included because it relies on an IDL pipeline to mask the bright emission and enhance the LSB features in the DECaLS images. This pipeline is in the process of being converted to Python, but as this is a work in progress it has not been included.

One of the key reasons why we did not attempt to make every value in the paper (including those that only appear in the text) fully reproducible is because at present this is not feasible. It is normal for astronomers to upload large data tables or large survey datasets to services such as Centre de Données astronomiques de Strasbourg (CDS) where they can be queried and accessed in an automated way, but this frequently does not cover all the required information. To take an explicit example, consider our Table 3. This table lists the SFRs and gas masses of the four core galaxies in HCG 16 along with the trivial result of the calculation of gas consumption timescale for each. The SFR data were taken from Lenkić et al. (2016) and the molecular mass data from Martinez-Badenes et al. (2012), while the H I mass values are from this work. In order to make this table completely reproducible (e.g. if a hypothetical mistake in one of those papers was to be found and corrected) all these data would need to be stored in a service like CDS and in addition we would also need an automated way to read the distance used (or the Hubble constant if distance is inferred directly from velocity). Our intention here is not to criticise these works, indeed their data tables are (mostly) available at the CDS, and we are also guilty of shortcomings in this respect. For example, someone wishing to automatically access the distance which we used for HCG 16 would face the same problem, while the B-band luminosities we list in Table 1 come from a private database compiled by the AMIGA team. Another issue is when the information provided by another paper is an equation; for example we use a scaling relation from Jones et al. (2018), but equations can also have mistakes and typographical errors, and there is no automated method to connect them between different works.

To overcome these shortcomings would require a standardised electronic method of reporting all this information so that all external information used in a paper can be fully clarified within a workflow. For data, services like CDS and the International Virtual Observatory Alliance are suitable solutions, but all data tables need to be included, even those which seem insignificant. However, for values such as the Hubble constant that a paper uses, the distances to objects (for papers that use a small sample of objects), or any other values which appear only in the text, there is currently no good solution of which we are aware. The same is true for equations. Ideally, such information would be included in some form of metadata that accompany a paper or its data.

At present making projects reproducible requires an enormous additional effort and it is unsurprising that it is so often neglected. However, once tools and standardised approaches have been established, the burden placed on individual scientists will be greatly decreased. There are many other tools and services available to help this process and we plan to continue to explore and test them in future works.

A.3. Workflow overview and final products

The workflow¹⁴ begins by constructing the necessary software environments as discussed above and then runs the following steps:

1. Download data: The raw VLA data, the CLEANing mask, the masks that we made using SlicerAstro to separate each galaxy and feature, and the GBT and HIPASS cube are all downloaded from the EUDAT B2SHARE service¹⁵ which provides a cloud repository for data¹⁶.

The reason the VLA data were re-hosted here is both because it was necessary to include the masks in addition to the raw data and because of the nature of the VLA archive it is not currently possible to access its data in an entirely automated manner.

2. Flagging and calibration: The exact steps (performed in CASA) to flag and calibrate the raw data are run.

3. Imaging: The imaging steps are run to generate the H I cubes with CASA. A simple threshold moment-zero map is also created and the masks for each separate feature and galaxy are used to generate mini-cubes and moments.

4. Masking: SoFiA is used to make two masks of the entire group (one with a lower threshold intended to included all real emission and one with a higher threshold to make the velocity field). The first mask is later used to generate the integrated spectrum. This step also produces a moment-zero and first-moment map. SoFiA is also used to make a mask of HCG 16 in the HIPASS cube.

The above steps generate all the final cubes and maps needed to produce the figures in the paper. For each figure in the paper we have provided a Python script in the form of a Jupyter¹⁷ notebook. These notebooks begin by listing the input data used in making each plot and explain which steps of the workflow those data dependent on and direct the user to the relevant scripts, should they wish to modify them and re-run the workflow. As we do not anticipate the majority of readers will want to modify our methodology, but may be interested in modifying or reusing our plots, these notebook and the final data products have also been provided (separate from the workflow) and may be accessed through the notebook server¹⁸.

At present the workflow can be downloaded and run by any user locally (or the notebook can be accessed in the cloud). We have also tested running the workflow in the cloud using a virtual machine hosted by the EGI Federated Cloud, which we accessed through the European Open Science Cloud (EOSC) catalogue, and containers on the SciServer service¹⁹. While there are definite advantages to working in these environments, for example on SciServer the software environment is already fixed for the end user which immediately simplifies the process of reproducibility, thus far our workflow has not been fully integrated with either of these services, and so is only available for local execution. In the near future, with the data volumes that upcoming surveys will produce, there will cease to be a choice between working locally and in the cloud, and we will continue to work on cloud-based implementations of reproducible workflows.

Appendix B: Integrated flux and profile discrepancies

Figure 4 shows the comparison of the VLA spectrum measured in this work and the spectrum we extracted from the relevant HIPASS cube. In addition to the slight excess of flux in the VLA spectrum (which is formally not possible) the profile shapes of the H I emission are quite different, suggesting more investigation is needed.

While regions that are brighter in HIPASS could be explained by the presence of diffuse gas that is resolved out by the VLA, the reverse does not apply. The HIPASS profile ends at ∼4100 km s⁻¹, whereas the VLA spectrum continues until ∼4300 km s⁻¹. Investigating this further we found faint negative regions in the HIPASS data, possibly indicating problems with baseline subtraction around this extended source. We also found that the disagreement between the two datasets is a function of position, with them agreeing well on a diagonal line stretching SE to NW, but diverging (in opposite senses in opposite directions) away from this line. As this is not a radial dependence, it is unlikely to result from the VLA primary beam correction. We also found a similar trend when comparing HIPASS with the previous reduction of the same VLA data (Borthakur et al. 2010), indicating the problem is also not due to the new reduction in CASA as opposed to AIPS, or MS-CLEAN versus standard CLEAN. Unfortunately we were unable to identify the exact cause of this disagreement.

Here we also note that the flux scale in our cube is slightly higher than the previous reduction published in Borthakur et al. (2010). These latter authors only show the VLA spectrum weighted to correspond to the GBT observation, finding that there is 5% less flux in the VLA spectrum than in the GBT spectrum, whereas we find 12% more. Part of this difference can be explained by the masking of the VLA data. The SoFiA-generated mask is designed to smooth over several scales in order to recover as much flux as possible and is thus more extended than more traditional masks. Applying our mask to the cube of Borthakur et al. (2010) results in a profile that falls between our data and the GBT data. Therefore, the mask is only a partial explanation for the difference in flux between our reduction of the VLA data and that of Borthakur et al. (2010). That work quotes the flux of the flux calibrator (3C48) as 16.5 Jy, while we use the CASA-based model (Perley & Butler 2013) which has a flux between 16.43 and 16.49 Jy across the bandwidth. This is a negligible difference that, if anything, would reduce rather than increase the flux we measure.

It might be suggested that using MS-CLEAN instead of standard CLEAN could have resulted in us recovering more flux. However, this is a misinterpretation of what MS-CLEAN does (discussed in Appendix C). In reality, the integrated flux in our MS-CLEAN image is lower than an equivalent reduction using standard CLEAN.

It is also worth noting that the flagging performed in this work and in Borthakur et al. (2010) is markedly different, but from the final image cube that we have available it is impossible to say exactly how. The header implies that almost none of the C array data were flagged previously and that only about 7% of the D array data were flagged. We flagged a total of 9% of the C array data and 37% of the D array data, either due to RFI, shadowing, or otherwise inconsistent baselines. This minimalist versus zealous flagging approach may contribute to the differences we find.

In conclusion we have not been able to fully resolve the small discrepancies between these datasets and we argue that without a complete workflow, as we have sought to create for this project, or completely repeating previous data reductions (which is both beyond the scope of this work and not guaranteed to be successful) it is not possible to reliably identify the exact cause.

Appendix C: Flux recovery with MS-CLEAN

The MS-CLEAN algorithm is sometimes regarded as being able to recover more flux because the CLEAN model (for MS-CLEAN) for an extended source will typically have more flux in its components than a standard CLEAN model based on the same data. In fact we find that our MS-CLEAN cube contains less flux than when we imaged using standard CLEAN, even though the MS-CLEAN model components contain more flux. This is because it is more able to CLEAN diffuse emission and can therefore remove the positive platform (and surrounding negative bowl) on which the entire group sits in a standard CLEAN image. Although in the case of HCG 16 this feature is quite subtle, it nevertheless appears in almost every channel where these is emission from the core group, and thus contributes significantly to the total flux of a standard CLEAN image. Therefore, it is more accurate to say that the MS-CLEAN algorithm gives a more robust measure of the total H I flux in the group and a more realistic representation of the morphology of the H I, but it is misleading to say that it recovers more flux.

Appendix D: Channel maps

The channel maps show the contours of H I emission in individual spectral channels in the region of the group core and around NGC 848 and the SE tail.

Fig. D.1. Channel maps showing the H I emission in the core group overlaid on the DECaLS r-band image. Contour levels: −1, 1, 2, 4, 8, 16 × 0.9 mJy per beam, ×1.8 × 1019 cm−2, or ×0.15 M⊙ pc−2. The contours are coloured (in order of increasing flux): black (dashed), black, blue, purple, red, and orange.

Fig. D.1. continued.

Fig. D.1. continued.

Fig. D.1. continued.

Fig. D.1. continued.

Fig. D.1. continued.

Fig. D.2. Channel maps showing the H I emission in the SE tail and NGC 848 overlaid on the DECaLS r-band image. Contour levels: −1, 1, 2, 4, 8, 16 × 0.9 mJy per beam, ×1.8 × 1019 cm−2, or ×0.15 M⊙ pc−2. The contours are coloured (in order of increasing flux): black (dashed), black, blue, purple, red, and orange. The noise level rises towards the South East as this it is near the edge of the VLA primary beam.

Fig. D.2. continued.

Fig. D.2. continued.

Fig. D.2. continued.

Fig. D.2. continued.

Fig. D.2. continued.

All Tables

All Figures

	Fig. 1. DECaLS grz colour image of HCG 16 with the member galaxies labelled.
In the text

Fig. 2. Left: moment-zero map (primary beam corrected) of the H I emission of HCG 16 calculated using a 3σrms mask in each channel (this mask is intended only for visual purposes, the whole analysis uses the SoFiA-generated mask shown in the right panel). The black ellipse in the lower right corner indicates the beam size (37.2″ × 30.3″ or 10.0 kpc × 8.1 kpc) and the small red cross shows the centre of the GBT pointing from Borthakur et al. (2010). The extended features that we separated from the galaxies are indicated by arrows, with the exception of the emission joining galaxies HCG 16c and d. The filled black circle and the white plus sign indicate the pointing centres of the VLA D and C array data respectively. Right: moment -zero contours (uncorrected for primary beam) overlaid on a DECaLS r-band image. In this case the map was generated using the SoFiA-generated mask described in Sect. 3.1, which includes more extended emission. The galaxies in the main band of H I emission going from top-right (NW) to bottom-left (SE) are: HCG 16b, a, c, d, NGC 848, and the single galaxy to the lower right of the group is PGC 8210. Contour levels: −2.45, 2.45, 9.80, 24.4, 49.0, 73.5, and 98.0 × 1019 cm−2, where 2.2 × 1019 cm−2 corresponds to the 3σ sensitivity in one channel. In order of increasing flux the contours are coloured: black (dashed), black, blue, purple, red, orange, and yellow.

In the text

Fig. 3. First-moment map of the H I emission calculated using the SoFiA-generated mask with a threshold of 5σrms and without dilation. The rainbow colour scale indicates the recessional (optical) velocity in km s−1. The black ellipse in the lower right indicates the beam size, the grey lines are isovelocity contours separated by 40 km s−1, and the small black star symbols indicate the locations of the optical centres of the galaxies in the group.

In the text

	Fig. 4. Integrated VLA spectrum of the entire group (thick black line) compared to the GBT spectrum of the core group (Borthakur et al. 2010, high spectral resolution green line), the VLA spectrum within the GBT beam area (thick black dashed line), and the spectrum extracted from the HIPASS cube (orange line). Vertical dashed lines indicate the optical redshifts of the member galaxies.
In the text

Fig. 5. Spectral profiles of each of the six galaxies detected in H I near the core group calculated from the separation of tidal features and galaxies performed in Sect. 3.1. The black vertical dashed lines show the optical redshifts of each galaxy as given in Table 1. The dot-dashed and dotted vertical lines show the central and extreme velocities, respectively, from the rotation curve measurements of Rubin et al. (1991) in red and Mendes de Oliveira et al. (1998) in blue. These measurements are only available for the four core galaxies; also Rubin et al. (1991) do not specify Vmax values for HCG 16b or c due to the peculiar shape of their rotation curves. The left panels have vertical scales that go to 11 mJy, whereas the scales on the right are a factor of approximately six higher (60 mJy).

In the text

	Fig. 6. H I spectral profiles of each of the eight separate tidal features in the group calculated from the separation of tidal features and galaxies performed in Sect. 3.1. Each row has a different vertical scale.
In the text

Fig. 7. Top: segmented position–velocity diagram showing the spatially and kinematically continuous H I emission spanning HCG 16. The vertical dashed lines show the locations of the nodes making up the segmented slice through the data cube. Galaxies are labelled in red and tidal features in blue. We note that the noise in this plot increases significantly near the left edge as this region is approaching the primary beam edge. Bottom, left panel: position–velocity slice as a segmented red line on top of the moment-zero map (using the standard 3σ threshold mask as in the left panel of Fig. 2) and right panel: same line but overlaid on the DECaLS r-band image.

In the text

	Fig. 8. Left: moment-zero map of HCG 16a. The white and black stars indicate the optical centres of HCG 16a and b respectively, and the black ellipse shows the beam. Right: first-moment contours for HCG 16a overlaid on the DECaLS grz image of HCG 16a and b. The green lines are separated by 20 km s−1.
In the text

Fig. 9. Left: moment-zero map of HCG 16b. The black stars indicate the optical centres of HCG 16a and b, and the black ellipse shows the beam. Right: first-moment contours for HCG 16b overlaid on the DECaLS grz image of HCG 16a and b. The green lines are separated by 20 km s−1. The signature of a small, consistent velocity gradient across the optical disc is visible, but the H I distribution is strongly off centre and disturbed.

In the text

	Fig. 10. Left: moment-zero map of HCG 16c. The white and black stars indicate the optical centres of HCG 16c and d respectively, and the black ellipse shows the beam. Right: first-moment contours for HCG 16c overlaid on the DECaLS grz image of HCG 16c and d. The green lines are separated by 20 km s−1.
In the text

Fig. 11. Left: moment-zero map of HCG 16d. The white and black stars indicate the optical centres of HCG 16d and c, respectively, and the black ellipse shows the beam. A central depression is clearly visible around the white star marking the centre of HCG 16d. Right: first-moment contours for HCG 16d overlaid on the DECaLS grz image of HCG 16d and c. The central region of the galaxy is removed as the H I absorption feature would contaminate the first-moment map there. The green lines are separated by 20 km s−1. While the velocity field is highly irregular, there is a slight gradient roughly aligned with the minor axis of the optical disc.

In the text

	Fig. 12. H I spectral profile of the centre of HCG 16d. The dashed vertical line shows the systemic velocity as determined from optical lines. The profile was extracted from the robust = 0 cube (higher resolution) by placing a double Gaussian window at the position of the optical centre of the galaxy with the shape and orientation of the synthesized beam.
In the text

	Fig. 13. Left: moment-zero map of NGC 848. The white star indicates the optical centre, and the black ellipse shows the beam. Right: first-moment contours for NGC 848 overlaid on the DECaLS grz image. The green lines are separated by 20 km s−1.
In the text

	Fig. 14. Left: moment-zero map of PGC 8210. The white star indicates the optical centre, and the black ellipse shows the beam. Right: first-moment contours for PGC 8210 overlaid on the DECaLS grz image. The green lines are separated by 20 km s−1.
In the text

	Fig. 15. Top: moment-zero maps of the NW, E, and S clumps (left to right) overlaid on the DECaLS grz image. Contour levels: 0.98, 1.47, 1.96, 2.45 × 1020 cm−2. Bottom: first-moment maps of the same features with the beam shown as a black ellipse in the corner.
In the text

Fig. 16. Top: enhanced DECaLS gr image of the core group with LSB features highlighted. The black regions around the galaxies are due to over-subtraction of the sky that occurs in the DECaLS pipeline. The grey regions are the masked high-surface-brightness emission. Bottom: as above but for the area surrounding NGC 848. This image also includes H I contours from the combination of the three channels where the NGC 848S tail is most pronounced. Levels: −1, 1, 2, 3, 4, 5 × 0.05 Jy km s−1 per beam, ×4.9 × 1019 cm−2, or ×0.39 M⊙ pc−2.

In the text

	Fig. D.1. Channel maps showing the H I emission in the core group overlaid on the DECaLS r-band image. Contour levels: −1, 1, 2, 4, 8, 16 × 0.9 mJy per beam, ×1.8 × 1019 cm−2, or ×0.15 M⊙ pc−2. The contours are coloured (in order of increasing flux): black (dashed), black, blue, purple, red, and orange.
In the text

	Fig. D.1. continued.
In the text

	Fig. D.1. continued.
In the text

	Fig. D.1. continued.
In the text

	Fig. D.1. continued.
In the text

	Fig. D.1. continued.
In the text

Fig. D.2. Channel maps showing the H I emission in the SE tail and NGC 848 overlaid on the DECaLS r-band image. Contour levels: −1, 1, 2, 4, 8, 16 × 0.9 mJy per beam, ×1.8 × 1019 cm−2, or ×0.15 M⊙ pc−2. The contours are coloured (in order of increasing flux): black (dashed), black, blue, purple, red, and orange. The noise level rises towards the South East as this it is near the edge of the VLA primary beam.

In the text

	Fig. D.2. continued.
In the text

	Fig. D.2. continued.
In the text

	Fig. D.2. continued.
In the text

	Fig. D.2. continued.
In the text

	Fig. D.2. continued.
In the text