Connection Between Galaxies and Hi in the Circumgalactic and Intergalactic Media:
Variation According to Galaxy Stellar Mass and Star-formation Activity

The standard picture of galaxy formation within the gravitational instability paradigm indicates that galaxy formation and evolution is closely linked to its surrounding gas called the circumgalactic medium (CGM) and intergalactic medium (IGM) (e.g., Rauch, 1998; Mo et al., 2010). The inflowing gas from the IGM provides the fuel for star formation in galaxies, and promotes the growth of galaxies and their central supermassive black hole (SMBH). As important as the inflow is the energetic feedback from massive stars and SMBHs which blows the gas away into the CGM and IGM. Therefore, determining the Mpc-scale distribution of gas as a function of time and space is quite important for understanding galaxy formation and evolution.

The connection between the CGM/IGM and galaxies has been studied using Ly$\alpha$ forest absorptions in quasar spectra. The most common method to clarify the CGM/IGM–galaxy connection is the cross-correlation analysis between Ly$\alpha$ forest absorption and galaxies (e.g., Adelberger et al., 2003, 2005; Chen et al., 2005; Ryan-Weber, 2006; Faucher-Giguère et al., 2008; Rakic et al., 2011, 2012; Rudie et al., 2012; Font-Ribera et al., 2013; Prochaska et al., 2013; Tejos et al., 2014; Bielby et al., 2017). Alternatively, comparisons of galaxy and neutral hydrogen (Hi) overdensities have also been discussed in the literature (Mukae et al., 2017, 2019; Mawatari et al., 2017). Specific high-density regions with abundant Hi gas and highly clustered galaxies have also been studied (e.g., Cai et al., 2016; Lee et al., 2016; Mawatari et al., 2017; Hayashino et al., 2019). All of the above studies have revealed that galaxy distribution correlates with the IGM up to tens of comoving Mpc scales.

Physical properties of Hi gas in the CGM or IGM have been studied in detail theoretically (e.g., Meiksin, 2009; Meiksin et al., 2014, 2017; Fumagalli et al., 2011; van de Voort & Schaye, 2012; van de Voort et al., 2012; Rahmati et al., 2015), aided by powerful cosmological hydrodynamic simulations such as EAGLE (Schaye et al., 2015), Illustris (Vogelsberger et al., 2014a, b; Genel et al., 2014; Sijacki et al., 2015), and IllustrisTNG (Villaescusa-Navarro et al., 2018; Nelson et al., 2019). Turner et al. (2017) presented the median Hi optical depth vs. line-of-sight or transverse distance around galaxies in the EAGLE simulation, and found that it is sensitive to dark matter halo mass. Sorini et al. (2018) compared the mean Ly$\alpha$ absorption profile as a function of the transverse distance from galaxies in both observations and simulations, and have shown a reasonable match between them beyond $2$ proper Mpc (pMpc), but significant differences at $0.02-2$ pMpc. The subsequent study by Sorini et al. (2020) has investigated the impact of feedback around mock quasars and found that the stellar feedback is the more dominant driver (rather than AGN, active galactic nuclei) to determine the average properties of the CGM.

The correlation between galaxies and IGM morphology has also been examined in the literature. Martizzi et al. (2019) have demonstrated that galaxies with lower stellar masses than the median are in voids and sheets of the IGM, whereas galaxies with higher stellar masses are more likely to be in filaments and knots of the IGM with higher gas densities. In addition, the correlation of mock Ly$\alpha$ forest absorption spectra with galaxy overdensity has been studied and compared with observations (e.g., Stark et al., 2015; Miller et al., 2019). Although these theoretical studies provide further evidence of a strong link between the CGM/IGM and galaxies, there are many aspects that are still unclear and our understanding is still insufficient.

