Understanding the Formation and Evolution of Dark Galaxies in a Simulated Universe

The IllustrisTNG (hereafter TNG) project (Nelson et al., 2018; Naiman et al., 2018; Marinacci et al., 2018; Springel et al., 2018; Pillepich et al., 2018a) conducts a suite of large-volume, cosmological simulations performed with the moving-mesh code, AREPO (Springel, 2010; Weinberger et al., 2020). As the successor to the original Illustris project (Vogelsberger et al., 2014a, b; Genel et al., 2014; Sijacki et al., 2015), the TNG simulations, which stands for ‘The Next Generation’, aim to model the formation and evolution of galaxies within the framework of the $\Lambda$CDM cosmology. The new sub-resolution physics and numerical improvements of TNG are described in Weinberger et al. (2017) and Pillepich et al. (2018b).

The TNG simulations are based on the $\Lambda$CDM cosmological parameters from Planck Collaboration et al. (2016): $\Omega_{\Lambda,0}=0.6911$, $\Omega_{m,0}=0.3089$, $\Omega_{b,0}=0.0486$, $\sigma_{8}=0.8159$, $n_{s}=0.9667$, and $H_{0}=100\,h\,\text{km}\,\text{s}^{-1}\,\text{Mpc}^{-1}$ with $h=0.6774$. The simulations incorporate self-gravity and hydrodynamics, as well as cosmic magnetic fields. Given the complex nature of baryonic physics and the wide range of scales involved, TNG employs various sub-resolution models to account for processes such as gas cooling, star formation, metal enrichment, and feedback mechanisms from supernova (SN) and active galactic nuclei (AGN). They also include the effects of cosmic reionization by activating a spatially uniform UV background radiation at $z\sim 6$. This cosmic reionization will be discussed in detail in Section 4.3, which turns out to be important for the formation of dark galaxies. The TNG project provides three realizations of simulations: TNG50, TNG100, and TNG300, with each number representing the side length of the simulation box in units of cMpc.

In this work, we use TNG50 (Pillepich et al., 2019; Nelson et al., 2019), as it provides the highest resolution for analyzing individual galaxies. We obtain the snapshots containing full particle data and group catalogs containing information on galaxy groups/clusters and individual galaxies at several selected redshifts covering $0\leq z\leq 11$. We also use the galaxy merger trees constructed with SUBLINK (Rodriguez-Gomez et al., 2015) and LHaloTree (Springel et al., 2005). Because most galaxies are first identified after $z\sim 11$ on their merger trees, we trace them until $z=11$.

A group catalog at each snapshot includes two types of groups: FoF (friends-of-friends) and SUBFIND groups. FoF groups, corresponding to galaxy groups/clusters, are initially identified using a friends-of-friends algorithm (Davis et al., 1985). This algorithm connects particles into a single FoF group if their separation is less than the linking length, denoted as $b$ times the mean inter-particle distance. Here, $b$ represents the percolation parameter with a specific value of 0.2.

For each FOF group, the SUBFIND algorithm (Springel et al., 2001; Dolag et al., 2009) is employed to identify gravitationally self-bound substructures. Such structures are identified in overdense regions in the FoF group, and then gravitationally unbound particles are removed. We use these SUBFIND groups (hereafter galaxies) that initially form as centrals in their own FoF groups and evolve on unbroken branches of merger trees.

We note that all measurements in this work are obtained from the SUBFIND group catalog. Therefore, when we use the term ‘halo’, we are specifically referring to the self-bound structures at the scale of individual galaxies, rather than larger group or cluster-scale halos.

We classify galaxies into two groups, luminous or dark galaxies, using their stellar-to-total mass ratios at redshift $z=0$ (see Figure 1). Luminous galaxies are defined as those with stellar masses greater than or equal to 1/10,000 of their total masses (i.e., $M_{*}/M_{\text{tot}}\geq 0.0001$), while dark galaxies have stellar masses below this threshold (i.e., $M_{*}/M_{\text{tot}}<0.0001$). Our choice of this mass ratio threshold is motivated by observational studies (e.g., Fig. 4 of van Dokkum et al., 2018), as we aim to focus on the objects that are darker than those typically observed. Minor adjustments to this threshold are unlikely to significantly impact our results.

