Unveiling the atomic hydrogen-halo mass relation via spectral stacking
Garima Chauhan, Claudia del P. Lagos, Adam R. H. Stevens, Matias Bravo, Jonghwan Rhee, Chris Power, Danail Obreschkow, Martin Meyer
MMNRAS , 1–21 (2020) Preprint 25 February 2021 Compiled using MNRAS L A TEX style file v3.0
Unveiling the atomic hydrogen–halo mass relation via spectral stacking
Garima Chauhan, , ★ Claudia del P. Lagos, , Adam R. H. Stevens, , Matías Bravo, Jonghwan Rhee, , Chris Power, , Danail Obreschkow, , Martin Meyer , International Centre for Radio Astronomy Research, The University of Western Australia, 35 Stirling Highway, Crawley, WA 6009, Australia ARC Centre of Excellence for All Sky Astrophysics in 3 Dimensions (ASTRO 3D)
Accepted XXX. Received YYY; in original form ZZZ
ABSTRACT
Measuring the Hi–halo mass scaling relation (HIHM) is fundamental to understanding the role of Hi in galaxy formationand its connection to structure formation. While direct measurements of the Hi mass in haloes are possible using Hi-spectralstacking, the reported shape of the relation depends on the techniques used to measure it (e.g. monotonically increasing withmass versus flat, mass-independent). Using a simulated Hi and optical survey produced with the Shark semi-analytic galaxyformation model, we investigate how well different observational techniques can recover the intrinsic, theoretically predicted,HIHM relation. We run a galaxy group finder and mimic the Hi stacking procedure adopted by different surveys and find wecan reproduce their observationally derived HIHM relation. However, none of the adopted techniques recover the underlyingHIHM relation predicted by the simulation. We find that systematic effects in halo mass estimates of galaxy groups modifythe inferred shape of the HIHM relation from the intrinsic one in the simulation, while contamination by interloping galaxies,not associated with the groups, contribute to the inferred Hi mass of a halo mass bin, when using large velocity windows forstacking. The effect of contamination is maximal at 𝑀 vir ∼ − . M (cid:12) . Stacking methods based on summing the Hi emissionspectra to infer the mean Hi mass of galaxies of different properties belonging to a group suffer minimal contamination but arestrongly limited by the use of optical counterparts, which miss the contribution of dwarf galaxies. Deep spectroscopic surveyswill provide significant improvements by going deeper while maintaining high spectroscopic completeness; for example, theWAVES survey will recover ∼
52% of the total Hi mass of the groups with 𝑀 vir ∼ M (cid:12) compared to ∼
21% in GAMA.
Key words: software: simulations – galaxies: haloes – galaxies: groups: general – radio lines: galaxies – galaxies: ISM
Since the establishment of the Λ cold dark matter model (hereafter Λ CDM) as the standard cosmological model, our understanding ofgalaxy formation and evolution is shaped by the idea that galaxiesform, evolve, and merge within host dark matter (DM) haloes. It iswell known that gravitational instability led to the growth of over-densities in the primordial matter density field, further leading gas todecouple from DM and dissipate, falling to the gravitational potentialcentre, and then cooling and forming stars, which eventually growinto galaxies (see Baugh 2006; Benson 2010). Thus, the growth, in-ternal properties, and spatial distribution of galaxies and DM haloesare closely linked (see review by Wechsler & Tinker 2018). This linkbetween galaxy properties and DM haloes has led to the establish-ment of some fundamental scaling relations, which characterise thedependence of the abundance on baryons to the host halo mass. Thestellar–halo mass relation (SMHM) is the most well-studied of theserelations (e.g. Behroozi et al. 2010; Moster et al. 2010). Our focus inthis paper is instead on the less-explored neutral atomic hydrogen–halo mass scaling relation (HIHM). Neutral atomic hydrogen (Hi)plays a fundamental role in our understanding of galaxy formationand evolution. Not only is it the raw fuel for galaxies, it also is an ★ Contact e-mail: [email protected] excellent tracer of galactic interactions and the underlying DM dis-tribution (e.g. Gunn & Gott 1972; Chung et al. 2009; Dénes et al.2016).Unlike the SMHM relation, the HIHM relation is much more com-plex, which we can see from the correlation between Hi mass andstellar mass in galaxies, which is characterised by a large scatter (e.g.Catinella et al. 2010; Brown et al. 2015, 2017), and thus, understand-ing the HIHM complexity is of particular interest. The form of theHIHM relation has been explored in various empirical studies, usingfor example, Hi clustering measurements (Padmanabhan & Refregier2017; Obuljen et al. 2019), which infer the HIHM relation by deter-mining the Hi associated with a halo by matching it with haloes thathave the same clustering strength. This requires a fair amount ofmodelling, as the shape of the clustering of DM haloes of a givenmass can significantly deviate from that of the Hi clustering. To solvethis, additional halo properties have been invoked, e.g. assembly age(see Guo et al. 2017).A more direct way of measuring the mean
Hi content of groupsis by employing Hi spectral stacking, whereby one co-adds the Hidetections and non-detections associated with an optically identi-fied group in an attempt to boost the signal-to-noise ratio, therebyrecovering a more statistically significant detection. This does notneed prior assumptions about relationships between halo propertiesand Hi as required when using Hi clustering. Recently, Guo et al. © a r X i v : . [ a s t r o - ph . GA ] F e b G. Chauhan et al. (2020) calculated the Hi content of groups identified in an opticalredshift survey by using an estimate of the halo radius to choosethe sky-projected area within which Hi emission was stacked. Guoet al. (2020) used SDSS (The Sloan Digital Sky Survey; York et al.2000) as their spectroscopic survey and ALFALFA (Arecibo LegacyFast ALFA; Giovanelli et al. 2005; Haynes et al. 2018) for their Histacking measurements. They first divided their galaxy groups (aswere presented in the Lim et al. 2017) according to the number ofmember galaxies, and then derived the HIHM relation for each sub-set. They found a monotonically increasing HIHM relation, whichstarts to plateau at 𝑀 vir ≥ . M (cid:12) . At those high halo masses,Guo et al. (2020) found that all selections in the number of membergalaxies produce the same relation (due to haloes of these massesalways having a large number of member galaxies). Guo et al. (2020)also found that as they select haloes to have a high occupation ofgalaxies (i.e. 𝑁 g ≥ 𝑁 g ≥ (cid:38) M (cid:12) ).Another way of estimating the Hi content of a group, employedby Rhee et al. (in prep), is through stacking the Hi spectra of indi-vidual galaxies identified as belonging to that group. Rhee et al. (inprep) used the GAMA (Galaxy And Mass Assembly; Driver et al.2011; Liske et al. 2015) survey as their spectroscopic survey and, theASKAP DINGO (Deep Investigation of Neutral Gas Origin; Meyer2009) early science observations for their Hi stacking. In this method,the spectra of galaxies are divided in host halo mass bins (based onthe group catalogue presented in Robotham et al. 2011), and thenstacked to measure the mean Hi mass of central/satellite galaxieswithin that halo mass bin. To then measure a total Hi mass for thehaloes in a bin, the mean Hi mass is simply multiplied by the numberof central/satellite galaxies belonging to groups in that halo mass bin,followed by summing up the total Hi mass of central/satellite galaxiesin a halo mass bin. Unlike Guo et al. (2020), Rhee et al. (in prep) finda non-monotonic HIHM relation, which dips at 𝑀 vir ∼ . M (cid:12) ,followed by an increasing Hi mass from 𝑀 vir (cid:38) . M (cid:12) .These are innovative new ways to further our understanding ofthe HIHM relation and provide valuable constraints for it. Although,novel in approach, these two techniques yield very different mea-surements of the HIHM relation. For example, there is a differenceof ∼ . 𝑀 vir ∼ . M (cid:12) between Guo et al. (2020) and Rhee et al. (in prep). It has so farbeen unclear whether this is simply due to the difference in theirmethodologies. We explicitly address this in this paper.The HIHM relation is of particular interest to theorists as well,and has been extensively investigated using various theoretical mod-els, ranging from semi-analytic models of galaxy formation (SAMs;e.g. Baugh et al. 2019; Spinelli et al. 2019; Chauhan et al. 2020) tohydrodynamical simulations (Villaescusa-Navarro et al. 2018). Sim-ulations have the power to investigate not only the average HIHMrelation, but also the scatter associated with it. Although the shapeand scatter of the HIHM relation is model-dependent (Chauhan et al.2020), it has been reported independently that the shape and scatterof the HIHM relation is heavily influenced by feedback from Ac-tive Galactic Nuclei (AGN) (see Baugh et al. 2019; Spinelli et al.2019; Chauhan et al. 2020). Spinelli et al. (2019) find that the HIHMrelation depends on the detailed assembly history of haloes, whichagrees with inferences based on Hi clustering studies in Guo et al.(2017). Chauhan et al. (2020) find that the scatter around the HIHMrelation can be constrained instead using halo spin and the ratio ofsubhalo mass to host halo mass.It is evident that progress has been made in our understandingof the HIHM relation on both observational and theoretical fronts, independently. On both of these fronts, some challenges need to beovercome. For example, Hi stacking is a very powerful tool for mak-ing Hi measurements of groups, but it heavily relies on group findersand halo mass estimates based on optical redshift surveys, whichmay sometimes provide inaccurate results due to survey limitations(Campbell et al. 2015; Bravo et al. 2020). Interestingly, simulationsand models, which do not suffer from observational limitations, showa range of different HIHM relations (see figure 2 in Chauhan et al.2020), which are significantly different from each other and dependon the models employed in the simulations. Here, we try to bridgethe gap between the observed and predicted HIHM relations, by pro-ducing a mock survey made using the state-of-the-art semi-analyticmodel of galaxy formation, Shark (Lagos et al. 2018; Chauhan et al.2019), and deriving mock-observed HIHM relations that mimic thetechniques of Guo et al. (2020) and Rhee et al. (in prep). Our mainaim is to see if we can reproduce the observed HIHM relation usingour mock survey and try to understand the limitations and advantagesof the stacking techniques used in the observational HIHM literature.As we will show, the HIHM relation that we derive from the mocksurvey is notably different from Shark’s intrinsic HIHM relation.We explore the origin of this difference, and test whether there is an“optimal-” observing method that can minimise it.The structure of this paper is as follows. Section 2 details theconstruction of the mock survey and mock group catalogue used forour analysis. Section 3 compares the halo mass estimates made forgroups by different methods against the intrinsic halo mass of groups.In Section 4, we delve into the systematic effects involved in Hi stack-ing measurements and compare different Hi stacking techniques withthe intrinsic prediction of Shark. In Section 5, we discuss the causesof the discrepancy seen between the intrinsic and the mock-observedHIHM relation. We draw our conclusions in Section 6. In order to make a fair comparison with the Hi stacking done ob-servationally, it is imperative to replicate stacking procedures on oursimulated “mock surveys". For this, it is necessary to build a simu-lated lightcone that replicates the limitations of the optical surveysthat serve as a base for the Hi stacking experiments. In the followingsubsections, we describe how these are built.
We use the semi-analytical model of galaxy formation ‘Shark’ (La-gos et al. 2018) to provide us with the simulated galaxies for ourHi stacking experiment. SAMs, such as Shark, use halo mergertrees, which are produced from a cosmological DM-only 𝑁 -bodysimulation, to follow the formation and evolution of galaxies by solv-ing a set of equations that describe exchange of mass, metals andangular momentum produced by a series of physical processes. Inthe following sections, we describe the DM-only 𝑁 -body simula-tion (surfs) that Shark is run on, followed by a description of thebaryon physics included in Shark. We end the section with a briefdescription of how the spectral energy distribution (SED) of galaxiesis computed in Shark. We use the surfs suite of DM-only 𝑁 -body simulations (Elahi et al.2018), which consists of 𝑁 -body simulations of differing volumes,from 40 to 210 ℎ − cMpc (comoving megaparsec) on a side, and MNRAS , 1–21 (2020)
I-Stacking: Systematics at play particle numbers ranging from ∼
130 million up to ∼ . Λ CDM Planck Collaboration et al. (2016)cosmology, which assumes total matter, baryon and dark energy den-sities of Ω 𝑚 = 0 . Ω 𝑏 = 0 . Ω Λ = 0 . ℎ = 0 . ℎ − cMpc, 1536 DM particleswith a mass of 2 . × ℎ − M (cid:12) and a softening length of 4 . ℎ − ckpc (comoving kiloparsec). surfs contains 200 snapshots for eachsimulation, with a typical time span between snapshots in the range6 −
80 Myr.Merger trees and halo catalogues were constructed using the phase-space finder VELOCIraptor (Elahi et al. 2019a; Cañas et al. 2019)and the halo merger tree code treefrog (Poulton et al. 2018; Elahiet al. 2019b), developed to work with VELOCIraptor. It has beenshown in Poulton et al. (2018) that treefrog + VELOCIraptor leadto well behaved merger trees, with orbits that are well reconstructed.We refer the readers to Lagos et al. (2018) for more details on how themerger trees and halo catalogues are constructed for Shark, and toElahi et al. (2018, 2019a); Cañas et al. (2019); Poulton et al. (2018)for more details on the VELOCIraptor and treefrog software. Shark is an open-source, flexible and highly modular SAM thatmodels the key physical processes of galaxy formation and evo-lution. These include (i) the collapse and merging of DM haloes;(ii) the accretion of gas onto haloes, which is governed by the DMaccretion rate; (iii) the shock heating and radiative cooling of gasinside DM haloes, leading to the formation of galactic discs via con-servation of specific angular momentum of the cooling gas; (iv) theformation of a multi-phase interstellar medium and subsequent starformation (SF) in galaxy discs; (v) the suppression of gas coolingdue to photo-ionisation; (vi) chemical enrichment of stars and gas;(vii) stellar feedback from evolving stellar populations; (viii) thegrowth of supermassive black holes (SMBH) via gas accretion andmerging with other SMBHs; (ix) heating by AGN; (x) galaxy mergersdriven by dynamical friction within common DM haloes, which cantrigger bursts of SF and the formation and/or growth of spheroids;and (xi) the collapse of globally unstable discs leading to bursts ofSF and the creation and/or growth of bulges.Shark also includes several different prescriptions for gas cooling,AGN feedback, stellar and photo-ionisation feedback, and SF. Sharkadopts a universal Chabrier (2003) IMF (Initial Mass Function).Shark uses these models to compute the exchange of mass, metals,and angular momentum between the key baryonic reservoirs in haloesand galaxies, which include hot and cold halo gas, the galactic stellarand gas discs and bulges, central black holes, as well as the ejectedgas component that tracks the baryons that have been expelled fromhaloes.The models and parameters used in this study are the Shark de-faults, as described in Lagos et al. (2018) and used in Chauhan et al.(2019); Chauhan et al. (2020) to study the Hi content of galaxies andhaloes. These have been calibrated to reproduce the 𝑧 = 0 , , and 2stellar mass functions; the 𝑧 = 0 black hole–bulge mass relation; andthe disc and bulge mass–size relations. This model also successfullyreproduces a range of observational results that are independent of https://github.com/icrar/VELOCIraptor-STF/ https://github.com/pelahi/TreeFrog https://github.com/ICRAR/shark those used in the calibration process. These include the total atomicand molecular hydrogen–stellar mass scaling relations at 𝑧 = 0; thecosmic star formation rate (SFR) density evolution up to 𝑧 ≈
4; thecosmic density evolution of the atomic and molecular hydrogen at 𝑧 ≤ The electromagnetic spectrum produced by the integrated contribu-tion of gas, dust and stars in galaxies, provides information aboutthe formation and evolution of the observed galaxies. This integratedelectromagnetic spectrum, also referred to as spectral energy distri-bution (SED), encompasses information of a galaxy’s stellar popu-lation, its interstellar medium (ISM) and the dust distribution (seeConroy 2013, for a review). As previously mentioned, for Hi stack-ing, the groups identified by optical surveys are matched with theirHi counterparts, to get an estimate of the Hi content of the groupcontained in a halo. In order to have a fair comparison betweensimulations and observations, we need to compute luminosities ofall Shark galaxies, which can then be used to apply optical surveylimits and group finding for these mock surveys.The process of computing the SEDs for galaxies in Shark isdetailed in Lagos et al. (2019), but here we briefly describe how isit done. In Shark, we have access to the star formation histories(SFH) and metallicity histories (ZH) for all the galaxies producedat every simulation snapshot. To produce the SEDs, we use twopackages: ProSpect and Viperfish . ProSpect (Robotham et al.2020) combines the GALEXev stellar synthesis libraries (Bruzual & https://github.com/asgr/ProSpect https://github.com/asgr/Viperfish MNRAS000
4; thecosmic density evolution of the atomic and molecular hydrogen at 𝑧 ≤ The electromagnetic spectrum produced by the integrated contribu-tion of gas, dust and stars in galaxies, provides information aboutthe formation and evolution of the observed galaxies. This integratedelectromagnetic spectrum, also referred to as spectral energy distri-bution (SED), encompasses information of a galaxy’s stellar popu-lation, its interstellar medium (ISM) and the dust distribution (seeConroy 2013, for a review). As previously mentioned, for Hi stack-ing, the groups identified by optical surveys are matched with theirHi counterparts, to get an estimate of the Hi content of the groupcontained in a halo. In order to have a fair comparison betweensimulations and observations, we need to compute luminosities ofall Shark galaxies, which can then be used to apply optical surveylimits and group finding for these mock surveys.The process of computing the SEDs for galaxies in Shark isdetailed in Lagos et al. (2019), but here we briefly describe how isit done. In Shark, we have access to the star formation histories(SFH) and metallicity histories (ZH) for all the galaxies producedat every simulation snapshot. To produce the SEDs, we use twopackages: ProSpect and Viperfish . ProSpect (Robotham et al.2020) combines the GALEXev stellar synthesis libraries (Bruzual & https://github.com/asgr/ProSpect https://github.com/asgr/Viperfish MNRAS000 , 1–21 (2020)
G. Chauhan et al.
Table 1.
Mock-survey parameters
Lightcone parameter Value
Area coverage 6900 deg Redshift range 0 − . 𝑟 -band magnitude ≤ 𝑀 ★ ≥ M (cid:12) Charlot 2003) and/or EMILES (Vazdekis et al. 2016) with the two-component dust attenuation model of Charlot & Fall (2000) and dustre-emission using the templates of Dale et al. (2014), which coversa rest-frame wavelength of up to 1 , 𝜇𝑚 . Viperfish is a wrapperfor ProSpect, which enables the use of Shark SFHs and ZHs togenerate a series of desired SEDs through target filters.An important novelty of the SED modelling in Shark is the waythe Charlot & Fall (2000) parameters are informed from galaxy prop-erties. Using a 3D radiative transfer analysis of galaxies in the EA-GLE hydrodynamical simulations (Schaye et al. 2015), Trayford et al.(2020) found that the resulting attenuation curves could be well fittedwith a parametrization ála Charlot & Fall (2000), in which the pa-rameters scale with the dust surface density and are independent ofredshift. Lagos et al. (2020) adopted this parametrisation, and henceevery galaxy has its own set of Charlot & Fall attenuation parametersbased on dust surface density. When combined with ProSpect andViperfish, Lagos et al. (2019) showed that Shark can reproducethe panchromatic emission of galaxies throughout cosmic time; mostnotably, Shark reproduces the number counts from the GALEX UVto the JCMT 850-microns bands, the redshift distribution of sub-millimetre galaxies, and the ALMA bands number counts. Bravoet al. (2020) showed that Shark reproduces the optical colour distri-bution of galaxies across a wide range of stellar masses and redshiftreasonably well, along with the fraction of passive galaxies as afunction of stellar mass. To produce mock catalogues that mimic the surveys used for theHi stacking experiments in the observations, we embed the galaxypopulation generated by Shark in a survey volume by applying theALFALFA survey’s angular and radial selection function.To construct the lightcone, we use the code stingray (Obreschkowet al, in prep; Chauhan et al. 2019), which is an extended version of thelightcone builder code used in Obreschkow et al. (2009). stingraytiles simulation boxes together to build a 3D field along the line-of-sight of an observer. The galaxies are drawn from simulationsnapshots corresponding to the closest lookback time, which for thisanalysis ranges over 𝑧 = 0 − .
1. To ensure that we have a largestatistical sample, we set the area of our lightcone to be ∼ containing all galaxies with 𝑀 ★ ≥ M (cid:12) . Once this lightcone isbuilt, we compute the SEDs of galaxies as described in Section 2.1.3,and apply in post-processing an apparent 𝑟 -band magnitude cut-offof 19 . 𝑟 -band AB magnitudes < . > 𝑧 ∼ .
25. The survey consists of five fields, amounting to a total sky area of 230 deg with ∼ , The concordance Λ CDM model predicts that all galaxies form andevolve in DM haloes. Galaxy groups are the observable tracer of DMhaloes offering direct insight into the physics that occur inside them.Thus, assigning galaxies to groups becomes imperative to further ourunderstanding of galaxy formation and evolution. The galaxies areassigned to groups by employing ‘group-finding algorithms’ whichaim to associate galaxies with common DM haloes.Many large spectroscopic surveys use group finders based on thefriends-of-friends (FoF) method, which identifies galaxy systems asmember galaxies that are linked by some adopted linkage criteria.There are differences in the way the halo masses of these groupsare later estimated, which we will elaborate on in Section 3. Thegroup-finding algorithm used for our mock survey is the same groupfinder that was used to construct the GAMA galaxy group catalogue(G C; Robotham et al. 2011). Here we give a brief description of theprocess of assigning galaxies to groups.(i)
Step 1:
A luminosity correction is applied to the galaxies that arefed to the group finder. As our mock survey is similar to GAMAin magnitude completeness, a luminosity correction is needed forsome Shark galaxies, whose luminosities are above the brightestgalaxies observed by GAMA. This luminosity correction is requiredas the group finder produces a luminosity function by integratingthe luminosities of galaxies ranging from minimum to the highestdetection limit (see Robotham et al. 2011; Bravo et al. 2020, fordetails), and having brighter than GAMA leads to an integrationerror.(ii)
Step 2:
A FoF algorithm is used to allocate the luminosity-correctedgalaxies into groups. Since we work in redshift space, separate linkinglengths need to be defined along the line of sight and in projection.These linking lengths are denoted by 𝑙 𝑧 and 𝑙 𝑝 , respectively. Thelinking lengths used by the group finder to define a group scaleas a function of the observed density contrast. This leads to linkingparameters that depend on both the position of the group and the faintmagnitude limit of the survey. These parameters therefore need to betuned to the survey completeness (see Equation 1-7 in Robothamet al. 2011).(iii) Step 3:
Once galaxies are assigned to a group, an initial estimate ofthe centre of the group is made by calculating the centre of luminosity.Then the process is iterated after removing the most distant galaxy,and recalculating the centre of luminosity, until only two galaxiesare left, of which the brightest is defined as the group central, andthe rest of the galaxies are flagged as satellites ranked in the order oftheir distance from the centre.For the rest of the paper, we use the ‘mock group catalogue’generated using the procedure described above for our mock-stackingresults. Our mock group catalogue recovers galaxy groups up to 𝑧 = 0 . MNRAS , 1–21 (2020)
I-Stacking: Systematics at play Figure 1.
Stellar mass ( left panel ) and Hi mass ( right panel ) distributions of the central galaxies in the full lightcone constructed using stingray, with thegalaxies generated with Shark. ‘Mock’ includes all the centrals present in the lightcone, with ‘GAMA (selection)’ and ‘SDSS (selection)’ corresponding to theGAMA-like, and the SDSS-like mocks (which include the corresponding flux magnitude cut and group finder definition of centrals). See Sections 2.2 and 2.3for full description of the mocks. an apparent 𝑟 -band magnitude of 17 .
77 mag. The latest group cata-logue, as used by Guo et al. (2020), contains galaxies from the SDSSDR7 Main Galaxy Sample (Albareti et al. 2017). To make our mockgroup catalogue more comparable to the SDSS group catalogue, weremove all galaxies that have an 𝑟 -band AB magnitude > .
77, andadapt the number of galaxies associated with each group accordingly.A caveat of using this approach is that the groups in the SDSS-likeselection will actually be better recovered than the recovery of groupsin the original SDSS catalogue, i.e. the groups recovered will be bet-ter matched to the underlying mock. This is due to us using a deepersurvey to define our initial groups, and then selecting galaxies thatare above the SDSS magnitude limit within those groups. We leavethe creation of a consistent group catalogue for the flux limit of SDSSto future work, as this implies re-calibrating to the observed 𝑟 -bandluminosity distribution as described in Bravo et al. (2020).To confirm that the group finder recovers the underlying baryondistribution, in Figure 1 we plot the distribution of the stellar andHi masses of the central galaxies present in the lightcone for all ofthe galaxies in Shark catalogue (using all the central galaxies inthe lightcone), comparing this against the same distributions afterapplying the FoF group finder and GAMA/SDSS flux limits. Wewould expect the high-mass end in all three criteria to match, as thegroup finder is expected to recover the most luminous centrals. InFigure 1, we define mock as the underlying mock survey describedin Section 2.2, containing all the centrals (as defined by Shark). GAMA-selection refers to the centrals that are defined by the groupfinder that has been run on the mock survey.
SDSS-selection includesall the centrals (as defined by the group finder) that are above theSDSS magnitude limit. A shift in the high-mass end of the stellar-mass distribution (left panel) is seen, with the GAMA and SDSSselection not matching with the lightcone distribution. This shift isa result of introducing luminosity corrections to the galaxies thatwere passed to the group finder, as described earlier. A consequenceof this correction (lowering Shark galaxies luminosities) is that thestellar masses of those galaxies are also decreased assuming a 1:1light-to-mass ratio(Bravo et al. 2020), resulting in the mismatch seenat the high stellar masses. Taking this into account, we can see thatthe SDSS selection is complete for 𝑀 star (cid:38) . M (cid:12) , while theGAMA selection matches the underlying lightcone distribution for 𝑀 star (cid:38) . M (cid:12) . This is to be expected, as the sensitivity of bothsurveys differs by approximately 2 mag. The high-mass end of the Hi mass distribution (right panel) of the centrals matches up in allthe three cases, confirming that the centrals with 𝑀 HI ≥ . M (cid:12) arerecovered irrespective of the survey magnitude limitations. Note thatthe correction to luminosities and stellar masses above is not appliedto the Hi masses ( as the Hi masses of galaxies are not determinedby the optical luminosity of the galaxy), hence the agreement hereis very good. We see the numbers declining at 𝑀 HI ≤ . M (cid:12) ,with the SDSS-like selection having a lower number of galaxies thanthe GAMA-like selection in the same Hi mass bin. This happensbecause at fixed 𝑀 HI , the GAMA limit has the capability of goinglower in stellar mass than the SDSS limit, capturing better thoselow-stellar-mass galaxies that are relatively gas-rich. One of the major possible systematic effects in measuring the HIHMrelation is the estimation of the halo mass of a galaxy group. Thisis easier in simulations than it is for observations. For example,in Shark, halo mass is determined from the VELOCIraptor halocatalogues, using well-defined quantities such as 𝑀 crit200 (hereafter 𝑀 vir ). In this section, we briefly describe how the halo masses areassigned by group finders and then compare them with the halomasses intrinsic to Shark. As stated earlier, it is possible to use certain properties of galaxiesin groups to obtain a mapping to properties of the underlying DMhalo.Therefore, once the galaxy groups have been defined, the nextstep is to assign properties to the underlying DM halo, namely itsmass. The way that halo masses are assigned to galaxy groups variesbased on the group finder and the preferences of the surveyor. Herewe describe the two ways that halo mass estimates were made forSDSS (abundance matching; Yang et al. 2005; Lim et al. 2017) andGAMA (dynamical mass estimation; Robotham et al. 2011).
MNRAS000
MNRAS000 , 1–21 (2020)
G. Chauhan et al.
It is expected that the total luminosity of the galaxies in a halo scaleswith the virial mass of the halo. Using this relation, we can matchup the number densities of galaxies and haloes, a process termed‘abundance matching’. Note that this process assumes no scatter inthe luminosity–halo relation. We use the luminosity–halo relationobtained from Lim et al. (2017) to estimate the halo masses of thegroups belonging to our SDSS-like mock to make a more one-to-onecomparison to the Hi stacking based on SDSS groups. Lim et al.(2017) developed their luminosity–halo relation by comparing theirfit from the mean luminosity–halo relation, derived from SDSS and2dfGRS (The 2dF Galaxy Redshift Survey; Colless et al. 2001) withthe results obtained from the EAGLE hydrodynamical simulation(Schaye et al. 2015; Crain et al. 2015).In order to assign halo masses to the galaxy groups in our cata-logue, we use the following relation for isolated centrals: (1)log 𝑀 h 𝑀 (cid:12) ℎ − = 10 .
595 + 4 . × − exp (cid:18) log 𝐿 c 𝐿 (cid:12) ℎ − . (cid:19) , where 𝑀 h and 𝐿 c refer to the halo mass and the luminosity of thecentral galaxy, respectively. In order to estimate the halo masses ofgalaxy groups, Lim et al. (2017) use the ‘GAP correction’ based onLu et al. (2016) but modified using EAGLE. The GAP correctionuses a combination of the luminosity of the central galaxy and thedifference in the luminosity to the 𝑛 th brightest galaxy in the group toestimate a correction to the halo mass of the galaxy group. We haveused the same procedure as Lim et al. (2017) (see their Equations9-10) to assign halo masses to our galaxy group using the SDSSbest-fitting parameters for the GAP correction (see Table 2 in Limet al. 2017). This means, our estimates are close to the SDSS groupcatalogue used for the Hi stacking in Guo et al. (2020). With the information gathered from the galaxy–galaxy links madeby the FoF algorithm, we can recover the group velocity dispersion( 𝜎 ) and radius ( 𝑅 ), which are key for estimating the dynamical massof the group. Galaxy groups are assumed to be virialized systems,in which case we expect the dynamical mass to scale as 𝑀 ∝ 𝜎 𝑅 .Robotham et al. (2011) use the dynamical mass estimate to assign ahalo mass to the galaxy group, using the following relation: 𝑀 dyn = 𝐴𝐺 ( 𝜎 dyn ) 𝑅 dyn , (2)where 𝑀 dyn is the dynamical mass of the system, 𝜎 dyn and 𝑅 dyn arethe velocity dispersion among the galaxies in the group and radiusof the FoF group, respectively. 𝐺 is the gravitational constant and 𝐴 is a scaling factor dependent on the number of galaxies in a groupand its redshift (see Equation 19 and Tables 2–4 in Robotham et al.2011).Unlike abundance matching, where the DM halo properties arereliant on the halo distribution obtained from Λ CDM predictions, thedynamical mass method depends more on the physically measurableproperties of galaxy groups to estimate the properties of their halo.This makes the dynamical mass method more physical in nature andgives us more information about the dynamics inside a halo. Notethat a velocity dispersion can only be strictly computed with morethan 2 galaxies. Hence, for isolated galaxies (i.e., those that are notassigned to any groups), abundance matching is the only way ofassigning them halo masses.
In order to reproduce the Hi-stacking results in our simulated light-cone, presented in Section 4.2, we compute halo masses as doneby the surveys (i.e. using abundance matching when comparing toSDSS-based stacking, and dynamical masses when comparing toGAMA-based stacking). To fully understand the systematic effectsintroduced by these different methods to measure the halo massesin the derived HIHM relation, we first investigate how well theseestimates reproduce the intrinsic halo masses of the simulation.In this section, we investigate the differences between the halomasses inferred from abundance matching and dynamical mass es-timates, to the intrinsic values in Shark. Here, we focus on galaxygroups with at least 2 or more members ( 𝑁 g ≥ 𝑁 g ≥
5, Robotham et al. 2011).Abundance matching, on the other hand, relies on the mass-to-lightratio from the selected galaxies, which is less sensitive to the groupmembership. It is hence easier to apply, regardless of the occupationof groups. With this knowledge at hand, we divide the groups intwo categories, all and high membership , with all containing groupswith two or more galaxies ( 𝑁 g ≥
2) and high membership containinggroups with at least five galaxies ( 𝑁 g ≥ centrals identified by the group finder back to the simulation boxand use the host halo mass associated with that galaxy as our intrinsicvirial mass for that galaxy group. Note that we do this regardless ofwhether that galaxy is classified as central or satellite in Shark.Figure 2 shows the distribution of halo masses assigned by thegroup finders (abundance matching and dynamical) against the truehalo masses (intrinsic) of the haloes containing those groups. As canbe seen from the left panels , group membership plays an importantrole in the accuracy of the dynamical mass estimate. For 𝑁 g ≥ 𝑁 g ≥
2, thoughwith a significant underestimation in the number of groups for 𝑀 vir ≥ . M (cid:12) . We plot the same distribution for high-membershipgroups (in the right panels ) and see a remarkable improvement in theshape of the dynamical mass estimate distribution, as it now closelyresembles the intrinsic distribution. The abundance-matching distri-bution follows the intrinsic distribution quite well (better than the all distribution), which is to be expected as the Lim et al. (2017) groupmasses were calibrated to reproduce a simulation (though, a differentsimulation and halo finder were used). Though the distributions aresimilar to each other, we do find discrepancies between the groupfinder allocations and the intrinsic distribution. The most prominentdiscrepancy is the underestimation of halo masses at the high-massend, which is seen in both the all and high-membership distributions,for both the abundance matching and the dynamical mass estimates.To investigate the effect of group membership further and to iden-tify where it falters, Figure 3 compares the two halo mass estimatemethods against the true halo masses in Shark, with the left panels comparing abundance matching and right panels looking at the dy-namical mass. The lower panels compare the running medians andpercentiles of the distribution. We find that abundance matching, ir- MNRAS , 1–21 (2020)
I-Stacking: Systematics at play Figure 2.
Halo mass distribution comparison between the Shark output (light grey) and the estimates made by abundance matching ( upper panel ) and dynamicalmass ( lower panel ). The left panel comprises all groups with 𝑁 g ≥
2, whereas the right panel shows the distribution of groups with at least five galaxies, 𝑁 g ≥
5. The distribution obtained for abundance matching is similar in shape to the intrinsic distribution for both group membership criteria at 𝑀 vir ≤ . M (cid:12) and 𝑀 vir ≤ M (cid:12) , for 𝑁 g ≥ 𝑁 g ≥
5, respectively. The number of groups at the high-mass end is significantly smaller for abundance matching thanShark. The distribution derived from dynamical mass is very different from the intrinsic distribution for 𝑁 g ≥ 𝑁 g ≤ 𝑁 g ≥ respective of member allocation, stays close to the 1:1 line, althoughthere is significant scatter around the relation. The same effect isreflected in the median relations plotted below ( lower-left panel ),where the purple and yellow points represent the median of all andhigh-membership groups, respectively, with the shaded region be-ing the 16 th –84 th percentile range of the distribution. As is evidentfrom the left panel of Figure 3, abundance matching underestimatesthe halo masses from 𝑀 vir (cid:38) M (cid:12) and 𝑀 vir (cid:38) . M (cid:12) , and over-estimates the halo masses for 𝑀 vir (cid:46) . M (cid:12) and 𝑀 vir (cid:46) M (cid:12) ,for all and high-membership groups, respectively. This underestima-tion explains the discrepancies seen in the upper panels of Figure 2.Given that the abundance matching is based on a different simulationwith a different halo mass function driven by the different assumedcosmology and halo mass definition, differences like these are to beexpected.For the dynamical mass estimates, we see a significantly largerscatter than the abundance matching comparison for all groups. Thisreverts for high-membership groups, where the scatter displayed bydynamical mass estimates is smaller than the one seen in abun-dance matching ( right panel ). We show the medians for the samedistribution ( right bottom panel in Figure 3), and find that the dy-namical masses overestimate the halo masses for 𝑀 vir (cid:46) . M (cid:12) and 𝑀 vir (cid:46) M (cid:12) , for all and high-membership groups, respec-tively. Though, it can also be seen that for higher halo masses ( 𝑀 vir ∼ . M (cid:12) ), where abundance matching underestimates the halomasses by ∼ . . ∼ . . centrals , satellites and orphans , which are assigned based on the merger trees and sub-halo catalogues used by Shark as a skeleton. We define centrals ( type = 0 ) to be the central galaxy of the most massive subhaloin the group. Satellite ( type = 1 ) galaxies are the galaxies that arehosted by the other subhaloes of the group. Orphans ( type = 2 ) MNRAS000
5, respectively. The number of groups at the high-mass end is significantly smaller for abundance matching thanShark. The distribution derived from dynamical mass is very different from the intrinsic distribution for 𝑁 g ≥ 𝑁 g ≤ 𝑁 g ≥ respective of member allocation, stays close to the 1:1 line, althoughthere is significant scatter around the relation. The same effect isreflected in the median relations plotted below ( lower-left panel ),where the purple and yellow points represent the median of all andhigh-membership groups, respectively, with the shaded region be-ing the 16 th –84 th percentile range of the distribution. As is evidentfrom the left panel of Figure 3, abundance matching underestimatesthe halo masses from 𝑀 vir (cid:38) M (cid:12) and 𝑀 vir (cid:38) . M (cid:12) , and over-estimates the halo masses for 𝑀 vir (cid:46) . M (cid:12) and 𝑀 vir (cid:46) M (cid:12) ,for all and high-membership groups, respectively. This underestima-tion explains the discrepancies seen in the upper panels of Figure 2.Given that the abundance matching is based on a different simulationwith a different halo mass function driven by the different assumedcosmology and halo mass definition, differences like these are to beexpected.For the dynamical mass estimates, we see a significantly largerscatter than the abundance matching comparison for all groups. Thisreverts for high-membership groups, where the scatter displayed bydynamical mass estimates is smaller than the one seen in abun-dance matching ( right panel ). We show the medians for the samedistribution ( right bottom panel in Figure 3), and find that the dy-namical masses overestimate the halo masses for 𝑀 vir (cid:46) . M (cid:12) and 𝑀 vir (cid:46) M (cid:12) , for all and high-membership groups, respec-tively. Though, it can also be seen that for higher halo masses ( 𝑀 vir ∼ . M (cid:12) ), where abundance matching underestimates the halomasses by ∼ . . ∼ . . centrals , satellites and orphans , which are assigned based on the merger trees and sub-halo catalogues used by Shark as a skeleton. We define centrals ( type = 0 ) to be the central galaxy of the most massive subhaloin the group. Satellite ( type = 1 ) galaxies are the galaxies that arehosted by the other subhaloes of the group. Orphans ( type = 2 ) MNRAS000 , 1–21 (2020)
G. Chauhan et al.
Figure 3.
Comparison of halo mass estimate methods with the intrinsic halo masses in Shark.
Upper panels:
Each point represents a galaxy group that hasbeen allocated a halo mass by abundance matching (left) and the dynamical mass method (right), coloured according to the number of group members.
Lowerpanels shows the median for both all (purple-diamonds) and high-membership (yellow-square) groups, for both abundance matching (left) and the dynamicalmass method (right), with the shaded region representing the 16 th -84 th percentile range of the distribution. The dotted line in all the panels represents the 1:1line. Abundance matching is seen to be closer to the 1:1 relation for 𝑀 vir (cid:46) M (cid:12) , but deviates thereafter. Dynamical masses, on the other hand, are not asclose to the 1:1 relation for lower halo masses, but tend to get closer as we move towards higher halo masses Figure 4.
Similar to Figure 3. Each point represents a galaxy group that has been allocated its halo mass by abundance matching ( left panel ) and the dynamicalmass method ( right panel ), coloured according to the intrinsic galaxy type (as defined by Shark) of the central of the galaxy group. In Shark, centrals, satelliteand orphan galaxies are termed type = 0,1 and , respectively. The dotted line in all the panels represents the 1:1 line.MNRAS , 1–21 (2020) I-Stacking: Systematics at play galaxies are the galaxies that cease to be tracked by VELOCIraptor(either because their number of particles drop below the thresholdrequired to consider a detection, or because it becomes indistinguish-able from the underlying density-velocity field), i.e. lack a subhaloentirely.Figure 4 compares the halo mass estimates from the two methodsagainst the true halo masses in Shark, in this case, coloured by theintrinsic galaxy type of the galaxies flagged as centrals by the groupfinder as they appear in Shark. During our analysis, we found that ∼ . The strategy we follow below is to get closer to the way Hi stacking isdone in observations step by step, so that we can critically analyse thesystematic and random effects introduced at each step of the process.We start by estimating the total Hi mass of the haloes defined by thegroup finder. We first obtain the HIHM relation of groups by usingthe true halo masses and then by using the estimated halo masses,using both abundance matching and dynamical masses.
Chauhan et al. (2020) derived a mean HIHM relation directly fromthe Shark outputs and found it to have a distinct shape. Unlikethe observed HIHM relation (see Guo et al. 2020), which was amonotonically increasing function, Chauhan et al. (2020) found a dipin their derived HIHM relation at 𝑀 vir ∼ M (cid:12) , which plateauedup to 𝑀 vir ∼ . M (cid:12) , followed by an increase towards highermasses (see figure 1 in Chauhan et al. 2020). The characteristicshape of the intrinsic HIHM relation is a result of the Hi masscontribution from the central and satellite populations of the haloes.Chauhan et al. (2020) found that central galaxy is the dominant Himass contributor for haloes with 𝑀 vir ≤ . M (cid:12) , with satellitestaking over thereafter. The dip in the HIHM relation seen at 𝑀 vir ∼ M (cid:12) is caused by the AGN feedback. One of the possible reasonsgiven for the discrepancy between the observed and intrinsic SharkHIHM relation was the uncertainty in group definitions around thathalo mass.In Figure 5, we compare the mean HIHM relation as obtainedfrom Shark (Shark-ref) with the mean HIHM relation obtainedwhen we use the groups defined by the group finder. Shark-refrefers to the intrinsic HIHM relation that was obtained for the haloes Figure 5.
Mean Hi content of groups identified by the group finder withina spherical aperture (as labelled) versus the intrinsic halo mass of the cor-responding halo in Shark containing the galaxy flagged as a central by thegroup finder. These apertures are based on the intrinsic halo masses of thegalaxy groups. It can be seen that as we move to larger aperture sizes, the Hicontent of the haloes increases, though the characteristic shape of the intrinsicrelation (maroon solid line) remains qualitatively the same. that have at least two subhaloes (dictated by the resolution of theDM-only simulation), irrespective of whether they are detectable byany survey or not, with the total Hi mass being the sum of all thegalaxies associated with the same host halo in the simulation box at 𝑧 = 0, regardless of their luminosity. We calculate the mean Hi massfor the groups in our mock group catalogue by assuming the central,as defined by the group finder, to be the true central of the host halo,and then calculating the Hi mass of all the galaxies that are withina spherical volume of 1 − 𝑅 vir ) of the hosthalo that galaxy belongs to. The virial mass assumed for this analysisis the true halo mass provided by Shark. From Figure 5 we find thatjust by using the group finder defined groups, we overestimate theHi mass in the region 𝑀 vir ∼ − M (cid:12) . However, even though theuncertainty in group definitions cause a change in the Hi content ofthe group, it does not change the qualitative shape of the relation.As we move to 𝑀 vir ≥ M (cid:12) , the total Hi mass associated withthe halo cannot be captured by an aperture of 1 𝑅 vir around the‘central’ galaxy (although, in this instance it can be captured with anaperture of 1 . 𝑅 vir ). One of the probable reasons for this departurefrom the underlying intrinsic HIHM relation is the misidentificationof satellite galaxies as centrals by the group finder. We find thatalmost 25 (9 .
5) per cent of satellites are misidentified as centrals for 𝑀 vir ∼ M (cid:12) (10 M (cid:12) ). If a satellite residing at the edge of thehalo is misidentified as a central, we consider the position of thatsatellite as the halo centre, and calculate the Hi mass of everything inthe vicinity of that satellite, using the 𝑅 vir of the host halo the satelliteresides in to define our sphere. This would lead us to miss many ofthe galaxies belonging to the original host halo. Another possiblereason for the Hi mass not being captured by 1 𝑅 vir is that the haloestend to be more elongated than spherical at high halo masses, making1 𝑅 vir quite small for the major axis of a massive halo (see figure 1in Cañas et al. 2020). Both of these cases are remedied by using alarger aperture, as demonstrated by using 2 𝑅 vir , which captures atotal Hi mass more in line with the VELOCIraptor haloes. Apartfrom these, there is also the possibility of the ‘true’ central of thathost halo being classified as a central in a different group. This willresult in the same host halo being used twice (or maybe more) forcalculating the Hi mass, which will end up being different each time. MNRAS000
5) per cent of satellites are misidentified as centrals for 𝑀 vir ∼ M (cid:12) (10 M (cid:12) ). If a satellite residing at the edge of thehalo is misidentified as a central, we consider the position of thatsatellite as the halo centre, and calculate the Hi mass of everything inthe vicinity of that satellite, using the 𝑅 vir of the host halo the satelliteresides in to define our sphere. This would lead us to miss many ofthe galaxies belonging to the original host halo. Another possiblereason for the Hi mass not being captured by 1 𝑅 vir is that the haloestend to be more elongated than spherical at high halo masses, making1 𝑅 vir quite small for the major axis of a massive halo (see figure 1in Cañas et al. 2020). Both of these cases are remedied by using alarger aperture, as demonstrated by using 2 𝑅 vir , which captures atotal Hi mass more in line with the VELOCIraptor haloes. Apartfrom these, there is also the possibility of the ‘true’ central of thathost halo being classified as a central in a different group. This willresult in the same host halo being used twice (or maybe more) forcalculating the Hi mass, which will end up being different each time. MNRAS000 , 1–21 (2020) G. Chauhan et al.
Figure 6.
Mean Hi content of groups identified by the group finder within a spherical aperture (as labelled) versus the halo mass estimates from abundancematching (left panel) and dynamical mass (right panel) methods for the corresponding halo in Shark containing the galaxy flagged as a central by the groupfinder. The apertures are also based on the halo masses estimations made for the galaxy groups by abundance matching (left) and the dynamical mass method(right). The departure from the characteristic shape seen when using the halo mass estimations alludes to the fact that correct halo mass estimation plays a majorrole in determining the shape of the mean HIHM relation.
Even though we find slight deviations from the intrinsic HIHMrelation when using group-finder-defined groups, the characteristicshape of the relation remains qualitatively the same. We then investi-gate the systematic effect that halo mass definition has on the derivedHIHM relation.In Figure 6, we repeat the same exercise as for Figure 5, but thistime using the halo masses estimated by abundance matching ( leftpanel ) and the dynamical mass method ( right panel ), which appearin the x-axis and are also used to calculate 𝑅 vir to estimate the totalHi mass associated with that galaxy group. It becomes clear that thehalo mass estimates causes a major deviation in the HIHM shapefrom the intrinsic one. By using the abundance-matching halo massestimates, we find that the characteristic shape of the relation is lost. Itbecomes a monotonically increasing relation, which underestimatesthe mean Hi mass for 𝑀 vir ≤ . M (cid:12) and 𝑀 vir ≥ . M (cid:12) ,and overestimates for the region in between. We do not recover theshape by using dynamical mass estimates either (see right panel) andfind larger deviations from the intrinsic shape than seen when usingabundance matching. The HIHM relation obtained by dynamicalmass estimates remains almost constant for 𝑀 vir ≤ . M (cid:12) , abovewhich it starts to increase.The median for all groups in the case of abundance matching (seelower-left panel in Figure 3), stays close to the true halo massesfor lower-mass haloes. As for the HIHM relation obtained usingdynamical mass, the higher Hi mass estimate at the low-mass endcould be a consequence of centrals being misidentified as satellitesresulting in a higher total Hi mass for a group in that region; wefound ∼
35 per cent of centrals in 𝑀 vir ∼ M (cid:12) groups are actuallysatellites. Figure 6 summarises the effect that uncertainties in the halomass estimates have on the mean HIHM relation, which becomesparamount for our analysis of the parameters involved in Hi stacking,covered in the following section. The general principle of Hi-stacking analysis is to identify the po-sitions and redshifts of galaxies using an external source cataloguecontaining optical redshifts, extract Hi information at these coordi-nates and then co-add the data. In this section, we describe differentHi stacking methods employed by surveys to recover the HIHM re- lation, along with analysing how tweaking their parameters changesthe derived HIHM relation. For the latter we make use of our mockgroup catalogue and survey. Our main aim here is to highlight thedifferences between these stacking techniques and how well theyrecover the underlying intrinsic HIHM relation.
Traditionally, Hi stacking has been done on individual galaxies, usingan optical catalogue to estimate the position and redshift of the targetgalaxy. One then extracts spectra associated with that area withina stacking velocity window and co-adds them, generally weightingeach spectrum by its RMS noise (Fabello et al. 2011).Guo et al. (2020) extended this technique to a galaxy group scale,which we will refer to as ‘group Hi stacking’. They use the SDSSgroup catalogue (Lim et al. 2017) as their optical spectroscopic cat-alogue and ALFALFA as their Hi survey. The ALFALFA (AreciboLegacy Fast ALFA) survey is a ‘blind’ Hi survey that has mappednearly 6900 deg in the Northern Hemisphere, with ∼ ,
500 directdetections (Giovanelli et al. 2005; Haynes et al. 2018) out to redshift 𝑧 = 0 .
06. The typical beam size is ∼ . (cid:48) × . (cid:48) with a velocityresolution of ∼
10 km s − . Guo et al. (2020) use 2 𝑅 as the haloradius as the projected aperture over which Hi stacking is done, withR being the radius containing a mean mass density that is 200times the mean density of the universe at a given redshift (not to beconfused with 𝑅 crit200 , which corresponds to the radius containing amean mass density that is 200 times the critical density of the Uni-verse). Using the software of Fabello et al. (2011), they extracted asingle spectrum for each group, which was then stacked with othergroup spectra in the same mass bin and, finally, fitted with a Gaussianprofile. The Hi mass was measured by integrating the signal within ∼ ± 𝜎 width of the Gaussian profile fitted to the stacked spectra,centred on the central galaxy of the group. This method, in princi-ple, should ensure that all the Hi associated with the target galaxygroups is captured. Figure 7 shows a schematic of how this stackingtechnique, which we refer to as “group Hi stacking", works. Notably,the Hi content of galaxies that were not detected in the parent opticalcatalogue would still contribute to the total Hi mass measured. Butbecause of the large stacking velocity window and projected aperture,some contamination is expected, as the Hi content of galaxies that MNRAS , 1–21 (2020)
I-Stacking: Systematics at play Figure 7.
Illustration of group Hi stacking, showing both the projected view(upper panel) and the orientation in redshift space (lower panel). The ma-roon ellipses symbolise the galaxies that are identified by the group finderas belonging to the group shown. The orange ellipses denote the galaxiesthat are missed by the group finder due to the magnitude limitations of theparent survey, which would include Hi gas in galaxies that are below themagnitude cut of the optical survey and the absence of Hi gas not in galaxiesin the group environment. The blue ellipse represents a galaxy that does notbelong to the group, but is still in the vicinity of the group, and would thuscontribute towards (and contaminate) the stack. The half maroon-half blueellipse represents a misidentified galaxy by the group finder, where the galaxyis not intrinsically a part of the common halo but is still assigned to the groupby the group finder. The Hi content associated with this particular group willcontain the Hi in all galaxies that are in the yellow shaded area. are not actually part of the group would still contribute to the stack.If the adopted apertures (both in projected and velocity space) arevery large, then more contamination is expected. While very smallapertures minimise such contamination, they come at the expense ofexcluding galaxies that belong to that group.In this section, we test if we can reproduce the observationallyderived HIHM relation using our mock survey, when we mimic theobservational procedure.Unlike observations, where they have to consider various obser-vational limitations, like noise and flux sensitivity, our mock surveydoes not suffer from either of those. For keeping our analysis sim-ple, we do not work with Hi emission spectra for mock-stacking ourgalaxies. Instead, we use the projected aperture (as used by observa-tions considered here) and a stacking velocity window, to define thevolume within which we sum the Hi mass of all the galaxies within.For approximating the ∼ ± 𝜎 width, we use Δ 𝑣 = ±
700 km s − as our stacking velocity window, which ends up being close to the ∼ ± 𝜎 width for most of the haloes (barring the most massive ones).This choice was made after discussions with the authors of Guo et al.(2020).To understand the effect of the stacking velocity window, we repeatthe same exercise, but this time adopting Δ 𝑣 = ±
350 km s − aroundthe group central. For both of these Hi stacking measurements, weuse an aperture of 2 𝑅 , which is the aperture used by Guo et al.(2020). The latter is based on the halo mass estimated via abundancematching (see Section 3.1.1). In Figure 8, we show the stacked resultsobtained for Δ 𝑣 = ±
700 km s − and ±
350 km s − , and compare themagainst the HIHM relation presented in Guo et al. (2020). We noticethat by using the Δ 𝑣 = ±
700 km s − for our stacking velocity window,we are able to recover the observed HIHM relation for groups thatcontain at least one galaxy ( 𝑁 g ≥
1; left panel) and for groupswith two or more members ( 𝑁 g ≥
2; right panel) quite closely. Aslight underestimation of Hi mass in the 𝑀 vir range 10 . − M (cid:12) is, however, seen. When we compare the derived HIHM relationobtained using Δ 𝑣 = ±
350 km s − as our stacking velocity window,we find an underestimation of the Hi mass for haloes with 𝑀 vir ≥ M (cid:12) and 𝑁 g ≥
1, and for the entire halo mass range for 𝑁 g ≥ 𝑀 vir ≥ M (cid:12) , which could be due to the stackingwindow not being large enough to encompass all the associatedHi sources of the group. Despite these deviations, we can say thatmimicking the stacking procedure of Guo et al. (2020) yields verygood agreement between simulations and their observations to within0 .
15 dex. This is not the conclusion we would have arrived at had weinstead simply compared the intrinsic relation to the observationalinferences of Guo et al. (2020), as shown in Chauhan et al. (2020).We notice that changing our stacking velocity window does make adifference in the amount of Hi that is measured for a certain group,but it is not a drastic change. An important conclusion here is thatwe are unable to recover the intrinsic HIHM shape by changing Δ 𝑣 ,and we generally see that the derived HIHM relation starts deviatingfrom the intrinsic one at 𝑀 vir ≥ M (cid:12) . The latter is due to thesystematic effects of halo mass definition as shown in Section 4.1.In Figure 9, we keep our stacking velocity window fixed at Δ 𝑣 = ± 𝑅 vir to 2 𝑅 vir , we see achange of ∼ . 𝑀 vir ≥ M (cid:12) ), forboth 𝑁 g ≥ ± 𝜎 window adopted by Guo et al.(2020) for the higher halo masses, where Δ 𝑣 = ±
700 km s − is toosmall compared to the velocity dispersion of the groups at that halomass, which can be as high as, ≈ − (for 𝑀 vir ≈ M (cid:12) ,which corresponds to the last halo mass bin). The velocity dispersionof a halo correlates strongly with the virial velocity. The halo massthat we consider here is obtained from abundance matching, whichwe use to estimate the virial velocity of that halo, thus, ending witha Δ 𝑣 = ± 𝑉 vir , where 𝑉 vir is the virial velocity of the halo. We aimto see if using a Δ 𝑣 that adapts to the halo mass of the galaxy group MNRAS000
700 km s − is toosmall compared to the velocity dispersion of the groups at that halomass, which can be as high as, ≈ − (for 𝑀 vir ≈ M (cid:12) ,which corresponds to the last halo mass bin). The velocity dispersionof a halo correlates strongly with the virial velocity. The halo massthat we consider here is obtained from abundance matching, whichwe use to estimate the virial velocity of that halo, thus, ending witha Δ 𝑣 = ± 𝑉 vir , where 𝑉 vir is the virial velocity of the halo. We aimto see if using a Δ 𝑣 that adapts to the halo mass of the galaxy group MNRAS000 , 1–21 (2020) G. Chauhan et al.
Figure 8.
Mean Hi content of haloes obtained by mock group stacking our mock group catalogue, over a fixed projected aperture, for groups with 𝑁 g ≥ 𝑁 g ≥ 𝑀 vir ≥ M (cid:12) for 𝑁 g ≥ 𝑁 g ≥ Figure 9.
Similar to Figures 8. Here we keep our stacking velocity window fixed at Δ 𝑣 = ±
700 km s − , and change the projected aperture for stacking aslabelled. The effect of changing the projected aperture is small at the low-mass end, but becomes prominent at higher halo masses. Figure 10.
Mean Hi content of haloes obtained by group Hi stacking of the groups identified by the group finder, over a fixed projected aperture, for groups with 𝑁 g ≥ 𝑁 g ≥ 𝑀 vir ≥ M (cid:12) is recovered for both 𝑁 g ≥ 𝑁 g ≥ , 1–21 (2020) I-Stacking: Systematics at play will help recover the intrinsic shape of the HIHM relation. For our‘adaptable stacking velocity window’ test, we use the same aperture(2 × 𝑅 ) as we used for our fixed stacking velocity window test. Wecan see that by using an adaptable window, we are able to recoverthe HIHM relation in the low-mass ( 𝑀 vir ≤ M (cid:12) ) and high-mass( 𝑀 vir ≥ . M (cid:12) ) end for 𝑁 g ≥
1. We are also able to recoverthe HIHM relation for the high-mass end for 𝑁 g ≥
2, even thoughwe overestimate the Hi content at the low-mass end. We do find,however, that the derived HIHM using these stacking parameters inthe mock catalogues is below the observationally derived one of Guoet al. (2020). If we instead use an adaptable stacking velocity windowof 1 . 𝑉 vir , which is closer to the ± 𝜎 window that was adopted forthe observed HIHM relation, we find that there is still a difference of ∼ . 𝑀 vir ∼ M (cid:12) . This could be hinting ata limitation of our simulation, where we are producing satellites thatare too Hi-rich in high-mass haloes. We discuss the limitation andeffect this has in Section 5.3. For the groups with 𝑀 vir ≤ M (cid:12) (see right panels in Figures 8, 9 and 10, which correspond the HIHMrelation for groups with 2 or more galaxy members), which generallyhave few occupants ( 𝑁 g ≤ 𝑀 vir ≥ M (cid:12) , correspondsto the area where the scatter between the halo mass estimates fromabundance matching against the intrinsic halo masses (see Figure 3)is largest. This can cause many smaller groups to be assigned a largerhalo mass than the one they reside in, leading to a bigger aperturesize, which in turn can lead to an increase in the Hi mass due to theinclusion of Hi sources not associated with the group. Even in theabsence of contamination, a systematically higher (lower) halo masswould lead to the Hi mass at fixed halo mass being lower (higher),leading to changes in the shape of the relation (as shown in Figure6). As opposed to the group Hi stacking, where whole group spectraare stacked, ‘individual Hi stacking’ involves stacking individualgalaxies’ Hi spectra that are part of the groups of interest. Figure11 shows a schematic of the individual Hi stacking process. All thegalaxies’ spectra are stacked individually, and then combined in thehalo mass bins of their respective groups, to determine the mean Himass of groups in that halo mass bin. Using this method ensuresthat the contamination from the Hi sources not part of the group (asrecognised by the group finder) is kept at its bare minimum, as thosesources would not be stacked. However, because to its reliance onthe optical group catalogue, the Hi content of galaxies that are part
Figure 11.
Illustration of individual-galaxy Hi stacking, showing both theprojected view (upper panel) and the orientation in redshift space (lowerpanel). The galaxies in the illustration are the same as shown in Figure 7.It can be seen that in the Individual Stacking technique, non-detections thatare part of the group are missed, thus making individual stacking technique away of providing a lower limit to the Hi content of the group, when comparedwith group Hi stacking. of the group but below the detection limit of the optical survey willbe missed.This method has been used by Rhee et al (in prep) to estimatethe Hi mass of galaxy groups in the G23 field of the GAMA sur-vey, observed using the ASKAP (Australian Square Kilometre ArrayPathfinder; Johnston et al. 2008) radio telescope as a part of the early-science phase of DINGO (Deep Investigation of Neutral Gas Origin;Meyer 2009). DINGO is a deep Hi survey that aims to observe theHi gas content of galaxies out to 𝑧 ∼ .
43 in the five GAMA regions.For their Hi stacking, Rhee et.al. (in prep) used the early-science datafor DINGO, which has a velocity resolution of ∼ − , with asynthesised beam size corresponding to ∼ (cid:48)(cid:48) × (cid:48)(cid:48) . The externalspectroscopic optical catalogue used was from Bellstedt et al. (2020),which is an update from the initial Robotham et al. (2011) catalogue,with the photometry performed using ProFound (Robotham et al.2018) for the GAMA G23 region. For the Hi stacking, Rhee et.al.(in prep) extract the Hi spectrum over the spatial pixels covering asquare sky region of ∼
49 kpc centred on the target galaxy position,over a velocity window of 300 km s − , calculating the flux densityfor each spectral channel using the method described in Shostak & https://dingo-survey.org/ MNRAS000
49 kpc centred on the target galaxy position,over a velocity window of 300 km s − , calculating the flux densityfor each spectral channel using the method described in Shostak & https://dingo-survey.org/ MNRAS000 , 1–21 (2020) G. Chauhan et al.
Figure 12.
Mean Hi content of haloes obtained by individual Hi stackingof the groups identified by the group finder, over a fixed aperture coveringa square sky region of 49 kpc, for groups with 𝑁 g ≥
2, and comparing itagainst the intrinsic HIHM relation obtained from Shark. The scatter pointsshow the observed HIHM relation as obtained from Rhee et. al. (in prep),by stacking the galaxy groups in G23 field. The observed, mock-stacked andintrinsic HIHM relations all agree for 10 . M (cid:12) ≤ 𝑀 vir ≤ M (cid:12) , at thelevel of observational uncertainties. Allen (1980). The extracted spectra are then co-added, weighted bythe RMS noise. In this analysis, we consider the Hi stacking resultsof Rhee et al (in prep) for the redshift range of 0 . < 𝑧 < . 𝑧 = 0 . 𝑁 g ≥
2) derived by Rheeet al. to our mock catalogue. For our measurements, we stack all thegalaxies in our group catalogue, using a projected aperture covering asquare sky region of 49 kpc and adding all the Hi sources around ourtarget galaxy in a velocity window of ±
300 km s − . We group thesegalaxies according to their corresponding group dynamical mass andcalculate the average Hi mass in each bin. Our mock-stacked andthe observed HIHM relations show similar features, in very goodagreement, and remain fairly constant over the entire range. Ourmock-stacking results are close to the Hi mass predictions made byour intrinsic HIHM relation for 10 M (cid:12) ≤ 𝑀 vir ≤ M (cid:12) . Theobserved HIHM points are close to the intrinsic HIHM relation inthat region. Surprisingly, the HIHM relation obtained from groupHi stacking (see Figure 8) starts diverging from the intrinsic HIHMrelation at 𝑀 vir ≥ M (cid:12) . This hints at the fact that group Histacking at these halo masses starts to include too many Hi sourcesthat are likely not part of the group. For the range, 10 M (cid:12) ≤ 𝑀 vir ≤ M (cid:12) , the magnitude limit of GAMA is sufficient to recoverall the relevant Hi sources associated with the group, as can beseen from the agreement between the mock galaxy stacking and theintrinsic HIHM relation. However, as we move to 𝑀 vir ≥ M (cid:12) ,the GAMA magnitude limit ceases to be deep enough to recoverall the Hi sources. For the groups with 𝑀 vir ≤ M (cid:12) , which aregenerally small groups ( 𝑁 g ≤ highlighting fundamental limitations of the current observationaltechniques.4.2.3 Central and satellite HI estimates In the previous sections, we have discussed the effect that groupstacking and individual stacking have on the total average Hi mass ofa halo. Here, we analyse how the central and satellite contributionsto the total Hi content of the halo are affected by the Hi stackingtechnique used.In Figure 13, we compare the Hi stacking results for centrals (leftpanel) and satellites (right panel) for 𝑁 g ≥
1, from Guo et al. (2020)with our mock group stacking results. We compare them against thecontribution of Hi mass by centrals and satellites to the intrinsicHIHM relation. For stacking the centrals in our mock catalogue, wefollow the approach adopted by Guo et al. (2020). We take a projectedaperture of 200 kpc around the central and sum the Hi masses ofall the galaxies that are within a velocity window of Δ 𝑣 = ± − . We find that we are able to reproduce the observed Hi masscontribution of centrals to the HIHM relation with our mock stackingquite well. In the region of 𝑀 vir ≤ M (cid:12) , we find that our mock-stacked and observed Hi values of the centrals, are comparable withthe intrinsic Hi mass contribution of centrals. This is not surprising,as the major contributor to the total Hi mass of haloes in this regionare isolated centrals. Any contribution from confused satellites inthis region will be small. As soon as we enter into the regime wherecentrals are typically no longer isolated (i.e. have satellites), 𝑀 vir ≥ M (cid:12) , we start seeing the effect of confusion. That is, thelarge aperture and wide velocity window allows some satellites tocontribute towards the Hi mass inferred for the centrals. The Hi massof centrals steadily increases as we go to higher halo masses ( 𝑀 vir ≥ M (cid:12) ), and the distance between the intrinsic central Hi mass andthe mock-stacked Hi mass of centrals keeps increasing (maximumbeing 0 . 𝑀 vir ≈ . M (cid:12) ), before coming closer againat 𝑀 vir ≈ M (cid:12) . The latter happens because the satellites closeto centrals at galaxy cluster scales are extremely gas-poor or evendevoid of gas. To provide a comparison of how much difference abigger stacking velocity window has on the total Hi mass of thecentral, we repeat the same exercise, but this time with Δ 𝑣 = ± − , and find that the Hi mass measured for centrals is slightlyhigher than with the Δ 𝑣 = ±
300 km s − counterpart.For estimating the satellite Hi contribution (right panel) to ourmock Hi-stacked groups, we subtract the Hi mass measured for thecentral (with a stacking velocity window of Δ 𝑣 = ±
300 km s − ) fromthe total Hi content of the halo, which we had measured for a groupstacking velocity window of Δ 𝑣 = ±
700 and 350 km s − , mimickingthe approach of Guo et al. (2020). We are able to reproduce theobserved Hi mass contribution of satellites by subtracting the centralHi contribution from the total Hi mass measured with a stackingwindow of Δ 𝑣 = ±
700 km s − , which is closer to the values Guo et al.(2020) used. Our Hi estimates for satellites derived by mock stackingshow a similar shape to the intrinsic Hi contribution of satellites,but display a different slope. Our mock-stacking results overestimatethe underlying satellite intrinsic Hi contribution for 𝑀 vir ≤ M (cid:12) , MNRAS , 1–21 (2020)
I-Stacking: Systematics at play Figure 13.
Central (left) and satellite (right) contributions to the total Hi content of haloes obtained by stacking as shown in previous figures. The results areshown for two stacking windows Δ 𝑣 = ±
300 km s − (initially used by Guo et al. 2020) and Δ 𝑣 = ±
700 km s − , for comparison. The black crosses are theobserved central and satellite contributions from Guo et al. (2020). Our mock-stacking results are able to reproduce the observed group-stacked HIHM relationfor centrals and satellites. Though, they only match with the intrinsic HIHM relation for 𝑀 vir ≤ M (cid:12) and 𝑀 vir ≥ M (cid:12) for central and satellite Hicontribution, respectively. and underestimate thereafter. The difference between the Hi massmeasured by observations and our mock-stacked results in the halomass bin of 𝑀 vir ≈ M (cid:12) is a physical effect. This is due to thelack of modelling of ISM stripping for Shark satellites. It is likelywe would find less Hi in Shark satellites if the SAM included ISMstripping from environmental effects (Stevens & Brown 2017). Thelack of ISM-stripping modelling is most significant for haloes with 𝑀 vir ≥ M (cid:12) . This is because generally the intra-group mediumis not dense enough as to ram pressure strip significant amounts ofthe ISM content of satellites (see Marasco et al. 2016, for an analysisof this in hydrodynamical simulations).Unlike the central Hi measurement for group stacking, we find thatwith individual stacking the observed measurements remain lowerthan the intrinsic central Hi contribution. Figure 14 compares the Hicontributions of centrals (left panel) and satellites (right panel) to thetotal Hi mass of the group for (i) observations as will be presentedin Rhee et al. (in prep), (ii) our mock individual stacking results and(iii) the intrinsic contributions as calculated from Shark-ref. Despitebeing able to reproduce the group Hi mass measurements for theindividual stacking observations via our mock-stacking technique,we predict higher Hi masses for centrals in the 10 M (cid:12) ≤ 𝑀 vir ≤ . M (cid:12) region. This difference is caused by the higher percentageof red centrals that are observed by Rhee et al. (in prep) in thatmass bin. Rhee et al. (in prep) stack the Hi spectra from red andblue galaxies, separately, measuring different mean Hi masses forthem and then using these to estimate the mean Hi mass of thecentrals in a halo mass bin. Red galaxies tend to be gas-poor. Evenwhen stacked, red centrals contribute relatively little Hi comparedto their blue counterparts. The fraction of red centrals in the regionof disagreement between mock-stacked results and observations ishigher than the blue centrals. This brings down the observed averageHi mass of centrals in that halo mass bin. As for our mock-stackedcentrals, we find that the region of 10 M (cid:12) ≤ 𝑀 vir ≤ . M (cid:12) is dominated by lowly star-forming galaxies (which might or mightnot appear red when observed). The average Hi mass associatedwith them is 𝑀 HI ≈ . × M (cid:12) . Compared to our values, Rhee etal. (in prep) find the average Hi mass of their red centrals to be 𝑀 HI ∼ . × M (cid:12) , and for their blue centrals to be 𝑀 HI ∼ . × M (cid:12) .Individual Hi stacking uses both a small aperture and velocity window, thereby minimising the effect of confused sources. In ouranalysis for mock individual Hi stacking, detections of a secondsource of Hi mass around the targeted galaxy (which could contributeto confusion) are negligible, with ≤ 𝑀 HI ≥ . M (cid:12) , the radius at which the density drops below 1M (cid:12) /pc islarger than 49 kpc (see Wang et al. 2016; Stevens et al. 2019b),thus a fixed projected aperture for galaxies at those Hi masses mightnot be enough to encompass all the Hi mass associated with them.Irrespective of the difference between observed and mock-stackedresults, the characteristic shape of the central HIHM relation is lost.Our mock-stacked centrals lose the bump seen in the intrinsic relationat 𝑀 vir ≈ M (cid:12) , though they start agreeing with the intrinsiccentral Hi mass from 𝑀 vir ≥ . M (cid:12) .In the right panel of Figure 14, we find that the observed Hi stack-ing measurements for satellites are very different from the intrinsicsatellite Hi mass contribution. Our mock Hi stacking results under-estimate the Hi contribution of satellites for 𝑀 vir ≤ M (cid:12) andoverestimate thereafter, compared with Rhee et al. (in prep). Un-like the observed stacking points, which show a weak decline fromlow masses up to 𝑀 vir ≤ . M (cid:12) and rise thereafter (barringthe last halo mass bin), our mock-stacking results show a steady in-crease in the satellite Hi contribution as we move towards higher halomasses. The satellite Hi mass contribution rises sharply after 𝑀 vir ≥ . M (cid:12) in our mock-stacking results, which is in tension with theobserved stacking results. At these high halo masses, the differencebetween the observations and mock-stacking results is ∼ MNRAS000
700 km s − , for comparison. The black crosses are theobserved central and satellite contributions from Guo et al. (2020). Our mock-stacking results are able to reproduce the observed group-stacked HIHM relationfor centrals and satellites. Though, they only match with the intrinsic HIHM relation for 𝑀 vir ≤ M (cid:12) and 𝑀 vir ≥ M (cid:12) for central and satellite Hicontribution, respectively. and underestimate thereafter. The difference between the Hi massmeasured by observations and our mock-stacked results in the halomass bin of 𝑀 vir ≈ M (cid:12) is a physical effect. This is due to thelack of modelling of ISM stripping for Shark satellites. It is likelywe would find less Hi in Shark satellites if the SAM included ISMstripping from environmental effects (Stevens & Brown 2017). Thelack of ISM-stripping modelling is most significant for haloes with 𝑀 vir ≥ M (cid:12) . This is because generally the intra-group mediumis not dense enough as to ram pressure strip significant amounts ofthe ISM content of satellites (see Marasco et al. 2016, for an analysisof this in hydrodynamical simulations).Unlike the central Hi measurement for group stacking, we find thatwith individual stacking the observed measurements remain lowerthan the intrinsic central Hi contribution. Figure 14 compares the Hicontributions of centrals (left panel) and satellites (right panel) to thetotal Hi mass of the group for (i) observations as will be presentedin Rhee et al. (in prep), (ii) our mock individual stacking results and(iii) the intrinsic contributions as calculated from Shark-ref. Despitebeing able to reproduce the group Hi mass measurements for theindividual stacking observations via our mock-stacking technique,we predict higher Hi masses for centrals in the 10 M (cid:12) ≤ 𝑀 vir ≤ . M (cid:12) region. This difference is caused by the higher percentageof red centrals that are observed by Rhee et al. (in prep) in thatmass bin. Rhee et al. (in prep) stack the Hi spectra from red andblue galaxies, separately, measuring different mean Hi masses forthem and then using these to estimate the mean Hi mass of thecentrals in a halo mass bin. Red galaxies tend to be gas-poor. Evenwhen stacked, red centrals contribute relatively little Hi comparedto their blue counterparts. The fraction of red centrals in the regionof disagreement between mock-stacked results and observations ishigher than the blue centrals. This brings down the observed averageHi mass of centrals in that halo mass bin. As for our mock-stackedcentrals, we find that the region of 10 M (cid:12) ≤ 𝑀 vir ≤ . M (cid:12) is dominated by lowly star-forming galaxies (which might or mightnot appear red when observed). The average Hi mass associatedwith them is 𝑀 HI ≈ . × M (cid:12) . Compared to our values, Rhee etal. (in prep) find the average Hi mass of their red centrals to be 𝑀 HI ∼ . × M (cid:12) , and for their blue centrals to be 𝑀 HI ∼ . × M (cid:12) .Individual Hi stacking uses both a small aperture and velocity window, thereby minimising the effect of confused sources. In ouranalysis for mock individual Hi stacking, detections of a secondsource of Hi mass around the targeted galaxy (which could contributeto confusion) are negligible, with ≤ 𝑀 HI ≥ . M (cid:12) , the radius at which the density drops below 1M (cid:12) /pc islarger than 49 kpc (see Wang et al. 2016; Stevens et al. 2019b),thus a fixed projected aperture for galaxies at those Hi masses mightnot be enough to encompass all the Hi mass associated with them.Irrespective of the difference between observed and mock-stackedresults, the characteristic shape of the central HIHM relation is lost.Our mock-stacked centrals lose the bump seen in the intrinsic relationat 𝑀 vir ≈ M (cid:12) , though they start agreeing with the intrinsiccentral Hi mass from 𝑀 vir ≥ . M (cid:12) .In the right panel of Figure 14, we find that the observed Hi stack-ing measurements for satellites are very different from the intrinsicsatellite Hi mass contribution. Our mock Hi stacking results under-estimate the Hi contribution of satellites for 𝑀 vir ≤ M (cid:12) andoverestimate thereafter, compared with Rhee et al. (in prep). Un-like the observed stacking points, which show a weak decline fromlow masses up to 𝑀 vir ≤ . M (cid:12) and rise thereafter (barringthe last halo mass bin), our mock-stacking results show a steady in-crease in the satellite Hi contribution as we move towards higher halomasses. The satellite Hi mass contribution rises sharply after 𝑀 vir ≥ . M (cid:12) in our mock-stacking results, which is in tension with theobserved stacking results. At these high halo masses, the differencebetween the observations and mock-stacking results is ∼ MNRAS000 , 1–21 (2020) G. Chauhan et al.
Figure 14.
Similar to Figure 13, the central (left) and satellite (right) contributions to the total Hi content of haloes as measured by individual Hi stacking. Theblack points with error-bars are the individual stacking results from Rhee et al. (in prep). Our mock-stacking results estimate higher Hi masses for group centralsin the halo mass region of 10 M (cid:12) ≤ 𝑀 vir ≤ M (cid:12) . the satellite Hi contribution in 𝑀 vir ≤ . M (cid:12) and underestimatethereafter. We suspect the higher satellite contribution in our mockHi stacking result at lower halo masses ( 𝑀 vir ≤ . M (cid:12) ) could be adirect consequence of centrals being grouped together causing themto occasionally be misidentified as satellites. This will result in ahigher Hi mass contribution from the misidentified galaxy, resultingin increasing the mean Hi contribution from satellites in that halomass bin. As described in the previous sections, despite using the same under-lying mock catalogue for our Hi stacking experiments, changing thestacking technique leads to a remarkable change in the Hi measuredfor galaxy groups. Our main aim in this section is to discuss whatcauses this change in the shape, and also assess the limitations of theHIHM relation derived using different stacking techniques.
The main difference between the intrinsic HIHM relation and the onemeasured via stacking is the group definitions. Intrinsic HIHM rela-tion uses the halo merger tree catalogues to identify all the subhaloesassociated with the host halo, and then calculates the Hi masses of thegalaxies residing in those subhaloes. As for our Hi-stacking experi-ment, we rely on the groups as defined by our FoF based group finder,which links galaxies together based on some adopted linkage crite-ria. Where an ideal group finder would result in a perfect mappingto a simulation halo catalogue, group finders have to work within thelimitations of observations in redshift space. According to Campbellet al. (2015), there are three major challenges faced by a group finder:(i) halo mass estimation, (ii) central/satellite designation error, and(iii) group-member allocation.Already discussed in detail in Section 3.2, comparing the halomass assigned by a group finder with the intrinsic halo mass of thesimulation results in a large scatter. Dynamical mass estimates showa particularly large scatter for 𝑀 vir ≤ M (cid:12) . This region is mainlydominated by groups with a low occupancy, mostly 𝑁 g ≤
5, for whichthe dynamical mass estimates are not reliable and are prone to errors.We demonstrate that for high-membership groups, dynamical massestimates show a significantly smaller scatter than the abundance matching estimates, and performs better at recovering the intrinsichalo mass. Although, abundance matching is more reliable in the low-halo-mass region, there is still a scatter of ∼ . M (cid:12) ≤ 𝑀 vir ≤ . M (cid:12) mass region, which is washed out due to theuncertainties inherent to the method. This is a large limitation of thestacking techniques discussed here, as it is exactly around this massregion that AGN feedback imprints its effect, and is also the regionthat shows the largest differences among simulations (see figure 2 inChauhan et al. 2020).The second challenge of the group finder is to identify the cen-tral galaxy of the galaxy group, which has paramount importancewhen estimating the halo mass for abundance matching. The bright-est member of a galaxy group is normally recognised as a centralgalaxy of the group, based on the idea that central galaxies growin mass by cannibalising their satellites (Dubinski 1998; Cooray &Milosavljević 2005). This way of identifying centrals works in a sta-tistical sense, but it has been shown by van den Bosch et al. (2005);Skibba et al. (2011) that in a significant fraction of dark matter haloesthe brightest group member is a satellite rather than a true central.This is usually the case in dynamically young groups, in which themost massive galaxies have not yet had time to merge. Campbell et al.(2015) quantify this fraction and show that the fraction of brightestgalaxies not being centrals ranges from ∼
10 per cent to ∼
30 per centfor haloes of masses 10 M (cid:12) and 10 M (cid:12) , respectively. Misidenti-fication of the central galaxy would affect the luminosity and thus,the halo mass estimates made by the abundance matching technique.This is highlighted by the fact that the scatter in the virial mass com-parison made for the abundance-matching allocation increased withincreasing halo mass (see Figure 4) – which also coincides with theincrease in the fraction of bright satellites. As dynamical mass is notdependent on the luminosities for their halo mass estimates, they donot suffer from it, but still tags a fraction of satellites as centrals.The final challenge of the group finder is to correctly partitiongalaxies into their groups. A group finder can fail at this by eitherassigning galaxies from different haloes to the same group or byassigning galaxies from the same halo to different groups. The formercase will lead to two true centrals being assigned to a single group,which will lead to one of them being recognised as a satellite andthus provide spurious results. For the latter case, a true satellite couldbe identified as a central, leading again, to errors. The effect of both MNRAS , 1–21 (2020)
I-Stacking: Systematics at play Figure 15.
Comparing the contamination from Hi sources not part of the group for group (left panel) and individual (right panel) Hi stacking measurementsfor high-membership groups. The red points are the observed HIHM relation from Guo et al. (2020) for 𝑁 g ≥ 𝑀 vir ≥ M (cid:12) , whereasindividual Hi stacking shows minimal contamination for 10 . M (cid:12) ≤ 𝑀 vir ≤ . M (cid:12) . these phenomena results in deviations in halo mass estimates, whichin turn would lead to uncertainty in the Hi mass measured by Histacking of those groups. In Section 4.2, we mentioned that group Hi stacking is expected tohave some contamination from the Hi content of galaxies not partof the group, due to those galaxies falling within the large projectedaperture and velocity window employed for the stacking. Individ-ual Hi stacking, by virtue of using a smaller aperture and stackingvelocity window, should not suffer from as much contamination.To quantify the contamination, we compare the Hi mass measuredvia Hi stacking with the total Hi mass contained in the “best matchhalo" counterpart of the group being stacked in our simulation box.We define a halo to be the best match counterpart of our galaxygroup by comparing the number of galaxies residing in that groupwith the number of GAMA/SDSS detectable galaxies residing in acommon halo. We use a ‘purity fraction’ (Robotham et al. 2011) forour analysis, which is defined as 𝑓 pf = ( 𝑁 galshared ) 𝑁 galgroup × 𝑁 galhalo , (3)where 𝑁 galshared , 𝑁 galgroup and 𝑁 galhalo is the number of GAMA or SDSS de-tectable galaxies that are shared between the intrinsic halo in Sharkand the FoF group in the mock group catalogue, the total numberof galaxies in the FoF mock group, and the total number of GAMAor SDSS detectable galaxies in the intrinsic Shark halo, respec-tively. We consider the halo with the highest purity fraction to bethe best match counterpart of our galaxy group. For example, if agalaxy group defined by the group finder, consists of five membergalaxies, where three of these galaxies share a halo that containsnine detectable members, while the other two share a halo with threedetectable member galaxies, the purity fraction for these two haloeswill be 3 / × / . / × / ≈ .
27, respectively. Thereforethe latter is considered the best match. Once we have identified abest-match halo, we compute the total Hi mass content of that haloby simply adding the Hi content of all galaxies belonging to it and use that as our ‘true’ Hi mass associated with the group. Any de-viation from this result (by stacking) is considered to be caused bycontamination from Hi sources not part of the halo. We also producean intrinsic HIHM line for comparison, where we consider the to-tal Hi mass of the haloes in our lightcone that have at least five ormore GAMA or SDSS detectable galaxies. For this comparison, wemock stack only the high-membership groups which have had theirbest match halo identified (94 per cent of the total high-membershipgroups).In Figure 15, we compare the contamination in the measured Himass of high-membership groups by both group (left panel) and indi-vidual (right panel) Hi stacking. We can see that the Guo et al. (2020)observed HIHM relation for high-membership groups is higher from 𝑀 vir ≥ . M (cid:12) onwards, relative to both the intrinsic and the best-match relations. The mock-stacked (group) measurements (see leftpanel) underestimate the Hi mass for 𝑀 vir ≤ M (cid:12) , agree for 𝑀 vir ∼ − . M (cid:12) and overestimate for 𝑀 vir ≥ . M (cid:12) , comparedto the observed HIHM relation. Irrespective of the agreement withthe observed HIHM relation, the HIHM relation derived from themock-stacked (group) Hi measurements starts suffering the effects ofcontamination from 𝑀 vir ≥ . M (cid:12) and 𝑀 vir ≥ . M (cid:12) for thefixed and adaptable velocity stacking window, respectively. The levelof contamination – that is, the difference between the Hi mass ofthe best-matched haloes (yellow triangles) and the derived Hi massvia group stacking – is largest at 10 . M (cid:12) ≤ 𝑀 vir ≤ . M (cid:12) ,according to our mock catalogue. In our analysis we find that, onaverage, the fraction of contaminants that are central galaxies fromother smaller haloes is ∼
81 per cent. These centrals have 10 M (cid:12) ≤ 𝑀 ★ ≤ . M (cid:12) and are hosted in haloes of 10 . M (cid:12) ≤ 𝑀 vir ≤ . M (cid:12) . The average stellar mass and halo mass weighted by theHi mass of the contaminant galaxies are 10 . M (cid:12) and 10 . M (cid:12) ,respectively. We find that these galaxies have a median Hi-to-stellarmass ratio (cid:39) . M (cid:12) ≤ 𝑀 vir ≤ . M (cid:12) ,signalling that the contamination is at its bare minimum when usingthis Hi stacking technique. This also signifies that the magnitude limitof GAMA (apparent 𝑟 -band magnitude ≥ .
5) is able to recover allthe prominent Hi sources in that halo mass region that contribute
MNRAS000
MNRAS000 , 1–21 (2020) G. Chauhan et al.
Figure 16.
Comparing the Hi mass contribution from the galaxies in the best-match haloes above the GAMA detection limit (black crosses), against theexpected Hi mass measurements from WAVES for the best-match haloes inour simulation for high-membership groups. It becomes evident that both theGAMA and WAVES detection limits are able to recover the major Hi sourcesfor 10 . M (cid:12) (cid:46) 𝑀 vir (cid:46) . M (cid:12) , and remain equivalent to the true mean Himass of the best-match haloes. For 𝑀 vir > . M (cid:12) , the expected Hi massfrom WAVES starts diverging from both the GAMA estimates and the trueHi mass of the best-match haloes, where it recovers more Hi than GAMA butnot all. An increase of 0 . 𝑀 vir ≥ . M (cid:12) . We use the virial mass from Shark asour halo mass in this plot. meaningfully to the total Hi content of groups. As we move towards 𝑀 vir ≥ . M (cid:12) , our mock-stacked (individual) Hi measurementsare lower than the total Hi mass of the corresponding best-matchhalo. This is due to GAMA magnitude limit not being sufficient todetect all the major Hi contributing galaxies, as the majority of them(from this analysis) appear to be under the detection limit of GAMA.During our analysis, we found that the median 𝑟 -band and 𝑧 -band magnitude of satellites residing in haloes of masses 𝑀 vir ∼ (10 . M (cid:12) , 10 M (cid:12) , 10 . M (cid:12) , 10 M (cid:12) ) is (20.65, 22.04, 22.72,22.83) mag and (20.18, 21.76, 22.46, 22.49), respectively. In Table2 we list the 5 th –95 th percentile distribution of the galaxies for 𝑟 -and 𝑧 -band magnitudes for different halo mass bins. We find thatfrom the total number of satellites belonging to the haloes of masses 𝑀 vir ∼ (10 . M (cid:12) , 10 M (cid:12) , 10 . M (cid:12) , 10 M (cid:12) ), GAMA is able todetect ∼ (45.4, 27.6, 10.7, 7.6) per cent of satellites, which make upfor (84.8, 68.8, 34.4, 21.5) of the total Hi mass of those haloes (seeTable 2). Even though the GAMA magnitude limit is sufficient fordetecting all the major Hi sources in the 10 M (cid:12) ≤ 𝑀 vir ≤ . M (cid:12) range, it is too bright to recover the major Hi contributing galaxiesfor 𝑀 vir ≥ . M (cid:12) . The satellite population also spreads about4 mag in both 𝑟 -band and 𝑧 -band magnitudes, though the majorityof them might not be the major Hi contributors. Nevertheless, theneed for a deeper spectroscopic survey arises to capture more of themissed faint satellites.The Wide-Area VISTA Extragalactic Survey (WAVES; Driveret al. 2019) is an upcoming spectroscopic survey that aims to map1200 deg of sky area up to redshift 𝑧 = 0 .
2, complete to an appar-ent 𝑧 -band magnitude of 21 . ∼ . 𝑀 vir ≥ . M (cid:12) . We find that with the deeper spectroscopic range ofWAVES, we will be able to detect ∼ (60, 46, 28, 25) per cent ofall satellites belonging to the haloes at 𝑀 vir ∼ (10 . M (cid:12) , 10 M (cid:12) ,10 . M (cid:12) , 10 M (cid:12) ), and recovering ∼ (88.5, 78.8, 54.2, 51.8) percent of the total Hi mass at those halo masses, according to our mockcatalogue (see Table 2). The difference between the intrinsic Hi masspredicted by Shark for 𝑀 vir ∼ . M (cid:12) , and the Hi mass measuredby GAMA and a survey of the depth of WAVES is 0 . . M (cid:12) ≤ 𝑀 vir ≤ . M (cid:12) . The final piece to solving the puzzle of the discrepancy between theintrinsic and observed HIHM relation lies in the modelling of gasstripping in satellites in Shark. In Shark, the model of “instanta-neous ram-pressure stripping" (as described in Lagos et al. 2014) isused, which assumes that as soon as galaxies become satellites, theirhalo gas is instantaneously stripped and transferred to the hot gas ofthe central galaxy. Thus, gas is only allowed to accrete onto the cen-tral galaxy in the halo and not onto the satellites. In addtion, satellitegalaxies are cut off cosmological accretion as soon as they becomesatellites, which prevents their halo gas from being replinished.Though the halo gas of the satellites is stripped in Shark, the coldgas in the discs of the galaxies is not stripped. This lack of modellingof “ISM stripping" does not affect the total Hi in haloes in the regimewhere the major Hi contributor is the central galaxy. As we move intothe regime of halo masses where satellites are the major Hi reservoir,we see the effects of ram pressure stripping. Stevens & Brown (2017);Stevens et al. (2019a) in their analysis, find that the ram pressurestripping will start affecting the total Hi mass of the haloes as soon aswe enter the regime where satellites dominate the Hi content of haloes(over centrals). The effect of ram pressure stripping continuouslyincreases with halo mass, with satellites residing in more massivehaloes displaying lower average gas content relative to those in lower-mass haloes. Other simulations, however, show that ram pressurestripping becomes effective only at 𝑀 vir ≥ M (cid:12) (Marasco et al.2016). We find that the hot gas stripping model employed by Shark issufficient to match observations well for satellites residing in haloesat 𝑀 vir ≤ M (cid:12) , as displayed by our intrinsic satellite Hi masscontribution matching with the observed satellite Hi mass in thegroup Hi stacking results (see Figure 13). Though, at higher halomasses, a higher Hi mass is seen in both mock-stacked as well theintrinsic HIHM relation for the last halo mass bin (corresponding to 𝑀 vir ≈ M (cid:12) ). This is when the tension between the mock stackedand observed Hi mass is seen. This shows that this is an importantarea that requires revision in the Shark model, which we leave forfuture work. MNRAS , 1–21 (2020)
I-Stacking: Systematics at play Halo Mass 𝑟 -band [mag] 𝑧 -band [mag] GAMA galaxies WAVES galaxies[M (cid:12) ] 5 𝑡ℎ 𝑡ℎ 𝑡ℎ 𝑡ℎ 𝑡ℎ 𝑡ℎ Satellites Hi mass Satellites Hi mass10 . . . . . . . .
4% 84 .
8% 60% 88 . . . . . . . .
6% 68 .
8% 46% 78 . . . . . . . . .
7% 34 .
4% 28% 54 . . . . . . . .
6% 21 .
5% 25% 51 . Table 2.
Table shows the 5 th –95 th percentile distribution of the galaxies of our mock-survey, along with the median value (50 th percentile), in 𝑟 - and 𝑧 -bandsresiding in different halo mass bins. The column ‘Satellites’ in the GAMA- and WAVES-galaxies section refer to the percentage of galaxies that are detectedby the respective surveys to the total number of galaxies residing in that halo. The ‘Hi mass’ column refers to the percentage Hi mass contribution of thesedetectable galaxies to the total Hi mass of the halo they reside in. We can see that with a deeper survey, such as WAVES, there is a significant increase in thedetection of galaxies which will provide better constrains for the HIHM relation. In this work, we have used the Shark semi-analytic galaxy formationmodel to create a mock survey with the area 6900 deg and redshift 𝑧 ≤ .
1. We have also produced a corresponding mock group cat-alogue for galaxy groups till redshift 𝑧 = 0 . • The correct estimation of halo masses for galaxy groups plays amajor role in defining the shape of the derived HIHM relation. Ir-respective of a group being well recovered or not, as soon as westart using the halo mass estimates from abundance matching or thedynamical mass method, we lose the characteristic intrinsic HIHMshape. Abundance-matching halo mass estimates are more reliablefor isolated centrals and small group ( 𝑁 𝑔 ≤
5) compared to thosederived by dynamical estimates, as they follow a 1:1 relation withthe virial masses of haloes in Shark, albeit with significant scat-ter ( ∼ . higher-membership groups, the dynamical massmethod provides a more reliable halo mass estimate and shows lessscatter than the abundance matching estimates. The difference of themedian halo mass estimates for high-membership groups betweenthe Shark intrinsic values and, those derived from dynamical andabundance-matching estimates is 0 . . • We find that by mimicking the Hi stacking procedure used by differ-ent surveys, we are able to reproduce all the observed HIHM relation.We also find that despite making changes to the Hi stacking projectedapertures and stacking windows, used by the surveys, we are unableto recover the intrinsic HIHM relation, which again points to theimportance of correct halo mass estimation. • We find that the group Hi stacking suffers from contamination from 𝑀 vir ≥ M (cid:12) , due to its reliance on large projected aperture andstacking velocity window. This contamination amounts up to 0 . 𝑁 g ≥ 𝑀 vir ≤ M (cid:12) , and is able to recover the total Hi mass associated withgroups in the halo mass range 𝑀 vir (cid:38) M (cid:12) . Individual Hi stackingshows minimal contamination and recovers the intrinsic total Hi massfor groups residing in 𝑀 vir ≈ − M (cid:12) . Due to the detection limitof GAMA, the Hi masses of groups thereafter are underestimated, asthe major Hi contributing galaxies lie below the detection limit. • We estimate that a deeper spectroscopic survey, such as WAVES,will be able to recover ∼ −
88 per cent of the total Hi mass ofthe haloes. This will lead to an improvement of 0 . − . 𝑀 vir ∼ . − M (cid:12) and reliable dynamical mass estimates (with 𝑁 g ≥ M (cid:12) . This improvement is likelysufficient to unveil the true shape of the HIHM relation at the criticalhalo mass range of 10 –10 . M (cid:12) , where differences are largestamong simulations and the effect of AGN feedback becomes mostprevalent. • We note that the Hi mass estimates from the intrinsic HIHM relationand the results of the mock Hi stacked HIHM relation for 𝑀 vir ≈ M (cid:12) are higher than their observed counterparts. This is likelya result of lack of modelling of ISM stripping in Shark satellites,which results in satellites being more Hi-rich in the 𝑀 vir ≥ M (cid:12) halo mass range than they need to be to recover observations.The current paucity of observational constraints on the shape,scatter and evolution of the HIHM relation is likely to change in thecoming decade, ultimately with the Square Kilometre Array (SKA;Abdalla & Rawlings 2005), but also with its pathfinders. With thecoming of Wide-Area VISTA Extragalactic Survey (WAVES; Driveret al. 2019), a new era of optical surveys will start, with surveysgoing as deep as an apparent 𝑧 -band magnitude of 21 . 𝑀 vir ∼ M (cid:12) , which we find in thiswork to be required to recover the true shape of the underlyingHIHM relation. The depth of these surveys will certainly lead toimprovements over the previous Hi and optical surveys; however,careful consideration of systematic effects such as those describedhere will be necessary to make measurements that can be robustlycompared with the simulation predictions. ACKNOWLEDGEMENTS
We would like to thank Hong Guo and Michael Jones for their con-structive comments, guidance and useful discussions. We also thankAaron Robotham, Rodrigo Tobar and Pascal Elahi for their contri-bution towards surfs and Shark, and Mark Boulton for his IT help.GC is funded by the MERAC Foundation, through the PostdoctoralResearch Award of CL, and the University of Western Australia.Parts of this research were carried out by the ARC Centre of Ex-cellence for All Sky Astrophysics in 3 Dimensions (ASTRO 3D),
MNRAS , 1–21 (2020) G. Chauhan et al. through project number CE170100013. CL is funded by ASTRO3D. ARHS acknowledges receipt of the Jim Buckee Fellowship atICRAR-UWA. DO is a recipient of an Australian Research CouncilFuture Fellowships (FT190100083) funded by the Australian Gov-ernment. MB acknowledges the support of the University of WesternAustralia through a Scholarship for International Research and AdHoc Postgraduate Scholarship. This work was supported by resourcesprovided by the Pawsey Supercomputing Centre with funding fromthe Australian Government and the Government of Western Australia.
DATA AVAILABILITY
The data that support the findings of this study are available uponrequest from the corresponding author, GC. The surfs simula-tions used in this work can be freely accessed from https://tinyurl.com/y4pvra87 (micro-surfs) and https://tinyurl.com/y6ql46d4 (medi-surfs).
REFERENCES
Abdalla F. B., Rawlings S., 2005, MNRAS, 360, 27Albareti F. D., et al., 2017, ApJS, 233, 25Amarantidis S., et al., 2019, MNRAS, 485, 2694Baugh C. M., 2006, Reports on Progress in Physics, 69Baugh C. M., et al., 2019, MNRAS, 483, 4922Behroozi P. S., Conroy C., Wechsler R. H., 2010, ApJ, 717, 379Bellstedt S., et al., 2020, MNRAS, 496, 3235Benson A. J., 2010, Physics ReportsBlitz L., Rosolowsky E., 2006, ApJBravo M., Lagos C. d. P., Robotham A. S. G., Bellstedt S., Obreschkow D.,2020, arXiv e-prints, p. arXiv:2003.11258Brown T., Catinella B., Cortese L., Kilborn V., Haynes M. P., Giovanelli R.,2015, MNRAS, 452, 2479Brown T., et al., 2017, MNRAS, 466, 1275Bruzual G., Charlot S., 2003, MNRAS, 344, 1000Cañas R., Elahi P. J., Welker C., del P Lagos C., Power C., Dubois Y., PichonC., 2019, MNRAS, 482, 2039Cañas R., Lagos C. d. P., Elahi P. J., Power C., Welker C., Dubois Y., PichonC., 2020, MNRAS, 494, 4314Campbell D., Van Den Bosch F. C., Hearin A., Padmanabhan N., Berlind A.,Mo H. J., Tinker J., Yang X., 2015, MNRAS, 452, 444Catinella B., et al., 2010, MNRAS, 403, 683Chabrier G., 2003, PASP, 115, 763Charlot S., Fall S. M., 2000, ApJ, 539, 718Chauhan G., Lagos C. d. P., Obreschkow D., Power C., Oman K., Elahi P. J.,2019, MNRAS, 488, 5898Chauhan G., Lagos C. d. P., Stevens A. R. H., Obreschkow D., Power C.,Meyer M., 2020, MNRAS,Chung A., van Gorkom J. H., Kenney J. D. P., Crowl H., Vollmer B., 2009,AJ, 138, 1741Colless M., et al., 2001, MNRAS, 328, 1039Conroy C., 2013, ARA&A, 51, 393Cooray A., Milosavljević M., 2005, ApJ, 627, L85Crain R. A., et al., 2015, MNRAS, 450, 1937Dale D. A., Helou G., Magdis G. E., Armus L., Díaz-Santos T., Shi Y., 2014,ApJ, 784, 83Davies L. J. M., et al., 2019, MNRAS, 483, 1881Dénes H., Kilborn V. A., Koribalski B. S., Wong O. I., 2016, MNRAS, 455,1294Driver S. P., et al., 2011, MNRAS, 413, 971Driver S. P., et al., 2019, The Messenger, 175, 46Dubinski J., 1998, ApJ, 502, 141Elahi P. J., Welker C., Power C., Lagos C. d. P., Robotham A. S. G., CañasR., Poulton R., 2018, MNRAS, 475, 5338 Elahi P. J., Cañas R., Poulton R. J. J., Tobar R. J., Willis J. S., Lagos C. d. P.,Power C., Robotham A. S. G., 2019a, Publ. Astron. Soc. Australia, 36,e021Elahi P. J., Poulton R. J. J., Tobar R. J., Cañas R., Lagos C. d. P., Power C.,Robotham A. S. G., 2019b, Publ. Astron. Soc. Australia, 36, e028Fabello S., Catinella B., Giovanelli R., Kauffmann G., Haynes M. P., HeckmanT. M., Schiminovich D., 2011, MNRAS, 411, 993Giovanelli R., et al., 2005, AJ, 130, 2598Gunn J. E., Gott III J. R., 1972, ApJ, 176, 1Guo H., Li C., Zheng Z., Mo H. J., Jing Y. P., Zu Y., Lim S. H., Xu H., 2017,ApJ, 846, 61Guo H., Jones M. G., Haynes M. P., Fu J., 2020, arXiv e-prints, p.arXiv:2004.04762Haynes M. P., et al., 2018, ApJ, 861, 49Johnston S., et al., 2008, Experimental Astronomy, 22, 151Lagos C. d. P., Davis T. A., Lacey C. G., Zwaan M. A., Baugh C. M.,Gonzalez-Perez V., Padilla N. D., 2014, MNRAS, 443, 1002Lagos C. d. P., Tobar R. J., Robotham A. S. G., Obreschkow D., MitchellP. D., Power C., Elahi P. J., 2018, MNRAS, 481, 3573Lagos C. d. P., et al., 2019, MNRAS, 489, 4196Lagos C. d. P., da Cunha E., Robotham A. S. G., Obreschkow D., ValentinoF., Fujimoto S., Magdis G. E., Tobar R., 2020, arXiv e-prints, p.arXiv:2007.09853Lim S. H., Mo H. J., Lu Y., Wang H., Yang X., 2017, MNRAS, 470, 2982Liske J., et al., 2015, MNRAS, 452, 2087Lu Y., et al., 2016, ApJ, 832, 39Marasco A., Crain R. A., Schaye J., Bahé Y. M., van der Hulst T., Theuns T.,Bower R. G., 2016, MNRAS, 461, 2630Meyer M., 2009, in Panoramic Radio Astronomy: Wide-field 1-2 GHz Re-search on Galaxy Evolution. p. 15 ( arXiv:0912.2167 )Moster B. P., Somerville R. S., Maulbetsch C., van den Bosch F. C., MacciòA. V., Naab T., Oser L., 2010, ApJ, 710, 903Obreschkow D., Klöckner H. R., Heywood I., Levrier F., Rawlings S., 2009,ApJ, 703, 1890Obuljen A., Alonso D., Villaescusa-Navarro F., Yoon I., Jones M., 2019,MNRAS, 486, 5124Padmanabhan H., Refregier A., 2017, MNRAS, 464, 4008Planck Collaboration et al., 2016, A&A, 594, A13Poulton R. J. J., Robotham A. S. G., Power C., Elahi P. J., 2018, Publ. Astron.Soc. Australia, 35Robotham A. S. G., et al., 2011, MNRAS, 416, 2640Robotham A. S. G., Davies L. J. M., Driver S. P., Koushan S., Taranu D. S.,Casura S., Liske J., 2018, MNRAS, 476, 3137Robotham A. S. G., Bellstedt S., Lagos C. d. P., Thorne J. E., Davies L. J.,Driver S. P., Bravo M., 2020, MNRAS, 495, 905Sargent M. T., et al., 2014, ApJ, 793, 19Schaye J., et al., 2015, MNRAS, 446, 521Shostak G. S., Allen R. J., 1980, A&A, 81, 167Skibba R. A., van den Bosch F. C., Yang X., More S., Mo H., Fontanot F.,2011, MNRAS, 410, 417Spinelli M., Zoldan A., Lucia G. D., Xie L., Viel M., 2019, The atomic Hy-drogen content of the post-reionization Universe ( arXiv:1909.02242 )Stevens A. R. H., Brown T., 2017, MNRAS, 471, 447Stevens A. R. H., et al., 2019a, MNRAS, 483, 5334Stevens A. R. H., Diemer B., Lagos C. d. P., Nelson D., Obreschkow D.,Wang J., Marinacci F., 2019b, MNRAS, 490, 96Trayford J. W., Lagos C. d. P., Robotham A. S. G., Obreschkow D., 2020,MNRAS, 491, 3937Vazdekis A., Koleva M., Ricciardelli E., Röck B., Falcón-Barroso J., 2016,MNRAS, 463, 3409Villaescusa-Navarro F., et al., 2018, ApJ, 866, 135Wang J., Koribalski B. S., Serra P., van der Hulst T., Roychowdhury S.,Kamphuis P., N. Chengalur J., 2016, MNRAS, 460, 2143Wechsler R. H., Tinker J. L., 2018, ARA&A, 56, 435Yang X., Mo H. J., van den Bosch F. C., Jing Y. P., 2005, MNRAS, 356, 1293York D. G., et al., 2000, AJ, 120, 1579van den Bosch F. C., Yang X., Mo H. J., Norberg P., 2005, MNRAS, 356,1233MNRAS , 1–21 (2020)
I-Stacking: Systematics at play This paper has been typeset from a TEX/L A TEX file prepared by the author. MNRAS000