Why do we need to study the distribution of Hi and its bias beyond the well-established linear theory of the $\Lambda$CDM model? One of the main reasons is that the bias between baryons, galaxies, halos, and dark matter is non-trivial, and it is of substantial interest for future projects of IGM tomography, such as the Subaru PFS, Euclid, DESI, and 21cm cosmology by SKA. Understanding nonlinear bias at intermediate and small scales is becoming increasingly more important for the BAO studies of precision cosmology. What we are trying to achieve in this paper is to dig deeper into this nonlinear bias between baryons (particularly Hi), galaxies, and dark matter, beyond the simple halo bias, and evaluate the impact of limited observational data on the cross-correlation function (CCF). Feedback effects from galaxies can alter the Hi distribution in the intergalactic space due to ionization, and it is crucial to use a cosmological hydrodynamic simulation that solves heating and cooling self-consistently with the impact of the UV background radiation field. As observations can only trace a certain phase of gas, it is unclear if the simple bias of linear theory is applicable to Hi. In fact, it has been shown that Hi has a complicated scale-dependent bias by many authors (e.g., Villaescusa-Navarro et al., 2018; Ando et al., 2019, 2020; Sinigaglia et al., 2020). It is thus crucial to examine the IGM Hi–galaxy connection more deeply, including the dependence on galaxy mass and galaxy population. Current observational data is still very limited to study the CCF, both in terms of the survey volume of Ly$\alpha$ forest absorptions and the number of galaxies with spectroscopic redshift (Momose et al., 2020; hereafter M21). In comparison, numerical simulations provide an alternative opportunity to study the IGM Hi–galaxy connection, and this paper presents such systematic investigations for the first time. Our results certainly include the effect of halo bias, but they also reflect the intricate nonlinear bias between Hi, galaxies, and dark matter, beyond the simple linear theory.

In order to unveil the gas distribution on Mpc-scales around galaxies, we systematically investigate it using numerical simulations. We aim to: 1) establish the methodology to evaluate the CGM/IGM–galaxy connection; 2) examine its dependency on galaxy properties; and 3) verify the cause of dependencies using GADGET3-Osaka cosmological hydrodynamic simulations (Shimizu et al., 2019). What is new in this paper is a systematic investigation of the IGM–galaxy connection depending on several galactic properties. We make subsamples based on galaxies’ stellar mass $M_{\star}$, halo mass $M_{\text{DH}}$, star-formation rate (SFR) and specific SFR (sSFR = SFR/$M_{\star}$), and compare the CCFs. Calculating the CCFs in three-dimensions taking various galaxy properties into account is also a novel aspect of this work. Studying the IGM–galaxy correlation with the same resolution as the latest observations and comparing the simulation results with them is clearly a worthwhile exercise to check the validity of the $\Lambda$CDM model because the application of linear theory to the IGM Hi–galaxy connection has not yet been fully investigated. In this study, we particularly focus on the statistical comparisons between simulations and observations. For a fair comparison, we use the same parameters as in our companion observational paper (M21). As we describe in more detail in M21, for observations, we use the CLAMATO (COSMOS Ly$\alpha$ Mapping And Tomography Observations) which is publicly-available Ly$\alpha$ forest 3D tomography data (Lee et al., 2014, 2016, 2018), and several other catalogs in the archives.

This paper is organized as follows. We introduce our numerical simulations in Section 2 and the methodology to examine the IGM–galaxy connection in Section 3. Results and discussion are presented in Sections 4 and 5, respectively. Finally, we give our summary in Section 6. In the appendix, we discuss the dependencies of our results on redshift width, redshift uncertainty, sample size, and cosmic variance. We note that “cosmic web” and “IGM” are used specifically for those traced by Hi gas unless otherwise specified in this paper. In addition, we mainly use $h^{-1}$ Mpc in comoving units in the following sections.

In this paper, we use the cosmological hydrodynamic simulations performed with GADGET3-Osaka (Shimizu et al., 2019), which is a modified version of the Tree-PM SPH code GADGET-3 (originally described in Springel, 2005). Some physical processes important for galaxy formation such as star formation, supernova (SN) feedback and chemical enrichment have been implemented and described in detail by Shimizu et al. (2019). Our simulations reproduce various observational results such as stellar mass function, SFR function, stellar-to-halo-mass ratio, and cosmic star formation history within observational uncertainties at $z\geq 2$ (Shimizu et al. 2020, in preparation).

