MUSEQuBES: Mapping the distribution of neutral hydrogen around low-redshift galaxies

We present a detailed study of cool, neutral gas traced by Lya around 4595 z<0.5 galaxies using stacks of background quasar spectra. The galaxies are selected from our MUSEQuBES low-z survey along with data from the literature. These galaxies, with a median stellar mass of log (M*/Msun)= 10.0, are probed by 184 background quasars giving rise to 5054 quasar-galaxy pairs. The median impact parameter is b = 1.5 pMpc (median b/Rvir=10.4) with 204 (419) quasar-galaxy pairs probing b/Rvir<1 (2). We find excess absorption out to at least ~ 15 Rvir transverse distance and ~ 600 km/s along the line of sight. We show that the median stacked profile for the full sample, dominated by the pairs with b>Rvir, can be explained by a galaxy-absorber two-point correlation function with r0 = 7.6 pMpc and gamma = -1.57. There are strong indications that the inner regions (<Rvir) of the rest equivalent width profile are better explained by a log-linear (or a Gaussian) relation whereas the outer regions are well described by a power-law, consistent with galaxy-absorber large-scale clustering. Using a sub-sample of 339 galaxies (442 quasar-galaxy pairs, median b/Rvir = 1.6) with star formation rate measurements, we find that the Lya absorption is significantly stronger for star-forming galaxies compared to passive galaxies, but only within the virial radius. The Lya absorption at b ~ Rvir for a redshift-controlled sample peaks at M* ~ 10^9 Msun~ (Mhalo ~ 10^11 Msun).


INTRODUCTION
It is now well established that the luminous parts of galaxies are surrounded by a large reservoir of low-density, diffuse gas called the circumgalactic medium (CGM).The CGM serves as a bridge between galaxies and the intergalactic medium (IGM).However, there is no consensus on the actual extent of this medium surrounding galaxies.
Absorption line spectroscopy of bright background sources such as quasars paved the way to study these elusive media, which are otherwise difficult to detect in emission.The advancement of sensitive space-based detectors such as the Cosmic Origin Spectrograph (COS) onboard the Hubble Space Telescope () has drastically increased our ability to study the diffuse gas around galaxies over the last two decades.Several surveys have established the connections between the CGM and many key galaxy properties using absorption line spectroscopy of background quasars (e.g., Tumlinson et al. 2017;Péroux & Howk 2020).Alongside the rapidly accumulating data from galaxy surveys to probe the CGM, the advent of new ★ E-mail: sayak18@iucaa.inphysical prescriptions for gas flows with new numerical methods and faster computers continues to pose fundamental questions about galaxy formation and evolution.Theoretical studies indeed suggest that the distribution of gas and metals surrounding galaxies is intricately linked to processes such as gas accretion, galactic winds, mergers, and stripping that dictate galactic evolution (e.g., Rahmati et al. 2016;van de Voort et al. 2019;Oppenheimer et al. 2020;Appleby et al. 2022;Mitchell & Schaye 2022).It is thus crucial to map the distributions of diffuse matter surrounding galaxies with a range of galaxy properties (e.g., stellar mass ( * ), star formation rate (SFR), specific-SFR (sSFR)) to gain insight into such physical processes.
The semi-empirical relation between the stellar mass to halo mass ratio ( * / halo ) as a function of  halo shows a peak at  halo ≈ 10 12 M ⊙ , suggesting that these halos are most efficient in converting baryons into stars (e.g., Conroy et al. 2006;Behroozi et al. 2013Behroozi et al. , 2019)).It is not well understood why halos of a certain mass have higher star formation efficiencies (SFE).Energy and momentum feedback due to supernovae (SNe) and active galactic nuclei (AGN) is thought to suppress the SFE for the low-mass and highmass halos, respectively (see e.g., Bower et al. 2006;Somerville et al. 2008;Crain et al. 2015) .Theoretical models suggest two distinct modes of gas accretion depending on the halo mass: (a) "hot mode" for galaxies with  halo > 10 12 M ⊙ and (b) "cold mode" for galaxies with  halo < 10 12 M ⊙ (see e.g., Kereš et al. 2005;van de Voort et al. 2011).The halos of masses ≈ 10 12 M ⊙ benefits from both the modes.Thus, the higher SFE for these halos may be related to the gas accretion process.It is therefore essential to probe how the cool, neutral gas reservoir surrounding galaxies changes as a function of stellar/halo mass.
The connection between galaxies and the cool, neutral gas surrounding them has long been a subject of research, both observational and theoretical.Chen et al. (2005) first observed a morphology dependent galaxy-Ly absorber cross-correlation function.They found that the cross-correlation function between Ly absorbers (with column density  (H i) > 10 14 cm −2 ) and absorption-line dominated galaxies within a projected distance of 1ℎ −1 Mpc is significantly lower compared to emission-line dominated galaxies (see also Chen & Mulchaey 2009).On the other hand, later CGM surveys such as COS-Halos did not find any significant difference in the H i content (rest-frame equivalent width (REW or   ) and column density ()) around star-forming and passive galaxies on CGM scales (typically within 0.5 vir ; see e.g., Thom et al. 2012;Tumlinson et al. 2013).Although a small fraction of the passive galaxies in their sample did not exhibit any detectable Ly absorption, overall the detected Ly absorption was similar around the star-forming and passive galaxies.
Earlier, Chen et al. (1998Chen et al. ( , 2001) ) noted that the   -profile of Ly absorption around galaxies depends on the -band and -band luminosities (proxies for the  * and recent SFR, respectively).Combining the observational data from the COS-Halos and COS-Dwarf (Bordoloi et al. 2014) surveys, Borthakur et al. (2016) found a strong anti-correlation between Ly   and sSFR.Additionally, they reported similar exponential scale lengths of the Ly   -profile for passive and star-forming galaxies.The larger dispersion in the radial distribution seen for passive galaxies has been interpreted as an indication of patchiness in their CGM.A positive correlation between Ly   and stellar mass was reported by Bordoloi et al. (2018).More recently, Wilde et al. (2021) presented the Ly covering fraction-and column density-profiles for three different bins of stellar mass.The covering fraction profiles were found to be very different for the high-mass ( * > 10 9.9 M ⊙ ) and low-mass ( * < 10 9.2 M ⊙ ) galaxies in their sample.The Ly covering fraction is as high as ≈ 90% within  vir for the high-mass subsample, but gradually declines to only ≈ 20% at / vir ≈ 6.The covering fraction for the low-mass sub-sample, on the other hand, remains roughly constant at ≈ 40% all the way out to ≈ 6 vir .
Although the CGM is thought to play a crucial role in driving galaxy evolution, there is no consensus on the extent of the CGM around galaxies.It is generally thought that the CGM extends out to the virial radius of a galaxy (Tumlinson et al. 2017).However, Wilde et al. (2021) recently argued that the CGM traced by neutral hydrogen extends beyond the virial radius (see also Wilde et al. 2023).In the study by Prochaska et al. (2011), the detection rate and   of Ly absorbers are found to be correlated with galaxies out to at least 1 pMpc, far beyond the virialized halos of the galaxies.They, however, suggested that the weak Ly absorbers arising at large impact parameters () may be unrelated to the gaseous halos around galaxies, and may be tracing the large-scale environments in which the galaxies are embedded.Tejos et al. (2014) studied the galaxy-Ly absorption correlation function and found that it is significantly different from the galaxy-galaxy autocorrelation function, and the difference is primarily driven by 'weak' H i absorbers ( (H i ) < 10 14 cm −2 ).They concluded that > 50% of weak absorbers are not correlated with galaxies, and hence, galaxies and these absorbers may not trace the same underlying dark matter distribution.Wakker et al. (2015) used the Ly absorption to study nearby galaxy filaments using /COS spectra of 24 background AGN.They observed a trend of increasing Ly equivalent width and line width with decreasing filament impact parameter.The Ly absorption detection rate is ≈ 80% within 500 pkpc (proper kpc) of galaxy filaments, but no absorption is seen at ≳ 2 pMpc.A study with a statistically significant number of quasargalaxy pairs with a wide range of impact parameters, such as the one presented here, is essential to probe the inner and outer regions of galaxy halos simultaneously.This, in turn, allows one to investigate the extent of the so-called CGM.
In order to build a statistically significant sample of background quasar-foreground galaxy pairs, we compiled data from several low- CGM surveys in the literature along with our own data obtained from the MUSE Quasar-fields Blind Emitters Survey (MUSEQuBES).MUSEQuBES is a dual MUSE program with 16 quasar fields to study the CGM of low- galaxies and 8 quasar fields to study the gaseous environments of high- Ly emitters (see Muzahid et al. 2020Muzahid et al. , 2021)).Integral field spectroscopy (IFS) with MUSE allows us to search for galaxies, particularly the continuum-faint ones, around background quasars more efficiently than with multi-object spectroscopy (MOS) and conventional long-slit spectroscopy.On the flip-side, MUSE has a relatively small (1 ′ × 1 ′ ) field of view (FoV) compared to a typical MOS spectrograph.Consequently, the probability of finding massive galaxies in the MUSE FoV is small, as they are rare.Therefore, combining IFS and MOS/long-slit observations of quasar fields provides an optimal way to probe the CGM and large-scale structures around galaxies with a wide range of stellar masses simultaneously.
Although MUSE can detect pure line emitters without any detectable continuum, we restricted our analysis to the continuumdetected galaxies in order to be able to estimate their stellar (halo) mass.Our deep observations with MUSE (2-10 hours of exposure time per field) enable us to obtain a galaxy sample of relatively low mass (median  * ∼ 10 8 M ⊙ ) compared to the existing CGM surveys in the literature.The MUSEQuBES galaxies along with the archival galaxy samples allow us to explore the uncharted territory of the CGM around low-mass, intermediate redshift galaxies (see orange points in Fig. 1).
Instead of the usual practice of identifying individual absorption features in quasar spectra and associating them with foreground galaxies, we used the spectral stacking method to map the H i gas, traced by the Ly absorption, in the CGM.The stack of Ly absorption at the rest-frame of galaxies provides statistical inference on the mean or median H i absorption around the galaxies without any prior knowledge regarding individual absorption systems.Although the presence of saturated absorbers complicates the possible inference of physical quantities such as column density, the strength of the absorption can be easily determined from the   (see e.g., Steidel et al. 2010).Recently, Ho et al. (2021) showed that the common approach of galaxy-absorber association based on LOS velocity cuts (e.g., ±300 km s −1 or ±500 km s −1 ) suffers from projection effects, which are more prominent at larger impact parameters.Spectral stacking enables us to be agnostic about the individual galaxy-absorber associations.Further, stacking a large galaxy sample can significantly improve the spectral signal-to-noise ratio (S/N).Stacking, however, erases the kinematic information of individual absorbers, which provides important insights on gas flow processes in galaxies (e.g., Bouché et al. 2013;Muzahid et al. 2015).
The connections between neutral gas and galaxies have been studied for high- galaxies as well (e.g., Steidel et al. 2010;Muzahid et al. 2021;Lofthouse et al. 2023).Rakic et al. (2012) produced 2D Ly optical depth maps as a function of impact parameter and LOS velocity separation for  ≈ 2.3 star-forming galaxies.Such maps show a clear evidence for redshift-space distortion along the LOS direction which cannot be fully attributed to redshift errors (Turner et al. 2014).Moreover, comparing with the eagle simulation (Schaye et al. 2015), Turner et al. (2017) found that infalling gas can account for the redshift-space distortion observed in the optical depth maps.Kinematic information in such maps are a useful tool to understand gas flow processes in galaxies (e.g., Chen et al. 2020).Here we present Ly optical depth maps for the first time for low- galaxies.
This paper is organized as follows.Section 2 summarizes the data sample used in this work.In Section 3, we provide the absorption data analysis.Section 4 presents the results of this work.Section 5 presents a discussion of the main results of our analysis followed by a summary in Section 6.Throughout this paper, we adopt a ΛCDM cosmology with Ω M = 0.3, Ω Λ = 0.7, and  0 = 70 km s −1 Mpc −1 .

