Compact groups of dwarf galaxies in TNG50: late hierarchical assembly and delayed stellar build-up in the low-mass regime Open Access

ABSTRACT

Compact groups of dwarf galaxies (CGDs) have been observed at low redshifts (z < 0.1) and are direct evidence of hierarchical assembly at low masses. To understand the formation of CGDs and the galaxy assembly in the low-mass regime, we search for analogues of compact (radius ≤100 kpc) groups of dwarfs (7 ≤ log [M_*/M_⊙] ≤ 9.5) in the IllustrisTNG highest resolution simulation. Our analysis shows that TNG50-1 can successfully produce CGDs at z = 0 with realistic total and stellar masses. We also find that the CGD number density decreases towards the present, especially at z ≲ 0.26, reaching |$n \approx 10^{-3.5} \ \rm cMpc^{-3}$| at z = 0. This prediction can be tested observationally with upcoming surveys targeting the faint end of the galaxy population and is essential to constrain galaxy evolution models in the dwarf regime. The majority of simulated groups at z ∼ 0 formed recently (⁠|$\lesssim 1.5 \ \rm Gyr$|⁠), and CGDs identified at z ≤ 0.5 commonly take more than 1 Gyr to merge completely, giving origin to low- to intermediate-mass (8 ≤ log [M_*/M_⊙] ≤ 10) normally star-forming galaxies at z = 0. We find that haloes hosting CGDs at z = 0 formed later when compared to haloes of similar mass, having lower stellar masses and higher total gas fractions. The simulations suggest that CGDs observed at z ∼ 0 arise from a late hierarchical assembly in the last ∼3 Gyr, producing rapid growth in total mass relative to stellar mass and creating dwarf groups with median halo masses of |$\sim 10^{11.3} \ \rm M_\odot$| and B-band mass-to-light ratios mostly in the range 10 ≲ M/L ≲ 100, in agreement with previous theoretical and observational studies.

1 INTRODUCTION

Galaxy formation models in a Lambda cold dark matter (ΛCDM) Universe predict hierarchical galaxy formation (White & Rees 1978; Blumenthal et al. 1984; White & Frenk 1991; De Lucia & Blaizot 2007; Frenk & White 2012; Primack 2012), where processes like mergers and interactions play a crucial role in the mass assembly of galaxies. Although the relative importance of these processes to the build-up of the galaxy mass varies across different mass regimes (Cattaneo et al. 2011), the hierarchical assembly is thought to occur at all scales, from dwarf to massive galaxies. Dwarf galaxies (⁠|$M_\ast \lesssim 10^9\,$| M_⊙) stand out as the most numerous systems across all cosmic time (Fontana et al. 2006; Grazian et al. 2015) and are often considered the building blocks of more massive galaxies. None the less, most studies primarily focus on the assembly of high-mass galaxies (Pérez-González et al. 2008; Ilbert et al. 2010; Rodriguez-Gomez et al. 2016), while the role of interactions and mergers in the evolution and mass assembly of low-mass galaxies remains unclear. Furthermore, comprehending the interplay between baryons and dark matter in this particular low-mass regime holds critical significance in evaluating the standard cosmological model (Sales, Wetzel & Fattahi 2022). Given the inherent challenge of directly measuring dark matter, this knowledge becomes indispensable for thoroughly validating the ΛCDM model through the observation of galaxies.

With large spectroscopic surveys, the number of dwarf galaxy systems discovered and studied in detail has been increasing (Stierwalt et al. 2015, 2017; Pearson et al. 2016, 2018; Besla et al. 2018; Paudel et al. 2018; Kado-Fong et al. 2020; Luber et al. 2022; Gao et al. 2023). Most recent studies focus on pairs of dwarfs since they are more frequent than groups of dwarf galaxies. Close pairs, which have projected separations |$R_{\rm sep} \le 50\,$|kpc, offer valuable insights into the impact of interactions on galaxy evolution in the low-mass regime. Stierwalt et al. (2015) identified 104 pairs of dwarf galaxies in the Sloan Digital Sky Survey (SDSS) at 0.005 < z < 0.07 and with pair member masses 7 < log (M_*/M_⊙) < 9.7. Their sample covers a wide range of environments, including isolated dwarf pairs that lie at distances greater than 1.5 Mpc from their closest massive neighbour (⁠|$M_\ast \gt 5 \times 10^9\, {\rm M}_\odot$|⁠). They found that the star formation rate (SFR) in close pair dwarfs is enhanced by a factor of 2.3 compared to a control sample of single dwarfs. The enhancement decreases with increasing R_sep, suggesting that interactions trigger the starbursts. They also find that the SFR enhancement is independent of the isolation of the pair relative to massive neighbours, and that is |$\sim 30~{{\ \rm per\ cent}}$| higher than that observed for massive galaxies in close pairs. These results suggest that interactions may play a more important role in the stellar mass assembly of dwarfs compared to that of massive systems. Detailed analysis of a dwarf galaxy pair by Privon et al. (2017) also reveals evidence of SFR enhancements. They found that the specific SFR (sSFR) of the system is more than an order of magnitude higher than that of non-interacting dwarfs in the same mass range. According to their work, in low-mass galaxy mergers, starbursts may be triggered by large-scale interstellar medium (ISM) compression rather than driven by the funnelling of gas to the nucleus and giving rise to a nuclear starburst. Therefore, the SFR enhancements are more distributed throughout the galaxies. These results and the findings from other recent studies (Kimbro et al. 2021; Martin et al. 2021) indicate that the interaction-driven star formation shows important hydrodynamical differences from what is observed in massive galaxy interactions.

Theoretical work by Martin et al. (2021) using the NewHorizon simulation revealed that enhancement of SFR in dwarf galaxies with |$M_\ast \lesssim 10^8\,$|M_⊙ due to interactions and mergers can be more than three times compared to single dwarfs at z = 1, while the increase is less pronounced in more massive galaxies (⁠|$M\ast \gtrsim 10^{9.5} \rm \ M_\odot$|⁠). These findings reinforce the differing roles of interactions and mergers in low- and high-mass galaxy evolution. Additionally, they highlight that, although mergers are rare among dwarf galaxies (1 or 2 with a ratio greater than 1:10 between z = 1 and 5), interactions between galaxies are more important than mergers in the mass assembly of dwarf galaxies. Together, these processes are responsible for |$\sim 20~{{\ \rm per\ cent}}$| of the stellar mass formed in dwarfs at low redshifts (z ≲ 0.5).

While interactions and mergers might seem to accelerate dwarf galaxy mass build-up and enhance their star formation (Subramanian, Mondal & Kalari 2024 ), they can also enhance gas turbulence, delaying star formation (Ellison, Catinella & Cortese 2018). Dwarf–dwarf tidal interactions can also remove the gas from the dwarf shallow potential and place it at large distances. This gas needs to be re-accreted for star formation, possibly leading to a delay in the stellar mass build-up of these systems (Pearson et al. 2016, 2018). Although gas in the outskirts can provide a long-lived supply to the system in isolated environments, this pre-processing of baryons in dwarf-only systems can affect the efficiency of gas stripping when dwarfs are accreted by a massive host, preventing gas re-accretion. A cosmological zoom-in N-body simulation predicts that this pre-processing is not negligible, with 30 per cent of galaxies with |$10^7-10^9\,$| M_⊙ undergoing pre-processing before entering a Milky Way-like halo (Wetzel, Deason & Garrison-Kimmel 2015). Moreover, recent work by Kado-Fong et al. (2023) suggests that dwarf–dwarf interactions can induce quenching in the low-mass galaxies, although only for a relatively short period of time (≲600 Myr).

In the local Universe, ensembles of dwarfs with three or more members are rare (Besla et al. 2018) and mostly observed as satellites of more massive galaxies, where they are highly susceptible to environmental effects. One of the first studies on these type of dwarf-only systems in the local Universe was carried out by Tully (1987, 1988), who identified structures with low luminosity density (⁠|$2.5 \times 10^8 \lt \rho _L \lt 2.5 \times 10^9 \ \rm L_\odot \,$| Mpc⁻³) within the local volume (⁠|$\lesssim 5\,$| Mpc) containing only faint galaxies. In later work, Tully et al. (2006) confirmed five previously identified associations and argued that these associations are gravitationally bound, which would imply high mass-luminosity ratios (300 ≲ M/L ≲ 1200). Additionally, these associations were found to have large diameters (⁠|$\sim 600\,$| kpc), resulting in distances between member galaxies that are considerably large when compared to the typical sizes of dwarf galaxies. Thus, the associations—characterized by large sizes—of dwarfs are not the ideal sites to study the effects of interactions and mergers on the mass build-up in the low-mass regime. For this purpose, the dwarf galaxies must be closer to each other, inhabiting smaller systems like compact groups of dwarfs.