Here, we briefly describe our simulations, which employ $N=2\times 512^{3}$ particles in a comoving volume of ($100\,h^{-1}$ ${\rm Mpc})^{3}$. The particle masses of dark matter and gas are $5.38\times 10^{8}\,{\rm M_{\odot}}$ and $1.00\times 10^{8}\,h^{-1}\,{\rm M_{\odot}}$, respectively. The gravitational softening length is set to be $8\,h^{-1}$ ${\rm kpc}$ in comoving units. For SPH particles, the smoothing length is allowed to be 10% of the softening length ($800\,h^{-1}$ ${\rm pc}$ in comoving units) in regions where baryons dominate over dark matter. This means that at $z=3$ the physical maximum resolution for gas may reach 200 pc (proper) in our simulations. In this study, we pay particular attention to the physical properties of gas at $r>100\,h^{-1}$ ${\rm kpc}$ from galaxies, and our resolution is sufficient for this study.

Star particles are generated from gas particles when a set of criteria are satisfied. Note that the mass of gas particles changes over time due to star formation and stellar feedback (by supernova and AGB stars). As described in detail in Shimizu et al. (2019), we employ the CELib package (Saitoh, 2017) which can treat the time and metallicity-dependent metal yields and energies from SNII, SNIa and AGB, allowing for more realistic chemical evolution of simulated galaxies. In order to identify the simulated galaxies, we run a friends-of-friends (FoF) group finder with a comoving linking length of $0.2$ in units of the mean particle separation to identify groups of dark matter particles as dark matter halos. We then identify gravitationally-bound groups of minimum 32 particles (dark matter + SPH + star) as substructures (subhalos) in each FoF group using the SUBFIND algorithm (Springel et al., 2001), which is a standard practice in many of the simulation works. We regard substructures that contain at least five star particles as our simulated galaxies. Moreover, we define the most massive galaxy in a halo as the central galaxy, and the rest as satellite galaxies. We also calculate the virial halo mass $(M_{\rm DH})$, which is defined by the total enclosed mass inside a sphere of 200 times the critical density of the Universe. This means that the member galaxies (central and satellite galaxies) in a dark matter halo have the same $M_{\rm DH}$, even though the substructures can have different subhalo masses calculated by SUBFIND.

Note that each gas (star) particle has some associated physical properties such as mass, SFR and metallicity. In our SPH simulation, these values are automatically computed following hydrodynamic and gravitational interactions, rather than given by hand, which is a fundamental difference from the semi-analytic models of galaxy formation. The properties (gas mass, stellar mass $M_{\star}$, and SFR) of a simulated galaxy are defined by summation of these quantities in each subhalo.

In order to directly compare our simulations and observations, we create the light-cone output of gas particles and galaxies by connecting $10$ simulation boxes of different redshifts following our previous work (Shimizu et al., 2012, 2014, 2016). The redshift range of our light-cone output is from $z\sim 1.8$ to $3.1$ which can cover the redshift range of recent Ly$\alpha$ absorption line surveys (e.g., CLAMATO; Lee et al., 2014, 2016, 2018) and future PFS Ly$\alpha$ absorption survey (Takada et al., 2014). We then randomly shift and rotate each simulation box so that the same objects do not appear multiple times on a single line-of-sight (LoS) at different epochs.

With this light-cone output, we calculate the Ly$\alpha$ optical depth ($\tau_{\rm Ly\alpha}$) along the LoS. First, we calculate the important physical quantities, $A_{\rm grid}(x)$, at each grid point $x$ along LoS, such as Hi density, LoS velocity and temperature as follows:

where $A_{j}$, $m_{j}$, $\rho_{j}$ and $h_{j}$ are the physical quantity of concern, gas particle mass, gas density, and smoothing length of $j$-th particle, respectively. $W$ is the SPH kernel function, and $r$ is the distance between LoS grid points and gas particles. For simplicity, the grid size ($dl$) is set to a constant value of $100\,h^{-1}$ kpc in comoving units which corresponds to 10 ${\rm km\,s^{-1}}$ in the velocity space at $z\simeq 2$. This is a higher resolution than any of the relevant Ly$\alpha$ observations. Then, we calculate the Ly$\alpha$ optical depth $\tau_{\rm Ly\alpha}(x)$ using these physical values at each grid point as follows:

