Observational Predictions for the Survival of Atomic Hydrogen in Simulated Fornax-like Galaxy Clusters

The presence of dense, neutral hydrogen clouds in the hot, diffuse intragroup and intracluster (IC) medium is an important clue to the physical processes controlling the survival of cold gas and sheds light on cosmological baryon flows in massive halos. Advances in numerical modeling and observational surveys mean that theory and observational comparisons are now possible. In this paper, we use the high-resolution TNG50 cosmological simulation to study the H i distribution in seven halos with masses similar to the Fornax galaxy cluster. Adopting observational sensitivities similar to the MeerKAT Fornax Survey (MFS), an ongoing H i survey that will probe to column densities of 1018 cm−2, we find that Fornax-like TNG50 halos have an extended distribution of neutral hydrogen clouds. Within 1 R vir, we predict the MFS will observe a total H i covering fraction of ∼12% (mean value) for 10 kpc pixels and 6% for 2 kpc pixels. If we restrict this to gas more than 10 half-mass radii from galaxies, the mean values only decrease mildly, to 10% (4%) for 10 (2) kpc pixels (albeit with significant halo-to-halo spread). Although there are large amounts of H i outside of galaxies, the gas seems to be associated with satellites, judging both by the visual inspection of projections and by comparison of the line of sight velocities of galaxies and IC H i.