Recently, Stierwalt et al. (2017) spectroscopically confirmed the existence of seven groups of dwarfs with projected group radius |$R_{\rm group} \lesssim 100\,$|kpc that lie at distances greater than |$1.5\,$|Mpc from their closest massive neighbour. Even though each one of these groups has between three to five members, these groups are more compact and brighter than the associations presented by Tully et al. (2006) and contain more galaxies than systems analysed by Stierwalt et al. (2015), Privon et al. (2017), and others. Stierwalt et al. (2017; hereafter S17) suggest that, given time, merging will turn some of these groups into intermediate-mass galaxies. Therefore, these compact (⁠|$R_{\rm group} \le 100\,$| kpc) groups of dwarfs (7 ≤ log [M_*/M_⊙] ≤ 9.5) are excellent sites to assess the aforementioned impacts of mergers and interactions in the baryon cycle and mass growth of dwarf galaxies. Furthermore, isolated compact groups of dwarf galaxies (CGDs) are excellent cosmic laboratories for studying the effects group environment on the evolution of dwarfs in the absence of massive (⁠|$M_{\ast } \gt 10^{10} \ \rm M_\odot$|⁠) galaxies, and to probe the hierarchical assembly in the low-mass regime. Although rare in the local Universe, these groups are expected to be more common at early epochs (Stierwalt et al. 2017). Therefore, their local analogues provide a unique window into the evolutionary processes that may operate more frequently on dwarf galaxies at higher redshifts (Kimbro et al. 2021). These systems also constitute small-scale structures (sizes <1 Mpc) within the broader cosmological framework, making their examination pivotal for gaining a deeper comprehension of matter clustering in domains where the standard cosmological model presents tensions with observations (Silk & Mamon 2012; Somerville & Davé 2015; Bullock & Boylan-Kolchin 2017; Sales, Wetzel & Fattahi 2022).

Previous studies indicate that ΛCDM can naturally account for the existence and properties of dwarf-only systems (Besla et al. 2018; Yaryura et al. 2020, 2023), such as associations and groups. However, it is worth noting that these important studies have certain limitations. They do not directly use fully hydrodynamical simulations with mass resolutions high enough to resolve low-mass galaxies (⁠|$M_{\ast } \lt 10^9 \ \rm M_\odot$|⁠). Additionally, they do not comprehensively explore the detailed mass assembly of different components within the low-mass regime represented by CGDs, nor conduct a thorough analysis of the formation and merging time-scales of dwarf groups. To complement these previous studies, in this work, we use the state-of-the-art TNG50-1 high-resolution cosmological simulation to investigate the properties, formation, and mass assembly of isolated CGDs within a cosmological context. We analyse their dynamical and halo properties and provide quantified predictions of the number-density statistics of these groups. We compare some observable physical quantities of simulated CGDs with observations at z ∼ 0 to address the consistency between the CGDs in TNG50-1 and the observational constraints. Our goal is to bring new insights into the baryon cycle of dwarf galaxies and the formation of low- to intermediate-mass galaxies through hierarchical merging.

This paper is organized as follows. In Section 2, we describe the simulations and observational data used in this work and the methodology adopted to identify and characterize the CGDs in the simulations and the observations. Our results are presented and discussed in Sections 3 and 4, respectively. Finally, in Section 5, we summarize the main findings and conclusions of our study. In Appendix  A, we provide extra figures, while in Appendix  B, we provide information about catalogues and supplementary material of this work.

2 DATA AND METHODS

2.1 IllustrisTNG cosmological simulations

In this work, we use data from the IllustrisTNG simulation suite (Nelson et al. 2019; Pillepich et al. 2019), which is a series of magneto-hydrodynamical simulations developed with the arepo code (Weinberger, Springel & Pakmor 2020) that incorporates a comprehensive galaxy sub-grid model. These models implement different AGN feedback modes, individual chemical element tracing, stellar evolution, and other relevant baryonic processes (Weinberger et al. 2017; Pillepich et al. 2018). The simulation self-consistently evolves the dynamics of dark matter, gas, stars, and magnetic fields within uniform periodic-boundary cubes (in the case of TNG50-1 run, a cube of side ∼51.7 cMpc). With many different cosmological simulations presented in the literature – with different resolutions, implemented physics and volume – a justification for why we choose one specific simulation is needed. TNG50-1 is one of the few hydrodynamic cosmological simulations able to resolve dwarf galaxies while simulating them in a larger cosmological context. Regarding the simulation resolution, the average initial stellar particle mass is 8.5 × 10⁴ M_⊙, the dark matter particle mass is 4.5 × 10⁵ M_⊙, while the Plummer equivalent gravitational softening lengths (ϵ) for both particles is 288 pc at z = 0. Thus, we choose TNG50-1 because of its spatial and mass resolution, enabling us to resolve groups of dwarf galaxies. In addition to that, among the modern cosmological simulations presented in recent years, the ones from the IllustrisTNG suite are among the simulations with the easiest access. This is an important feature to consider when carrying on work based on state-of-the-art simulations with limited computational resources. For example, a single snapshot of the TNG50-1 run occupies a disc space greater than 1000 GB; thus, it is not practical to download all the simulation data for local analysis. To facilitate data exploration, the IllustrisTNG team provides access to a server with a JupyterLab interface¹, which is an important difference relative to other public simulation suites.

To have a complete picture of the formation of CGDs, ideally, one would create multiple cosmological simulations with varying initial conditions – to avoid spurious results from stochastic variations (Genel et al. 2019; Borrow et al. 2023) – and varying cosmological parameters – to make predictions robust against changes in cosmology – however, this would be extremely expensive and is out of the scope of this work. Instead, we use IllustrisTNG suite of cosmological simulations rather than trying to find all the multiple initial conditions that could generate CGDs and doing individual simulations from scratch. The dwarf groups found in this work are direct by-products of the complex model employed in TNG simulations and, thus, can serve as tests for the model. Moreover, the CGDs found in TNG50-1 have evolved within a larger cosmological environment rather than complete isolation, which is essential if we want to study the large-scale environment of these objects. However, by studying a cosmological simulation that was not tailored to explore all the possible formation pathways of dwarf galaxy groups, we are subjected to a sampling bias. In other words, the simulated groups studied in this work are generated only from the limited set of initial conditions – e.g. density matter field – in TNG50-1 that can generate CGDs. Thus, the properties of dwarf groups that we present here may not represent all possible CGDs that are allowed to exist in a cosmological simulation within the ΛCDM framework. None the less, for the purpose of this work, we assume that TNG50-1 is homogeneously sampling the space of initial conditions that create CGDs.

Throughout this work, when we mention stellar masses, star formation rates, and other properties of individual simulated galaxies, we refer to these quantities computed in the specific aperture of 2 half-mass radius, which is commonly used in the literature. The r-band magnitudes of galaxies in the simulation are absolute magnitudes computed from the total luminosity of all stellar particles in the SDSS r-band, and dust attenuation is not taken into account. Since the extinction for these galaxies is small, it should not strongly affect our results and conclusions from the comparison with observations. Throughout the text, we use the terms ‘subhalo’ and ‘galaxy’ interchangeably since subhaloes in the simulation are supposed to represent galaxies in the real Universe. To be consistent with the TNG50 simulation, in this work we adopt the standard ΛCDM cosmology with Ω_m,0 = 0.3089, Ω_Λ,0 = 0.6911, Ω_b,0 = 0.0486, and H₀ = 67.74 km s⁻¹ Mpc⁻¹, from Planck Collaboration XIII (2016).

2.2 Nearest neighbours and environment

To avoid imposing a direct selection bias in the halo mass of the CGDs, we selected the groups based on the distance between galaxies. This approach has the advantage of being more easily reproduced in observations since halo mass is not directly observable. For each target galaxy in the simulation box, we computed the distances to the nth nearest neighbours and stored them in catalogues together with neighbours’ respective subhalo identifiers. In TNG50-1, we adopted a minimum stellar mass of |$M_{\rm \ast , min} = 10^7 \ \rm M_\odot$| for both the target and the neighbour galaxies because we are not interested in galaxies with stellar masses close to the resolution limit of the simulation. Besides, their observational counterparts are very difficult to detect at large distances and restricted to the local volume. Using information about distances, stellar masses, and star formation of nearest neighbours, we created the sample of simulated CGDs analysed in this work (see Section 2.3 for details on the selection).

We also use the catalogues described above to compute the densities of the environment around the groups. To analyse the large-scale environment, we define a volumetric density contrast following the approach by Muldrew et al. (2012), computed within a sphere of radius 1 Mpc centred on the galaxy of interest:

where N_g is the number of galaxies within the sphere around the target galaxy, and |$\bar{N}_{g}$| is the average number of galaxies that would be expected in the spherical aperture if galaxies were distributed uniformly throughout the simulation volume. To estimate |$\bar{N}_{g}$|⁠, we redistribute all the galaxies in new positions and sample the simulation volume with 10⁶ spheres. The rearrangement of the galaxies is done by sampling new positions from uniform distributions – one for each of the three coordinate components – spanning the values in the interval [0, 51.67] Mpc. The positions of the sphere centres are sampled in the same manner.

Alternatively, we also compute the volumetric dark matter density (ρ_1 Mpc) around the haloes of interest in the TNG50-1. We use snapshot data from the full simulation box and sum the masses of all dark matter particles in the simulation within 1 Mpc from the position of interest. Given that the computation of ρ_1 Mpc is expensive, we only do it for the halo samples analysed in Section 4.2, not all haloes in the simulation run.

2.3 Sample of simulated CGDs