Galaxy sample from the MUSEQuBES survey
The low- part of the MUSEQuBES survey targeted 16  ≈ 0.5 − 1.5 UV-bright quasar fields using VLT/MUSE.The MUSE observations were conducted between September 2014 and April 2017 (ESO programmes 094.A-0131, 095.A-0200, 096.A-0222, 097.A-0089 and 099.A-0159; PI: Schaye), with a total exposure time of 62.75 h.The quasars were selected solely based on the availability of high / FUV spectra obtained with /COS. and MUSE observation details of these 16 quasar fields are tabulated in Table 1.Each MUSE observation block of 1 h was split into 4 × 900 s exposures, which were rotated by 90 • and offset by a small ≈ 1 − 5 ′′ shift from each other.The data reduction is performed using the standard MUSE data reduction pipeline (v1.2;Weilbacher et al. 2020), adopting the default (recommended) set of parameters.A few additional reduction procedures using the CubEx package (Cantalupo et al. 2018) were carried out to improve flat-fielding and sky-subtraction using the CubeFix and CubeSharp routines, respectively.The procedure is detailed in Borisova et al. (2016), and more recently in Muzahid et al. (2021).
The effective seeing per field, corresponding to the full width at half maximum (FWHM) of a 2D Gaussian profile fitted to a point source at  = 7000 Å in the reduced and combined data cube, varies between 0.54 ′′ and 1.20 ′′ , but is typically 0.7 ′′ − 0.8 ′′ (see Table 1).With the MUSE field-of-view (FoV) of 1 ′ ×1 ′ centered on the quasar, we are able to observe a region of 480 pkpc × 480 pkpc around the QSO at  ≈ 1 (110 pkpc × 110 pkpc at  ≈ 0.1).The field is spatially sampled by a grid of 0.2 ′′ ×0.2 ′′ pixels.All MUSE observations were carried out using the standard wavelength range of 4750-9350 Å, sampled by 1.25 Å spectral pixels.The resolving power ranges from  ≈ 1800 at  = 5000 Å to  ≈ 3500 at  = 9000 Å, corresponding to a FWHM of 167 km s −1 to 86 km s −1 , respectively.
The details of galaxy identifications and galaxy property measurements will be presented in a future work.Here, we briefly outline the main steps.First, we run the Source Extractor (SExtractor; Bertin, E. & Arnouts, S. 1996) on the MUSE white-light images using a detection threshold of 1 per pixel (DETECT_TRESH = 1) and requiring a minimum number of neighbouring pixels above the threshold of 3 (DETECT_MINAREA = 3).The 1D spectra of the continuum-detected objects, extracted from the MUSE cubes using the SExtractor-generated segmentation maps, are then inspected by a modified version of MARZ (Hinton et al. 2016) to determine their redshifts based on the spectral features.The redshifts are further refined using a modified version of the code PLATEFIT (Brinchmann et al. 2004) by fitting Gaussian profiles to the available emission and absorption line features.Note that the wavelengths in the MUSE data cubes are given in air.We applied appropriate corrections while determining the galaxy redshifts.
The H or [O ii] (when H is not covered or not detected at > 3) line fluxes returned by PLATEFIT are used to estimate the SFRs of the galaxies using the relation from Kennicutt (1998) or from Kewley et al. (2004, for O ii luminosities) adjusted for the Chabrier (2003) initial mass function (IMF).The H emission-line flux is corrected for dust extinction using the flux ratio of H and H lines.By comparing H/H to its intrinsic value of 2.85, corresponding to Case B recombination at a temperature of  ∼ 10 4 K and electron densities of  e ∼ 10 2 − 10 4 cm −3 (Osterbrock & Ferland 2006), we derive a correction for the H flux, assuming a Cardelli et al. (1989) reddening curve.For galaxies with H coverage beyond the MUSE spectra, we use H to calculate the SFR, under the condition that we can correct the line flux for dust extinction using the H/H ratio.We require that both H and H be detected with / > 3. We then convert the corrected H flux into the H flux to obtain the SFR, making use of the known intrinsic ratio between the H and H fluxes.Dust correction is not performed for galaxies with O ii based SFR measurements.
The stellar masses of galaxies are estimated using stellar population synthesis (SPS) code FAST (Kriek et al. 2009), which fits SPS templates to a set of photometric flux values.Owing to the lack of ancillary photometric data of the quasar fields in our sample, we constructed 11 pseudo-filters with a width of 400 Å each and spanning the wavelength range from 4800-9200 Å, after masking the prominent emission lines.We calculated the filter flux by convolving the 1D galaxy spectrum-the same one as used for the redshift determination with MARZ-with a boxcar function centered at  = 5000, 5400, . . .Å.The halo mass and virial radius, defined as the mass and radius of a spherical region within which the mean mass density is 200 times the critical density of the universe, are estimated from stellar mass and redshift measurements of galaxies using the abundance matching relation from Moster et al. (2013) In total, we have 475 galaxies detected in the 16 MUSE fields.Note that these are galaxies detected 3000 km s −1 blueward of the corresponding quasar redshifts.The galaxies have redshifts ranging from 0.05 − 1.4, and stellar masses and SFRs ranging from 10 6.0 to 10 11.8 M ⊙ and 10 −3 to 10 2.7 M ⊙ yr −1 with the median values of 10 8.9 M ⊙ and 10 −0.7 M ⊙ yr −1 , respectively.The galaxies have impact parameters from the corresponding background quasars in the range 10-320 pkpc with a median value of 150 pkpc (median / vir = 1.7).

Galaxy sample from the literature
Besides the MUSEQuBES galaxies, we combined galaxy samples from six different CGM surveys from the literature; namely COS-Halos (Tumlinson et al. 2013), Liang & Chen (2014), COS-Dwarf (Bordoloi et al. 2014), COS-Gass (Borthakur et al. 2015), Johnson et al. (2015) (hereafter Johnson+15), and Keeney et al. (2018) (hereafter Keeney+18).Brief summaries of these surveys are presented in Appendix A1-A6.Before merging the galaxy catalogs, we confirmed that a given galaxy observed in different surveys was not counted multiple times.To eliminate repetition, we ensured that there are no two galaxies within 1 ′′ spatially and within 500 km s −1 along LOS distance with each other.In cases of multiple occurrences, we count them only once (60 such cases).This provided us with a large sample Note-a This is the S/N at  = 1350Å.Due to a Lyman limit system at  = 0.390, there is no flux at  < 1280 Å.
of ∼9000 galaxies.A significant fraction of galaxies in this sample come from Keeney+18.The redshift and stellar mass distributions of the galaxies in the combined sample are shown in Fig. A1.

Quasar spectra and continuum fitting
Galaxies from the MUSEQuBES survey were detected in 16 fields centered on 16 UV-bright quasars.As mentioned earlier, these quasars were chosen solely based on their FUV brightness and the availability of high / COS spectra.The COS spectra of the quasars that are part of the six studies from the literature are also available in the  public archive.In total, we obtained 190 COS spectra in reduced form from the  Spectroscopic Legacy Archive (HSLA; Peeples et al. 2017).COS has a resolving power of  ≈ 18, 000 (FWHM ≈ 18 km s −1 ).Among the 190 quasars, 135 are observed with both the G130M and G160M gratings covering 1150-1800 Å.The remaining 55 quasars are observed with the G130M grating only, which has a spectral coverage of 1150-1450 Å.The G130M and G160M grating spectra of a given quasar are spliced together at the wavelength where the average / per pixel becomes equal (typically around 1424 Å).The quasar redshifts range from 0.03-1.88 with a median of 0.43.
These 190 quasar spectra have a wide range of /, with some quasars having / as low (high) as 1-2 (85-90) per resolution element.The median / per resolution element in the G130M for the 190 spectra is ≈ 10, computed around 1250 Å.The median / for the G160M grating is ≈ 7, computed near 1650 Å.Note that we did not impose any / cut while selecting the quasars.This is because only 19 (39) spectra have / < 5 in the G130M (G160M) grating, covering redshifted Ly absorption for only 72 (25) galaxies.We verified that excluding these low-/ spectra does not have a significant effect on our results.
We performed continuum fitting for all of these 190 spectra using a custom-made, semi-automated python routine which finds the absorption-free regions in a spectrum in an automated way.First, each spectrum is divided into several chunks and fluxes are sigma clipped with an upper bound of 50 (to take into account the emission lines as well).The lower bounds of sigma clipping are varied between 2.5, 2.1, 1.6, and 1.3  for regions with median / (per pixel) of > 10, 3−10, 2−3, and < 2, respectively.These numbers were chosen after many trials on randomly selected spectra from our sample and by visually inspecting the quality of the continuum model.Next, we increased the chunk size (or decreased the number of chunks) in iterations and repeated the clippings, until the number of clipped pixels converges, or the number of chunks is reduced to 2, whichever happens first.Then, we used spline fitting for the absorption-free pixels to estimate the continuum.The number of 'knots' around the known emission lines of Ly C iv N v and O vi are increased (if present in the spectra) to take care of larger curvatures.