INTRODUCTION
The formation and evolution of galaxies is now understood to be strongly linked with the diffuse gas filling their dark matter halos.Depending on galaxy mass, this gas is called the circumgalactic medium (CGM) when the gas resides in the halos of galaxies of the Milky Way mass or lower, or the intracluster medium (ICM) for gas living in more massive cluster-sized halos (M tot ≳ 10 14 M ⊙ ).Various physical processes, such as gas accretion from the intergalactic medium, and feedback driven by stars and AGN, occur in this gaseous halo.These processes drive flows which regulate the rate at which gas cools on to the galaxy itself, controlling the amount of mass in the interstellar medium and hence the rate of star formation itself, including possibly quenching (see avi.chaturvedi@aip.deTumlinson et al. 2017, for a review).Therefore, understanding the physical processes in the CGM or ICM is crucial to building a comprehensive picture of galaxy evolution.We highlight here the important question of the origin and survival of cold gas in the hot ICM, by studying the phase-space distribution of cold gas in massive halos at low redshift, with a particular focus on the observability of this gas in ongoing and future HI surveys.
For low-mass, star-forming galaxies, it is expected that their halo can host an abundant amount of cold gas that may be the fuel for future star formation.The halos of massive galaxies, on the other hand, may lack cold gas, resulting in the lack of any recent star formation activity (Gauthier & Chen 2011).In addition, current observational (Chen et al. 2018;Berg et al. 2019;Zahedy et al. 2019) and simulation studies (Davé et al. 2020;Rahmati et al. 2015;Nelson et al. 2020) suggest that at intermediate redshift (0.3 ≤ z ≤ 0.8), massive ha-los have a significant amount of cold gas.For example, Chen et al. (2018); Zahedy et al. (2019) studied luminous red galaxies (LRG) at redshift z ∼ 0.21-0.55 and found that these galaxies host a high column density of the cold gas tracers HI and MgII.Other similar studies point toward the same conclusion that LRGs host an abundant amount of cold gas (Zhu et al. 2014;Lan & Mo 2018;Anand et al. 2021Anand et al. , 2022)).At low-redshift, simulation studies of Milky Way-like galaxies show that their CGM also host substantial amounts of cold gas (van de Voort et al. 2019).
Recently, simulations have made improved predictions for the cold gas distribution in halos that reproduce many of the observed galactic properties.For example, Nelson et al. (2020) have shown that cold gas in the halos of LRGs can be attributed to the thermal instability triggered by local density perturbations.They suggest that these perturbations are related to gas stripped from infalling galaxies via tidal interactions, or ram pressure stripping.Performing a comparative study between cosmological and idealized simulations (individual galaxy halo simulations), Fielding et al. (2020) have also shown that cold gas extends to the virial radius for Milky Waymass halos.They also suggest that non-spherical accretion and satellite galaxies contribute to the cold gas content in the outer halos.Previously, Villaescusa-Navarro et al. (2018) performed a detailed study using the Illustris TNG100 simulation, investigating the HI abundance and clustering properties in halos for z ≤ 5.They showed that HI density profiles are sensitive to various processes such as AGN feedback and tidal stripping.For massive halos, they found that HI is mostly concentrated in their satellite galaxies, whereas for small halos it is concentrated in the central galaxy.
In contrast to intermediate redshifts, the study of cold gas in the ICM of massive halos (M * ≥ 10 11 M ⊙ ) in the local Universe (z ∼ 0) is limited to a few studies.Nonetheless, using the HI 21 cm emission line, radio observations have demonstrated the abundant existence of cold neutral atomic gas around early-type galaxies (E/S0) galaxies (Serra et al. 2012(Serra et al. , 2013b;;Young et al. 2014;Serra et al. 2013a).However, these observations are limited to tens of kpc around the targeted galaxies and are typically not sensitive to HI column densities below ≤ 10 19 cm −2 , and therefore do not provide a comprehensive picture of the cold gas in the ICM.
The ongoing observations from the MeerKAT Fornax survey (hereafter MFS, Serra et al. 2019;Serra et al. 2023), a radio continuum and line survey of the Fornax cluster, provides an excellent opportunity to study the HI gas in great detail in the nearby Fornax ICM (d ∼ 20 Mpc).Photometric and spectroscopic studies (Cantiello et al. 2020;Chaturvedi et al. 2022) have shown that the Fornax cluster mass assembly is still ongoing, making it an interesting target to study.MFS is dedicated to studying the HI distribution and kinematics within the Fornax environment.The HI column density sensitivity of the MFS ranges from ∼5×10 19 cm −2 at a spatial resolution of ∼10 arcsec (∼1 kpc at Fornax distance) down to 10 18 cm −2 at ∼100 arcsec (∼10 kpc at Fornax distance).With a mosaic area of the 12 square degrees, MFS will detect HI in the Fornax intracluster (hereafter IC) region -which in this paper refers to the region within the massive dark matter halo and outside satellite galaxies.
The high-resolution TNG50 cosmological simulations (Nelson et al. 2019;Pillepich et al. 2019) provide a median spatial resolution of ∼100 parsec and its validation of cold gas (neutral and molecular hydrogen) against observational work (Popping et al. 2019;Diemer et al. 2019) makes TNG50 an ideal framework to explore the cold gas distribution in Fornax-like halos.This also provides a chance to forecast the upcoming MFS survey results and test the simulations against the MFS observations.In this work, we use the TNG50 simulations (Nelson et al. 2019;Pillepich et al. 2019) and adopt the observing criteria of the MFS to study the HI content in the TNG50 halos similar to the Fornax galaxy cluster.We also study the HI distribution in these halos and their IC region.We calculate the HI covering fraction for these halos and predict the expected observed MFS HI covering fraction.
In addition, both the spatial and velocity distribution of HI gas in the ICM of clusters and groups can be used to gain insight into the origin and survival of this cold gas.If the gas is correlated in both position and velocity with satellites, we can argue that the HI is likely either stripped from satellites or cooling is induced by satellites.However, if cold gas is not correlated with satellite galaxy positions or velocities, we might argue that either cold gas formation is related to the central galaxy or that cold gas survives in the ICM long enough to become virialized (e.g., Voit et al. 2017;Rohr et al. 2023).In this paper we take the first step of making these spatial and velocity maps, leaving gas particle tracking to future work.
The paper is organized as follows: In section 2, we briefly introduce the TNG50 simulation and present the methodology for calculating the HI covering fraction.Section 3 presents our results about the HI distribution and its covering fraction.In section 4, we compare our results to current simulation (Section 4.1) and observational (Section 4.2) studies, and discuss the likely origin of the cold gas in the intracluster medium of Fornax-like halos (Section 4.3).Section 5 presents the summary of the work.

SIMULATION AND METHODOLOGY
This section briefly introduces the TNG50 simulation that we use for our analysis as well as our criteria for selecting Fornax-like halos.In addition, we present our methodology for calculating the HI column density and HI covering fraction of the selected halos in TNG50 in order to compare to observational surveys.

The TNG simulations
For our study, we use the TNG50 simulation (Nelson et al. 2019;Pillepich et al. 2019), the highest resolution simulation of the IllustrisTNG cosmological magnetohydrodynamical (MHD) simulation suite (Nelson et al. 2018;Springel et al. 2018;Pillepich et al. 2018;Marinacci et al. 2018;Naiman et al. 2018).The IllustrisTNG project is a set of large cosmological simulations that include a variety of galaxy formation physics including AGN feedback.The model has been designed to match a wide range of observational constraints (Pillepich et al. 2018;Springel et al. 2018) and was carried out with the moving mesh code AREPO (Springel 2010).The AREPO code solves the coupled evolution of dark matter, gas, stars, and black holes under the influence of selfgravity and ideal MHD.Developed with the key motivation to study galaxy formation physics and understand the growth of cosmic structure physics, the IllustrisTNG project uses three distinct simulation box sizes.
TNG50 was carried out with a box size of 51.7 Mpc per side with 2160 3 gas and dark matter cells, resulting in a baryon mass resolution of 8.4 × 10 4 M ⊙ .In particular, we used TNG50-1 (hereafter TNG50), the highest resolution of the three variants run, which provides a median spatial resolution of ∼100 pc.We analyze this run, although the other larger boxes (TNG100 and TNG300) contain a larger number of Fornax cluster-sized objects, because we need high spatial resolution to study the interaction of cold high density columns and the hot group medium.TNG50 adopts initial conditions and cosmological parameters consistent with the Planck Collaboration et al. (2016) cosmology with h = 0.68, Ω b = 0.05, Ω m = 0.31, Ω λ = 0.69, and σ 8 = 0.82 and assuming a flat universe governed by a Λ cold dark matter (ΛCDM) cosmology.

Fornax-like Halo selection
In TNG50 a galaxy cluster and groups of galaxies are referred to as halo or FOF (hereafter referred as halo), identified through the friends-of-friends algorithm (Davis et al. 1985).Within each halo, the SUBFIND algorithm (Springel et al. 2001) identifies the subhalos including the primary (central) galaxy and other satellite galaxies (hereafter referred as satellite).To find halos similar to the Fornax cluster in TNG50 at snapshot 99 (redshift z = 0), we applied a virial mass selection criterion analogous to the Fornax cluster mass (M 200 ∼ 5 × 10 13 M ⊙ , adopted from Drinkwater et al. 2001), namely the mass range of 10 13.5 ≤ M 200 ≤ 10 14 M ⊙ , where M 200 is defined as the mass enclosed within a virial radius R 200 equal to 200 times the critical density of the Universe.With this condition, we find a total of seven halos.
For these halos, we measured the stellar velocity dispersion of their central galaxy and found that this value is quite close to that of NGC1399, the central galaxy of the Fornax cluster.Except for halo IDs 4 and 9, all other halo central galaxy stellar velocity dispersions fall within ± 20 % to the stellar velocity dispersion value of 315 kms of NGC1399 (Vaughan 2019).In addition to this, the velocity dispersion of all the subhalos within these halos agrees to within 15% of the observed Fornax cluster members (giants and dwarf galaxies) mean velocity dispersion value of 374 ±25 km/s (Drinkwater et al. 2001).A previous study of the hot X-ray emitting medium of galaxy groups in the TNG50 cosmological simulations (Truong et al. 2020) showed a good match to observations.More recently, a new set of zoom-in cosmological simulations using the same model 'TNGcluster' (Truong et al. 2023) have demonstrated good agreement on a larger sample.These successes in reproducing observable properties of galaxy clusters indicates that our seven halos are reasonable matches to the Fornax cluster.From here onward, we refer to these halos as Fornax-like halos.In Table 1, we list the physical properties of these halos.

Atomic HI content
To determine the HI mass of gas cells in TNG50 Fornax-like halos, we use the Popping et al. (2019) molecular hydrogen fraction (H 2 ) catalogue, previously calculated for the TNG simulations.In this work we use their fiducial recipe, which is based on the work by Gnedin & Kravtsov (2011).Gnedin & Kravtsov (2011) performed detailed simulations including non-equilibrium chemistry and simplified 3D on-the-fly radiative transfer calculations.Based on these simulations, the authors presented fitting formulae for the H 2 fraction of neutral gas as a function of the dust-to-gas ratio of the gas, the impinging UV radiation field, and surface density of the neutral gas.Popping et al. (2019) assume that the dust-to-gas ratio scales with the metallicity of the neutral gas, that the local UV radiation field scales with the SFR of the gas cell with an additional contribution from the ionising UV background field and they calculate the gas surface density of a gas cell by multiplying its density by the Jeans length of the cell.A detailed description of the implementation of the Gnedin & Kravtsov (2011) fitting formulae within the TNG simulation suite is presented in Popping et al. 2019, (see their Section 2).

Halo HI covering fraction
To understand and quantify the HI distribution in the halos, we measure their HI covering fraction in different column densities bins (hereafter denoted as N HI ), adopting a range of N HI > 10 18 , 10 19 , and 10 20 cm −2 .We first measured the HI column density, by performing a two-dimensional binning of all HI gas cells along the three projected axis X, Y and Z (spatial position), regardless of the velocity space and adopted a pixel size of 2 kpc, similar to the MFS spatial resolution limit.This assumes that gas particle sizes are smaller than 2 kpc, which is generally true for the TNG50 halos gas cells.We checked that larger gas cells have negligible contribution to the HI mass and hence to HI column density.In Section 3.2, we present the HI covering fraction results as a function of velocity.
We measure the HI covering fraction in two ways, similar to Rahmati et al. (2015), as follows: Cumulative HI covering fraction, hereafter denoted as f HI (R), is defined as the fraction of surface area covered by the binned pixels having column density higher than a given N HI value within a radius R divided by the total area of pixels within radius R. f HI (R) is expressed as: Here A N HI is the single pixel area with a column density equal or higher than the given N HI value and is summed over N such pixels in a given area of radius R divided by the total area of pixels in radius R.
Differential HI covering fraction, hereafter denoted as f HI (∆R), is defined similarly to be cumulative HI covering fraction, except here we consider only the covering fraction within the radial bin defined between radius R j and R j+1 and f HI (∆R) is expressed as: We also separately measure the total HI covering fraction and the HI covering fraction in the intracluster (IC) regions of the halos, that is, the regions well away from identified galaxies.For the IC measurement, we first select the satellite galaxies within a halo with a stellar mass ≥ 10 8.5 M ⊙ and having at least 10 gas cells (varying this to 1 gas particle has no effect on our results).Then, we remove the HI gas cells associated with these satellite (including the central galaxy) using the SUB-FIND algorithm of TNG50.This removal of HI gas is done out to ten times the stellar half mass radius (denoted as R 1/2 * ) of a satellite.In TNG50, the SUBFIND algorithm identifies all the gas cells that are gravitationally bound to a specific satellite.After removing the gas cells of a satellite identified with the SUBFIND procedure out to 10×R 1/2 * , we measure the HI covering fraction in the same way as defined earlier.The R 1/2 * radii of satellite vary from a few kpc to tens of kpc, and with our adopted radial limit (10×R 1/2 * ), we make sure that we exclude all gas cells that are within the domain of an individual galaxy.However, the IC measurement may contain the very extended tidal/stripped gas tails originating from the individual galaxies.This section presents our results showing the HI distribution and covering fraction in the TNG50 halos and their IC regions.For our study, we considered gas cells gravitationally bound to a halo, including gas cells extending to an average 1.5 R vir radii of halos.Figure 1 shows the HI distribution in TNG50 halo 6 (left panel) and in its IC region (right panel).Three pink circles indicate the virial radius drawn at 0.5, 1.0 and 1.5 R vir .These images demonstrates both the patchy nature of HI in the simulated clusters, as well as it's distributed nature.The IC image (right) emphasizes that much of the HI (at least by area) is not immediately connected to galaxies -that is, it is at least 10 stellar half-mass radii from any galaxy in the simulation.
Figure 2 shows the HI distributions for all other Fornax-like halos (top two rows) and in their IC regions (bottom two rows).The first visual impression we get from these plots is that the large-scale distribution of HI extends beyond 0.5 virial radii (∼350 kpc) for all the halos.Otherwise, these halos demonstrate diverse HI distributions with significant variations in the amount of HI mass (Table 1).The diverse and extended HI distribution in these halos could potentially be related to the merger/accretion history of these halos or the ac-tivity of supermassive black holes (Zinger et al. 2020).We can also see the streams or filamentary structures connecting the central and satellite galaxies.We notice, in addition, a large number of small HI regions with relatively high column densities, which are possibly not related to any satellite.We refer to these as clouds.We caution that, unlike Nelson et al. (2019), we have not performed Voronoi tessellation over the gas cells to identify these clouds-like structures and it merely represents HI clumps around the satellite galaxies and in the intra-cluster regions of these halos.For the halos HI map, we see that the centers of the satellite are dominated by HI column densities N HI between 10 20 to 10 21 cm −2 .Looking at these maps, it is quite clear that HI clouds extend out to the virial radius (corresponding to ∼ 700 kpc) covering the IC regions.Within the inner region of each halo, around ∼ 0.25 R vir , a large fraction of the HI gas cells are associated with the central galaxy of the halo.In the IC region, the observed HI structures primarily have column densities N HI < 10 20 cm −2 , whereas the structures beyond 0.5 virial radii lie mainly at 10 18 < N HI /cm −2 < 10 19 .We present the results of the HI covering fraction of halos and in the IC regions in subsections 3.1 and 3.2, respectively.
Log NHI cm -2    3, but for the HI distribution in the intracluster (IC) region of Fornax-like halos HI maps (Figure 2, bottom two rows) measured along the arbitrarily chosen z-axis for a pixel size of 2 kpc.We obtain the IC HI by removing the gas cells gravitationally bound to all galaxies within 10× their stellar half mass radii.Note the y-axis ranges differ from Figure 3.

HI covering fraction profiles
We used the HI projected maps as shown in Figure 1 to measure the HI covering fraction as discussed in section 2. We measured the HI covering fractions of halos for three projections, along the X, Y, and Z directions, using a pixel scale of 2 kpc.
In this section, we present the HI covering fraction for the full halo maps (top two rows in Figure 2) and in the next section, we discuss the covering fraction of IC regions.In Figure 3 we show the cumulative and differential covering fraction profiles measured along the (arbitrarily chosen) z-axis in the left and right panel, respectively.For these panels, we include all of the HI gas in the covering fraction calculation, whether or not it is within the central or a satellite galaxy.We do this because, in a blind HI survey, the satellite galaxies may not be identified, so the HI would be measured globally.By definition, the innermost point of the cumulative and differential covering fractions are the same, then the differential covering fraction begins decreasing more steeply than the cumulative covering fraction.The first, middle, and bottom rows in both panels show the covering fraction for HI column density for N HI ≥ 10 18 , 10 19 and 10 20 cm −2 respectively.The thin lines indi-cate the individual halos, and the thick lines mark the average value.
We focus mainly on the N HI bins of 10 18 and 10 19 cm −2 , which are the optimal range for studying the HI distribution at ∼kpc scales for the MFS.We find that regardless of the projection axis, the average covering fraction for the N HI ≥ 10 18 cm −2 bin remains between 10-15% at 0.5 R vir and drops to 6-10% at 1 R vir .With increasing column density, the covering fraction decreases, such that for the N HI ≥ 10 19 cm −2 bin the covering fraction drops to 5-10% at 0.5 R vir and to less then 5% at 1 R vir .The differential covering fraction at 0.5 virial radius is between 5-10% and drops to less than 5% at 1 R vir .
These covering fractions quantify our visual impressions.Although all the halos have some HI gas within 0.5 R vir , it is distributed non-uniformly in small structures that look like filaments or clouds.In Figure 3 we verify that the covering fraction of these structures is low, even when including column densities down to 10 18 cm −2 .These structures and clouds could potentially be associated with the central galaxy and satellite galaxies, but from Figure 2 we see that the HI is clearly extended well beyond the satellite stellar radii.However, this gas might have been stripped from satellites to form part of the IC region, which we discuss in the next section and in Section 4.3.

Intra-cluster HI covering fraction
The bottom two rows of Figure 2 show the HI distribution in the IC regions of the halos.The majority of pixels with column density N HI ≥ 10 18 cm −2 extend into the IC regions (the cumulative HI covering fraction drops by about 30% at 1 R vir ).It is primarily pixels with high column density HI, N HI > 10 20 cm −2 , which are removed when only including the IC gas cells in our column density measurement (the cumulative HI covering fraction drops by about 70% at 1 R vir ).
Similar to Figure 3, we show the HI covering profiles for the IC regions in Figure 4.For the N HI bins 10 18 and 10 19 cm −2 , the cumulative covering fraction at 0.5 R vir is between 5-10 % and drops to less than 5% at 1 R vir .
We checked the differences in the HI covering fraction along the different projected directions and find that, on average, the HI covering fraction remains the same and varies only by a few percent when changing the projected axis.Figure 5 shows the cumulative HI covering fraction at 0.5 R vir (upper panel) and 1 R vir (lower panel) between the Fornax-like halo and their IC region in three different projected axes for a column density N HI ≥ 10 18 cm −2 .The open blue (green) star, circle and square in- dicate the covering fractions at 0.5 R vir (1 R vir ) along the X, Y and Z axis projections, respectively.The open black (red) star, circle and square indicate the average values measured at 0.5 R vir (1 R vir ).On average, the IC HI covering fraction is between 70-80% of the total HI covering fraction.This suggests that a large fraction, around ∼ 75% (by covering fraction) of the low-column density HI gas that is distributed throughout these massive halos is well outside the satellite galaxies.We also note that some pixels that had a high column density when including the satellite gas have lower column densities when only including IC gas.This can be seen as some red pixels in the top panels of Figure 2 turning into green pixels in the bottom panels.
Depending on the HI column sensitivity, MFS will map the HI distribution in Fornax at different spatial resolutions, varying from 1 to 10 kpc.We created similar HI maps, as shown in Figure 2, with different pixel  In Figure 6, we show the average halo and IC cumulative HI covering fraction of the 7 halos at 1 R vir radius for different pixel scales.As anticipated, when degrading the pixel size of HI maps, the HI covering fraction increases and, for N HI ≥ 10 18 cm −2 , reaches 12% at 10 kpc pixel scale for the full halo and is around 10% for the IC halo map.These are then our resolution-matched predicted values for the N HI ≥ 10 18 cm −2 covering fraction that will be observed by MFS.
Finally, in order to investigate whether we can see any signature of diffuse HI gas in the inter-galactic medium or cosmic filamentary structures around these mock Fornax-like halos, we made HI maps extending out to 3 R vir , which we show in Appendix B. Visually inspecting these maps, we do not find diffuse HI in the inter-galactic medium of halos.We notice that outside 1.5 R vir , there are several infalling smaller satellite galaxies, particularly for Halos 1, 2 and 4, but they contribute a negligible amount in the HI covering fraction profile.In particular, we computed HI covering fractions for these maps and found similar results as for the maps shown in Figure 2.

HI covering fraction in velocity space
We also measured the cumulative HI covering fraction in the line-of-sight velocity space, adopting a velocity bin size of ∼100 kms −1 within a range of -700 to 700 kms −1 .This allows us to connect the HI gas at different column densities in velocity space.A correlation of HI gas in velocity space with the satellite velocity distribution can give us hints about its possible stripping origin from satellite galaxies.Figure 7 shows the HI distribution in velocity space for Halo 6 (left panel) and in its IC region (right panel).Similar to Figure 1, we created these velocity space maps (Sec.2.4) by performing twodimensional binning in the phase space and taking the mean velocities in bins corresponding to the X, Y and Z axis, with a pixel scale of 2 kpc.In the velocity maps, we consider the halo velocity as the zero velocity point.Remaining HI velocity maps of other halos and their IC regions are shown in Figure 8.In these maps, we can see the distribution of larger satellite galaxies having different velocities then the surrounding lower density HI (first two rows in Figure 8).We used these maps to measure the HI covering fraction in velocity space for these halos which is shown in Figure 9.The top, middle, and bottom panels show cumulative HI covering fraction for N HI column densities of 10 18 , 10 19 and 10 20 cm −2 with a velocity bin size of 100 kms −1 .The thin lines indicate the individual halos, and the thick continuous and dashed lines mark the average value of halo and IC region respectively.Vertical dashed gray line marks the halo velocity which is taken as the zero velocity point.
We find that, for a velocity bin size of 100 kms −1 , on average the cumulative HI covering fraction is less than 4 % (2 %) for the column density N HI ≥ 10 18 (10 19 ) cm −2 and significantly drops to less than 1% for N HI ≥ 10 20 cm −2 .The halo and IC HI velocity covering fraction looks bimodal for most of the halos (except halo 4 and 9), as also shown by the average HI covering fraction (continuous and dashed blue lines in Figure 10).On a halo-by-halo comparison, we have verified that the IC HI does not show a more Gaussian velocity distribution than the total HI, as would be expected by virialized, ballistic gas clouds.
We can also use our HI velocity covering fraction analysis to learn about the origin of the HI gas in the IC region.To do this, we looked for a possible correlation between the HI velocity covering fraction and the satellite velocity distribution.In Figure 10, we show the velocity distribution of all satellite galaxies (all identified satellites with at least 100 cells and stellar mass ≥ 10 8.5 M ⊙ ), as a function of the number of galaxies (orange leftmost y-axis of each panel) and their stellar mass (black left y-axis of each panel).The grey bars in Figure 10 indicate the summed stellar mass of galaxies binned within velocity ranges of 100 kms −1 .The orange histogram represents the number counts of all satellite galaxies in each velocity bin.Additionally, we overlay the velocity covering fraction (right y-axis of each panel) for log N HI ≥ 18 cm −2 in blue (total and IC HI in solid and dashed lines, respectively).For most of the halos, we find that the velocity covering fraction for both the halo and IC regions follows the velocity distribution of satellite galaxies.This points towards a scenario where the IC HI gas originated from satellite galaxies.We performed a similar analysis considering only the HI-rich galaxies (satellite having HI gas cells outside 10 R 1/2 * ) and found similar correlation between the stellar mass distribution and HI velocity covering fraction, except that the stellar mass distribution is less centrally peaked (as we would expect for gas rich galaxies that have likely recently entered the cluster).
We show both the number distribution and the stellar mass distribution of satellites because both could effect the total HI mass brought into the cluster.The stellar mass distribution of satellite galaxies may even be more likely to predict where HI will be found as, for example, a single 10 10 M ⊙ galaxy is likely to fall into the cluster with more HI than two 10 9 M ⊙ galaxies.In Figure 10, we note that besides the agreement with the number distribution of satellites, the HI velocity covering fraction also correlates with peaks in the velocity distribution of the stellar masses of all satellite galaxies (grey histogram).Indeed, in Halos 3 and 7, where the satellite stellar mass distribution shows multiple peaks, the HI velocity distribution seems to follow the stellar mass more closely than the number of satellites.This general agreement between the satellite and IC HI velocity distributions suggests that IC HI originates from satellite galaxies, and the possibly stronger relation between the satellite stellar mass and HI velocity gives a hint that the IC HI may even originate from more massive galaxies (although we stress that this does not necessarily imply the gas is stripped from massive disks but may fall in as part of the galaxy's circumgalactic gas).
Although a detailed investigation is outside the scope of this paper, we note that it is expected that some processes occurring in clusters affect the gas dynamics, and therefore some offset between the HI and satellite velocity distributions is unsurprising.For example, a past merger can cause long-lived motions in the ICM (Vaezzadeh et al. 2022), or black hole activity could affect gas dynamics (Weinberger et al. 2017).

DISCUSSION
In this section, we discuss our results and compare the HI covering fraction to available observations and other cosmological simulation studies.

Comparison to other cosmological simulations
We begin by comparing our measured cumulative HI covering fraction with the available HI covering fraction  from the studies of Nelson et al. (2020).They measured the abundance of cold gas in TNG50 halos for massive halos with mass > 10 11 M ⊙ at intermediate redshift z ∼ 0.5.Although we have only 7 TNG50 halos, and apart from the evolution of the halos from redshift z∼0.5 to 0, the simulations and HI model we used are the same as Nelson et al. (2020); therefore, the HI covering fraction should be of the same order.Our measured HI covering fraction for N HI agrees well with the Nelson et al. (2020) measured values.For a column density of N HI > 10 17 cm −2 , we find a covering fraction around 70 ± 15 % at 10 kpc, dropping to 30 ± 15% at 100 kpc, and for N HI > 10 20 cm −2 , at 100 kpc, the covering fraction is roughly 10%, similar to the findings of Nelson et al. (2020).Rahmati et al. (2015) used the EAGLE simulation to study the HI distribution around high-redshift massive galaxies.They found a strong evolutionary trend in the HI covering fraction within the virial radius with redshift.For an averaged HI column density in between 10 17.3 < N HI /cm −2 < 10 21 , the HI covering fraction drops from 70 % at z = 4 to 10 % at z =1.The HI content of galaxies in the EAGLE cosmological simulations was also investigated (Marasco et al. 2016;Crain et al. 2017), finding that the highest resolution simulations reproduced the HI masses of galaxies as well as their clustering.In addition, in dense group and cluster environments, they found that ram pressure stripping was the primary HI mass removal process but that galaxy interactions also played a role.
Studying the spatial distribution and ionisation state of cold gas in the CGM, Faerman & Werk (2023) performed semi-analytical modelling of cold gas in the CGM of low-redshift star-forming galaxies.Assuming that cold clouds in the CGM are in local pressure equilibrium with the warm/hot phase, they reported that cold gas can be found out to 0.6 R vir or beyond.Although we examine more massive halos (10 13.5 ≤ M 200 ≤ 10 14 M ⊙ ) compared to Faerman & Werk (2023), we also find that the CGM of Fornax-like halos normally shows a spatially extended distribution of cold gas clouds out to more than 0.5 R vir .
Previously, van de Voort et al. ( 2019) have shown that standard mass refinement and a high spatial resolution of a few kpc scale can significantly change the inferred HI column density.Studying zoom-in simulations of a Milky-way mass galaxy within the virial radius, they found that the HI covering fraction of N HI ≤ 10 19 cm −2 at 150 kpc is almost doubled from 18% to 30% when increasing the spatial resolution of the CGM.Although the simulation setup of van de Voort et al. ( 2019) and TNG50 is different, we compare our findings to theirs based on the similar resolutions (∼1 kpc).We find that the cumulative HI covering fraction within 150 kpc is around 25%, quite close to the van de Voort et al. ( 2019) result.
Figure 10.HI covering fraction profiles of Fornax-like halos in velocity space for log NHI ≥ 18 cm −2 gas over-plotted on the number counts of all satellite galaxies (orange histogram) and the stellar mass velocity distribution (grey histogram).The profiles and histograms are measured in 100 km/s bins.The leftmost y-axis (orange) of each subplot represents the number count of satellite galaxies, while the second axis (black) represents the log stellar mass of satellite galaxies.Solid and dashed blue lines indicate the velocity covering fraction of the halo and IC regions, respectively (right y-axis of each panel).
In addition to the van de Voort et al. ( 2019) work, several papers have focused on studying CGM properties using higher spatial and mass resolution simulations.For example, the FOGGIE group (Peeples et al. 2019) has used the cosmological code Enzo Bryan et al. (2014) to carry out a set of simulations with high resolution in the CGM, finding that a great deal of small-scale structure emerged in the multiphase gas, producing many more small clouds.However, they found that the HI covering fractions did not change significantly.Similarly, Hummels et al. (2019) studied the simulated CGM with an enhanced halo resolution technique in the TEM-PEST simulation, again based on ENZO cosmological zoom simulations.They found that increasing the spatial resolution resulted in increasing cool gas content in the CGM.With an enhanced spatial resolution in the CGM, they found that observed HI content and column density increases in the CGM.In a similar vein Suresh et al. (2019) explored the CGM of a 10 12 M ⊙ mass halo with a super-Lagrangian zoom in method, reaching up to ∼ 95 pc resolution in CGM.They reported that enhanced resolution results in an increased amount of cold gas in the CGM and this increase in the cold gas also results in a small increase in the HI covering fraction.For a column density of N HI ≥ 10 19 cm −2 , they found a HI covering fraction value of around 10 % at one virial radius, which agrees well with our measured covering fraction of 12 % at one virial radius.More recently, Ramesh & Nelson (2024) re-simulated a sample of Milky Way galaxies at z ∼ 0 from TNG50 simulation with the super Lagrangian refinement method.Going down to a scale of 75 pc, they also reported that the abundance of cold gas clouds increases with enhanced resolution but did not find a large change in the covering fraction.

Detection of HI clouds in the Intracluster (IC) region through observational work
For all the TNG50 Fornax-like halos, we found that the HI clouds can be found beyond 0.5 R vir , which corresponds to an average physical scale of around 350 kpc.
Although only a few, there are observational studies showing the existence of remote HI clouds associated with galaxies.An important example is the detection of HI clouds in the inter-galactic medium of the galaxy group HCG44, where the HI clouds extend to more than ∼ 300 kpc (Serra et al. 2013b).Another example is the case of NGC 4532, in the Virgo cluster, where the HI tail of the galaxy with some discrete clouds extends to 500 kpc and constitutes around 10% of the total HI mass (Koopmann 2007).
A number of observational properties of the multiphase nature of the CGM in groups and clusters have been explored through absorption lines of background quasars with the HST-Cosmic Origins Spectrograph (COS).Studying a sample of low redshift luminous red galaxies (LRG) with metal absorption lines such MgII, CIII, and SiIII, Werk et al. (2013) measured the occurrence of cool metal enriched CGM and reported a MgII ion covering fraction (down to very low column densities) of 0.5 within 160 kpc radius.Using a sample of 16 LRG at z ∼ 0.4 observed with the HST/COS, (Chen et al. 2018;Zahedy et al. 2019) found a high HI covering fraction for column density N HI > 10 17.2 cm −2 of about 0.44 within 160 kpc impact parameter.In comparison to this, we measured an HI covering fraction of around 0.75 for column density N HI > 10 17.2 cm −2 at 160 kpc.Studying the CGM of a sample of 21 massive galaxies at z ∼ 0.5, Berg et al. (2019) measured an HI covering fraction of column density for N HI > 10 17.2 cm −2 within the virial radius of 15 %, which closely agrees with our average measured HI covering fraction of 12 %.Finally, Emerick et al. (2015) compared HI covering fractions within the virial radius in Virgo-like clusters between simulations and observations, finding values consistent with those found here.
Most recently, using the MeerKAT observations of the Fornax A subgroup, Kleiner et al. (2021) reported the detection of HI clouds at ∼220 kpc from NGC 1316, the central galaxy of Fornax A. Another study done with the MeerkAT telescope have detected a large extended HI cloud, extending ∼ 400 kpc in proximity to a large galaxy group at a redshift of z ∼ 0.03 (Józsa et al. 2022).Although observational detections of HI clouds in the IC regions around massive galaxies are few and rare, our and other simulation work like Rahmati et al. (2015); van de Voort et al. (2019); Nelson et al. (2020) strongly suggest the existence of dense small HI clumps within the ICM.For our work we find that within the IC region, HI tends to have a column density log N HI ∼ 19 cm −2 or less, and current observations are mostly not that sensitive yet.A strong test of cosmological simulations will be to compare our predicted HI covering fractions with MFS observations.If simulations overpredict the covering fraction of cold gas, some combination of i) excess cold gas removal from satellites, ii) excess cold gas added to the ICM from feedback or filamentary accretion, and iii) suppressed heating of cold gas in the ICM are likely at play.On the other hand, if simulations underpredict the HI covering fraction, some combinations of these effects are leading to too little cold gas.