The observational studies on groups of dwarf galaxies are mostly restricted to low redshifts (z ≲ 0.1), but with simulations, we can analyse these objects at higher redshifts. Therefore, we create samples of CGDs for all snapshots over the redshift interval of 0 ≤ z ≤ 0.5.

We are interested in compact groups, meaning we must define a group size measurement and use it as a selection criterion. We then define the group radius (r_group) as the radius of the smallest sphere that encloses all member galaxies. With the coordinates of galaxies, we compute r_group and group centres, using the Python package miniball², which has a recursive implementation of Welzl’s algorithm (Welzl 1991) to solve the smallest bounding sphere problem. In this sense, the group centre is simply the centre of the smallest sphere found to enclose all galaxies of the group, not necessarily the barycentre of the group. To avoid selecting groups suffering the influence of close galaxies, we select only isolated groups, which means that the closest galaxy must be at least 200 kpc away from the group. The distance from the group to neighbour galaxies is obtained using the nearest neighbour catalogue described in Section 2.2 and the coordinates of the group centre. We then generate catalogues of groups by iterative searching for groups of N star-forming dwarf galaxies that are enclosed by a sphere with radius |$\lt \rm 100 \ kpc$|⁠. The search is iterative in the sense that we first look for groups with three galaxies, then four galaxies, and so on, looking for groups up to a maximum of 10 star-forming dwarf galaxies. The criteria described above might lead to the selection of systems that are not exactly representative of the CGDs we want to study; thus, further sample cleaning is needed. For example, the initial catalogue of groups can contain a CGD with N members that are a substructure of other CGD with N + 1 members, and the first step is to remove groups within groups. We want to study CGDs as isolated structures, so we eliminate CGDs that pass the isolation criteria described above but are still part of a larger halo, such as a more massive group or a cluster.

In summary, to select a CGD in the simulation, we use the following criteria:

• (Compact) The group radius must be smaller than 100 kpc;

• (Dwarfs) All members of the group must have |$M_{\rm \ast } \le 10^{9.5} \ \rm M_\odot$|⁠;

• (Star-forming) All members must have |${\rm sSFR} \ge 10^{-11} \ \rm yr^{-1}$|⁠;

• (Rich) The group must have 3–10 members;

• (Isolated) Absence of neighbour galaxies within 200 kpc from the group centre;

• (Not substructure) One of the group members must be the primary³ galaxy of a halo;

The coordinates, stellar masses, and star formation rates used in sample selection were obtained from the TNG subhalo catalogues. Subhaloes with |$M_{\rm \ast } \lt 10^7 \ \rm M_\odot$| can be close to or inside CGDs selected in our sample. However, we do not consider these objects in our group selection because they may be at the limit of what we can call a ‘resolved’ galaxy in the simulation (number of particles N_particle ≲ 200) and because observations of dwarfs in this mass regime are very limited.

We want to study the properties of the haloes that host CGDs in simulations, and in TNG, the haloes are identified using a standard friends-of-friends (FoF) algorithm (Davis et al. 1985) with linking length b = 0.2 which is run over the dark matter particles. For simplicity, we only analyse groups that inhabit a single halo, in opposition to groups composed of more than one halo. This cut removes less than 5 per cent of the groups initially found at z = 0, so it does not introduce a strong bias in our final sample of simulated CGDs at z = 0. In total, 38 CGDs (∼97 % triplets and ∼3 % quartets) are found in the last snapshot of TNG50-1 – which corresponds to z = 0 – and in Section 3.2, we show examples of such groups. The total number of CGDs found in other snapshots is presented in Table 1 of Section 3.1.

2.4 Halo control sample at z = 0

To characterize the haloes of CGDs at z = 0, we can compare their properties with the overall halo population. To control the effect of halo mass in other halo properties, we must build a control sample of haloes with a similar distribution of this variable. We adopt as halo mass the total mass of the halo enclosed in a sphere whose mean density is 200 times the mean density of the Universe at the time the halo is considered (M_200,m).

To create the halo control sample at z = 0, we apply the propensity score matching (PSM) technique (Rosembaum & Rubin 1983), using the matchit⁴ package (Ho et al. 2011), developed in R language. This method allows us to select from a sample of haloes not hosting CGDs a control sample in which the distribution of the control variable (i.e. M_200,m) is as close as possible to that of those haloes hosting CGDs (see Fig. 1). For the matching, we adopted the Mahalanobis distance (Mahalanobis 2018) approach and the nearest-neighbour method. The matching process yielded a sample of 190 control haloes for the 38 haloes hosting CGDs, i.e. five control haloes for each CGD. We select more than one control halo for each CGD to have better statistics for the properties of the overall halo population in the same mass interval. We tested different ratios between treat and control samples (1:1, 1:3, and 1:10) and found that our conclusions remain the same.

2.5 Merger trees and cosmic evolution

To trace the evolution of the CGDs through different simulation snapshots, we used the subhalo merger trees constructed with the sublink algorithm (Rodriguez-Gomez et al. 2015). The algorithm links progenitor subhaloes to their unique descendants in future snapshots based on the particles that they have in common and a merit function that takes into account the binding energy of the descendants’ particles. With the merger trees, we can obtain past properties of the groups identified at z = 0 and properties of descendants identified in snapshots at z > 0. We trace progenitors of galaxies using the main progenitor branch of their merger trees, while their descendants are identified using main descendant branches.

The merger trees are used to compute relevant time-scales for the evolution of the groups, which are outlined as follows. First, we define t₁₀₀ as the maximum look-back time in the simulation when a CGD has r_group ≤ 100 kpc and all its members have |$M_{\ast } \ge 10^7 \ \rm M_\odot$|⁠, interpreted as the moment that the group formed and all dwarf galaxies began to interact in a compact region of space. We also use a formation time-scale based on halo assignment, t_FoF, defined as the maximum look-back time in the simulation in which all group members are subhaloes of the same FoF group, that is, all dwarfs are hosted by the same halo. Third, we define t_50,halo as the look-back time when the halo hosting a CGD reached half of its current virial mass at present, that is M_200,c(z_50,halo) = M_200,c(z = 0)/2, where z_50,halo is the correspondent redshift of t_50,halo. The values of t_50,halo are obtained from the values of the scale factor, a_form, given in the supplementary catalogue of halo structure from TNG50-1 (Anbajagane, Evrard & Farahi 2022). We also quantify the coalescence time-scale of CGDs, i.e. the amount of time that the galaxies in a group take to merge into one single galaxy. It is defined simply as the difference between the look-back time in which a group was identified (t_obs) and the maximum look-back time in which all group members have the same descendant subhalo (t_coalesc).

We also use the merger trees to trace the internal dynamical evolution of the groups until their final state at z = 0. Analysing the trees of all group members simultaneously, we determined the fate of their respective groups as one of the following:

• ‘Steady’ – all original group members remain the same at z = 0 and r_group did not grow significantly, i.e. |$r_{\rm group}(z=0)\lt 200\,$| kpc, which corresponds to two times the maximum radius adopted to select the CGDs;

• ‘Coalesced’ – all original members of the group merged into a single galaxy;

• ‘Dispersed’ – at least one original member of the group is not present anymore or r_group grew significantly, i.e. |$r_{\rm group}(z=0)\gt 200\,$| kpc;

• ‘Altered’– some members of the groups merged, or a new galaxy is part of the group;

The four classes described above cover all the possible final states of CGDs in the simulation.

2.6 Projected quantities

To make fair comparisons between the physical properties from simulations and observations, it is convenient to project the simulated groups in two dimensions. To do this, we project the positions and velocity vectors of galaxies in the simulations on to 150 randomly oriented 2D planes. For each plane on to which a group is projected, we compute the quantities of interest (e.g. radius, velocity dispersions, etc.) and the 25th, 50th, and 75th percentiles of their respective distributions. In this manner, we do not privilege any specific orientation of the CGDs, and we consider variations on the projected quantities that can arise purely from differences in orientations.

To obtain group masses, we use a virial mass estimator based on projected velocities and positions of galaxies. We use the projected mass estimator described in Heisler, Tremaine & Bahcall (1985), given as follows:

where R_p,i is the projected distance between an ith group member and the geometric centre of the group; Δv_i is the difference between the line-of-sight velocity of ith group member and the mean velocity of the group members; N is the number of group members; G is the gravitational constant; f_pm = 20/π and α = 1.5 are constants as adopted in Stierwalt et al. (2017).

With the multiple projections on simulations, we are able to access how the group mass estimator adopted is sensible to specific configurations of CGDs. As shown in Fig. 5, Section 3.2.1, the variation in the values of M_PME,group can be large (greater than 0.5 dex). For each simulated CGD, we adopt the 75th (25th) percentiles of the distribution of masses estimated from the 150 projections as the upper (lower) uncertainties in M_PME,group. However, for the observational data, we have a single projection, so we compute the mean interquartile ranges (IQRs) of M_PME,group for all simulated CGDs and assume the uncertainty to be symmetric and equal to half of the mean of IQRs. Here, we consider that the projection effects on the M_PME,group estimates in observations are the dominant source of error. These effects are similar to those in simulated systems.