Quasar-galaxy pairs
Using the 190 quasars in the background of ≈ 9000 galaxies we constructed quasar-galaxy pairs using the following two conditions: (1) The projected distance of a galaxy from the corresponding quasar sightline,  < 3000 pkpc.(2) The line of sight (LOS) separation of a galaxy from the background quasar is > 2000 km s −1 .The large impact parameter cut-off of 3000 pkpc was chosen to sample the large-scale environments of the galaxies.The limiting LOS separation of 2000 km s −1 from the quasar was chosen to avoid the quasar's proximity zones and associated absorbers arising from the quasars (e.g., Muzahid et al. 2013).Using these criteria, we obtained a total of 6020 quasar-galaxy pairs with 5310 unique galaxies.Some of the galaxies in this work are probed by multiple quasars.For example, the COS-GASS and COS-Dwarf surveys are restricted to impact parameters of 250 pkpc and 150 pkpc, respectively.However, a fraction of these galaxies is also probed by quasars at a transverse distance of ≈ 1 pMpc.These provide us with more skewers at larger impact parameters, and hence are included in our studies (see blue squares and red crosses at ≈ 1 pMpc in top-left panel of Fig. 1).We obtained the redshift, stellar mass, virial radius, and SFR (when available) for all the galaxies from the literature, except for the MUSEQuBES galaxies.There are only 634 galaxies, constituting 882 quasar-galaxy pairs, for which SFRs are known (including upper limits).Boogaard et al. (2018) (at  = 0.1) and its 1 spread.The data points in all of these panels are color-coded by the survey from which a given galaxy (or a quasar-galaxy pair) is drawn.Noted that Keeney et al. (2018) and Johnson et al. (2015) did not report SFRs.It is worthwhile to note that the MUSEQuBES data points (orange circles) make up for the uncharted territory of the parameter space covered by the previous CGM surveys (low-mass and higher redshift) and increase the low-mass galaxy sample size significantly.Note that the redshifted Ly absorption cannot be observed with COS for all of these galaxies as the spectral coverage of COS allows a maximum observable redshift of 0.48 (0.19) for Ly for the G160M (G130M) grating.Below, we describe our scheme to obtain the galaxies contributing to the Ly stacks.
First, we masked the following regions of the COS spectra: • 1212-1220 Å -to exclude the geocoronal Ly emission and Galactic Ly absorption.
• ±0.5 Å around known strong Galactic absorption lines (e.g., N i, Si ii).At this point, we note that while ±0.5 Å does not exclude all high-velocity clouds (HVCs), we verified that a broader mask of ±1.0 Å does not change any of our conclusions.
The only galaxies for which the redshifted Ly wavelengths (±0.05 Å) fall within the spectral coverage of the corresponding quasar spectra are considered here.The above conditions led to a total of 5054 galaxy-quasar pairs, with 184 background quasars probing 4595 foreground galaxies.The subsample of galaxies with measured SFR (including upper limits) is reduced to 339 probed by 157 background quasars, giving rise to a total of 442 quasar-galaxy pairs.We note here that a significant fraction (152/442) of them are arising from our MUSEQuBES survey.The SFR is measured (not a limit) for 289 galaxies, constituting 389 galaxy-quasar pairs.

Properties of the galaxies contributing to Ly𝛼 stacks
In Fig. 1 we plot the different galaxy properties against one another.The top-left and top-right panels show the impact parameter and normalized impact parameter (/ vir ) of different quasar-galaxy pairs plotted against the redshift of the galaxies.The bottom-left panel shows the stellar mass of the galaxies plotted against redshift.The bottom right panel shows the SFR of the sub-sample of galaxies plotted against stellar mass.The different markers indicate the studies from which the galaxies are drawn.The median (68% range)  and / vir for the 5054 quasar-galaxy pairs used in this study are 1.5 pMpc (0.5-2.5 pMpc) and 10.4 (3.8-17), respectively.There are 204 (419) quasar-galaxy pairs having / vir < 1 (< 2).The median (68% range) log 10 ( * /M ⊙ ) and  for the 4595 unique galaxies are 10 (9.1-10.6)and 0.13 (0.07-0.19), respectively.
The presence of correlations between different galaxy properties is clearly visible in Fig. 1, particularly between the stellar mass and redshift.The Spearman rank correlation test results are summarized in Table 2 for the full sample and the MUSEQuBES galaxies separately.The redshift and stellar mass are more tightly correlated for the full sample (  = 0.48) as compared to the MUSEQuBES galaxies (  = 0.27).This is due to the fact that most of the galaxy surveys in the literature are magnitude limited.
For the MUSEQuBES survey, the FoV of MUSE sets the upper limit on the impact parameter at a given redshift, which in turn leads to a tight correlation between redshift and impact parameter for the MUSEQuBES galaxies (  = 0.52).However, the combination of different surveys eliminates this strong correlation for the full sample (  = 0.24).The mild anti-correlation seen between normalized impact parameter and redshift is a direct consequence of choosing a fixed upper bound on the impact parameter ( max = 3 pMpc). 1At this point, we emphasize that the galaxies observed in MUSEQuBES 1 Since the galaxies at higher redshifts tend to have higher stellar masses survey (orange circles) populate a region of parameter space ( most prominent in  −  * plot of Fig. 1) that previous CGM surveys did not cover.
The distributions of individual galaxy properties (SFR, redshift, stellar mass, impact parameter, and normalized impact parameter from left to right) are shown in Fig. 2. The distributions corresponding to the complete sample are shown in blue, while the subsample with SFR measurements is shown in red.Among the surveys from the literature used in this work, Keeney+18 and Johnson+15 did not report SFRs for their galaxies.The SFR of the subsample with the rest of the 339 galaxies varies from 10 −1.82 -1 M ⊙ yr −1 (68%) with a median of 0.1 M ⊙ yr −1 . 2 The range and median redshift for the galaxy subsample with SFR measurements are consistent with the full galaxy sample.The median stellar mass of the galaxy subsample with SFR measurements (log 10 ( * /M ⊙ ) = 9.3) is lower than for the full sample (log 10 ( * /M ⊙ ) = 10.0).The typical impact parameter of this subsample is also lower than for the full sample (median  ≈ 160 pkpc or / vir ≈ 1.6 compared to ≈ 1.5 pMpc or / vir ≈ 10.4 for the full sample).

ABSORPTION DATA ANALYSIS
We study the distribution of neutral hydrogen around the galaxies in our sample by analyzing the Ly absorption signal.Instead of fitting individual Ly absorption lines associated with each galaxy, a statistical approach is adopted.We stack the quasar spectra after shifting them to the rest-frames of the corresponding foreground galaxies.Here, we focus on the median-stacked Ly absorption.However, we verified that mean-stacked Ly absorption produces consistent conclusions.Spectral stacking is an efficient technique to statistically study the CGM without having to go through rigorous absorption line analysis such as identification, deblending, and profile fitting of individual lines.This method has been successfully used to analyze CGM absorption signals in previous studies (see e.g., Steidel et al. 2010;Rakic et al. 2012;Turner et al. 2014;Chen et al. 2020;Muzahid et al. 2021).

Analysis with rest-frame equivalent width
To obtain the stacked spectra, each normalized quasar spectrum is shifted to the rest-frames of each of the foreground galaxies with  < 3 pMpc.The median stacked spectrum is then generated by calculating the median flux in line-of-sight velocity bins of 40 km s −1 .We have verified that our results are insensitive to the bin size.We have not applied any weighting for the stack, all the galaxy-quasar (and hence higher halo masses and larger virial radii), / vir values decrease with redshift, leading to the observed anti-correlation. 2 The upper limits on the SFR are regarded as measured values for the median and 68% measurements.
pairs are treated equally because they provide different and independent probes of the CGM.
To briefly summarize our analysis procedure: a) we first fit a global continuum to each quasar spectrum and normalize the spectrum by this continuum (see Section 2.3).b) For every foreground galaxy, we select a region of ±1500 km s −1 around the galaxy's redshifted Ly wavelength of the normalized spectrum of the corresponding background quasar.c) We obtain the median normalized flux in 40 km s −1 wide velocity bins for a set of quasar-galaxy pairs.This is referred to as the observed median stacked flux profile.d) Next, a local pseudo-continuum is estimated and subtracted out from the observed median stacked spectra to account for the suppression of the overall continuum below unity due to uncorrelated absorbers.The median flux within a LOS velocity window of ±2000 km s −1 for a stack of random redshifts is used to determine the pseudocontinuum (see Section 4 for details).The offset of the observed pseudo-continuum from unity is added to the median stacked spectra to obtain the continuum-subtracted median stacked spectra.
Finally, e) the Ly   s are measured from direct integration of this continuum subtracted median stacked spectra using LOS velocity windows of ±300 km s −1 (  ,300 ) and ±500 km s −1 (  ,500 ) from the line center.These velocity windows are commonly used in the literature.The errors on   measurements are obtained from 1000 bootstrap realizations of the galaxy sample 3 .We confirmed that convergence is reached in bootstrap distribution for ≳ 200 realizations.

Analysis with pixel optical depth
Owing to the large dynamic range compared to the normalized flux, the pixel optical depth () is another useful quantity that has been used to study the connections between galaxies and gas around them (e.g., Rakic et al. 2012;Turner et al. 2014).We obtained pixel optical depth from the continuum normalized flux  as We set a flag value of  = 10 −6 for  > 1.For heavily saturated pixels with  ≤ 0 or  < Δ where Δ is the error in flux, we set the flag value to  = 104 .The flag values are chosen arbitrarily small and large to ensure that they do not affect the measured median.Similar to our analysis with flux stacking, we masked the aforementioned spectral regions before obtaining the optical depths.
To produce 2D optical depth maps, we stack the Ly optical depth of the absorption associated with galaxies in bins of impact parameter (transverse distance from galaxies) the same way we do as flux stacking.Instead of integrating the median stacked flux profile over the LOS velocity to obtain rest-frame equivalent width, the median optical depth is color-coded as a function of LOS Hubble distance along the y-axis (by converting the LOS velocity to a distance assuming pure Hubble flow at the median redshift of the sample) and transverse distance along the x-axis.These optical depth plots retain information about Ly absorption along the LOS direction along with the transverse direction, thus providing insight into the average kinematics of Ly absorption of a galaxy sample alongside the absorption strength.
3 Each bootstrap realization produces a stack of  quasar-galaxy pairs from a sample of  quasar-galaxy pairs, but with replacement.