Dark galaxies are further divided into two groups: star-poor galaxies, which have a positive stellar mass, and starless galaxies, which have a stellar mass of zero. It should be clarified that we still call starless objects as galaxies, even though they do not have any visible stars. Similarly, we use the terms of ‘luminous galaxies’ or ‘star-rich galaxies’ when referring to the galaxies not classified as dark galaxies depending on the context: i.e. ‘luminous’ vs. ‘dark’, ‘star-rich’ vs. ‘star-poor’ and ‘starless’.

Figure 1 shows our classification of galaxies based on their stellar masses (x-axis) and total masses (y-axis). A white dashed line corresponds to the threshold where $M_{\star}/M_{\text{tot}}=0.0001$, which separates dark galaxies from luminous galaxies. In the TNG50 group catalog, there is a total of 5,686,764 galaxies. Of these, only 116,200 galaxies are classified as luminous galaxies, while the majority of galaxies are classified as dark galaxies in our definition. The dark galaxies are composed of 47,270 star-poor galaxies and 5,524,844 starless galaxies.

Figure 2 represents the mass distribution of the dark matter component for each sample. The total dark matter mass distribution shows an abundance of galaxies with lower-mass dark matter halo and fewer galaxies with higher-mass dark matter halo. This trend is consistent with the initial power spectrum predicted by the $\Lambda$CDM cosmology. Luminous galaxies are prevalent at the high-mass end, while dark galaxies, particularly starless galaxies, are more abundant at the lower-mass end. The large number of starless galaxies in the low-mass range would be partly because of simulation resolution coupled with the star formation recipe. For example, the star particles do not form below a mass resolution in simulations, which can contribute to the increase in the number of starless galaxies. Additionally, the specific choice of star formation recipe can further impact the modeling of star formation in low-mass galaxies (Federrath & Klessen, 2012; Nuñez-Castiñeyra et al., 2021).

To ensure a fair comparison between luminous and dark galaxies, we select only the galaxies with dark matter mass in the range of $1.0\times 10^{9}\,h^{-1}\text{M}_{\odot}\leq M_{\text{DM}}\leq\,3.0\times 10^{9}\,h^{-1}\text{M}_{\odot}$. The magenta region in Figure 2 corresponds to this dark matter mass selection criterion. We adopt this criterion for two reasons. First, at this mass range, the three types of galaxies (star-rich, star-poor, and starless galaxies) have the most overlap in frequency, indicating similar sample sizes. Second, the selection allows us to investigate variations in baryonic content among galaxies by fixing the dark matter mass of the galaxies at $z=0$.

To summarize, our final galaxy sample includes 14,206 star-rich galaxies, 19,245 star-poor galaxies, and 14,318 starless galaxies. When we explore other redshifts, we continue tracking these galaxies across time. Hence, in this work, the type of galaxy is always referenced to their state at $z=0$.

We first investigate the gas properties of both luminous and dark galaxies. Because gas is the fuel of star formation, this investigation may enable us to find hints as to why dark galaxies fail to form stars.

Using the merger tree of TNG50, we trace the evolutionary history of the key properties of galaxies, i.e., their mass, size, and density. Figure 5 shows the evolutionary histories of these properties. Here, we present the plots for dark matter and baryonic components (gas and star) separately. Each curve on the plot represents the median value with the shaded region indicating the 1-$\sigma$ scatter. We note that the curves and 1-$\sigma$ scatters are measured from the galaxies with positive stellar or gas mass.

To explore the spatial distributions of luminous and dark galaxies, we analyze their local densities ($\rho_{20}$) as well as their distances to the nearest filament ($D_{\text{skel}}$).