To facilitate the comparison of our results with literature, we also computed velocity dispersions and group radius as defined in other studies addressing dwarf systems. As a measure of group physical size, we computed the projected inertial radius (R_I) analogously to Tully et al. (2006) and Yaryura et al. (2020) (hereafter Y20):

where r_i is the 2D projected distance of a galaxy from the system centroid, and the sum for each group is computed over all N members. We also compute a measure of the groups’ velocity dispersion using an analogous expression to that of Y20:

where v_i is the difference between the line-of-sight velocity of a galaxy and the group mean. As a complementary measure to σ_Y, we also compute the velocity dispersion of groups using the gapper estimator (Wainer & Thissen 1976; Beers, Flynn & Gebhardt 1990). This robust scale estimator is based on the gaps of an order statistics, x_i, xi + 1,..., xn, with gaps defined by

and approximately Gaussian weights given by

The estimator of scale then is given by

where x_i is the ith element of an ordered set of differences between the line-of-sight velocities of the galaxies and the group mean. We found σ_Y values that are strongly correlated with σ_WT values, with both estimators returning similar velocity dispersions, with a relative difference of |$\sim 13~{{\ \rm per\ cent}}$| in both simulated and observed CGDs. When computing the velocity dispersions for observations, the scaling factor |$1/(1 + \bar{z})$| multiplies the right-hand side of equations (4) and (7), where |$\bar{z}$| is the mean redshift of the group members.

2.7 Observational sample at z ∼ 0

To compare the simulation results with observations, we built a sample of CGDs using the spectroscopic catalogue from Sloan Digital Sky Survey Data Release 18 (Almeida et al. 2023; SDSS-DR18) in the redshift interval 0.01 ≤ z ≤ 0.1. The search for CGDs was made using the CasJobs asynchronous data base tool for long, complex queries in the SDSS Catalog Archive Server (CAS).⁵ We selected groups with projected radii R_group < 100 kpc and imposed the limits of |$M_\ast \le 9.5 \ \rm M_\odot$| and M_r ≥ −20 for the stellar masses and r-band absolute magnitudes of group members, respectively. The upper limit in stellar mass is compatible with the limit we adopted when searching CGDs in the simulation. On the other hand, the lower limit in r-band absolute magnitude guarantees that galaxies in the observational sample will not be brighter than the brightest galaxies in simulated CGDs, and we apply this limit to make comparisons between the simulations and observations more meaningful. As is shown in Fig. 2, all dwarf galaxies inhabiting CGDs in TNG50-1 have absolute magnitudes M_r fainter than −20, justifying the adopted magnitude limit mentioned above. The additional limit in magnitude also accounts for the fact that not all galaxies in the SDSS have estimates of stellar mass, so we can select them using M_r.

The criteria described above to select CGDs in observations are similar to those employed in TNG50-1, but some differences are important to mention. The projected radius cannot exceed the real group radius in three dimensions. However, due to projection effects, observed CGDs can have a projected radius that is below the upper limit used in the selection (R_group < 100 kpc) but a real radius that is actually larger. Additionally, chance alignments can significantly contaminate observational samples of compact groups (Taverna et al. 2022). Thus, to search for real compact groups, we also impose a maximum difference between line-of-sight velocity (Δv) of group members of 300 km s⁻¹.

In the spectroscopic catalogue of SDSS, some sources are not individual galaxies but clumps of larger galaxies. Thus, after we built the group catalogues using CasJobs queries, we visually inspected all CGD candidates and removed false group members from the sample, in which a clump of a larger galaxy was interpreted as an individual galaxy. At the end of the sample selection process, the observational sample of CGDs was composed of 58 groups (∼93 % triplets and ∼7 % quartets). A table with their properties is provided as supplementary material for this work (see Appendix  B).

Once we had a clean sample of CGDs from SDSS data, we cross-matched the SDSS CGDs with data from the Arecibo Legacy Fast ALFA (ALFALFA) survey (Haynes et al. 2011, 2018) to obtain information about the atomic neutral gas content of observed groups. The ALFALFA extragalactic neutral hydrogen (H i) source catalogue provides estimated masses for individual detections based on the 21-cm line flux density. Thus, we obtain the H i gas masses (M_H i) of CGDs by summing the individual masses of all detections within 1.5R_group and |$\Delta v \lt 500\, {\rm km\, s}^{-1}$| from the centre of each group. Since the SDSS spectroscopic catalogue has only a partial overlap with the ALFALFA extragalactic catalogue, we do not have information about M_H i for all groups in our observational sample. More precisely, we only present results about the neutral hydrogen content of 20 CGDs in SDSS (see Section 3.2.1).

We consider as CGD candidates all groups having at least three members with spectroscopic redshifts, which introduces a limit in the magnitude of the faintest members. Thus, throughout the work, we always use the magnitude of the third brightest group member (⁠|$M_{r,\rm 3rd}$|⁠) as a reference when comparing the faintest members in simulated and observed groups. It is possible that fainter objects are part of the CGDs in the observational sample, but we are only considering group members with spectroscopic redshifts because we want to avoid contamination from chance alignments and analyse only groups of galaxies that are truly close in real space.

To compare some properties of simulations and observations more directly, we created sub-samples of simulated and observed CGDs. We match groups in simulations with groups in observations with the PSM technique described in Section 2.4 using as control variables their total absolute magnitudes (⁠|$M_{r,\rm tot}$|⁠) and third brightest member magnitudes, both in r band. We use both magnitudes to constrain, at the same time, the luminosity of the faintest group members and of the whole group. Due to the different distributions of |$M_{r,\rm tot}$| and |$M_{r,\rm 3rd}$| in both simulations and observations, only 9 CGDs are paired for each sample.

The magnitudes of galaxies selected in SDSS are Petrosian magnitudes corrected by extinction. The luminosity distance used to compute the absolute magnitudes is obtained using the FlatLambdaCDM class of astropy cosmology package⁶ and the parameters mentioned at the end of Section 2.1. We also apply K-correction to the magnitudes using the code kcorrect (Blanton & Roweis 2007). Magnitudes in the B band are computed using the following Lupton (2005) transformation⁷:

where M_g and M_r are absolute magnitudes in the g and r band, respectively. To convert the magnitudes to luminosities, we consider M_B = 5.44 for the Sun (Willmer 2018).

2.8 General caveats

It is important to note that the simulated systems analysed in this work are taken as analogues of objects found in observations, not perfect representations. Cosmological simulations assume recipes to model physical phenomena that occur in scales smaller than their particle resolution; these coarse-grained models may not be appropriate for all situations (Crain & van de Voort 2023), resulting in galaxies with spurious properties. If we assume that the recipes work well as effective models, then the coarse-grained physics may not be a problem. In other words, we can expect the cosmological simulations to fail in reproducing some properties of galaxies or galaxy groups. However, we can also expect that the simulations will give us useful insights about the formation and evolution of galaxy populations, even when the objects simulated are not perfect representations of observed galaxies. In addition, here we are focusing on the analysis of IllustrisTNG simulations, which adopt a specific code, cosmology, and sub-grid physics. Other simulations may present different results about the nature and properties of CGDs. Thus, further analysis is needed in order to address this dependence on simulation parameters.

Regarding the number counts statistic of CGDs, we expect the dominant effect in our observational sample to be the limiting magnitude that severely limits the observed volume for dwarf galaxies. Due to our choice for high-resolution simulations, we also have a limited simulated volume, which may influence our results. Simulations with similar resolution to TNG50-1 but large volumes (volume |$\gt 100^3 \ \rm Mpc^3$|⁠) are not available at the time, thus, the dependence of our results with simulation volume may be tested only in the future. Additionally, cosmic variance may also play a role in our results since the volume we are sampling – in simulations and observations – is relatively small.

3 RESULTS

3.1 Space number density

After applying the selection criteria described in Section 2.3, we obtained catalogues of simulated groups for different snapshots corresponding to different redshifts. In Table 1, we show basic information about each catalogue of CGDs. Catalogues with properties of simulated CGDs will be provided as supplementary material (see Appendix  B).

The space number density is an important quantity to be evaluated when studying populations of extragalactic sources. Thus, we estimate the comoving number density of CGDs in TNG50-1 and interpret it as a prediction of the cosmological simulation. In the top panel of Fig. 3, we show the evolution of the comoving number density of CGDs. To estimate the errors for the number densities found in TNG50, we create a ‘super-box’ of size ∼310 cMpc by randomly rotating 6³ identical TNG50 boxes. Then, we assemble them in a cube to create a larger box, 216 times the size of TNG50. We sample the super-box with 10⁵ cubes of the same size as TNG50, computing the number density of CGDs. The errors shown in the top panel of Fig. 3 are the 16th and 84th percentiles of the log (n) distribution obtained with this sampling. The number density estimated by counting the number of CGDs in the TNG50-1 standard box (shown as cyan circles in the figure) coincides with the median value of our sampling in all snapshots, indicating that the method is working as expected.

Our results for the number density at z = 0 are in agreement with previous theoretical work from Y20, who found |$n \approx 10^{-3.4} \ \rm Mpc^{-3}$| for the systems they define as ‘Groups’, with three members and adopting similar mass limits in the stellar mass of group members as those used in this work. As can be seen in Fig. 3, the density of groups decreases towards the present, with a more pronounced decline at z ≤ 0.26. From the values of the number densities in Fig. 3, we can also see that CGDs are scarce at redshifts z < 0.1, which is expected from ΛCDM cosmologies and is in agreement with previous studies which suggest that groups and associations of dwarf galaxies are rare (Stierwalt et al. 2017; Yaryura et al. 2020).