Possible Origin Scenario
A detailed study of the origin or production of the large amount of HI gas in TNG50 Fornax-like halos is beyond the scope of this work.However, we speculate here on the possible source of the HI gas we do find.Nelson et al. (2020) studied the cold gas distribution in TNG50 massive halos at intermediate redshift z ∼ 0.5.Using Lagrangian tracer analysis, they argued that cold gas in TNG50 halos is related to gas that is removed from the halos in infalling satellites.These gas clouds can later stimulate the cooling process leading to a significant amount of cold gas.
Most recently, using the TNG50 simulation Ramesh & Nelson (2024) studied the cold gas clouds in the CGM of Milky-Way like galaxies.They reported that these high density gas clouds show clustering behaviour and this over-density increases around satellite galaxies, suggesting that a fraction of these clouds originate from ram-pressure stripping.This suggests that, not only for massive halos (like our Fornax-like halos), but even for smaller galaxies, the satellite gas stripping scenario can be important in producing HI clouds in the CGM.We however do not discount the possibility that a fraction of these clouds, whether around the satellite galaxies or isolated, may originate from the in-situ condensation of hot halo CGM gas in the cluster environment, or may be associated with outflows from the central galaxy (Fraternali & Binney 2006).Studying the formation mechanism of high velocity clouds around the Milky-Way like disk galaxies, Binney et al. (2009);Fraternali et al. (2015) have suggested that condensation of hot CGM gas can produce the HI clouds.It will be interesting to learn which is the effective mechanism for producing these HI clouds in the CGM.Possibly all the mechanisms: a) satellite stripping, b) thermal instability in CGM, or c) the feedback from the galaxy play roles in the formation of these clouds and could be explored by characterising these cloud properties in phase space (spatial and velocity) such as was done in Ramesh & Nelson (2024).
In addition to the analysis in Section 3.3, which indicates that the IC velocity covering fraction is associated with the satellite galaxies' velocity distribution, we checked if we could see any correlation between the satellite galaxies' number with the total IC HI mass. Figure 11 shows the halo intracluster HI mass as function of satellite galaxy number.In Figure 11, red stars denote all of the satellite galaxies and open squares mark the satellites having HI gas cells outside their 10 × R 1/2 * .With only 7 halos, its hard to quantify any relation, but we find that the halo IC HI mass increases with increasing number of satellite galaxies.The IC HI mass correlates more steeply with satellites having HI and shows less scatter.This further suggests that the HI mass in the IC regions of Fornax-like clusters could be associated with the stripped HI gas from satellite galaxies, similar to the findings of Nelson et al. (2020).We find a similar level of correlation between the total stellar mass in satellites and the halo IC HI mass.