RESULTS
The red (blue) absorption profile in the left panel of Fig. 3 shows the median (mean) stacked normalized flux for Ly as a function of LOS velocity for all galaxies in our sample.The pseudo-continua of the mean and median stacked spectra are lower than the actual continuum level, which is unity for individual normalized spectra.Random absorbers with uniform velocity distribution with respect to a galaxy redshift cannot produce any coherent absorption line.They will only suppress the overall continuum below unity, leading to a pseudo-continuum.
In order to model the pseudo-continua we determined the median and mean Ly flux around random redshifts.The random redshifts are chosen from a uniform distribution between the lowest and highest galaxy redshifts in each quasar field.The number of random redshifts is kept the same as the number of galaxies in each field. 4 The average of median (mean) flux within LOS velocity of ±2000 km s −1 is used as the median (mean) pseudo-continuum.The 68% confidence intervals of the median/mean flux, shown by the greyshaded regions, are calculated from 1000 bootstrap realizations of the complete sample.
The observed flux profiles are conventionally normalized by the pseudo-continua before measuring the   s (see e.g., Steidel et al. 2010;Prochaska et al. 2013).Instead of normalizing by the pseudocontinuum, we integrate the observed median stacked flux spectrum over the velocity window of ±300 km s −1 (±500 km s −1 ) and subtract out the contribution stemming from the pseudo-continuum within the same velocity window to obtain the   s.This ensures that the   values do not depend on the pseudo-continuum level, which is determined by stochastic absorption.We found that due to the small values of flux decrement, owing to the lower density of the low- Ly forest, the conventional choice of continuum normalization and our adopted choice of continuum subtraction produce consistent results.
The observed   ,300 and   ,500 for the median stack for the full sample are 0.051 +0.003 −0.003 Å and 0.061 +0.003 −0.003 Å, respectively.For the mean stack,   ,300 and   ,500 are 0.120 +0.005 −0.005 Å and 0.149 +0.006 −0.006 Å, respectively.The quoted errors on the   s are obtained from the 68% confidence intervals of the   distribution of 1000 bootstrap realizations.The measured   s indicate that Ly absorption is detected around our galaxy sample with > 99% confidence interval for both the median and mean stacks, yet the absorption strength is very weak.The significant difference between the mean and median stack is due to skewed flux distribution.This difference is discussed in detail in Muzahid et al. (2021).Unlike their work, the distribution of pixel flux remains left-skewed for a velocity window of ±100 km s −1 , giving rise to a higher median flux value for both the pseudo-continuum and absorption centroid compared to the mean, resulting in a smaller   .For all the subsequent analyses, we used the median stacked spectra.
In the right panel of Fig. 3, we show the 2D median optical depth map for the complete sample.At the low optical depth end, the color scale saturates to a dark-blue color which represents the median Lya optical depth of random regions (median Ly optical depth within ±2000 km s −1 of 4595 random redshifts).A more detailed analysis of the optical depth map is presented in Section 4.2.An excess in optical depth compared to random regions is clearly evident out to ≈ 8 pMpc   The red and blue points denote rest equivalent widths calculated within velocity windows of ±500 km s −1 and ±300 km s −1 , respectively, and are plotted against the median  (/ vir ) of each bin.The error bar on the equivalent width is the 68% confidence interval obtained from 1000 bootstrap samples.The error on  (/ vir ) represents the 68% percentile of the  (/ vir ) distribution in each bin.The median stellar masses and redshifts in each bin are listed in the legends with different markers at the bottom of the plots.The solid blue and red lines in both these plots are the best-fitting power-laws to the blue and red data points, respectively.These best-fit single component power-law relations are provided in the plots.
LOS Hubble distance (≈ 600 km s −1 assuming pure Hubble flow at the median redshift of 0.1), consistent with the width of the median stacked flux in the left panel.The elongation of the excess optical depth along the LOS direction is the manifestation of redshift space distortion.

Ly𝛼 rest frame equivalent width profile
One of our primary goals is to measure the H i Ly equivalent width profile to understand how the cool, neutral gas is distributed in and around galaxies.With this motivation, we divided the impact parameter range of our sample (6-3000 pkpc) into 7 bins and produced stacks of Ly absorption.In order to sample both the inner and outer regions of the CGM, the first two bins are confined within 6-50 pkpc and 50-100 pkpc, respectively.The remaining  range is divided into 5 logarithmic bins.The rest-frame Ly equivalent widths measured from the median stacked spectra are plotted against the median of the impact parameter bins in the left panel of Fig. 4. The red points denote   ,500 and the blue points denote   ,300 .The error bars along the y-axis represent 68% confidence intervals of the median   distribution obtained from 1000 bootstrap realizations.The error bars along the x-axis represent 68% confidence intervals of the impact parameter distribution in each bin.The median stellar mass and redshift of each impact parameter bin are indicated by the legends at the bottom of the plots.
We have a large dynamic range in stellar mass for our galaxy sample.Empirical studies at low redshift as well as theoretical models have shown that the halo mass, as well as the halo size, increases with the stellar mass of galaxy.We, therefore, generated stacked Ly equivalent width profiles as a function of normalized impact parameter (/ vir ; see the right panel of Fig. 4).Instead of dividing the whole / vir range in logarithmic bins, we created 3 bins for / vir < 1 with roughly equal number of galaxies in each, and divided the remaining / vir range (upto 25  vir ) in 4 logarithmic bins.The median stellar mass of each bin are indicated in the legends.
(3) A fixed pivot near the middle of the impact parameter distribution at 400 pkpc is used for the   -profile as a function of the impact parameter to reduce the correlation of errors in slope and normalization.No such pivot is needed for the normalized impact parameter because log 10 (b/R vir ) = 0 is near the middle of the / vir distribution.Note that the power law indices in Eq. 2 and Eq. 3 are consistent with each other within 1, and that an excess Ly absorption around galaxy redshifts is detected with > 99% confidence interval out to ≈ 2 pMpc (or equivalently ≈ 15 vir ) in the transverse direction.
Although a single power-law can adequately explain the  profile, it is to be noted that the second point in the left panel of Fig. 4 is barely consistent with the model.In Section 5.1 we revisit the   -profile for the full sample where we explored the possibility of a power-law + log-linear (or Gaussian) model for the profile.All the measurements related to the   -profiles presented in this work are given as "online only" tables (see Appendix B).

Variation with redshift
The galaxies contributing to the Ly stack span a redshift range of 0.01 − 0.48, corresponding to ≈ 5 Gyr of cosmic time.Since this is almost 35% of the age of the universe, it is interesting to investigate whether the Ly   ,500 -profile evolves with time.
Fig. 5 shows the redshift evolution of the median Ly   ,500profile.First, we split the galaxy sample into three tertiles of the redshift distribution, indicated by the blue (low-, median  = 0.05), black (intermediate-, median  = 0.12), and red (high-, median  = 0.18) colors.Each subsample is then further divided into 4 bins of impact parameters and normalized impact parameters.The first two normalized impact parameter bins are within the virial radius with roughly equal number of galaxies in each.The last two bins outside the virial radius are logarithmic with a bin size of 1.0 dex, 0.8 dex, and 0.7 dex for the low-, intermediate-, and high-mass galaxies, respectively (see the left panel of Fig. 5).For the impact parameter, the first bin extends to 100 pkpc.The rest of the impact parameter range is divided into 4 logarithmic bins with a bin size of ≈0.5 dex for all three mass bins (right panel of Fig. 5).The median   ,500 measured from the stacked spectra are plotted against the median values for / vir and  in each bin.Median stellar masses of each bin are listed in the legends below the plot.
In the left panel of Fig 5, the low- points appear to lie beneath the intermediate-and high- points for  <  vir .This trend is reversed for / vir > 1, where the low- points lie well above the other points.A similar trend is also seen in the right panel of the figure, with the   ,500 -profile for the low- sample showing a shallower slope compared to the intermediate-and high- samples.This is supported by the best-fitting power-law relations indicated in the plots.The power-law indices for the intermediate-and high- bins are consistent with each other in both the left and right panels.However, the low- sample shows a significantly different (shallower) power-law slope.We point out that except for the second / vir bin, the median stellar mass of the low- sample is almost an order of magnitude lower compared to the high- sample.This is not unexpected, since we already noticed a strong trend between redshift and stellar mass in Section 2.5 (see Table 2).As such, the low- subsample has predominantly low-mass galaxies and that may be one of the reasons for the apparent difference.Further, the environment of these galaxies can give rise to the apparent redshift evolution.To mitigate the effects of this apparent redshift evolution, we control the redshift for further analysis whenever necessary.

Variation with stellar mass
The galaxies in our sample span a large range in stellar mass.The dependence of Ly absorption strength on stellar mass at a given (normalized) impact parameter can shed light on the structural variation and self-similarity of the CGM.In order to investigate the possible mass dependence of the Ly   ,500 -profile, we proceed with a similar binning procedure as described in Section 4.1.1.
Fig. 6 shows the dependence of the Ly   ,500 -profile on stellar mass as a function of / vir and  in the left and right panels, respectively.The blue, black, and red colors represent, respectively, the low-(median log 10 ( * /M ⊙ ) = 8.8), intermediate-(median log 10 ( * /M ⊙ ) = 9.8), and high-mass (median log 10 ( * /M ⊙ ) = 10.5)galaxy samples.In the left panel of Fig. 6, the low-mass points are consistently above the intermediate-and high-mass points for  >  vir , indicating stronger Ly absorption around low-mass galaxies outside the virial radii.However, no clear trend is seen for  <  vir .Overall, the Ly   ,500 -profile is significantly shallower for the low-mass sample (slope = −0.55 ± 0.07) compared to the intermediate-mass (slope = −0.77± 0.07) and high-mass (slope = −0.89± 0.06) samples.
In the right panel of Fig. 6, the low-mass galaxies show suppressed Ly absorption compared to their high-and intermediate-mass counterparts for a given impact parameter for the two inner-most bins.However, no significant difference in   is seen among the three different mass samples for  > 500 kpc.Here we also notice a shallower slope for the low-mass sample, but the power-law indices are not as significantly different as in the left panel.
The presence of a strong correlation between  and  * in our complete sample (Table 2) is reflected in the median redshifts of the three mass bins.The redshifts in the legends show a general trend that the median  is higher for the high mass subsample at a given  or / vir bin (with the exception of the second  and / vir bin).However, by controlling the redshifts of the galaxy samples, we confirmed that the trend of   ,500 with  * is not driven by the underlying  −  * correlation (see Section 5.3).

Variation with SFR
The SFR of a galaxy is intimately related to the availability of cool gas in the ISM which is fuelled by the CGM.Star formation driven outflows expel metal-rich gas from galaxies to the CGM and/or IGM.
A fraction of the metal-enriched gas may be recycled back to the galaxies aiding in further star formation.Therefore, connecting the CGM properties and SFR is imperative to understand the role of diffuse circumgalactic gas in galaxy evolution.
We used the 442 quasar-galaxy pairs comprising of 339 galaxies with measured SFR (including upper limits) to investigate the effects  of the SFR on the Ly profile.The left and right panels of Fig. 7 show the dependence on SFR of the   ,500 -profile plotted against / vir and , respectively.The data points in red and blue open squares represent the   measurements for the low-and high-SFR bins.The median SFRs of the galaxies in individual bins are tabulated in Appendix B (see Tables S7 and S8).
The 442 pairs with SFR measurements are divided into two SFR bins based on the median SFR.The galaxies with measured SFR and upper limits below the median SFR are included in the low-SFR bin (median 10 −1.6 M ⊙ yr −1 .Upper limits are treated as detection to compute the median).No galaxies with upper limits on SFR are included in the high-SFR bin (median 10 −0.4 M ⊙ yr −1 ).Galaxies in each SFR bin are further divided into two  or / vir bins.For / vir , the inner bin contains galaxies with / vir < 1 and the outer bin contains galaxies with / vir ≥ 1.For impact parameter , we adopted a separating value corresponding to the 33 percentile of the  distribution (119 pkpc and 107 pkpc for low-and high-SFR bin, respectively).In this way, the outer bin will have twice as many galaxies as the inner bin.This binning strategy is adopted to increase the signal-to-noise ratio of the stacked spectra for galaxies with large impact parameters for which the Ly absorption signal is intrinsically weak.The high-SFR galaxies show significantly stronger Ly absorption for / vir < 1 as indicated by the blue triangle in the left panel of Fig. 7.A similar trend is also seen in the right panel for  ≲ 100 pkpc.However, there is no clear difference in Ly absorption between the two SFR bins outside the virial radius (or at > 100 pkpc).Although we have not explicitly controlled for redshift, the median redshift in a given  or / vir bin is similar for the high-and low-SFR galaxies (the values are tabulated in Appendix Table S7 and S8, along with the 68% confidence intervals).
A similar exercise only with the MUSEQuBES galaxies gives rise to consistent results (open crosses with lighter shades).However, the larger error bars, owing to the smaller numbers of galaxies contributing to the stacks, reduce the statistical significance of the results.The suppression of Ly   outside virial radius for the full sample compared to MUSEQuBES only sample (for both high-and low-SFR bins) in the left panel of Fig. 7 is due to the presence of more high impact parameter quasar-galaxy pairs in the full sample.This is evident from the median and the 68% confidence interval of the / vir and  distributions shown in Fig. 7.The enhanced Ly absorption for the full sample within ≈ 100 pkpc is due to the enhanced median SFR of the contributing galaxies in this bin compared to the MUSE-QuBES only galaxy sample (log 10 (SFR/M ⊙ yr −1 ) = −0.2 for the full sample as opposed to −0.7 for the MUSEQuBES only sample (See Table S7 in Appendix B)).
A similar analysis only for the MUSEQuBES galaxies is not carried out in Section 4.1.1& Section 4.1.2owing to the inadequate number of pairs.