In the bottom panel of Fig. 3, we show the number density as a function of a minimum absolute magnitude, in this case, the r-band magnitude of the third brightest group member (⁠|$M_{r \rm ,3rd}$|⁠). For most simulated CGDs, |$M_{r \rm ,3rd}$| is the magnitude of the faintest galaxy in the group. We see that the number of groups decreases monotonically with |$M_{r \rm ,3rd}$|⁠, with a more pronounced decline for |$M_{r \rm , 3rd} \lesssim -15$|⁠, meaning that groups with fainter |$M_{r \rm ,3rd}$| are relatively more common. An interpretation for this result is that we should observe more groups that contain fainter galaxies (⁠|$M_{r \rm , 3rd} \gtrsim -15$|⁠), implying that some observed dwarf galaxy pairs may have fainter companions that are not detected in current spectroscopic surveys. However, the size of the simulation box of TNG50-1 is relatively small (L_box = 35 cMpc h⁻¹), meaning there may be an insufficient sampling of the bright end of luminosity functions estimated from this simulation. Ideally, our result for TNG50-1 should be compared with other simulations with the same resolution – to resolve dwarf galaxies – but with a much larger volume to check for variations in the number density of CGDs. Nevertheless, it is relevant to study if TNG50 correctly estimates the number of CGDs that we should find in a flux-limited survey, as we discuss in Section 4.1.

In the bottom panel of Fig. 3, we show how the number density of groups as a function of |$M_{r, \rm 3rd}$| varies with redshift. The largest difference is between the number densities at z = 0 and z = 0.26, as is expected from the steeper log (n) decline in this redshift interval on the top panel of the same figure. Except for the bright end (⁠|$M_{r \rm ,3rd} \lesssim -16$|⁠), the difference in number density does not seem to depend strongly on the magnitude of the third brightest member, with the median of differences being Δlog (n) ≈ 0.4. This shows that the number of groups decreases by more than half from z = 0.26 until the present, almost uniformly for all groups with |$M_{r \rm ,3rd} \gtrsim -16$|⁠. Thus, at least in the TNG50-1 simulation, the formation and destruction rates of CGDs do not seem to be strongly affected by the luminosity of the faintest group members.

3.2 CGDs at z = 0

In what follows, we focus on the properties of the 38 CGDs found in the latest snapshot of TNG50-1, corresponding to z = 0. In Fig. 4, we show examples of such CGDs. As one can see in the different panels, the stellar component of the groups resembles images of observations of CGDs in the optical (e.g. Stierwalt et al. 2017). As for the gas distribution in these groups, it is clear that it can be very scattered and complex, in agreement with HI observations of dwarf pairs by Pearson et al. (2016, 2018); Luber et al. (2022).

Given that dwarf galaxies are challenging to detect at large distances, the observational samples of CGDs are mostly restricted to the local Universe. Therefore, predictions for the properties of CGDs found at earlier times in the simulation cannot be properly constrained by current observations. Nonetheless, we explore the predictions of TNG50-1 about CGDs at higher redshifts (0 < z ≤ 0.5) in Sections 3.3 and 4.3.

3.2.1 Observable properties

One goal of this work is to use simulations to understand the formation of CGDs found in observations. But first, we investigate whether the general properties of CGDs in simulations are in accordance with what is expected from observations. To achieve this, we use two observational samples as reference for the properties of CGDs. The first observational sample is built in this work (see Section 2.7), and the second observational sample is from S17. Neither are supposed to be complete in volume/magnitude; thus, the observational constraints that they impose are also restricted to a small luminosity/mass range. None the less, we can use the ranges of the properties of the groups in the observations as weak constraints for what we can expect cosmological simulations to reproduce.

As shown in the top panel of Fig. 5, we have a diagram region where both simulated and observed CGDs overlap within the error bars. Specially groups from S17, lie within the 95  per cent contour of the simulated groups. Regarding the sub-samples in observations and simulations, which were matched by magnitude (see Section 2.7), we see that they are almost all within the limits of the 95  per cent contours of each other parent samples. However, TNG50-1 produces many low-mass groups with no counterparts in the observational sample. The incompleteness of our observational sample partially causes this disagreement, but the possibility that the simulation is producing spurious low-mass CGDs also cannot be excluded. Unfortunately, the current data is not enough to confidently discard any of the hypotheses. Additionally, the simulation does not produce CGDs as massive as the most massive CGDs in observations, which is more challenging to reconcile with observations but may be due to the limited simulation volume of TNG50-1. Since more massive groups will be relatively less numerous, they may be under-sampled in the cosmological volume of TNG50-1 (13.8 × 10⁴ cMpc³). On the other hand, in our observational sample, the more massive CGDs (⁠|$M_{\rm PME,group} \gt 10^{12} \ \rm M_\odot$|⁠) have a median absolute magnitude |$M_{r,\rm 3rd}=-17.8$|⁠, which, with a magnitude limit of 17.77 in the r-band, corresponds to a maximum detection redshift of z ≈ 0.0287. Considering approximately 7500 |$\rm deg^2$| as the observed sky area of the SDSS Legacy spectroscopic sample (Abazajian et al. 2009), the cosmological volume observed up to z = 0.0287 would be more than 12 times the volume of TNG50-1. Thus, the bright end of our observational sample is selected within a larger volume than the simulation.

Regarding the maximum difference in stellar mass of group members (Δlog (M_*)_max), we see that the most simulated CGDs lie within the limits imposed by S17. Still, the agreement between simulated CGDs and observed CGDs is limited. This means that the analogue systems only partially capture the stellar mass ranges within observed groups. However, we stress that the values of Δlog (M_*)_max may be directly affected by the stellar mass or magnitude selection criteria.

As we can see in Fig. 5, for almost the entire range of total masses, the observations lie in the same regions as the simulations. This is an important result because group total masses are never used as a selection criterion to find CGDs in the simulation. We only look for associations of dwarf galaxies by looking at their distances to each other. Thus, the similarity of the group masses in simulations and observations indicates that, at least in terms of mass, the analogue CGDs in the simulation reproduce the scale of these systems in observations.

Regarding the gas content of the groups, if we consider H I mass fractions (Fig. 6), we see hints that the simulated CGDs may have lower H I fractions at fixed stellar mass. However, this is only a first approximation of the total HI mass fraction in groups, and the observational samples are even smaller than those shown in Fig. 5, thus, further work and larger samples are needed to draw more robust conclusions.

3.2.2 Host halo properties

In this work, we consider the simulated groups as possible models of the real CGDs. Thus, it is relevant to analyse the basic properties of the haloes that host CGDs in simulations to access quantities that can be difficult to measure with observations. By physically characterizing the simulated CGDs, we can obtain insights about possible formation scenarios for their observed counterparts and how they fit in the larger context of galaxy formation and evolution. It is expected that properties of haloes scale with their mass (Yang, Mo & van den Bosch 2008), so it is essential to consider this scaling during analysis. Thus, we compare the properties of haloes hosting CGDs with control haloes having similar masses and not hosting a CGD. In Fig. 7, we show scatter plots of halo mass (M_200,m) versus four different halo properties for the control sample and the CGD hosts sample. The properties are defined as total stellar mass (M_*,halo), which is all stellar mass in the halo; total gas fraction (f_gas,halo), which is the fraction M_gas,halo/(M_gas,halo + M_*,halo); 50  per cent-mass age (t_50,halo), which is the look-back time when the halo reached half of its current virial mass, that is M_200,c(z = z_50,halo) = M_200,c(z = 0)/2; and a baryonic mass fraction (f_b,halo), which is the fraction (M_gas,halo + M_*,halo)/M_200,m.

The median differences in total stellar mass are of the order of −0.4 dex, with CGDs having less than half of the stellar mass present in the control haloes. Regarding how gas-rich the CGD haloes are, we see that they tend to have |$\sim 10~{{\ \rm per\ cent}}$| more gas than other haloes of similar mass, with this difference being slightly higher for CGDs in the most massive bin (11.45 ≲ M_200,c ≲ 11.8). Additionally, analysing the difference in the halo ages, we see that, on average, CGD host haloes assembled ∼4 Gyr later than other haloes of similar mass, with a smaller difference (∼2.5 Gyr) for the most massive bin. By analysing the fraction of baryons in both halo samples, we see that CGD hosts have similar amounts of baryons to control haloes, indicating that that the differences in halo properties mentioned so far are not a by-product of differences in the baryonic mass fraction.

In summary, we found that CGDs tend to inhabit haloes with less stellar mass and higher gas fractions than the control sample. The haloes that host the groups also assembled half of their virial mass later when compared to haloes of similar mass. The tendencies on halo properties described so far are significant over all the halo mass range probed and suggest that CGDs inhabit haloes with a different assembly history from the overall population of haloes within the same mass regime at z = 0. We explored this possibility further in the results presented in Section 3.2.3 and discuss them in Section 4.2.