SUMMARY AND CONCLUSION
In this paper, using the publicly available TNG50 simulation data, we have studied the distribution of HI gas in halos similar to the Fornax galaxy cluster.Adopting the MeerKAT Fornax survey (MFS) observational conditions, we have measured the HI covering fraction of the halos with a mass of 10 13.5 <= M 200 <= 10 14 M ⊙ .The following points summarise our findings and conclusions: 1. Atomic hydrogen in TNG50 Fornax-like halos shows a wide spatial distribution, appearing as clouds and filamentary structures (Figures 1 and  2).HI is non-uniformly distributed and extends in patches well beyond 0.5 virial radii of the central galaxy.On a physical scale, this corresponds to ∼ 350 kpc.
2. Using our HI covering fraction measurements, we find that individual Fornax-like halos in TNG50 show a wide scatter in the measured HI covering fraction ranging from 3% to 15% at 1 R vir (Figure 5).We predict the upcoming MFS should observe a total HI covering fraction of ∼ 25% at 0.5 virial radii and ∼ 12% at 1 R vir (Figure 6) at N HI ≥ 10 18 cm −2 (spatial resolution ∼ 10 kpc).For intracluster regions, this values drops to ∼ 20% at 0.5 virial radii and ∼ 9% at 1 R vir .
3. Intracluster (IC) regions (i.e. more than 10 stellar half-mass radii from identified galaxies) in Fornaxlike halos hold a substantial fraction of the HI.When using the N HI ≥ 10 18 cm −2 contour, the IC HI covering fraction at 1 R vir (spatial resolution ∼ 10 kpc) corresponds to around 75% of the total HI covering fraction (Figure 5).
4. The HI velocity covering fraction for the Fornaxlike halos (both in total and in the IC regions only) shows a broad velocity distribution that is not generally Gaussian, indicating that HI is not virialized in the halos (Figure 9).The HI velocity covering fraction for both halo and IC largely follows the velocity distribution of satellite galaxies, suggesting that IC HI is associated with satellite galaxies (Figure 10).
5. We find that halo HI intracluster mass increases with increasing number of satellite galaxies and shows an even stronger correlation with the satellites having HI presence in their outskirts (Figure 11).This also suggests that HI in the IC regions is associated with the stripped gas of satellite galaxies, similar to the results of Nelson et al. (2020).
With this work, we have demonstrated, based on TNG50 simulation data, that HI cold gas is predicted to co-exist and survive in the hot intracluster medium for Fornax-like clusters.Based on HI maps of Fornax-likes halos in TNG50, we expect MFS to find extended HI well beyond the satellites in the halo, but should generally follow the large-scale satellite distribution, both on the sky and in velocity space.This is also reflected in the asymmetry of the HI in velocity space -while the average distribution of the seven halos studied is symmetric, any individual halo can show strong asymmetries.We plan to perform a future follow-up study to pinpoint the origin of these HI clouds, whether they are possibly stripped or formed in situ in the cluster environment.It will be illuminating to see what MFS will observe within the Fornax cluster.With its higher sensitivity, there remains a good chance that the MeerKAT telescope can provide observational support and constraints for current and future simulation work on the multiphase nature of halo gas.In our paper we use pixels that are 2 kpc on a side in order to match the MFS resolution at HI column densities of about 10 19 cm −2 .However, the resolution at HI column densities of 10 18 cm −2 is about 10 kpc.In Figure 6 we show how the covering fraction at this column density increases with pixel size, and in this Appendix Figure 12, we show HI maps made with 10 kpc pixels so the reader can directly compare with the higher resolution maps shown in Figures 2  and 1.
B. TNG50 HALOS MAPS AT OUT TO 3 R VIR For our study we considered only the gas cells gravitationally bound to a halo (Sec.2), which generally includes cells out to an average of 1 to 1.5 R vir .Here in this Appendix Figure 13, we show TNG50 Fornax-like halos maps out to 3 R vir , demonstrating that there are several infalling satellite galaxies, particularly for Halos 1, 2 and 4. Their contribution to the HI covering profile are minimal.These maps also confirm that we do not detect any signature of diffuse HI (not closely associated with satellite galaxies) in the outer IGM of halos).