Variation with sSFR
Most of the galaxies with SFR estimates in our sample are scattered around the star-forming main sequence relation (see Fig. 1).The sSFR is thus a better diagnostic of star formation activity which is independent of stellar mass.
We investigate the dependence on sSFR by dividing the sample into star-forming and passive subsamples using a threshold sSFR of 10 −11 yr −1 .Each of these subsamples is again divided into two  and / vir bins following the same binning procedure as used in Section 4.1.3.The left panel of Fig. 8 shows the median   ,500profile for star-forming and passive galaxies with blue and red open boxes, respectively, plotted against / vir .The right panel of Fig 8 shows the same but plotted against .While the median   -profiles for the two subsamples do not show any significant difference when plotted against  in the right panel, the left panel reveals significantly stronger Ly absorption for the star-forming subsample for  <  vir .Nonetheless, no significant difference is seen outside the virial radius.Similar to Fig. 7, the redshifts of star-forming and passive galaxies in a given  or / vir bin are consistent (see Tables S9 and S10 With open stars of lighter shades, we show the Ly   ,500 -profiles for the two sSFR subsamples using only the MUSEQuBES galaxies.Consistent with the full sample, a significantly stronger Ly absorption is observed for the star-forming galaxies for  <  vir .
Moreover, only a marginal difference between the star-forming and passive subsamples is seen outside the virial radius.
The enhanced Ly absorption for the star-forming galaxies outside virial radius for the MUSEQuBES sample compared to the full sample in Fig. 8 can be attributed to the smaller values of / vir and .The apparent inconsistency in results between the full and MUSEQuBES samples of passive galaxies in the smallest impact parameter bins of Fig. 8 is likely due to small number statistics.Only 6 and 7 galaxies contribute to the stack for the passive MUSE-QuBES subsample in the two impact parameter bins (10 and 3 in the two normalized impact parameter bins).The 68% confidence interval was not possible to compute for the small sample size in some MUSEQuBES bins (labelled with ±0.00 range in the Tables in Appendix B).

Optical depth maps
In order to inspect the correlation between galaxies and cool, neutral gas around them simultaneously along the transverse direction and along the line of sight direction, we produced 2D median optical depth maps following the procedure described in Section 3.2.
The map for the complete sample is shown in the right panel of Fig. 3.The first LOS Hubble distance and impact parameter bin is constructed within 0.1 pMpc.Next, the bin size is 0.17 dex and 0.29 dex for the LOS Hubble distance and impact parameter, respectively.The LOS velocity around a galaxy is converted to LOS Hubble distance assuming pure Hubble flow at the galaxy redshift.The negative and positive LOS velocity differences are further merged to increase the / in each bin.The optical depth for the random region is generated following the strategy described in Section 4. The median optical depth of all pixels from the random redshift stack represents the median random optical depth, which is set as the minimum of the optical depth color scale in the map.A Gaussian filter with half of the bin size is used to smooth the raw optical depth map.
The right panel of Fig. 3 reveals enhanced Ly optical depth compared to random regions out to 600 km s −1 or 8 pMpc along the LOS direction.The enhanced optical depth is observed out to 2 pMpc along the transverse direction.The apparent elongation of the excess optical depth in the LOS direction is reminiscent of the "fingers of god" effect and is reported earlier in the literature (Rakic et al. 2012;Turner et al. 2014).This is owing to the peculiar motions of the infalling and/or outflowing gas rather than redshift uncertainties (see e.g., Rakic et al. 2013;Turner et al. 2017;Chen et al. 2020).
Next, we divided our sample into two stellar mass bins with log 10 ( * /M ⊙ ) = 9.9 as the separating value (i.e, the median stellar mass of the full quasar-galaxy pair sample).The median redshifts for the high-(median 10 10.4 M ⊙ ) and low-mass (median 10 9.2 M ⊙ ) galaxy subsamples are 0.145 and 0.088, respectively.However, we did not control the redshift distribution here, as this reduces the number of available pairs leading to insufficient /.The OD maps for the high-and low-mass galaxy samples are shown in the left and right panels of Fig. 9, respectively.The median optical depth of random region for the high-and low-mass galaxy samples are 10 −2.34 and 10 −2.29 , respectively.A common color scale is chosen for the maps which runs from the minimum of these two.We saturate the color scale at the maximum median optical depth of the two samples.The map for the high-mass galaxies shows a significantly stronger signal in the innermost transverse and LOS Hubble distance bin as compared to the low-mass sample.The excess optical depth around high-mass galaxies is also found to be more extended along the LOS direction compared to the low-mass counterparts due to higher pe- culiar velocities, as predicted by simulations (see e.g., Kim & Croft 2008;Rakic et al. 2013;Turner et al. 2017).

DISCUSSION
In this section, we discuss the main findings from Section 4.

Comparison of the Ly𝛼 𝑊 𝑟 -profile with the literature
We found a monotonically decreasing trend of the median   (Ly) with both impact parameter () and normalized impact parameter (/ vir ) (Fig. 4).Such a trend has been reported in the literature with individual Ly   measurements in the CGM of low- galaxies (see e.g., Prochaska et al. 2011;Borthakur et al. 2015) or from   of mean stacked Ly absorption (see e.g., Liang & Chen 2014).However, here we reported the presence of excess cool, neutral gas surrounding galaxies out to 2 pMpc, equivalently ≈ 15 vir .A single power-law with a slope of ≈ −0.75 can adequately represent the median-stacked Ly equivalent width profile (Eq. 2 & Eq. 3).The single-component power-law index is in good agreement with Prochaska et al. (2011), obtained from individual Ly   measurements for sub- * and dwarf galaxies at similar redshifts ( < 0.48) for  < 1 pMpc.Using an -test we confirmed that a two-component power-law model is unnecessary.However, this does not necessarily discard the possible presence of a secondary component other than a power-law.Previously, Borthakur et al. (2016) combined individual Ly   measurements in the CGM of galaxies from the COS-Halos and COS-GASS surveys and obtained a log-linear   ,500 -profile for / vir ≲ 2 with a slope of 0.387 ± 0.103.The power-law nature of the   ,500 −profile at large distances likely arises from the galaxyabsorber correlation which is well described by a power-law (Tejos et al. 2014).While going from smaller to larger impact parameters, a transition from a dark matter halo dominated environment to the regime dominated by halo-halo clustering has been observed in the gas surface density profile (see Zhu et al. 2014).Recently, Wilde et al. (2023) used a Gaussian 1-halo term (arising from the CGM) along with a power-law 2-halo contribution (which is essentially the contribution due to galaxy-absorber clustering) to explain the observed Ly covering fraction profile.The transition point between the 1-halo and 2-halo terms is defined as the extent of the CGM in their work.
In order to look for a similar transition region between 1-halo and 2-halo contributions, we fit a log-linear and a Gaussian profile along with the default power-law model to the observed   ,500 −profile plotted against / vir .The two models are as follows: ,500 = 10 −  (/ vir ) +  (/ vir )  (4) and, The left and right panels of Fig. 10 show the best-fit models with Eq. 4 and Eq. 5, respectively.The parameters of the best-fit models are summarized in Table 3.The slope for the log-linear relation we obtain is somewhat steeper than found by Borthakur et al. (2016) but consistent with 2.An −test rejects the null hypothesis that a single component power-law is a better representation of the   ,500 −profile compared to either of the two-component models with a -value of > 0.96.The 1-halo term, represented by the Gaussian or log-linear function, contributes negligibly to the profile at / vir ≳ 1 which closely resembles the CGM-scale ( CGM ) in the recent study of Wilde et al. (2023).
Using a single sightline (towards 3C 273), Morris et al. (1993) found that Ly absorbers indeed cluster around galaxies, but with a smaller amplitude than galaxy-galaxy clustering.Chen et al. (2005); Chen & Mulchaey (2009) found that strong Ly absorbers cluster around emission-line dominated galaxies with a clustering amplitude  comparable to the autocorrelation amplitude of emission-line galaxies.They also concluded that weak Ly absorbers cluster weakly around galaxies.Tejos et al. (2014) further showed that ≈ 50% of weak Ly lines are correlated with galaxies at large transverse distances.The low column density Ly absorbers can be attributed to the gas present in filament-like structures.Consistent with this picture, Wakker et al. (2015) found a trend of increasing Ly equivalent width and line width with decreasing filament impact parameter.The Ly detection rate is ≈ 80% within 500 pkpc of galaxy filament.Recently, Bouma et al. (2021) studied the relation between Ly absorbers and nearby galaxy filaments.They found an excess incidence rate (d/d) near filaments and a somewhat shallower slope for the column density distribution function compared to the general population of Ly absorbers at  ≈ 0. They also noted that the strongest Ly absorbers are preferentially detected near galaxies or filament axes, albeit with a significant scatter.The observed excess Ly   far outside the virial radii of galaxies can, in part, come from the CGM of other galaxies ('intra-halo' gas) residing in the filaments.The 'inter-halo' gas tracing the underlying density fields of the filaments can also contribute.The   ,500 −profile in our fit is dominated by a single component power-law outside the virial radius.Such a powerlaw is widely usually used in the context of galaxy-galaxy clustering.However, the power-law component is regarded as the 2-halo term in our work.Although, strictly speaking, the gas constituting the 2-halo component is expected to have an 'intra-halo' origin, we emphasize that this can account for gas in large-scale structures or filaments that are not part of the CGM of any (detectable) galaxy.The connection between Ly absorbers and galaxies, over a scale of > 500 ℎ −1 kpc, is investigated by Davé et al. (1999) using hydrodynamical simulations.They obtained a power-law index of −0.71 for the   ,500 -profile around low redshift galaxies, which is in good agreement with our single component power law fit.To investigate the origin of the   ,500 −  correlation in their simulations, they determined the power-law index of the   -profile for the diffuse gas phase with densities and temperatures typical of the IGM.The overall power-law index and the index for the diffuse gas phase were found to be statistically indistinguishable.Based on this finding, they argued that the trend between   and impact parameter can arise from the clustering of gas and galaxies over large-scales.In the next section, we will examine whether the stacked Ly profile can be explained by the large-scale clustering of gas and galaxies.

The role of galaxy-absorber clustering
The stacked spectral profile for the complete sample is mainly shaped by the velocity distribution of the Ly absorbers with respect to the host galaxies.If the width of the velocity distribution is mainly determined by the random redshift errors, a Gaussian is a reasonable choice to model the stacked profile.However, the left-top panel of Fig. 11 shows that a single-component Gaussian cannot explain the observed median stacked flux profile of Ly absorption for the full sample.Two distinct Gaussian components with  values of 114.9 ± 6.5 and 341.7 ± 26.0 km s −1 are required to adequately explain the entire profile.The narrower and stronger Gaussian component likely arises from the Ly absorbers originating in galactic halos (CGM).The origin of the shallow and wider Gaussian component, however, can be more complex.It can arise from weak and broad individual Ly absorption associated with galactic halos (i.e., so-called broad Ly absorbers (BLAs) buried in the noise of individual spectra).Another possible origin is the narrow and loosely correlated Ly absorbers likely arising from the large-scale structures tracing the same overdensities as the galaxies.
The large-scale clustering between galaxies and Ly absorbers has been characterized in the literature (see e.g., Chen et al. 2005;Tejos et al. 2014;Wilde et al. 2021;Borthakur 2022).The galaxy-absorber two-point correlation function can be parameterized as: where,  0 and  are the scale-radius and power-law slope, respectively.The 3D distance  = √︃  2 +  2 ∥ .The correlation function  (,  ∥ ) is essentially a measure of the excess number of absorbers as compared to random regions, hence it can be written as where,  obs and  rand are the number of observed and random absorbers at a given  and  ∥ .At large impact parameters of ∼ 1 pMpc, the absorbing gas is unlikely to arise from individual galactic halos.As discussed in Section 5.1, the contributing absorbers most likely reside in large-scale structures tracing the same overdensities as galaxies.At this length scale, it may be safe to assume a similar distribution of absorption strength for the uncorrelated Ly absorbers and for those absorbers which are correlated with the galaxies.The excess Ly absorption we see at large impact parameter can be interpreted as the increased number of absorbers around galaxies due to clustering, and not necessarily because of increased strength.With this assumption, the two-point correlation function can be directly related to the observed median stacked flux.
The optical depth can be written as  = ⟨⟩, where ⟨⟩ and  are the average strength and the number of contributing absorbers, respectively.For the median stack, the median optical depth at a given LOS velocity bin   can be written as The assumption of similar absorption strength distributions for the random and clustered absorbers implies ⟨⟩ med,obs ≈ ⟨⟩ med,rand .Thus, Note that the observed median stacked flux can be written as,  med,obs =  − (  med,obs −  med,rand + med,rand ) .
Combining Eq. 8 and Eq. 9 and using  −  med,rand =  med,rand , we obtain . (10) Here we dropped the superscript  from  med,rand , as the optical depth of the random region does not depend on the LOS separation from galaxy redshift.We fix  med to the median impact parameter of our galaxy sample of 1.5 pMpc.The value of  med,rand is obtained by converting the median flux at random regions (essentially the pseudo-continuum of the median stack) to optical depth.To allow for a small velocity offset of the centroid of the stacked spectrum from 0 km s −1 , we introduce a parameter  in the two-point correlation function such that: where, we wave used the 3D distance  = √︃  2 + ( ∥ − ) 2 in Eq. 6.The offset parameter  is a purely mathematical construct to account for the velocity offset that can arise due to the finite LOS velocity bin size of 40 km s −1 .
Fitting the observed median stacked rest-frame Ly absorption profile with Eq. 10 and Eq.11 yields best fit  0 and  values of 7.6 pMpc and −1.57, respectively.The best-fit5  value is 0.24±0.06pMpc, corresponding to 18 ± 4 km s −1 which is well within the LOS velocity bin size used in our analysis.The best-fitting model profile is shown in the left-bottom panel of Fig. 11 with the red solid line.The corner plot in the right panel of Fig. 11 shows the 1-and 2dimensional projections of the posterior distribution for the model parameters.
The reduced- 2 of 1.9 obtained for this fit is larger than the reduced- 2 of 0.6 obtained for a double-component Gaussian fit.However, a fit with this simple two-point correlation function provides a significant improvement compared to a single component Gaussian fit (reduced- 2 = 2.7).This indicates that the broad component in the stacked profile for large impact parameters can be partly explained by a simple two-point correlation function arising from the large-scale clustering between galaxies and Ly absorbers.
At this point, we emphasize that previous studies on clustering analyses measured the  0 and  from the correlation function projected along the line of sight, thus eliminating the effect of redshift space distortions.The correlation length-scale and power-law index obtained in our analysis however suffer from redshift space distortions since the two-point correlation function is not projected along the line of sight.The evidence for redshift space distortion is also clear from the optical depth maps (see Figs. 3 & 9).The correlation length of 7.6 pMpc found in this work is considerably larger compared to the projected galaxy -H i absorber correlation length of  0 = 1.6ℎ −1 70 Mpc found by Tejos et al. (2014).This is also true even when we compare the  0 by fixing  = −1.44 as reported in Tejos et al. (2014)).Recently, Wilde et al. (2021) presented projected correlation lengths of 3.4ℎ −1 68 Mpc and 6.2ℎ −1 68 Mpc for high ( * > 10 9.9 M ⊙ ) and low ( * < 10 9.2 M ⊙ ) mass galaxies at low redshift.The power-law index for the projected correlation function varies from −1.2 to −1.9 for low-to high-mass galaxies in their work.Nonetheless, this exercise demonstrates that the presence of the significant non-Gaussian wings in the median stack can be interpreted as a consequence of weak, correlated absorbers likely arising from the large-scale structures tracing the same overdensities as the galaxies rather than warm/hot, widespread gas associated with individual halos.
Another possible origin of the broader component could be the environment of the galaxies.We have employed a simple friendsof-friends (FoF) algorithm with a projected separation of 10 ′ and LOS separation of 500 km s −1 to identify overdense galactic environments.We have defined 'groups' as the overdense region of more than four galaxies identified by the FoF algorithm with the aforementioned linking lengths.A power-law galaxy-galaxy clustering with a slope of −1.8 and correlation length of 3ℎ −1 Mpc is expected to give rise to ≈ 4 galaxies in the region of interest (at  ≈ 0.1) due to clustering (see Cherrey et al. 2023, for the formalism).The adopted number of five or more galaxies within the linking lengths thus suggests an overdense region.At this point, we emphasize that these overdensities are not necessarily virialized.In Fig. 12, we show the stacked Ly absorption profile for isolated and group galaxies with blue and teal histograms, with the shaded region denoting the corresponding 68% confidence interval obtained from 1000 bootstrap realizations.The median stacked Ly absorption from the group galaxies showed symmetric secondary peaks at around ≈ ±500 km s −1 , but no such feature was present in the stacked spectrum of isolated galaxies.While we acknowledge that it is challenging to characterize the environment of galaxies observed from different surveys with varying depths and completeness, the presence of the secondary peaks in the stacked spectrum for the group galaxies is consistent with the findings of Muzahid et al. (2021) for a homogeneous sample of  ≈ 3 galaxies.Moreover, we obtained a consistent result when we considered the galaxies with > 90% completeness from the Keeney+18 survey.A detailed analysis of the effects of the galaxy environment on the CGM will be presented elsewhere.However, we emphasize that the secondary peaks observed in the stacked spectrum of group galaxies are very different compared to the smooth wings of the observed spectrum, and they cannot be explained by a power-law cross-correlation function.The observed broad component in the stack at large impact parameters is more likely the manifestation of the galaxy-absorber correlation function.