From the results presented so far, we see that the CGD host haloes are a special population because they have specific properties that deviate from the overall halo population at the same redshift and mass interval. However, the sample of CGDs in TNG50-1 is small (N < 50) and may be biased. Borrow et al. (2023) draws attention to the stochastic effects on the properties of individual galaxies, and they show that re-simulations of a galaxy can result in differences up to 25 per cent in stellar mass. However, median values for scaling relations should not be significantly affected. Thus, the predictions for populations of galaxies are more robust against stochastic variations than predictions for individual galaxies. We can apply this logic to the population of CGDs found in this work and argue that the overall median values of halo properties can possibly be robust against the stochastic effects of a single simulation run (TNG50-1).

3.2.3 Assembly history

To fully characterize the CGDs found at z = 0 in TNG50-1, we can go beyond their current properties and analyse their mass assembly history (MAH), readily available in a cosmological simulation. To understand the mass assembly of CGDs, it is useful to compare them with the control sample of haloes with similar mass. Thus, in Fig. 8, we present the mass growth of control and CGD haloes over the whole time span of the cosmological simulation.

As we can see in the normalized total MAH, the haloes hosting CGDs got their masses later when compared to the haloes in the control sample. In other words, CGDs, on average, assembled relatively more of their mass in the last few Gyr. As is shown in the gas MAH, the CGD and control haloes also only reached similar values of M_gas,halo recently. The median gas accretion rate of CGDs was higher for several Gyr, and this is illustrated by the steeper slope of the CGDs gas MAH since z ∼ 0.5. By analysing the evolution of gas distribution in simulated CGDs (see videos in Supplementary Material), it becomes clear that the accretion of gas in these systems happens hierarchically, that is, from the merger of similar low-mass haloes that only recently assembled into a single structure. This would be contrasted with the accretion of a single, more massive halo continuously receiving gas from its surroundings and the cosmic web. This crucial difference in the accretion mode could explain why the total and gas MAHs look steeper than in the control sample for z ≲ 0.5. Although the total and gas masses were accreted at higher rates and reached similar values for control and CGD haloes at z = 0, the picture differs for the stellar masses. As we can see in the stellar MAH, on average, CGDs have always had fewer stars than other haloes of similar mass at z = 0. The results presented in this subsection indicate that the stellar and gas fractions of CGD haloes differs from that of other haloes with similar mass because of the particular assembly history of these compact groups. More specifically, the halo properties of CGDs found at z = 0 seem to be the result of a late hierarchical assembly of low-mass (⁠|$M_{\rm 200,m} \lesssim 10^{10} \ \rm M_\odot$|⁠) haloes.

We can also study the history of the CGDs to quantify how long these groups can exist in the simulation. In Fig. 9, we show a histogram of group formation time indicators derived from the analysis of merger trees (see Section 2.5), which can be adopted as criteria to determine the moment a group can be considered to be ‘formed’. In practice, both time-scale indicators presented in Fig. 9 are challenging or even impossible to infer in observations, but through this work, we are going to consider the t₁₀₀ estimator as the closest to what we can call an age of the CGDs. We choose it because it is defined more straightforwardly than t_FoF and, in principle, can be more easily inferred from observations of the kinematics of galaxies within the groups. As is shown in the histogram of t₁₀₀, the majority of CGDs formed compact associations of dwarf galaxies in the last 2 Gyrs, but they also can be much older, with two of the groups being older than 10 Gyr. The results for t_FoF are quite different, with most groups forming between 2 and 4 Gyr ago, where the median for all the CGDs is ∼3 Gyr. These look-back times coincide with the epochs in which the CGDs hosts present steeper slopes in their normalized total MAH (see Fig. 8). Thus, t_FoF may capture the epoch of the hierarchical assembly that builds up CGDs haloes in the simulation. While t₁₀₀ may better trace the time since which galaxies start to interact more strongly inside the groups.

3.2.4 Position in the large-scale structure

Groups of galaxies in the real Universe and simulations are not formed in random positions of space, they are usually found close to matter clustering in the large-scale structure (LSS). Thus, it is relevant to explore whether low-mass groups – such as CGDs – follow this trend and if they are situated in specific regions of the LSS. The advantage of doing this kind of analysis in simulations is that we can access all positions in real space, precisely computing local densities and distances between galaxies in the simulation. For this brief exploration of the location of CGDs in the LSS, we use a proxy of the underlying matter density, the local galaxy density contrast (δ_1 Mpc), described in Section 2.2. This proxy of total matter density has the advantage of being accessible also in observations, especially with the advent of modern facilities, which will map the large-scale structure of the Universe through spectroscopic surveys, e.g. DESI (DESI Collaboratio) and Euclid (Euclid Collaboration 2024).

As is qualitatively shown in the left and centre panels of Fig. 10, the CGDs tend to be situated away from the centres of voids and knots in the cosmic web, differently from the control halo sample. When we look at the distribution of the adopted large-scale environment proxy – right panel of Fig. 10 – we see that CGDs are concentrated in a narrow range of local density, while the haloes in the control sample are distributed over a wider range of δ_1 Mpc values. Although both CGDs and control haloes avoid the highest densities, only the CGD hosts avoid the lowest density extremes in the simulation volume, being primarily found in regions with densities slightly above the 25th percentile. The avoidance of high-density environments of the cosmic web (e.g. knots) can be expected because low-mass haloes (⁠|$M_{\rm 200,m} \lt 10^{12} \ \rm M_\odot$|⁠) may be easily accreted into the larger haloes that inhabit these environments, becoming subhaloes of them or being disrupted. However, the fact that CGDs also avoid low-density environments (e.g. voids) is less straightforward. That CGDs do not inhabit the lowest density environments at present implies the possibility that these groups can only form above a certain threshold of local matter density, while control haloes do not have this restriction. Given that the control haloes are paired by mass only at z = 0, the differences in the large-scale environment at z = 0 may be a consequence of past local over-density and/or different halo assembly histories (see Sections 3.2.3 and 4.2).

3.3 Coalescence fractions and time-scales

One way to understand how the population of CGDs evolves with redshift is to look at their fate in the last snapshot of the simulation. We categorize the final states of CGDs found at 0 < z ≤ 0.5 using the merger trees of their group members and divide them into four classes: Steady, Altered, Coalesced, and Dispersed (see Section 2.5). In Fig. 11, we show the fraction of each class as a function of redshift. The fractions are computed for each group catalogue in a given redshift, that is, for each z, the sum of the fractions of the four classes is 1. In this way, we can see the fractions as probabilities of a CGD at z having a certain final state at z = 0.

The most obvious trend that we see in Fig. 11 is that the number of ‘Steady’ groups increases rapidly as we approach z = 0, and the number of ‘Coalesced’ groups increases slowly with redshift. Both trends are related to the coalescence time-scale of CGDs. For example, ∼30 per cent of the CGDs found in z = 0.5 coalesced into a single galaxy, while only ∼2 per cent remained unaltered (‘Steady’) up to z = 0. Meanwhile, the fraction of dispersed groups remains below 15 per cent for most of the redshift interval presented. From this, we can state some main results: first, even though it is a minority, the simulation predicts that a tiny fraction of CGDs can survive for more than 5 Gyr with all its original group members; second, the most common final state for CGDs identified at z > 0.2 is to be altered, that is, accrete new galaxies (massive or not) in the group and/or have mergers between original members of the group; third, the dismantling of groups is not common, meaning that a small fraction of them are transient unbound structures. These fractions presented in Fig. 11 provide important information to infer how much time we can expect the CGDs in observations to serve as isolated laboratories to study environmental effects in dwarf-only contexts.

Looking more specifically at the class of coalesced groups, we can determine the coalescence time-scales of CGDs and provide typical time-scales for the maximum interaction time of the dwarf galaxies in these environments. As shown in Fig. 12, most of the CGDs found at z ≤ 0.5 take more than 1 Gyr to completely merge into a single galaxy, with a typical time-scale of ∼3 Gyr. This is a reasonable amount of time in terms of galaxy evolution, allowing for the dwarf galaxies in CGDs to have multiple episodes of star formation that may be triggered due to close encounters or fly-bys (Martin et al. 2021).

We also investigate if the time it takes for a CGD to merge fully has any dependence on the radius of the group at the time it is identified. As shown in Fig. 12, there is no clear trend between the coalescence time and r_group over the whole range of group radius. Yet, we can see a tentative trend for |$r_{\rm group} \lesssim 40 \ \rm kpc$|⁠, where the larger the radius, the longer it takes for the group to coalesce.

4 DISCUSSION

4.1 CGDs in TNG50-1