Figure 1 .
Figure 1.HI distribution in the TNG50 halo 6 at redshift z=0, colour-coded with the HI column density.The maps are made using a pixel scale of 2 kpc and are shown projected along the (arbitrarily chosen) z-axis.The three pink circles indicate the viral radius of the halos marked at 0.5, 1.0 and 1.5 times Rvir.Dark blue color in the maps indicate the HI column density lower than NHI= 10 16 cm −2 .Left: Full halo HI distribution, right: HI distribution in the intracluster (IC) region (i.e. after removing HI within 10 stellar half-mass radii of all galaxies).

<Figure 2 .
Figure 2. HI maps in the Fornax-like halos (first two rows) and in their IC regions (bottom two rows) made with a pixel size of 2 kpc (projected along the arbitrarily chosen z-axis).The three pink circles indicate the viral radius of the halos marked at 0.5, 1.0 and 1.5 times Rvir.

Figure 3 .
Figure 3. HI cumulative (left panel) and differential (right panel) covering fraction profiles of Fornax-like halos HI maps (Figure 2, first two rows) measured along the arbitrarily chosen z-axis for a pixel size of 2 kpc.The first, middle, and bottom rows show the covering fraction for HI column densities log NHI ≥ 18, 19 and 20 cm −2 respectively.The thin lines indicate the individual halos, and the thick lines mark the average value.The vertical dashed lines indicate the average virial radii of the halos.