The role of stellar mass
The observed   ,500 -profile is shallower for the low-mass galaxies than for their high-mass counterpart (see Fig. 6).This owes to the fact that the high-mass galaxies show stronger Ly absorption in the inner regions ( < 0.5 vir ) whereas the low-mass galaxy sample shows relatively stronger Ly absorption outside the virial radius.We note however that we did not control for the redshifts of the different mass bins.Consequently, the median redshift of the high-mass bin is higher compared to the low-mass bin., measured within 2 vir , as a function of stellar/halo mass is plotted with black symbols.All galaxies with 0.25 < / vir < 2 are selected and split into 3 stellar mass bins with log 10 (  * /M ⊙ ) < 8, > 9.5 and 8 <log 10 (  * /M ⊙ ) < 9.5 to produce the Ly stacks.The magenta triangle represents the median   ,500 for stacked Ly absorption in low redshift ( < 0.5) cluster outskirts (/ 500 < 1.5) by Mishra et al. (2024).The x-error bars show the 68% range of  * in each bin, y-error bars are 68% confidence intervals of the median   .The median / vir and  (along with 68% confidence intervals) of each bin are given in the legends.The overlayed red dotted line showing the star-formation efficiency (see text) corresponds to the right y-axis.
In the left panel of Fig 13, we show the Ly   measurements for three mass bins within and outside the virial radius.We first divided the galaxy sample into two bins with / vir > 1 and ≤ 1. Galaxies in each / vir bin are split by log 10 ( * /M ⊙ )< 9, > 10, and 9 <log 10 ( * /M ⊙ )< 10.In each normalized impact parameter bin, we now controlled the redshifts of the galaxies so that the redshift distributions of the low-, intermediate-, and high-mass samples are similar. 6The median log 10 ( * /M ⊙ ) and median redshifts of different bins are tabulated in Appendix Table S11.Outside the virial radius, the Ly   is significantly higher for the low-mass bin compared to the intermediate-and high-mass bins.Inside the virial radius, we notice a decline in   for the high-mass bin compared to the low-and intermediate-mass bins, but consistent within the 1 uncertainties.The error bars are similar to Fig. 6.The data points inside the virial radius are plotted with small offsets in the abscissa for clarity.Measurements for the MUSEQuBES galaxies for a similar binning strategy are shown with a lighter shade and open star symbols.Similar to Fig. 7, the enhanced Ly absorption outside virial radius for MUSEQuBES galaxies compared to the full sample can be attributed to a lower median / vir compared to the full sample.The trend for the MUSEQuBES galaxies is consistent with the full sample, although not significant due to a smaller sample size.
The presence of significant Ly absorption outside the virial radius can arise from the extended CGM and/or the large-scale structures around galaxies.Cosmological hydrodynamical simulations by 6 We chose equal number of galaxies within some narrow Δ for the three mass bins.This ensures that the redshift distributions (and the median values) of the three mass bins are very similar.van de Voort & Schaye (2012) showed that the "cold gas", which can be traced by Ly, is mainly residing in filamentary structures around galaxies.Moreover, the geometry of such cosmic filaments is much more widespread (up to 4  vir ) for low-mass halos as compared to high-mass halos (see their figure 1), leading to higher gas covering fraction for low-mass galaxies outside virial radius.Johnson et al. (2017) reported the presence of Ly absorption well beyond the virial radius of dwarf galaxies.Recently, the observations of Wilde et al. (2021) also showed that the Ly covering fraction remains roughly constant at ≈ 40% out to 1 pMpc from galaxies with log 10 ( * /M ⊙ )<9.2.The covering fraction for the high-mass galaxies, however, declines rapidly with impact parameter.Although their quoted 1 errors have significant overlaps in different mass bins, this may hint at the steeper slope for the   -profile for the high-mass galaxies.
A possible reason for the enhanced Ly absorption outside the virial radius for low-mass galaxies could be the environment.However, different completeness of the different surveys hinders a proper investigation of environmental effects for our complete sample.To circumvent this, we selected a subsample of galaxies from Keeney+18 with total completeness >90% (see their Table 8).In that subsample, isolated galaxies are searched for using a simple 3D FoF algorithm with a linking separation of 10 ′ 7 and a linking velocity of 500 km s −1 .We found enhanced Ly absorption around low-mass galaxies outside the virial radius (median / vir ≈ 6) even for this well-chosen sample of isolated galaxies.Therefore, we conclude that the observed mass dependence cannot be directly attributed to the environment.However, if a significant fraction of these low-mass galaxies is satellite, then their  vir measurements can be underestimated.This can lead to the apparent excess 'outside' the virial radius of the low-mass galaxies.Indeed, the convergence of the   ,500 profile for the three mass bins at large  (see the right panel of Fig. 6) implies that there is no excess Ly absorption at large proper distances around low-mass galaxies.
The underlying gas distribution in cosmic filaments can explain the origin of excess Ly   far outside the virial radius for low-mass galaxies compared to high-mass galaxies when plotted against the normalized impact parameter.At a similar normalized impact parameter, massive galaxies trace a region further away in the filament compared to less massive or dwarf galaxies.This can cause the apparent enhancement of Ly   outside virial radius when plotted against normalized impact parameter.However, the underlying gas distribution tracing the overdensity in the filament will produce similar absorption for both low-and high-mass galaxies at large physical distances, consistent with our observation.In section 5.1, we suggested that the power-law 2-halo term outside the virial radius can account for the gas in the filaments of 'inter-halo' origin as well.
A marginal suppression in the Ly absorption inside the virial radius for high-mass galaxies compared to the low-/intermediatemass bin is observed when the redshift is controlled.To further explore this with larger galaxy samples, we selected all galaxies with 0.25 < / vir < 2.0 in three mass bins with log 10 ( * /M ⊙ ) < 8, 8 < log 10 ( * /M ⊙ ) < 9.5, and log 10 ( * /M ⊙ ) > 9.5 to generate Ly stacks.The lower cut in / vir is motivated by the lack of lowmass galaxies below this limit in our galaxy sample.Thus, a cut like this roughly brings the / vir distributions of the three mass bins on a similar footing.Although we did not explicitly control for the redshift, the median redshifts are comparable for the three mass bins (this is also true for the 68% ranges).
The right panel of Fig. 13 shows the median Ly   ,500 as a function of stellar mass (top axis) and halo mass (bottom axis).To convert the stellar masses to halo masses, we used the sub-halo abundance matching relation of Moster et al. (2013).The median   ,500 for our highest mass bin shows a clear suppression compared to the intermediate mass bin but is comparable to the low-mass bin.The (pink) data point is obtained from Mishra et al. (2024) which represents the median Ly   ,500 measured around  ≈ 0.17 clusters using background quasars with impact parameters within 1.5 200 .The Ly   ,500 measured for galaxy clusters is significantly lower compared to the   ,500 we obtained for the highest mass bin of our sample.
The dotted red line in the plot (right y-axis) shows the starformation efficiency (SFE) as a function of halo mass.The SFE is defined as SFE ≡  * / baryon =  * /( halo *   ), where   (≡ Ω b Ω M ) is the mean cosmic baryon fraction.As mentioned earlier, we used the stellar mass -halo mass (SMHM) relation from Moster et al. (2013) to obtain the SFE curve.Both quantities show a peak at some halo mass (∼ 10 11 M ⊙ for the Ly bearing gas and ∼ 10 12 M ⊙ for the SFE) and decline at both the low-and high-mass ends.Such a similarity may suggests that the cool neutral gas content of the CGM has direct consequences for the efficiency of star formation inside galaxies.However, we note that the comparison of a relative quantity such as the SFE with an absolute quantity like the   should be done with caution.
Although the connection of the Ly-bearing gas with the SFE is not straightforward, one particularly important aspect of the right-hand panel of Fig. 13 is the increase of   with decrease in  halo from 10 14 to 10 11  ⊙ .Although the SFE peaks at around ≈ 10 12  ⊙ and then decreases with decreasing halo mass, this is not reflected in the Ly bearing gas.The lower SFE of galaxies with  halo ≈ 10 11  ⊙ compared to ≈ 10 12  ⊙ halo mass galaxies, despite being more gas-rich, may be indicative of a larger gas depletion time for these galaxies.Conversely, despite a decreasing Ly   , the increased SFE of galaxies with  halo ≈ 10 12  ⊙ compared to ≈ 10 11  ⊙ halo mass may suggest that the circumgalactic gas is depleted before the star formation inside galaxies.The recent simulation by Appleby et al. (2022) showed that the column density distribution function (CDDF) of H i absorbers around green valley galaxies resembles that around quenched galaxies more closely than that of the star-forming galaxies.They interpreted this as a sign of CGM depletion before the quenching of the host galaxy (see also, Davies et al. 2020;Oppenheimer et al. 2020).This is in agreement with our observational findings.
Finally, the suppression of cool, neutral gas around high-mass galaxies can also be a consequence of higher virial temperatures leading to a higher degree of ionization of the CGM.However, the virial temperature corresponding to the intermediate-mass bin is already above the temperature required to collisionally ionize the neutral hydrogen.Simulations predict a transition from almost no virialized gas to substantial virialized gas at  halo ∼ 10 12 M ⊙ -which is the transition between cold to hot mode accretion (see e.g., Kereš et al. 2005) .The observed suppression of Ly absorption around massive galaxies can be a possible indication of this transition.Using a sample of Milky Way-type galaxies in the IllustrisTNG simulation, Ramesh et al. ( 2022) reported that the CGM of high mass galaxies exhibits less H i gas compared to low-mass galaxies.They attributed this effect to the kinetic mode feedback by SMBH residing in more massive galaxies which can be responsible for sweeping up of the cool gas in the CGM.This is consistent with our observations.
In passing, we note that our results contrast with Bordoloi et al. (2018), who reported a positive trend between Ly   and  * using a sample of 85 galaxies with  < 160 kpc and 8 ≤ log 10 ( * /M ⊙ ) ≤ 11.6.This apparent disagreement may be due to their small sample size, and the fact that they did not take into account the upper limits in their analysis.