In Section 3.1, we presented the number density of simulated CGDs and showed how it evolved within the redshift interval of 0 ≤ z ≤ 0.5. This is an important result for the statistics of CGDs, being relevant to predicting the number of observed CGDs, and possibly bringing useful insights for improvement in the numerical models. Considering this, we use the z = 0 CGD number density estimated in the TNG50-1 to predict the number of CGDs that may be found, given different r-band magnitude limits and a fixed sky coverage⁸ of 7500 |$\rm deg^2$|⁠. In Fig. 13, we show the predicted number of CGDs (N_groups) that should be found in observations according to the results we describe in Section 3.1. For the comparison with our observational sample, we considered only absolute magnitudes |$M_{r \rm , 3rd}$| where the maximum redshift of detection (z_max) is greater than 0.01, to have a reasonable sampling volume. Since our observational sample is constructed from the spectroscopic catalogue of the SDSS Legacy Survey, our reference magnitude limit for comparison in Fig. 13 is m_r = 17.77 (shown in blue). This causes the confrontation with the observations in Fig. 13 to be limited to a short magnitude range; however, it is still helpful to check whether, even in this range, we see some disagreement. By construction, the observational sample of CGDs in this work gives only a lower limit to the estimated number of CGDs in the local Universe, and we see that the predictions from TNG50-1 are roughly in agreement with such estimates, at least in the same order of magnitude. Still, it is essential to have larger samples of dwarf galaxy groups complete to fainter magnitudes (M_r > −16), to more strongly constrain current models of galaxy formation and evolution in the low-mass regime. Current (Darragh-Ford et al. 2023; Luo et al. 2023) and future deeper surveys – such as those conducted by Vera Rubin Observatory (Ivezić et al. 2019), DESI (DESI Collaboration 2024), and Euclid (Euclid Collaboration 2024) – will be essential for this goal. For example, the DESI Bright Galaxy Survey (Hahn et al. 2023) will have a magnitude limit of M_r < 19.5 for its ‘Bright’ target sample, allowing for the spectroscopic confirmation of fainter and further CGDs not present in our observational sample (see orange line on Fig. 13). Even without spectra, deep photometric surveys will show low surface brightness tidal features resulting from interactions, which are important signals to infer the proximity of dwarfs in pairs or groups.

The number density of CGDs can be viewed as a prediction from the TNG50-1 cosmological simulation, which may be helpful to test it as a model for galaxy formation and evolution. If the number density of CGDs is above the values expected from observations, it may indicate that this kind of group is created too easily in the cosmological simulation. On the other hand, if the number density is below the values expected from observations, then the models employed in IllustrisTNG may suppress the formation of CGDs. None the less, if simulation and observations agree within the error bars, it does not necessarily mean that the physical model of the simulation is completely correct; that is, it reproduces observations for the right reasons. Still, it is a good indicator in favour of it, as long as the number densities of other populations of galaxies and groups are also correctly reproduced. The analysis we make in this work is only part of a bigger series of ongoing comparisons that the community may carry in the literature (Jackson et al. 2020; Oppenheimer et al. 2020; Donnari et al. 2021; Goddy et al. 2023; to cite a few), and we hope the predictions presented here will be useful in future observational studies exploring the faint end of galaxy groups.

Regarding the properties of CGDs compared to observations at z ∼ 0, we see that even without directly selecting the CGDs by group mass, we could obtain reasonable agreement for the mass scales of the simulated groups (see Fig. 5). This means that, at least in terms of stellar and total masses, TNG50-1 can produce acceptable analogues for a restricted interval of these parameters. The fact that the simulation can predict the existence of a population of CGDs without them being explicitly put in the simulation is a positive point for the model. Thus, it is reasonable to try to use the simulation to understand the formation and nature of this kind of group.

The masses that we find for the CGDs in TNG50-1 are comparable to the theoretical results from Y20, where they find virial masses in the range |$\sim 10^{10}$| to |$\sim 10^{11.7} \ \rm M_\odot$| for groups of dwarfs⁹ in which all members inhabit the same halo. These groups described by Y20 also have similar mass-to-light ratios in the B band (see their Figs 4 and 7) to the majority of CGDs we found in TNG50-1 (10 ≲ M/L ≲ 100, see Fig. 15). Additionally, the typical inertial radius (R_I) that Yaryura et al. (2023) find for dwarf groups are of the same order of magnitude of R_I values that we find for the CGDs in TNG50-1 (see their Fig. 2 and our Appendix  B).

Our work independently demonstrates that compact groups of dwarf galaxies can emerge naturally in a ΛCDM framework, which evolves not only dark matter but also takes into account the baryonic processes and self-consistent hydrodynamical evolution of gas.

4.2 Formation of CGDs at z = 0

By comparing the host haloes of CGDs with control haloes of similar mass, we are trying to obtain insights into specific properties that CGD host haloes must have. In other words, we can explore what differentiates the haloes of CGDs at z = 0 from the general population of haloes at the same mass scale and redshift. As we present in Sections 3.2.2 and 3.2.3, the host haloes of CGDs have lower amounts of stellar mass, have higher total gas fractions and assemble all their mass later when compared with haloes of similar mass. Since we systematically find lower stellar masses in the host haloes of CGDs at present, it is important to understand if CGD hosts always deviate from the overall halo population through cosmic time. In Fig. 14, we see that, on average, the main progenitor haloes of CGDs deviate noticeably from the expected evolution for their masses at z ≳ 0.26. As expected, this redshift interval is compatible with the epochs at which the CGD were first hosted by a single halo, the look-back time traced by t_FoF (see Section 3.2.3).

From the results presented so far, we can interpret that the host haloes of CGDs at z = 0, on average, only reached their halo masses at present-day due to recent – last ∼3 Gyr – hierarchical merging, which increased their masses quickly. However, the stellar component did not grow at the same rate as the total mass, making CGD haloes more deficient in stellar mass. If the simulated CGDs are good representations of their observational counterparts, then we should expect both samples to have similar mass-to-light ratios. In Fig. 15, we present a diagram of the group masses versus luminosities in the B-band (L_{B, group}) and a histogram of mass-to-light ratios in the same band. It is evident from this diagram that most CGDs in the simulation have mass-to-light ratios spanning 10 to 100, with a median value of M/L ∼ 38 (25th and 75th M/L percentiles being ∼24 and ∼89, respectively), comparable to the median of groups from S17 (M/L ∼ 50). The observational sample of CGDs presents similar values of mass-to-light ratio to the simulated groups, with a median of M/L ∼ 42 (25th and 75th M/L percentiles being ∼14 and ∼83, respectively). Thus, the mass-to-light ratios estimated for observed groups could possibly be explained by a recent hierarchical formation, driving rapid growth in total mass with a delayed and relatively inefficient build-up in stellar mass. Additionally, by visually comparing the gas distributions of the CGDs and control haloes, we see that the gas can be greatly scattered in the groups by the interaction between the members. This may create an even larger delay in the stellar assembly of haloes hosting CGDs as compared to the control haloes, which typically have a single dominant galaxy accreting gas and smaller satellites, thus concentrating the gas into a single region.

In the simulation large-scale context, the CGDs at z = 0 can be viewed as structures that formed recently, and their formation can be interpreted as the direct product of a late hierarchical mass assembly of haloes inhabiting regions of the cosmic web with low-to-intermediate density – away from the centres of either knots or voids. A recent study by Martizzi et al. (2020) analyses the correlation between galaxy stellar masses, halo masses, and large-scale environment in the TNG100 simulation. The authors show that for a fixed halo mass, galaxies with stellar masses lower than the median are more likely to be found in voids and sheets, while galaxies with stellar masses higher than the median are more likely to be found in filaments and knots. Naturally, it is tempting to ask if CGD haloes had a different assembly history due to differences in the primordial large-scale environment of these haloes and other haloes of similar mass. To briefly explore this hypothesis, we show in Fig. 16 the difference in the dark matter density at z = 4 around the progenitor haloes of CGDs and control sample.

It is possible to see that at least for the present-day haloes with |$M_{\rm 200,m} \gtrsim 10^{11.1} \ \rm M_\odot$|⁠, there was a difference in the local dark matter density at z = 4 between CGD and control progenitor haloes, with the groups progenitors inhabiting regions slightly less dense (⁠|$\sim 15~{{\ \rm per\ cent}}$| lower than control). If, at a fixed cosmic time, less massive haloes have a higher probability of forming in less dense regions of the cosmic web, then the difference in the assembly of CGDs when compared to control haloes could be partially explained by a difference in the past large-scale environment.

4.3 Group coalescence and descendants

Dwarf galaxies inside the CGDs may suffer from cumulative environmental effects, so it is important to quantify the lifespan of these groups. The ages of most CGDs found at z = 0 in TNG50-1 indicate that these systems can survive for more than one billion years without the merging of their members (see Fig. 9). This means dwarf galaxies inhabiting CGDs may suffer dwarf-only environmental effects for a non-negligible amount of time. The typical coalescence time-scale of simulated CGDs that we found at z > 0 is even larger (>1 Gyr), implying that the dwarf galaxies may suffer even longer with tidal forces and gas shocks inside the groups. In their work, S17 extrapolates the merger time-scale of groups with massive galaxies (Barnes 1989) to predict the time until the observed CGDs would merge and suggests a time-scale of ≲1 Gyr. On the other hand, the CGDs that we found at z = 0 in the simulation exist as compact groups for time-scales longer than 1 Gyr. As is shown in Fig. 12, the typical time-scale for the coalescence of CGDs in the simulation is even longer, being |$\sim 3 \ \rm Gyr$| for groups found at z ≤ 0.5. This suggests that the extrapolation of results found for massive groups may not be adequate to predict properties of groups – such as their merger time-scales – in the low-mass regime.