Figure 4 .
Figure4.The same as Figure3, but for the HI distribution in the intracluster (IC) region of Fornax-like halos HI maps (Figure2, bottom two rows) measured along the arbitrarily chosen z-axis for a pixel size of 2 kpc.We obtain the IC HI by removing the gas cells gravitationally bound to all galaxies within 10× their stellar half mass radii.Note the y-axis ranges differ from Figure3.

Figure 5 .
Figure 5. Halo and intra-cluster cumulative HI covering fractions at 0.5 Rvir (upper panel) and 1 Rvir (lower panel) for three projected axes.Open star, circle and square indicate the X, Y and Z projections, respectively.The open black (red) star, circle and square indicate the average values.

Figure 6 .
Figure 6.Average halo and intra-cluster cumulative HI covering fractions of seven halos at 1 Rvir for three projected axes corresponding to different pixel scale size of HI maps.Numbers in the plot mark the pixel size in kpc used in creating the HI maps.
sizes ranging from 2 to 10 kpc.Examples of HI maps made with 10 kpc pixel scale are shown in Appendix A.

Figure 7 .
Figure 7. HI distribution in the TNG50 halo 6 at redshift z=0, colour-coded with the HI mean velocity (projected along the arbitrarily chosen Z axis) for a pixel size of 2 kpc.Left: Full halo HI distribution, right: HI distribution in the IC region.