The role of SFR
Using a subsample of galaxies with SFR estimates, we studied the impact of SFR and sSFR on the distribution of cool, neutral gas in and around galaxies.The left panel of Fig. 7 shows a clear indication of enhanced Ly absorption within the virial radius for galaxies with high SFR.However, outside the virial radius, no dependence on the SFR is seen.Similarly, in the right panel, a strong SFR dependence is seen only within 100 pkpc.Although we did not explicitly control the redshifts in this analysis, the median redshifts of the compared bins are very similar.The measurements obtained for only the MUSE-QuBES galaxies (shown with green and magenta points) are also consistent with the complete sample.
No statistically significant trend is reported between Ly   and SFR for the galaxies in the COS-Halos survey (Thom et al. 2012;Tumlinson et al. 2013).In fact, Thom et al. (2012) argued that the CGM of passive galaxies is equally rich in H i. However, Borthakur et al. (2015, COS-GASS survey), found a marginal correlation of individual Ly absorption strength with SFR using 45 low-redshift galaxies, with the majority of the sightlines passing within the virial radius of the galaxies.Combining the observations of the COS-HALOS survey along with COS-GASS, Borthakur et al. (2016) found a stronger correlation between SFR and impact parameter-corrected Ly equivalent width, defined as the ratio of observed Ly equivalent width and that predicted by the best-fit model for the   -profile for the entire sample.Their combined sample is also limited to sightlines within ≈ 1.4  vir .This is consistent with the results of our analysis using spectral stacking.A similar trend is also seen for  ≈ 3.3 Ly emitters (LAEs; see Muzahid et al. 2021).However, the LAE sightlines probe the CGM at  >  vir , where we did not see any trend with SFR.
The trend between Ly absorption around galaxies and SFR can arise from the two following scenarios.First, star-formation driven outflows can entrain and deposit cool gas in the CGM.We recall that the strong SFR dependence is only seen for impact parameters ≲ 100 pkpc.The gas traced by Ly at large galactocentric distances does not show any dependence on SFR.This suggests that in this scenario the outflows are only effective in determining the gas distribution out to ≈ 100 pkpc.Second, the availability of more cool gas in the CGM can lead to higher SFRs in galaxies.Although it is not straightforward to relate the gas in the CGM to star-formation activity inside galaxies (since the gas in the CGM has to go through different physical processes via which it can cool and become molecular before forming stars), the second scenario is predicted by simulations (see e.g., Davies et al. 2019Davies et al. , 2020)).
Finally, we note that the median stellar mass of the high-SFR bin is higher than for the low-SFR bin.This is not surprising since the majority of the galaxies in our sample follow the main sequence relation (see Fig. 1).However, the observed trend with SFR is unlikely due to the difference in stellar masses for the two SFR bins, since a similar trend is also seen for sSFR, and again only within the virial radius (see the left panel of Fig. 8).The fact that we see enhanced Ly absorption inside the virial radius for star-forming galaxies compared to quenched galaxies (sSFR< 10 −11 yr −1 ) is consistent with the findings of Johnson et al. (2015), who found an enhanced covering fraction inside the virial radius of late-type galaxies compared to early-type counterparts.

Implications of the optical depth maps
We produced 2D optical depth (OD) maps for the complete sample as well as for two sub-samples with high-mass and low-mass galaxies.In the ideal case of an isotropic distribution of gas around galaxies, circularly symmetric OD maps should be observed.The departure from circular symmetry in the maps indicates the presence of redshift space distortions.This is clear from the enhanced optical depth along the LOS direction at small impact parameters (Fig. 3 and Fig. 9).However, at larger impact parameters, this effect is less pronounced.Such a distortion along the line of sight direction is also observed by Tejos et al. (2014).They found that the extent of the distortion was consistent with the redshift uncertainty of their galaxy sample.We observed the excess OD up to ≳ 1 pMpc LOS Hubble distance for the lowest impact parameter bin.Considering a redshift uncertainty of ≈ 50 km s −1 for our galaxy sample, the length scale of the distortion would be ∼ 0.7 pMpc at the median redshift of  = 0.1.Hence, the observed elongation of the OD map along the LOS direction unlikely to be entirely due to redshift uncertainty of our galaxy sample.Turner et al. (2014) reported a similar redshift space distortion for a sample of LBGs at  ≈ 2.3.Comparing with the eagle simulation, Turner et al. (2017) argued that infalling gas is responsible for the redshift space distortion.A similar comparison of our observed optical depth maps with simulations could shed light on the gas kinematics around low-redshift galaxies.
The strong excess Ly absorption within ≈ 100 pkpc transverse and ≈ 1 pMpc LOS Hubble distance around high-mass galaxies is consistent with our analysis of the   -profile (right panel of Fig. 6).This is a possible indication of sightlines passing through denser, cool gas clouds originating in high-mass galaxies from galactic processes (e.g., extended galactic disk, clouds from galactic fountains).An increasing strength of Ly absorption with stellar mass (and hence halo mass) in the vicinity of galaxies has been predicted by simulations (see e.g., Rakic et al. 2013;Turner et al. 2017).
The OD maps are qualitatively in agreement with the maps produced by Chen et al. (2020) for high-redshift ( ≈ 2) galaxies.However, they find that the total optical depth in the redshifted region of the optical depth map is larger by > 50% than that of the blueshifted side (for 70 < (kpc) < 150).Producing the OD map without folding the negative and positive LOS velocities we confirmed that this asymmetry is not pronounced in our case.The total optical depth in the redshifted region is only ≈ 7% higher compared to the blueshifted region.The "least implausible" explanation for this asymmetry between redshifted and blueshifted flux according to Chen et al. (2020) was that the Ly emission from galaxy halos contaminates the absorption features.This can be non-negligible if a foreground galaxy is being probed by a background galaxy spectra, as was the case in their study.However, our galaxy sample is probed by much brighter background quasars, so this possible source of contamination can be neglected for all practical purposes.The lack of any significant asymmetry in the optical depth maps found in this work may be a direct consequence of the aforementioned reason.