The local density indicates how crowded galaxy neighborhoods are. For a given galaxy, we calculate the local density by using surrounding dark matter particles as tracers. Because there are a huge number of dark matter particles in the full snapshot data ($\#=10,077,696,000$), we choose a smaller subset by random resampling. This subset contains 0.001% of the total dark matter particles. With this resampled dataset, we compute the local dark matter density using the following equation:

Here, $n$ is the number of neighboring dark matter particles, while $W(|\mathbb{r}_{i}|,h_{sml})$ is the kernel function to smooth out the distribution of dark matter particles. Additionally, $\mathbb{r}_{i}$ denotes the distance between a given galaxy and the $i$-th dark matter particle, and $h_{\text{sml}}$ represents the smoothing length, corresponding to the distance to a $n$-th dark matter particle. We choose $n=20$, because the distance to the 20th nearest dark matter particle corresponds to a few Mpc, which is suited for the consideration of large-scale structural properties (Park et al., 2007; Muldrew et al., 2012). Furthermore, we normalize the $\rho_{20}$ value by estimating the mean mass density of dark matter over the entire simulation box, denoted as $\bar{\rho}$.

We also compute the distance of each galaxy to its nearest filament. To extract the filamentary structures from the resampled distribution of dark matter particles, we use the Discrete Persistent Structure Extractor code (DisPerSE, Sousbie, 2011; Sousbie et al., 2011). We adopt a minimum persistence level of 6$\sigma$, ensuring that the selected filamentary structures are distinct from random fluctuations at a confidence level of 6$\sigma$. For each galaxy, we identify the closest filament and measure the distance to it, denoted as $D_{\text{skel}}$.

Figure 6 represents the spatial distributions of galaxies in terms of these two parameters: $\rho_{20}/\bar{\rho}$ (x-axis) and $D_{\text{skel}}$ (y-axis), which indicate isotropic and anisotropic environments, respectively (e.g., Park et al., 2007; Kraljic et al., 2018).

The left panel of Figure 6 shows the spatial distribution of star-rich, star-poor, and starless galaxies at the present epoch ($z=0$). Individual galaxies are denoted by dots, while contours illustrate their distributions. Star-rich, star-poor, and starless galaxies correspond to black, blue, and red contours, respectively. We realize that the center of black contours is located towards the lower right side, indicating that star-rich galaxies tend to be closer to filaments and located in denser regions. In contrast, the centers of red and blue contours are located towards the upper left side, indicating star-poor and starless galaxies tend to be far from filaments and located in less dense regions.

The right panel of Figure 6 mirrors the left panel but shows the distribution in the early universe at $z=11$, corresponding to when the universe is $\sim 0.4$ Gyr old. The contours appear closer together because large-scale structures had not yet fully developed at this point. This leads to a small difference in the distribution of the contours among the samples. Despite the close proximity of the contours, there exists an offset between lumionus and dark galaxies, which is confirmed by the Kolmogorov-Smirnov (K-S) test (see the p-values in each histogram panel). This finding implies that luminous and dark galaxies originate from slightly different locations, and these disparate origins could impact their future evolution, such as influencing gas accretion or mass growth rates.

The spin of dark matter halo can be influenced by the tidal torque from its surrounding environment, which in turn can have implications for gas accretion. Therefore, it could be an important component in the evolution of galaxies. Here, we examine the spin parameter proposed by Bullock et al. (2001), which is defined as

where $\mathbf{J}$ is the angular momentum of a galaxy halo, $M$ and $R$ are the virial mass and virial radius, respectively, and $V$ is the virial circular velocity, given by $V=\sqrt{GM/R}$ with the gravitational constant $G$.

The virial radius of a halo can be described as the point where the enclosed mass density meets the criteria in the following equation:

Here, $\rho_{c}$ is the redshift-dependent critical density of the universe, defined as

with the Hubble parameter $H^{2}(z)=h^{2}[\Omega_{m,0}(1+z)^{3}+\Omega_{\Lambda,0}]$. In the context of a flat universe (Bryan & Norman, 1998), $\Delta_{c}=18\pi^{2}+82x-39x^{2}$ where $x=\Omega(z)-1$, and $\Omega(z)$ is the matter content at redshift $z$ defined as

Once measuring the virial radius, we can easily obtain the virial mass as follows:

Figure 7 displays histograms of the spin parameters measured with dark matter particles at different redshifts. In the bottom-most panel, it is evident that dark galaxies have larger spin parameters than luminous galaxies at $z=0$. In the top-most panel (corresponding to the highest redshift $z=11$), the distributions of spin parameters between luminous and dark galaxies overlap significantly. However, the small p-values derived from the K-S test (shown in each panel) indicate that we can reject the null hypothesis suggesting that spin parameters of luminous and dark galaxies are drawn from the same underlying distribution. Therefore, our finding indicates that dark galaxies have slightly larger spin parameters than those of luminous galaxies.

We also find a decreasing trend in the spin parameters over time for both luminous and dark galaxies. Notably, the spin parameters of luminous galaxies diminish more significantly after cosmic noon ($z<2$). This decline in the spin of star-rich galaxies can be attributed to the transfer of angular momentum through mergers and interactions with other galaxies. A more detailed discussion about this phenomenon will be provided in Section 4.2.

The distinct properties of present-day luminous and dark galaxies originate from their early environments, such as their birthplace and the amount of available gas for star formation (See Figures 4 and 6). Our findings suggest that galaxies formed in less dense environments with a scarcity of star-forming gas are more likely to become dark galaxies. In contrast, galaxies formed in denser environments with sufficient star-forming gas are more prone to develop into luminous galaxies.

Moreover, we also find a subtle disparity in the initial spins of dark galaxies, which tend to be slightly larger than those of luminous galaxies even though the differences in spin are not significant (see Figure 7). This is consistent with the findings reported in the previous studies (Jimenez et al., 1997; Kim & Lee, 2013; Jimenez & Heavens, 2020).

As galaxies evolve, the differences between luminous and dark galaxies become more pronounced. We use the merger history data (Rodriguez-Gomez et al., 2017; Eisert et al., 2023), and find that luminous galaxies in denser regions have experienced a larger number of mergers than dark galaxies in sparser regions. Here, it should be noted that the merger history data is constructed for the cases when galaxies have non-zero stellar mass. Therefore, galaxies with no star particles (i.e., starless galaxies in this study) are not properly considered in this analysis.

Figure 8 shows the total number of mergers across the cosmic history. It highlights that star-rich galaxies have experienced numerous mergers with the systems containing stars, but star-poor galaxies have largely remained unaffected. As noted above, the plot shows the cases when the mergers occur between galaxies with stars. If there are galaxies that may have merged with starless systems, the data in this plot do not reflect such cases.

We also use the merger tree to investigate the satellite fraction ($f_{\text{sat}}$) as a function of time for luminous and dark galaxies. The satellite fraction represents the proportion of luminous or dark galaxies that are satellites of other, more massive galaxies. This is computed by dividing the number of satellite luminous or dark galaxies by the total number of galaxies.

Figure 9 shows that the satellite fraction of star-rich galaxies has gradually increased up to approximately 0.5 by $z=0$. This indicates that half of them have become members of galaxy clusters or groups via frequent merging and infalling to the large system. In contrast, the satellite fractions of star-poor and starless galaxies never exceed 0.1 during the entirety of cosmic history, implying that about 90% of them continue to be isolated. This result provides further evidence supporting the idea that the disparity between luminous and dark galaxies becomes more significant as time goes on.

When galaxies undergo mergers and interactions with other galaxies, they exchange their angular momentum. The process involves the transfer of orbital angular momentum from less massive galaxies to large central galaxies. As a result, the spin parameter of satellites may decrease, contributing to a reduction in their sizes. Moreover, tidal stripping, as galaxies fall into clusters, can remove both gas and dark matter components (Moore et al., 1996; Gnedin, 2003; Aguerri & González-García, 2009; Smith et al., 2016; Niemiec et al., 2019). This phenomenon further contributes to a decrease in their sizes. As evidenced in Figure 8, such events tend to be more common among luminous galaxies, resulting in the reduction of sizes for their gas and dark matter components; this is illustrated in Figure 5.

Our results emphasize the important role of cosmic reionization in the evolution of dark galaxies. By heating up gas components and expelling them from galaxies, cosmic reionization can boost the differentiation between luminous and dark galaxies (Thoul & Weinberg, 1996; Hambrick et al., 2011).

Galaxies that are significantly influenced by cosmic reionization undergo a reduction in their gas reservoirs, which are essential for star formation. Consequently, these galaxies struggle to form new stars and evolve into dark galaxies. On the other hand, galaxies that withstand the impacts of cosmic reionization and maintain their gas supply continue their star formation activities. These galaxies eventually evolve into luminous galaxies (See Figures 4 & 5).

One of the factors that can make a galaxy survive from cosmic reionization is its dark matter mass before this epoch. As shown in the dark matter mass evolution history in Figure 5, luminous galaxies clearly contain more dark matter masses than their dark galaxy counterparts just before $z=6$. This abundance of dark matter in luminous galaxies will provide a deeper gravitational potential, which can reduce the gas loss during the cosmic reionization phase.

The formation time of galaxies could be another factor that can make a galaxy survive from cosmic reionization. Figure 10 shows the time when galaxies are identified as SUBFIND groups; i.e., a group should have at least 20 dark matter and star particles. It reveals that the majority of galaxies are identified before cosmic reionization. However, star-rich galaxies are generally identified earlier than star-poor and starless galaxies. Therefore, dark galaxies could not have enough time to be stabilized, which can make them more vulnerable to the effects of cosmic reionization from outside. As previously mentioned in Section 3.2.1, when it comes to formation time, star-rich galaxies require 1.2 Gyr to reach half of their final dark matter mass, while star-poor and starless galaxies need 2.2 Gyr and 3.5 Gyr, respectively. This further supports the idea that luminous galaxies form earlier so that they could be stabilized to survive from cosmic reionization effect.

However, it is important to be aware of the limitations of simulations in modeling cosmic reionization. In the TNG simulation, cosmic reionization starts to operate at $z=6$ based on the assumption of a uniform UV background radiation field as a heating source. Only the gas cells with a density greater than 0.001 $\rm cm^{-3}$ are able to cool via self-shielding. Gas cells with a density lower than this threshold experience heating up to $\rm 10^{4}\,K$ without cooling from self-shielding. Therefore, in this simple prescription, the effect of ionization could be exaggerated in low-density environments where the density of ionizing photons is also low (Rosdahl et al., 2018). Additionally, it could suppress gas infall onto low-mass galaxies by photoheating the intergalactic gas and raising its pressure (Shapiro et al., 1994; Gnedin, 2000; Hoeft et al., 2006).

There may be other factors that we did not explore in detail, but could still be important for understanding the origin and fate of dark galaxies in the cosmic landscape. We briefly list them here and leave them for future studies.

Baryonic feedback: Baryonic feedback could have an impact on the formation of dark galaxies. For instance, the supernova (SN) feedback mechanism injects energy into gas, affecting both their kinematics and thermodynamics. This energy injection can lead to the suppression of star formation and potentially play a role in the formation of dark galaxies.

Ram pressure stripping: When a galaxy falls into the center of a galaxy cluster, it experiences strong ram pressure due to the surrounding gas (Gunn & Gott, 1972; Hwang et al., 2018). This pressure can remove the gas component of a galaxy, which in turn hinders star formation and may lead to the birth of dark galaxies. According to our findings, roughly 10% of star-poor galaxies are satellites within galaxy clusters or groups. We suggest that some dark galaxies in relatively dense regions might have formed through this process. These galaxies, which experienced ram pressure as satellites, could have different origins compared to dark galaxies in lower-density environments.

Dark matter halo concentration: Although we have examined the density of dark matter halos in Section 5, it is important to note that this density is calculated within the half-mass radius. To gain further insights into the characteristics of dark matter halos, it is essential to investigate their concentration, as they can vary even when the integrated densities appear similar.

Dark matter halo shape: The shape of the dark matter halo itself could be one of the other factors impacting the formation of dark galaxies. While luminous galaxies seem less perturbed by cosmic reionization than dark galaxies, the disparity might stem from differences in halo configurations. Although both types of galaxies have similar halo masses, dark galaxies might possess more triaxial halo shapes than luminous galaxies, making them more prone to losing their gas.