Figure 8 .
Figure 8. HI distribution in the TNG50 halos (top two rows) and in their IC (bottom two rows) at redshift z=0, colour-coded with the HI velocity (projected along the arbitrarily chosen Z axis) for a pixel size of 2 kpc.

Figure 9 .
Figure 9. HI cumulative covering fraction profiles of Fornax-like halos in the velocity space.The top, middle, and bottom panels show the covering fraction for HI column densities log NHI ≥ 18, 19 and 20 cm −2 , respectively.The thin lines indicate the individual halos, and the thick lines mark the average value.The thick dashed lines mark the average IC HI cumulative covering fraction in velocity space.

Figure 11 .
Figure 11.HI intra-cluster (IC) mass as a function of number of satellites for our seven Fornax-like halos.
. would like to thank Abhijeet Anand and Roland Sazacks for helpful suggestions and discussions regarding this work.A.C. acknowledges the financial and computing support from Simons Foundation and thanks the Flatiron Institute pre-doctoral fellowship program through which this research work was carried out.GLB acknowledges support from the NSF (AST-2108470, AC-CESS), a NASA TCAN award, and the Simons Foundation through the Learning the Universe Collaboration.APPENDIX A. TNG50 HALOS MAPS AT 10 KPC PIXEL SIZE

Log NHI cm - 2 <Figure 12 .<Figure 13 .
Figure 12.Example of HI maps made with 10 kpc pixel scale, showing HI distribution in the TNG50 halo 3, 6 and 7 at redshift z=0, colour-coded with the HI column density.These maps are shown projected along the (arbitrarily chosen) z-axis.The three pink circles indicate the viral radius of the halos marked at 0.5, 1.0 and 1.5 times Rvir.Dark blue color in the maps indicate the HI column density lower than NHI= 10 16 cm −2 .Left: Full halo HI distribution, right: HI distribution in the intracluster (IC) region.