Possible caveats
In this study, we combined galaxy samples from six different surveys from the literature (non-IFS) along with our MUSEQuBES galaxies (IFS) leading to 5054 QSO-galaxy pairs.While such a large number of quasar-galaxy pairs is critical to probe any difference in the Ly   -profile as a function of galaxy properties, it inevitably introduces heterogeneity in the galaxy sample.The different archival surveys have different depths and spectroscopic completeness, and probe a wide range of redshifts ( ≈ 0.01 − 0.48 corresponding to 5 Gyr of cosmic time).Further, surveys such as COS-Halos and Liang & Chen (2014) were designed to study the CGM of 'isolated' galaxies in contrast with the Keeney+18 survey which includes a substantial fraction of group galaxies in it.Note that because of the lack of SFR measurements, the Keeney+18 sample does not contribute to our analysis of the SFR and sSFR dependence.Hence for the SFR and sSFR dependence, the contributing galaxies are mostly 'isolated', and environmental effects are not expected to be important.
We took two measures to minimize the effects of the heterogeneous galaxy sample: (1) Most of our analyses are based on the Ly  profile as a function of impact parameter and normalized impact parameter (/ vir ).This naturally takes care of the different distances from the galaxy centre when the dependence on a given galaxy parameter is investigated.(2) Whenever required, we controlled the redshift of the galaxies in a given  and / vir bin (see e.g., Fig. 13).
Finally, another form of bias can be introduced as a consequence of correlated galaxy properties (e.g., between  and  * and between  * and SFR).As pointed out earlier, the merging of different samples helps us to mitigate such intertwined correlations (see Fig. 1).We have not explored the effect of galaxy orientation with respect to the quasar sightlines in this work, as it is beyond the scope of this paper.

SUMMARY
In this study, we used 4595 galaxies with median  of 0.1 (68% range of 0.07 -0.19) from the MUSEQuBES and archival CGM

Figure 2 .
Figure 2. (a) SFR, (b) redshift, (c) stellar mass, (d) impact parameter, and (e) normalized impact parameter probability density distribution functions for the 4595 galaxies (339 with SFR measurements) in our complete sample.The galaxies with SFR measurements (including upper limits) are shown in red.

Figure 3 .
Figure 3. Left: The median (red) and mean (blue) stacked Ly absorption profile for the full sample in bins of 40 km s −1 .The mean and median pseudocontinua, shown by the blue and red dashed lines respectively, are obtained from stacks of random redshifts.The shaded regions indicate 68% confidence intervals of the mean and median flux distributions obtained from 1000 bootstrap realizations of the complete sample.Right: Median Ly optical depth around galaxies as a function of impact parameter and (absolute) LOS Hubble distance (left y-axis) and LOS velocity (right y-axis).The first bin is confined within 0.1 pMpc, onwards the bin size is 0.17 dex and 0.29 dex for LOS distance and impact parameter, respectively.The map has been smoothed with a Gaussian kernel with standard deviation of half of the bin size.The minimum of the optical depth color scale is set to the value obtained for the stack around random redshifts.The elongation of the signal along the LOS direction indicates the presence of redshift space distortions.

Figure 4 .
Figure 4. Median Ly rest equivalent width as a function of physical impact parameter (left) and normalized impact parameter (right).The red and blue points denote rest equivalent widths calculated within velocity windows of ±500 km s −1 and ±300 km s −1 , respectively, and are plotted against the median  (/ vir ) of each bin.The error bar on the equivalent width is the 68% confidence interval obtained from 1000 bootstrap samples.The error on  (/ vir ) represents the 68% percentile of the  (/ vir ) distribution in each bin.The median stellar masses and redshifts in each bin are listed in the legends with different markers at the bottom of the plots.The solid blue and red lines in both these plots are the best-fitting power-laws to the blue and red data points, respectively.These best-fit single component power-law relations are provided in the plots.

Figure 5 .
Figure5.Dependence of the median Ly   ,500 -profile on redshift.The galaxy sample is divided into three tertiles of the redshift distribution indicated by the blue, black, and red colors which correspond to median redshifts of 0.05, 0.12 and 0.18 as shown in the legends along with the 68% confidence intervals.The data points represent median   ,500 measurements plotted against median / vir (left) and median  (right) for each redshift bin.The median stellar masses of galaxies in each bin are listed in the legends below the plot.The error bars are similar to Fig.4

Figure 6 .
Figure 6.Dependence of the median Ly   ,500 -profile on stellar mass.The galaxy sample is divided into three tertiles of the stellar mass distribution indicated by the blue, black, and red colors.Median redshifts of galaxies in each bin are listed in the legends below the plot.The other details are similar to Fig. 5.The low-mass galaxy sample (blue) shows a significantly shallower slope as compared to the intermediate-(black) and high-mass (red) samples, particularly in the left panel.

Figure 7 .
Figure 7. Dependence of the median Ly   ,500 -profile on SFR as a function of / vir (left) and  (right).Only two SFR bins and two  (and / vir ) sub-bins are created for stacking due to the reduced number of galaxies with SFR estimates.The red and blue open squares denote the   ,500 measurements for the lowand high-SFR bins, respectively.The median SFR and redshift (along with the 68% confidence interval) of each bin are tabulated in appendix (S7 and S8).Upper limits are treated as detections to compute the median SFR.The open stars in lighter shades indicate measurements for the MUSEQuBES subsample only.Error bars are similar to Fig. 4. A strong SFR dependence of the circumgalactic Ly absorption is seen only within the virial radius.

Figure 8 .
Figure 8. Same as Fig. 7 but for sSFR.A strong dependence of the median Ly   ,500 profile on the sSFR is seen within the virial radius in the left panel.

Figure 9 .
Figure9.2D median optical depth maps of Ly absorption around high-mass galaxies (left, median log 10 (  * /M ⊙ ) = 10.4) and low-mass galaxies (right, median log 10 (  * /M ⊙ ) = 9.2).The two galaxy samples are separated at log 10 (  * /M ⊙ ) = 9.9 , which is the median stellar mass of our sample.The color scale is common for both maps.The median Ly optical depth is plotted as a function of impact parameter and LOS Hubble distance (LOS velocity in the right y-axis).Similar to Fig.3, positive and negative velocity sides are folded to increase the / ratio.The bin size and smoothing procedure are similar to Fig.3.The high-mass sample shows a significantly stronger absorption signal in the inner parts of the map.

Figure 10 .
Figure 10.Modeling the observed Ly   ,500 -profile with two components.The data points show   ,500 plotted against / vir from this work.We use a log-linear (left, orange dashed) or a gaussian (right, orange dashed) and a power-law component (green dashed, both panels) as shown by the legends.The solid red line is the combined best-fit model in both panels.The single-component power-law from Fig.4is overlayed in both panels with grey dotted lines.

Figure 11 .
Figure 11.Left-top: The observed median stacked Ly absorption profile for the complete sample (black histogram).The blue solid line indicates the two-component Gaussian fit to the observed profile (extrapolated region shown with dotted line), with the individual components shown in red and green dashed lines.The standard deviations of the individual Gaussian components are indicated in the legends.The grey error bars are 68% confidence intervals of the median flux in each bin coming from 1000 bootstrap realizations.Left-bottom: Observed profile fit with the 2-point correlation function given in Eq. 11 shown by the red solid line (extrapolated region shown with red dashed line).The shaded region shows 1 scatter around the best-fit relation.Right: Corner plots showing one and two-dimensional projections of the posterior probability distributions of the parameters used to fit the median stacked profile.

Figure 12 .
Figure12.The median stacked Ly absorption profiles for the isolated and group subsamples shown with magenta and teal histograms, respectively.The shaded region represents the 68% confidence interval of median flux distributions obtained from 1000 bootstrap realizations.

Figure A1 .
Figure A1.Probability density distribution of stellar mass (left) and redshift (right) of galaxies in the complete sample, shown with stacked histograms.The stellar mass and redshift distribution of galaxies from different CGM surveys are shown with histograms of different colors and symbols.

Table 1 .
Overview of the 16 MUSEQuBES fields.From left to right, the columns show the quasar name, right ascension (J2000), declination (J2000), redshift,  band magnitude, exposure time ( exp ) for the COS/G130M grating, / per resolution element at  = 1250 Å, exposure time for the COS/G160M grating, / per resolution element at  = 1650 Å, HST programme ID of the COS observations, exposure time for the MUSE observations, and the effective seeing measured in MUSE cubes at =7000 Å.
The horizontal black dashed lines represent the median  and / vir for the complete sample.Bottom-left: Stellar mass plotted against redshift of the 4595 galaxies.The black dashed lines indicate the median log 10 (  * /M ⊙ ) and .Bottom-right: SFR plotted against stellar mass for the subsample of 339 galaxies with SFR measurements.The downward arrows indicate upper limits.The black solid line and the shaded region represent the star-forming main sequence relation obtained by

Table 2 .
Correlation between galaxy properties  Spearman rank correlation coefficient.The -values are < 10 −10 for all cases.

Table 3 .
Best-fit parameters of multi-component fit to the   ,500 −profile Left:The effect of  * on the   ,500 -profile.The redshift is controlled for galaxies in a given / vir bin.We first divided the galaxy sample into two bins with / vir > 1 and < 1. Galaxies in each / vir bin are split by log 10 (  * /M ⊙ ) < 9, > 10, and 9 <log 10 (  * /M ⊙ ) < 10.Redshift-controlled galaxies are selected in each  * sample for a given / vir bin to produce the Ly stacks.  ,500 of median stacks are shown by blue, black, and red open squares for low-, intermediate-, and high-mass bins as a function of / vir .Measurements for the MUSEQuBES subsample only are shown with lighter shades and different markers (one bin is made by combining intermediate-and high-mass galaxies due to the lack of log 10 (  * /M ⊙ ) > 10 galaxies in MUSEQuBES sample).The data points inside the virial radius are plotted with offsets in the abscissa for clarity.The error bars are similar to Fig 6.Right: The Ly   ,500