where $e$, $m_{e}$, $c$, $f_{ij}$, $n_{\rm HI}$, and $x_{j}$ are the electron charge, electron mass, speed of light, absorption oscillator strength, Hi number density, and $j$-th grid point location, respectively. $\phi$ is the Voigt profile, and we use the fitting formula of Tasitsiomi (2006) without direct integration. In this study, after making the high resolution LoS data, we reduce our resolution by coarse-graining the grid size to match the observations. 1024 $(=32^{2})$ LoSs are drawn with regularly spaced intervals. The mean separation of each LoS is 3.3 $h^{-1}$ Mpc which is similar to the CLAMATO survey.

Finally, we note that the AGN feedback is not considered in the simulation that we use for this paper. We note that the impact of AGN feedback on Hi distribution is still under debate. For example, Sorini et al. (2020) have found little impact of AGN feedback on the surrounding CGM in general, and that SN feedback has more dominant effects. Besides, statistical results of the Hi–galaxy connection obtained from the same simulations as this study have shown good agreement with current observations (Nagamine et al., 2020). Some observations have found the proximity effect around quasars (e.g., Mukae et al., 2019; Momose et al., 2020), which requires full radiative transfer calculations to compare with simulations. Therefore, we defer the study of AGN feedback using radiation transfer to our future work.

One of the main purposes of this study is to compare with the CLAMATO, which is a 3D tomography data of Ly$\alpha$ forest transmission fluctuation ($\delta_{\text{F}}$: see the following definition) over $2.05<z<2.55$ in $0.157$ deg² of the COSMOS field (Scoville et al., 2007; Lee et al., 2016, 2018). The CLAMATO consists of $60\times 48\times 876$ pixels corresponding to $30\times 24\times 436$ $h^{-1}$ Mpc cubic with a pixel size of $0.5$ $h^{-1}$ Mpc. The average separation of background galaxies for measuring Ly$\alpha$ absorption are [$2.61$, $3.18$] $h^{-1}$ Mpc in [RA, DEC] directions, and the separation in LoS-direction is $2.35h^{-1}$ Mpc at $z\sim 2.3$. (See M21 for more detailed comparison with CLAMATO.)

To produce a similar data cube of $\delta_{\text{F}}$ from our LoS data, we first evaluate the Ly$\alpha$ transmission fluctuation $\delta_{\text{F}}$ in each LoS pixel by

where $F(x)=\exp{[-\tau_{\rm Ly\alpha}(x)]}$ is the Ly$\alpha$ flux transmission, and $\langle F_{z}(x)\rangle$ is the cosmic mean transmission. We adopt the following value derived by Faucher-Giguère et al. (2008):

because it is used in CLAMATO. Additionally, we also use the same setup used in CLAMATO with Hubble constant $h=0.7$ and redshift coverage of $2.05\leq\,z\,\leq 2.55$.

In the following two subsections, we present the methods for two analyses: (I) cross-correlation, and (II) overdensity analysis.

We systematically investigate the connection between galaxies and CGM/IGM, particularly traced by Ly$\alpha$ forest absorption. In this study, we use cosmological hydrodynamic simulation (Shimizu et al., 2019; Nagamine et al., 2020) and demonstrate the CGM/IGM–galaxy connection using two methods: one is the cross-correlation analysis, and the other is the overdensity analysis proposed by Mukae et al. (2017). Using our simulation, we also calculate CCFs of relative gas density (both total and Hi) around galaxies. All parameters for our analyses are chosen to match the observations presented in M21. The main results of this paper are summarized below.