More than half of the groups found at z = 0 already existed more than 1 Gyr ago, meaning that the galaxies in these groups may have been interacting for long enough to suffer the effects of the group environment. In the dwarf regime, interactions which do not end in a merger can be more important than mergers for the amount of stellar mass formed (Martin et al. 2021), so to understand how these interactions affect the star formation history (SFH) of dwarf galaxies that inhabit compact groups, we are also conducting a detailed study of the CGDs members SFH’s.

Considering the CGDs as an intermediate stage for the formation of more massive systems, it is relevant to characterize the properties of the galaxy that form from the complete coalescence of these groups. In Fig. 17, we present the position of the descendants of all coalesced groups at z ≤ 0.5 in the M_*–sSFR diagram. We also plot reference lines to represent the MS of star formation at z = 0 and z = 0.5. The MS at each redshift was obtained from the linear interpolation between the median sSFR in different stellar mass bins. The bins are defined within |$7 \le \log (M_\ast /\rm M_\odot) \le 10$|⁠, with ∼750 data points each. As we can see in the figure, the coalescence of CGDs gives origin to normally star-forming galaxies (SFGs) in the mass range |$8 \lesssim \log (M_\ast /\rm M_\odot) \lesssim 10$|⁠. These objects are a minority in the population of SFGs at z = 0; however, CGDs are expected to be more common at high-redshift, so they can possibly be an important evolutionary stage in the formation of low- to intermediate-mass galaxies at earlier epochs.

Furthermore, galaxies that recently formed from the coalescence of CGDs may differ from other galaxies of similar mass in terms of their gas reservoir. The interaction of the dwarf galaxies during group coalescence may have similar effects on the diffuse gas as in dwarf pairs, where stripping and scattering of the gas creates a large gas reservoir that can fuel future star formation in the absence of a more massive galaxy (Pearson et al. 2016). The fact that all merger remnants of CGDs in the simulation lie on the star-forming main sequence means that star formation can be maintained for a long time after the group coalescence. This is in line with results found for dwarf-dwarf mergers and their remnants in other studies (Amorisco & Evans 2012; Pearson et al. 2018; Kado-Fong et al. 2020), which point to a scenario where the merger remnants have a gas supply to sustain star formation continuously.

Regarding the structure of the galaxies formed from group coalescence, we found that discs of gas with active star formation are easily formed short (less than 2 Gyr) after the complete merging process. This finding is valid for a subset of galaxies with |$M_{\ast } \lt 10^{10} \ \rm M_\odot$| living in gas-rich haloes, but according to results from Sotillo-Ramos et al. (2022), this easiness to form discs may also be valid for higher masses since they find that Milky Way/Andromeda analogues in TNG50 can easily reform their stellar discs after a major merger given that there is sufficient gas available.

5 CONCLUSIONS

In this work, we search for analogues of CGDs in the IllustrisTNG cosmological simulations and compare their properties with some expectations from observations. We analyse the properties of analogue CGDs, investigating how these groups can form in a full cosmological context and what type of galaxies will result from their coalescence. We can summarize our conclusions as follows:

• The TNG50-1 simulation can reproduce isolated compact groups of dwarf galaxies in the Local Universe (z < 0.1), with the majority of these groups being triplets (Fig. 4) and having a median halo mass of |$10^{11.3} \ \rm M_\odot$| (Fig. 1).

• The comoving number density of CGDs in TNG50-1 at z = 0 is |$n \approx 10^{-3.56} \ \rm cMpc^{-3}$|⁠, and increases with redshift (Fig. 3), effectively doubling the number of CGDs by z = 0.26.

• The predicted number of CGDs observed in a survey like the SDSS Legacy Survey partially agrees with the lower limits imposed by the observational sample used in this work. However, the confrontation between simulation and observations is limited to a small range of magnitudes (⁠|$-17.5 \lesssim M_{r \rm ,3rd} \lesssim -16.5$|⁠) due to the SDSS sample detection limit and the TNG50-1 volume (Fig. 13).

• The total group masses (Fig. 5) and mass-to-light ratios (10 ≤ M/L ≤ 100; Fig. 15) of the simulated CGD are mostly compatible with observations, but the most massive observed groups are missing in the simulation, probably due to the limited volume of TNG50-1.

• Host haloes of CGDs at z = 0 have lower stellar masses, higher total gas fractions and were assembled later compared to other haloes of similar mass at the same epoch (Fig. 7).

• Simulated CGDs found at z = 0 have median formation times of ∼1.6 Gyr (Fig. 9), with ∼20 per cent of the groups being already formed for more than 3 Gyr.

• In the large-scale structure at z = 0, CGDs avoid regions of high matter density, such as the meeting point of filaments (knots), and regions with the lowest matter density, such as voids (Fig. 10).

• About 14 per cent of all unique CGDs identified since z ≤ 0.5 merged into a single galaxy, and their median coalescence time-scale is ∼3 Gyr (Fig. 9).

• All the CGDs that coalesced at z ≤ 0.5 formed star-forming galaxies with 8 < log (M_*/M_⊙) < 10 at z = 0 (Fig. 17). Part of these galaxies have intermediate stellar mass (9.5 ≤ log (M_*/M_⊙) ≤ 10) and represent 2.7 per cent of the galaxies within this mass regime in the latest snapshot of TNG50-1.

In summary, we showed that the IllustrisTNG highest resolution simulation can create CGDs with total masses comparable to observations. We also found that the simulated CGDs at z = 0 are formed from a late hierarchical assembly of haloes in the low-mass regime and that CGDs at z < 0.5 can commonly take more than 1 Gyr to fully merge in the simulations. Investigating the formation of galaxy groups in the low-mass regime is important to understand the environmental effect of dwarf groups on galaxies even before they infall in denser environments and also to obtain insights into the hierarchical assembly of low-mass haloes that may be common at the first stages of galaxy formation. In this manner, we hope that our work can help connect the models of galaxy formation with present and future observations in the still vastly unexplored low-luminosity regime.

ACKNOWLEDGEMENTS

The authors thank the referee for comments that led to an improved version of the manuscript. RFF thanks the support of Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq). MT thanks the support of CNPq (process #312541/2021-0). We acknowledge the use of SDSS data (http://www.sdss.org/collaboration/credits.html), TOPCAT Table/VOTable Processing Software (Taylor 2005, http://www.star.bris.ac.uk/mbt/topcat/), and R language and environment for statistical computing (R Core Team 2015). The IllustrisTNG simulations were undertaken with compute time awarded by the Gauss Centre for Supercomputing (GCS) under GCS Large-Scale Projects GCS-ILLU and GCS-DWAR on the GCS share of the supercomputer Hazel Hen at the High Performance Computing Center Stuttgart (HLRS), as well as on the machines of the Max Planck Computing and Data Facility (MPCDF) in Garching, Germany.

DATA AVAILABILITY

The IllustrisTNG simulations are publicly available at www.tng-project.org/data (Nelson et al. 2019). Nearest neighbour catalogues for IllustrisTNG are available at http://www.tng-project.org/floresfreitas24. Catalogues of simulated and observed CGDs analysed in this work are provided as supplementary material. Additional data directly related to this publication are available on request from the corresponding author.

Footnotes

References

APPENDIX A: CUMULATIVE DISTRIBUTIONS OF HALO PROPERTIES

In Section 3.2.2, we presented results to show the difference between the properties of haloes hosting CGDs and haloes on a control sample paired by mass. For completeness, in Fig. A1, we show the empirical cumulative distribution of the halo properties plotted in Fig. 7. In the first row of Fig. A1, we show the distributions of halo mass, and it is clear that both samples have similar distributions of mass even in the mass bins. As seen in the second, third, and fourth rows, the shape of the cumulative distributions is different for all mass bins. We quantified the difference between the distributions through KS-tests discussed in Section 3.2.2, and we found the same results using the Anderson–Darling k-sample test, which is more sensitive to the tails of distributions.

APPENDIX B: SUPPLEMENTARY MATERIAL AND NEAREST NEIGHBOUR CATALOGUE

We provide videos and catalogues as supplementary material for this work. The material includes videos showing the evolution of gas and stellar density of CGDs. Catalogues of simulated and observed CGDs analysed in this work will also be provided as machine-readable tables. Here, we also present details of the nearest neighbour catalogue, which will be provided as supplementary data for the IllustrisTNG simulations. The catalogue is constructed by measuring Euclidean distances between galaxies in the simulation box. The minimum stellar mass (⁠|$M_{\ast \rm ,min}^{\rm nn}$|⁠) for target galaxies and neighbours depends on the resolution of the simulation run, with |$M_{\ast \rm ,min}^{\rm nn} = 10^7 \ \rm M_\odot$| for TNG50-1, |$M_{\ast \rm ,min}^{\rm nn} = 10^{8.2} \ \rm M_\odot$| for TNG100-1, and |$M_{\ast \rm ,min}^{\rm nn} = 10^{9.1} \ \rm M_\odot$| for TNG300-1. For each galaxy, the catalogue stores information about its ten nearest neighbours within 15 Mpc and the total number of neighbours within a fixed spherical aperture centred on the target galaxy. Only subhaloes with subhalo_flag = 1 are used to build the catalogue.