1. 1.

   We calculate CCFs between Ly$\alpha$ forest transmission fluctuation ($\delta_{\text{F}}$) and galaxies as shown in Figure 1. The CCF obtained from all galaxies reproduce the one from LBGs in the literature (Adelberger et al., 2005; Bielby et al., 2017). Further investigations based on subsamples divided by $M_{\star}$, $M_{\text{DH}}$, SFR, and sSFR of simulated galaxies show following trends and variations in the CCF. For the $M_{\star}$ and $M_{\text{DH}}$ subsamples, we find a clear trend that the CCF signal becomes stronger with increasing galaxy masses. We also confirm that the turnover radius of CCFs becomes smaller with increasing mass (see Figure 1–2), indicating stronger clustering of massive galaxies around density peaks. Additionally, from the best-fit parameters of power-law fitting for the CCFs, clustering length $r_{0}$ and slope $\gamma$ are found to become longer and steeper with increasing galaxy masses in the IGM regime ($r\geq 1$ $h^{-1}$ Mpc: see also Figure 2). For the SFR samples, we find that they tend to have stronger signals with increasing SFR. Such a trend for SFR samples is also linked to the mass dependence of CCF, because $M_{\star}$ and SFR is almost linearly related with each other through star-formation main sequence of galaxies. However, we do not identify clear trends in the CCF of sSFR samples.

2. 2.

   We measure CCFs between gas density fluctuation ($\delta_{\rho_{\text{gas}}}$ and $\delta_{\rho_{\text{{\sc Hi}}}}$) and galaxies in Figures 3 and 4. Overall trends of CCFs are similar to those of CCFs in $\delta_{\text{F}}$ except for SFR–(v) and sSFR–(iii). It indicates that the variation in $\delta_{\text{F}}$ CCF reflects different relative gas densities around galaxies; i.e., galaxies with higher mass and SFR generally reside in higher density gas, and vice versa. For the SFR–(v) and sSFR–(iii) subsamples, we find the highest CCF signal at $r=0.3-2$ $h^{-1}$ Mpc among all SFR and sSFR samples. Because the two subsamples have a mild bimodal halo mass distribution with two peaks at $M_{\text{DH}}\sim 10^{11.3}$ M_⊙ and $M_{\text{DH}}\sim 10^{12.5}$ M_⊙, their highest CCF signals are probably due to high-density regions where massive host halos reside. We suggest that the observed variation in the CCF is caused by the dependence of gas density (both total and Hi) around galaxies.

3. 3.

   Our overdensity analysis between galaxy overdensity $\Sigma_{\text{gal}}$ and mean IGM fluctuation $\langle\delta_{\text{F}}\rangle$ is presented in Figure 5. We statistically identify anti-correlations from all subsamples of $M_{\star}$–11, $M_{\star}$–10, $M_{\star}$–9 and ALL. In addition, we also find that their slopes are decreasing with increasing $M_{\star}$, although within the error. It suggests that galaxies in the $M_{\star}$–11 subsample are more strongly correlated with higher density gas than those in $M_{\star}$–9 in terms of their spatial distribution.

4. 4.

   Considering all of our results together, we conclude that the mass, particularly the dark matter halo mass, is the most sensitive parameter to determine the Mpc-scale gas-density environment around galaxies. Galaxies in massive halos tend to be clustered in higher density regions of the cosmic web, resulting in a CCF with a higher amplitude, greater $r_{0}$, steeper $\gamma$, and shallower anti-correlation between $\langle\delta_{\text{F}}\rangle$ and $\Sigma_{\text{gal}}$ at $r\geq 1$ $h^{-1}$ Mpc.

Overall, our analyses confirm the strong connection between galaxies, dark matter halos, and IGM, providing further support for the gravitational instability paradigm of galaxy formation within the concordance $\Lambda$CDM model. Future observations of CCF studies between galaxies, Hi, and metals will provide useful information on the interaction between them and the details of feedback mechanisms which is important for the theory of galaxy formation and evolution such as galactic wind and associated ejection of metals into IGM.

By comparing the results of Figures 1 & 2 against predictions of linear perturbation theory, we can infer the mean bias parameters of Ly$\alpha$ forest and galaxies relative to the underlying dark matter density field (e.g., Croft et al., 2016; Bielby et al., 2017; Kakiichi et al., 2018; Meyer et al., 2019). In addition we can perform cross-checks by computing such bias parameters directly from our simulation output, and compare with those inferred from linear theory framework. Such bias parameters will further provide additional checks against the gravitational instability paradigm, and we plan to carry out such analyses in our further work.