Characterizing cool, neutral gas and ionized metals in the outskirts of low-z galaxy clusters

We present the first statistical detection of cool, neutral gas in the outskirts of low-z galaxy clusters using a sample of 3191 z $\approx$0.2 background quasar - foreground cluster pairs with a median cluster mass of $\sim 10^{14.2}$ M_sun by cross-matching the Hubble Spectroscopic Legacy Archive quasar catalog with optically- and SZ-selected cluster catalogs. We detect significant Lya, marginal CIV, but no OVI absorption in the median stacked spectra with rest-frame equivalent widths (REWs) of 0.043$\pm$0.006A, 0.020$\pm$0.007A, and<0.006A (3$\sigma$) for our sample with a median impact parameter ($\rho_{cl}$) of $\approx$5 Mpc (median $\rho_{cl}$/$R_{500}$ $\approx$7 ). The Lya REW shows a declining trend with increasing $\rho_{cl}$ ($\rho_{cl}$ / $R_{500}$) which is well explained by a power-law with a slope of -0.74 (-0.60). The covering fractions measured for Lya, CIV and OVI in cluster outskirts are significantly lower compared to the circumgalatic medium (CGM). We find that the CGM of galaxies residing in cluster outskirts is considerably deficient in neutral gas compared to their field counterparts. This effect is more pronounced for galaxies that are closer to cluster centers or that are in massive clusters. We argue that the cool gas detected in cluster outskirts arises from the circumgalactic gas stripped from cluster galaxies and to large-scale filaments feeding the clusters with cool gas.


INTRODUCTION
According to the standard structure formation model, galaxy clusters, being the most massive structure in the universe, mark the nodal points in the cosmic web where several filamentary strands intersect. The inflow of cosmic matter via these filaments feeds the growth of galaxy clusters. The infalling gas in the intracluster medium (ICM i.e., < 500 1 ), having been shock-heated to a very high temperature of ∼ 10 7−8 K (Davé et al. 1999;Voit 2005), is bright enough to be directly mapped in X-ray observations (Urban et al. 2014;Simionescu et al. 2015;Biffi et al. 2018). This hot phase of gas in the ICM is thought to account for 80% of the baryonic content in the clusters. On the other hand, clusters are constantly growing and evolving in their outskirts (> 500 ) as a result of a succession of galaxy mergers and the accretion of infalling gas from the intergalactic medium (IGM). Interestingly, despite the fact that a significant portion of the IGM is in the cool/warm phase ( ∼ 1 Radius within which the mean mass density of a cluster is 500 times the critical density of the universe. Similarly, 200 corresponds to the radius at which the mean mass density of a cluster is 200 times the critical density of the universe. 10 4−5 K; Davé et al. 2010;Kravtsov & Borgani 2012), the nature of this gas phase remains poorly understood in the outskirts of clusters. Given the lack of sensitive X-ray diagnostics for directly probing this gas in emission, absorption line spectroscopy of UV-bright background quasars can be leveraged as an ideal alternative to study this otherwise invisible yet crucial phase of the cluster outskirts.
There are a handful of studies focusing on the distribution of neutral gas traced by Ly but using limited numbers of quasar sightlines, primarily using Hubble Space Telescope Cosmic Origin Spectrograph ( /COS) spectra of background quasars (e.g., Yoon et al. 2012;Tejos et al. 2016;Muzahid et al. 2017;Yoon & Putman 2017;Burchett et al. 2018), but see Lanzetta et al. (1996); Tripp et al. (1998); Miller et al. (2002) for pre-COS studies. Using 23 quasar sightlines in the background of the Virgo cluster, Yoon et al. (2012) found that Ly absorbers avoid the hot ICM and are more abundant in the outskirts. Based on the concomitant occurrence of these absorbers with H i-emitting substructures in their study, the authors posited that the warm gas is tracing the large-scale structure (LSS). In a subsequent study of 29 and 8 quasar sightlines passing through the Virgo and Coma clusters respectively, Yoon & Putman (2017) concluded that Ly absorbers are more prevalent between 1−2 vir 2 distance from the Virgo cluster center, while no such trend is observed for the Coma cluster. Burchett et al. (2018) also studied the outskirts (0.2−2.4 R 200 ) of 5 X-ray-selected clusters. They reported seven (H i) > 10 13 cm −2 Ly absorbers in 4/5 clusters. Detection of significantly stronger H i absorbers (i.e., (H i) > 10 16.5 cm −2 ) is reported in the outskirts (1.6−4.7 500 ) of 3/3 SZ-selected clusters studied by Muzahid et al. (2017).
Recently, Mishra & Muzahid (2022), using a large sample of ≈ 80, 000 quasar−cluster pairs at ≈ 0.5 from the Sloan Digital Sky Survey (SDSS), reported a detection of significant Mg ii (7 ) and marginal Fe ii (3 ) absorption in the mean and median stacked spectra of the quasars. From the density and metallicity constraints in their study, the authors suggested that the absorption signal is likely originating from the stripped gas from infalling galaxies. The authors further argued that the stripping process can be effective up to a distance of ≈ 2.4 Mpc (≈ 3.6 500 ) from the clusters. Additionally, Anand et al. (2022) reported detection of Mg ii absorption in the outskirts of clusters from the Dark Energy Spectroscopic Instrument (DESI) survey. They reported a covering fraction (CF) of 1-5% within 500 for W (2796) > 0.4 Å, and concluded that the Mg ii absorption likely stems from stripped interstellar medium (ISM) and/or satellite galaxies, based on the lack of correlations between the absorbers and properties of nearest cluster galaxies within 200 .
The nature of neutral gas in the cluster outskirts has also been explored in a handful of theoretical studies. For example, Emerick et al. (2015), using hydrodynamical simulations, examined the distribution of Ly absorbers around a Virgo-like and a Coma-like cluster. The authors found that the majority of their fast-moving low column density Ly absorbers in the outskirts have filamentary origin. In the inner region, however, the increased column density and metallicity of the Ly absorbers imply that gas from galaxies has been stripped away. To investigate the interplay between the ICM, cluster outskirts, and circumgalactic medium (CGM) of cluster galaxies, Butsky et al. (2019), using a high-resolution hydrodynamical simulation, mapped the distribution of cool and warm gas around a ∼10 14 M ⊙ cluster. They found that, compared to the ICM, the cluster outskirts are more multiphased and richer in cool/warm gas. Due to the inadequate mixing of the ICM gas with the stripped gas from galaxies, the metallicity of the cool/warm gas phase in the outskirts shows a huge scatter compared to the metallicity of the hot ICM gas. In addition, they found the signature of stripping of the CGM of cluster galaxies out to ≈ 4 200 .
Cluster outskirts are ideal test-beds to study the environmental effects of cluster galaxies. It is well established both from observation and simulations that galaxies in dense environments, such as groups and clusters, are gas-deficient with elliptical morphology and redder colors than their field counterparts of comparable stellar mass (Davies & Lewis 1973;van den Bosch et al. 2008;Wetzel et al. 2012;Bahé et al. 2013;Fossati et al. 2017;Davies et al. 2019;Hough et al. 2023;Kim et al. 2023;Rohr et al. 2023). However, the effects of cluster environments on the CGM of galaxies, which is relatively loosely bound to galaxies as compared to the disc and the ISM, are not well explored observationally with statistically significant samples. Yoon & Putman (2013) and Burchett et al. (2018) reported a lower CF of Ly absorbers in the CGM of cluster galaxies compared to the field galaxies up to an impact parameter of 500 kpc from the galaxies. This trend has even been observed on group scales, with galaxies in groups showing an intermediate CF between cluster and field galaxies (Burchett et al. 2018). Interestingly, the simulation study by Bahé et al. (2013) has revealed that the ram pressure exerted by the extended halos of clusters can efficiently strip the hot gas from the infalling galaxies up to ≈ 5 200 . The authors found that other processes such as "overshooting" (Gill et al. 2005) and"pre-processing" (McGee et al. 2009) also play a role in stripping, depending on the mass of the infalling galaxy, its distance from the cluster core, and the cluster mass.
Here we present the first systematic study to probe and characterize the cool, neutral gas and ionized metals in the outskirts of low-clusters using a statistically significant sample of clusters. The purpose of this study is twofold: (i) to map the distribution of cool gas and metals surrounding clusters out to 10 500 , and (ii) to investigate the environmental effects on the CGM of cluster galaxies. We employ spectral stacking analysis to determine the average absorption strengths of the lines of interest. The major advantage of using spectral stacking is that this approach does not rely on any linking velocity to associate an absorber with a cluster, but at the expense of gas kinematics. However, in order to determine the covering fraction, we do use a linking velocity of ±500 kms −1 (as generally used in the literature) for the visually identified absorbers.
The paper is structured as follows: In Section 2, we discuss the quasar and cluster samples used in this study, followed by the construction details of the quasar-cluster pairs and the continuum normalisation of the HST/COS spectra. The key results are provided in Section 3. In Section 4, we discuss the nature and origin of H i gas and metals detected in this work followed by the conclusion in Section 5. We used a flat ΛCDM cosmology with 0 = 71 km s −1 Mpc −1 , Ω M = 0.3, and Ω Λ = 0.7.

Cluster sample
To build a statistically significant sample of quasar−cluster pairs suitable for detecting the Ly and UV metal absorption lines, we use seven cluster catalogs presented in (i) Bleem et al. (2015, hereafter BL15), (ii) Wen & Han (2015, hereafter WH15), (iii) Wen et al. (2018, hereafter WHF18), (iv) Bleem et al. (2020, hereafter BL20), (v) Huang et al. (2020, hereafter H20), (vi) Hilton et al. (2021, hereafter H21), and (vii) Zou et al. (2021, hereafter Z21). As different algorithms are used to identify clusters in these catalogs, their identification criteria have varying limits on redshift, mass, and signal-to-noise ratio (SNR) for which their identification algorithms are complete. Consequently, we have only considered clusters from each catalog that meet their respective completeness criteria. The details on mass, redshift range covered, and completeness criteria of these catalogs are given in Appendix A1. In addition, we restrict our sample to clusters from these catalogs that have a known spectroscopic redshift for the brightest member galaxy.
Merging these seven cluster catalogs resulted in a total of 247,844 spectroscopically confirmed clusters. Similar to Mishra & Muzahid (2022), to eliminate repeated entries with somewhat different redshifts and sky positions within these clusters, we flag the clusters that met the following criteria: (i) two clusters with a velocity offset < 1000 kms −1 and (ii) physical separation is less than the sum of their 500 values 3 . We consider the most massive one in our analysis if two or more clusters satisfied these two conditions.

Quasar sample
We use the UV quasar catalog from the Hubble Spectroscopic Legacy Archive (HSLA) survey Data Release 2 (Peeples et al. 2017). This catalog contains spectra of 799 quasars observed with /COS. We select 583/799 quasars whose spectra are available in median resolution ( ∼ 18,000) G130M and/or G160M gratings. The spectra obtained with G130M and G160M cover a useful spectral range between 1050−1450Å and 1400−1800Å respectively ,which are suitable for detecting the redshifted Ly and UV metal absorption features. The SNR per resolution element of these quasars varies across a range of 4 to 19 with a median SNR of 8.

Quasar-cluster pairs
We cross-match the catalog of 583 background quasars with the 186,378 galaxy clusters with spectroscopic redshifts. We impose two selection criteria for our initial quasar-cluster pairs: (i) The line-of-sight (LOS) velocity offset between the quasar and cluster redshifts is > 5000 kms −1 in order to minimize possible contamination due to absorption intrinsic to the background quasar and/or the quasar host-galaxy (see e.g., Muzahid et al. 2013). (ii) The projected separation between the foreground cluster and the background quasar at the redshift of the cluster ( cl ) should be < 10 500 , suitable to probe the cluster outskirts beyond several virial radii.
Our search after imposing the above-mentioned conditions yielded a total of 3230 quasar-cluster pairs with 2808 'unique' clusters and 431 'unique' quasars. 4 In our subsequent analysis, we found broad absorption line (BAL) features associated with the 4 Many of the quasars are probing multiple clusters while some clusters are probed by more than one quasar. quasar emission in 7 quasar sight-lines, accounting for a total of 21 quasar-cluster pairs. In addition, one quasar sight-line accounting for 18 quasar-cluster pairs was unusable due to a Lyman Limit System (LLS) at z = 0.812. We exclude these 39 quasar-cluster pairs from our sample. Therefore, our final sample consists of 3191 quasar-cluster pairs with 2785 and 423 unique clusters and quasars, respectively.
In Fig. 1, we show the scatter plot of cluster mass ( 500 ) vs redshift ( cl ) for the 2785 clusters (left panel) and cl / 500 vs cl for the 3191 quasar-cluster pairs (right panel). The cluster redshifts range from 0.01 to 0.76 with a median of 0.19. The 500 values range from 0.2 − 12.9 × 10 14 M ⊙ with a median of 1.3 × 10 14 M ⊙ . The normalized cluster impact parameter ( cl / 500 ) of our quasarcluster pairs ranges from 0.1 − 10.0 by design, with a median value of 7.0. Finally, the median cluster impact parameter of the quasarcluster pairs in our sample is 4.8 Mpc.

Continuum normalization
191 of the 423 quasars have both G130M and G160M grating spectra, while 198 (34) quasars are observed only with the G130M (G160M) grating. The G130M and G160M grating spectra of a given quasar are simply joined together from the wavelength at which the average SNR per pixel from the two gratings becomes roughly equal. Before continuum normalization, we exclude the spectral region from 1210−1220 Å and 1301−1307 Å to avoid the geocoronal Ly and O i emissions, respectively. In addition, we exclude the ±100 kms −1 region around the known strong galactic absorption lines. We bin each quasar spectrum by 3 pixels using the Python routine SpectRes (Carnall 2017). For continuum normalization of individual quasar spectra, we adopt a similar approach outlined in Mishra & Muzahid (2022).
Briefly, we first smooth each spectrum by 111 pixels to generate the pseudo-continuum level. This pseudo-continuum is then used to normalize each spectrum. We perform iterative boxed sigmaclipping with asymmetric sigma levels on the pseudo-continuumnormalized spectrum to ensure efficient clipping of absorption features while retaining the residual emission lines. The number of boxes is chosen based on the presence/absence of strong emission lines, while the asymmetric sigma levels are determined based on the median SNR of each spectrum. The iterative clipping is performed until the residual spectrum is sufficiently free from any absorption features. The spectrum is then interpolated linearly over the clipped spectral regions. We finally fit a spline over this interpolated spectrum. The knots of the spline are optimized based on the presence of prominent emission lines (such as the Ly , O vi, Ly , and C iv) and the SNR of each spectrum. This is done to take into account the larger curvature near the emission lines while preventing over-fitting in noisy spectral regions.

Spectral stacking of the full sample
To detect and analyze diffuse, cool/warm ( ∼ 10 4−6 K) gas in the outskirts of the clusters, we employ the spectral stacking technique (see Mishra & Muzahid 2022). For each quasar-cluster pair, we first shift the quasar spectrum to the rest-frame of the cluster by dividing the observed wavelengths by (1 + cl ). We only consider spectral regions between 1135−1450Å and 1400−1790Å for the G130M and G160M gratings, respectively, in the observed frame of each quasar. Whenever present, the spectral region blue-ward of the Lyman continuum break of a Lyman limit system is excluded. We only select those quasar-cluster pairs for our stacking analysis where the quasar spectra contribute at least 5 pixels with SNR per pixel > 1 within a velocity window of ±500 kms −1 centered on a given line (i.e., Ly , O vi 1031, and C iv 1548). We then calculate the median fluxes in bins of 50 kms −1 . We confirm that our results are not sensitive to the bin size used by varying the bin size to 100, 150, and 200 kms −1 .
The median composite spectra of Ly , O vi 1031, and C iv 1548 are shown in Fig. 2. The number of quasar-cluster pairs contributing to the Ly , O vi, and C iv stacks are 2146 (1828 unique clusters and 404 unique quasars), 1765 (1686 unique clusters and 349 unique quasars), and 688 (594 unique clusters and 201 unique quasars), respectively. To estimate the rest-frame equivalent width (REW), the stacked spectra are further normalized by the corresponding pseudo-continua estimated by fitting first-order polynomials to the line-free regions. As can be seen from the bottom panel of Fig. 2, the Ly line is detected at more than 99% confidence level (CL) in the stacked spectrum with REW, 1215 = 0.043 ± 0.006 Å. The uncertainty in the 1215 is estimated by quadratically adding the statistical uncertainty calculated from the stacks of 200 bootstrap realizations 5 of the 2146 quasar−cluster pairs and the uncertainty in the pseudo-continuum placement. In the linear part of the curve of growth (COG), the measured 1215 corresponds to a column density of 10 12.9 cm −2 . We fit a single component Gaussian to the stacked Ly absorption. The best-fit Gaussian has a velocity centroid of 13 ± 48 kms −1 and a velocity dispersion ( v, H i ) of 355 ± 60 kms −1 . The errors in the velocity centroid and velocity dispersion values are determined from the standard deviations of these values obtained from the fitting of the absorption profiles corresponding to the 200 bootstrap realizations.
Next, we detect marginal C iv absorption at CL of more than 95% with 1548 of 0.020 ± 0.007 Å, despite the fact that the number of quasar−cluster pairs contributing to the C iv stack profile is 5 We confirm that increasing the realizations to 500 or 1000 does not alter the results or conclusions drawn in this study. significantly smaller (i.e., 688). In the linear part of the COG, the 1548 corresponds to a column density of 10 12.7 cm −2 . We fit the C iv doublet with a single Gaussian component. The velocity width of the two doublet Gaussians are kept the same but the amplitudes are allowed to vary. This yield a velocity centroid and velocity dispersion of −9 ± 84 kms −1 and 145 ± 90 kms −1 , respectively. The fit suggests that the C iv absorption is partially saturated. While the velocity centroid is consistent with 0 kms −1 , the velocity dispersion is ≈ 2 times narrower than that of the Ly , suggesting that the metal-bearing gas correlates over smaller velocity scale compared to the Ly .
Finally, no significant O vi absorption is detected in the high SNR composite spectrum shown in the top panel of Fig. 2. We obtained a 3 upper limit on the O vi REW ( 1031 ) of 0.006 Å, assuming that the undetected line is spread over 31 pixels corresponding to the full-width at tenth maximum (FWTM) of the Ly line (≈ 4.29× v,H i ≈ 1523 kms −1 ). In the linear part of the COG, the 1031 upper limit corresponds to a column density of < 10 12.7 cm −2 . We also put constraints on the upper limits of the REWs of other ions such as C ii, C iii, Si ii, Si iii, and Si iv, based on the non-detection of signals in the composite spectra. Table 1 summarizes the results obtained from the stacking analysis.

The Ly equivalent width-profiles
To construct the Ly REW-profile, we split the sample of 2146 quasar−cluster pairs contributing to the Ly stack into five bins of cl and cl / 500 using binsize of 1.5Mpc and 2 respectively. Table A1 summarises the details of the measurements performed on the Ly subsample stacks. Fig. 3 shows the Ly REW-profiles as a function of cl (left panel) and cl / 500 (right panel). A declining trend in 1215 with increasing cl and cl / 500 is evident in Fig. 3. As can be seen from Fig. 3, a single power-law can explain the REW-profiles adequately with power-law slope of −0.74±0.17 for the 1215cl profile (left panel) and −0.60±0.15 for the 1215cl / 500 profile (right panel). We do not find any correlation between the line width ( v, H i ) of the Ly absorption signal with the cl or cl / 500 . Besides, we investigate the dependence of 1215 on redshift and cluster mass by dividing the sample into three redshift and mass bins while maintaining a similar distribution of cl / 500 in all three bins of redshift and mass. We find no significant dependence of 1215 on redshift. However, we find a tentative positive correlation between 1215 and cluster mass, but the data points in the three mass bins are consistent within 1 . Given the small number of quasar−cluster pairs for C iv stack and the absence of significant signal in the O vi stack, we could not perform similar analysis for the C iv and O vi stacks.

Covering fraction analysis
In this section we determine the covering fractions (CFs) of the commonly detected UV transitions such as the Ly , C iv, and O vi in the outskirts of the galaxy clusters by visual inspection of the relevant parts of the quasar spectra. The CF of a given transition is defined as: Here, det is the number of quasar-cluster pairs for which the given line is detected with a REW ( ) more than a threshold equivalent width ( th ), and tot is the total number of quasar−cluster pairs Relative Los Velocity ( 10 3 km s 1 ) Normalized flux Number of quasar-cluster pairs contributing to the stacked spectrum, median cluster redshift and rest-frame equivalent width (a 3 upper limit for O vi) estimated from the stacked profile within ± 500 kms −1 and ± 300 kms −1 for Ly and C iv respectively around the line centroid is indicated in each panel. The best-fit single Gaussian component for the Ly and C iv doublet absorption profiles are shown in red and the corresponding velocity dispersion from the fits are also indicated.  Notes -(1) Name of the species around which stack profile is generated.
(2) Number of quasar-cluster pairs. (3), (4), (5), (6), and (7) median values of cluster redshift, 500 , 500 , cl , and cl / 500 respectively. The values in the parenthesis for the parameters listed from (3) − (7) indicate the 16 and 84 percentiles of the parameter distribution. (8) REWs measured within ±500 kms −1 and ±300 kms −1 for Ly and C iv respectively and upper limits on the REWs for other ions are estimated from the median stacked spectra. (9) & (10) Velocity dispersion and Velocity centroids obtained from Gaussian fitting as explained in Section 3.1 in the median stacks of Ly and C iv. for which the quasar spectra are sensitive to detect the th (i.e., the 5 limiting equivalent width, lim , is lower than the th ). For each transition, we calculated the lim over a line-free region within ±500 kms −1 of the cluster redshifts using equation 6 of Hellsten et al. (1998).
For visual identification of Ly absorption systems, we start with the 2146 quasar-cluster pairs used for Ly stacking (see Table 1). To confirm the presence of Ly , we only select the pairs for which at least Ly , among the higher order Lyman-series lines, is covered. We point out that the requirement of the presence of Ly and/or higher-order lines, will miss out the weak, stand-alone Ly absorbers for which the higher order lines, including Ly , is too weak to be detected. Nonetheless, as identifying all/most of the absorption lines in so many quasar spectra is beyond the scope of this paper, we adopt this approach to avoid false positives. We do not impose any constraints on the detection significance for the Ly or higher-order Lyman-series lines. This reduces our sample of 2146 quasar−cluster pairs to 1160 pairs in the foreground of 295 quasars. We only search a ±500 kms −1 spectral region 6 around the cluster redshifts for identifying Ly absorption lines. We flag our detected Ly absorption systems in two categories (Flag-1 and Flag-2) based on the confidence level of the detection. Flag-1 represents the high confidence systems where the Ly absorption is accompanied by at least one higher-order Lyman series line and at least one metal line (e.g., C ii, C iii, C iv, Si ii, Si iii, Si iv, and O vi). The shape of the apparent column density profiles (ACD; using equation 8 from Savage & Sembach 1991) of the Lyman-series lines are also checked for consistency for this class. For Flag-2, the Ly systems satisfy one of the following two conditions: (1) have only associated Ly absorption with consistent column density profiles (2)  are not fully consistent with each other, indicating blends. In total, we identify 241 Ly systems, with 128 Flag-1 and 113 Flag-2 absorbers (in 218 unique clusters towards 118 unique quasars). Taking into account both Flag-1 and -2 systems (only Flag-1 systems), the CF of Ly is 0.21 +0.01 −0.01 (0.11 +0.01 −0.01 ) for th = 0.1 Å, at the median cl / 500 ≈ 7.1.
Among the 1160 quasar−cluster pairs probed by 295 quasars, we find that 128 clusters probed by 58 quasars have redshifts consistent within ±500 kms −1 of each other, but at different impact parameters. To avoid association of an absorption system with multiple clusters at varying impact parameters, we remove these 128 clusters and estimate the CF of the Ly for the remaining 1032 quasar−cluster pairs. The CF of this sample is 0.21 +0.02 −0.01 (0.12 +0.01 −0.01 ) for Flag-1 and -2 systems (only Flag-1 systems) for th = 0.1 Å which is consistent with the CF estimates for the full sample of 1160 quasar−cluster pairs.
Out of the 688 quasar−cluster pairs used for C iv stacking, we only consider the 649 pairs for visual inspection of C iv absorption systems, for which the quasar spectra cover the C iv 1548,1550 lines simultaneously. Similar to Ly , we only use the regions within ±500 kms −1 around the cluster redshifts. We also assign two flags to the identified C iv absorption systems based on the detection confidence. Flag-1 is assigned to the C iv absorption systems that have two lines of the C iv doublet with consistent ACD profiles, in addition to at least one Lyman series line or one other metal line. Systems containing only the two lines of the C iv doublet with consistent ACD profiles or the two lines of the C iv doublet with slight mismatching ACD profiles in addition to the presence of other lines are marked as Flag-2. We identify 49 C iv systems with 38 Flag-1 and 11 Flag-2 categories in the outskirts of 45 unique clusters (34 unique quasars). The CF of C iv including Flag-1 and -2 systems (only Flag-1 systems) with th = 0.05 Å is 0.10 +0.03 −0.04 (0.09 +0.04 −0.03 ) at the median cl / 500 ≈ 6.1. Of the 649 quasar−cluster pairs examined for C iv, we find  47 clusters close in redshift within ±500 kms −1 of each other towards 22 quasars but at different impact parameters. Similar to Ly , after excluding these 47 clusters, the C iv CF of the remaining 602 quasar−cluster pairs with Flag-1 and -2 absorption systems is 0.12 +0.03 −0.04 for th = 0.05 Å which is fully consistent with the CF of the full sample of 649 quasar−cluster pairs. Likewise, for the visual identification of O vi absorption systems, from among the 1765 quasar−cluster pairs used in O vi stacking, we only choose 1553 pairs that have simultaneous spectral coverage of the O vi 1031,1037 lines within ±500 kms −1 of the clusters' redshifts. The identified O vi systems are classified into Flag-1 and Flag-2 categories using the same conventions used for C iv. We identify 114 O vi 1031 systems (78 Flag-1 and 36 Flag-2) in the outskirts of 114 unique clusters towards 82 unique quasars. The CF of O vi with th = 0.05 Å is 0.10 +0.01 −0.02 (0.07 +0.01 −0.02 ) for Flag-1 and Flag-2 systems (only Flag-1 systems) at a median cl / 500 ≈ 6.9.
Of the 1553 quasar-cluster pairs searched for O vi, 120 clusters probed by 52 quasars have comparable redshifts within ±500 kms −1 , but different impact parameters. Similar, to Ly and C iv, after excluding these 120 clusters, the CF of the O vi for the remaining 1433 quasar−cluster pairs is 0.11 +0.01 −0.02 for Flag-1 and -2 systems with W th = 0.05 Å which is in agreement with the CF of the full sample of 1553 quasar−cluster pairs. Table A2 summarizes the results of the covering fraction analysis. In Fig. 4, the covering fraction profiles (CF-profiles) as functions of cl (top panel) and cl / 500 (bottom panel) for Ly , C iv, and O vi absorption systems are shown in the left, middle, and right panels, respectively, for three different threshold rest equivalent widths ( th ). The CF-profiles in Fig. 4 are based on the samples that include both Flag-1 and Flag-2 absorption systems. The corresponding CF-profiles only for the Flag-1 systems are shown in Fig. A1. As evident from Fig. 4, we find a moderate decreasing trend in CF with increasing cl and cl / 500 for the Ly and O vi while no trend is seen for the C iv CF-profile. The covering fraction decreases with decreasing the sensitivity limit (i.e. increasing the th ). As discussed in Section 4.2 and 4.3, the Ly , C iv, and O vi covering fractions measured in the outskirts of low-galaxy clusters are significantly lower compared to the measurements in the CGM of galaxies.

Probing the CGM of cluster galaxies
To explore the influence of cluster environments on the CGM, for each cluster, we search for galaxies within impact parameters ( gal ) of < 300 kpc from quasars and with spectroscopic redshifts within ±1000 kms −1 (Δ ) of the cluster redshift. Note that these galaxies in the outskirts of clusters may or may not be member galaxies, but it is foreseen that they will eventually be part of the cluster. Hereafter, we will call them 'cluster galaxies' for simplicity.
To identify the cluster galaxies, we use the catalog 7 by the Max Planck Institute for Astrophysics -Johns Hopkins University (MPA-JHU) team. This catalog provides information on the star formation rate (SFR) and stellar mass of approximately 1.8 million galaxies with redshift < 0.33 from SDSS DR8. The details of the stellar mass and SFR determination are explained in Brinchmann et al. (2004) and Kauffmann et al. (2003). We use the All Cluster Galaxies SFR-matched Cluster Galaxies Non-cluster Galaxies  Notes -(1) Galaxy sample.
CasJobs SDSS SkyServer 8 to obtain the galaxy properties from the catalog. We only select galaxies with reliable spectroscopic properties with "RELIABLE! = 0" flags. We only consider galaxies for which both stellar mass and SFR information are avail-8 https://skyserver.sdss.org/casjobs/ able in the catalog. For 94 clusters in the foreground of 83 unique quasars, we find 179 galaxies with gal < 300 kpc and |Δv| < 1000 kms −1 . We refer to this sample of 179 cluster galaxies with spectroscopic redshifts as the CLCGM1000 sample. The median redshift and stellar mass of the galaxies in the CLCGM1000 sample are gal ≈ 0.04 and log ( ★ / ⊙ ) ≈ 9.9 respectively. Similarly the median star formation rate (SFR) and specific SFR (sSFR) are log [SFR/M ⊙ yr −1 ] ≈ −0.9 and log [sSFR/yr −1 ] ≈ −10.76, respectively.
Next, we divide the CLCGM1000 sample into two bins of normalized clustocentric impact parameter ( cl / 500 ) and two bins of cluster mass ( 500 ) based on the corresponding median values. The median cluster masses of the two cl / 500 -bins are consistent with each other within 1 . The median cl / 500 values of the two mass-bins are also consistent within 1 . We generate median stacked Ly profiles in the rest-frame of the galaxies contributing to each bins of cl / 500 and 500 . The Ly REWs obtained from the stacking as a function of cl / 500 and 500 are shown with blue stars in the left and right panels of Fig. 5, respectively. Table 2 presents the summary of the measurements obtained from this analysis. It is evident from the left panel of Fig. 5 that the CGM of galaxies is relatively gas-poor when they are closer to clusters ( cl / 500 < 4). In addition, the right panel shows that the CGM of galaxies is relatively gas-poor when they are in the outskirts of massive clusters ( 500 > 1.2 × 10 14 M ⊙ ). The measured Δ 90 9 values suggest that the differences in REWs in the comparing bins are not owing to the significantly different line widths but due to the difference in optical depths. If at all, the Δ 90 values are somewhat higher for the bins showing lower REWs.
Even though the median stellar mass and SFR of the galaxies are consistent within 0.3 dex for the two bins of cl / 500 and 500 , a two-sample Kolmogorov-Smirnov (KS) test suggests a marginal difference in the SFR distribution of the galaxies in the two bins of cl / 500 and 500 . 10 A strong correlation between Ly REW and SFR within the virial radius of (non-cluster) galaxies is recently reported by Dutta et al. (2023). We thus generated SFR-matched subsamples with a tolerance of 0.1 dex for the above exercises. The measurements corresponding to the SFR-matched subsamples are shown by the open magenta circles in Fig. 5. The figure shows that even after matching the SFR of the galaxies, the trends persist. Consequently, it can be inferred that the trends are governed by external environments rather than internal galactic processes.
Further, to compare the REW of Ly absorption obtained for the cluster galaxies with that of non-cluster galaxies, we derive a sample of 58 non-cluster galaxies from Dutta et al. (2023) that are matched in SFR, ★ , galactic impact parameter ( gal ), and redshift of the cluster galaxy sample. The REW of Ly estimated from the median stacked profile in the rest-frame of the non-cluster galaxies is 1215 = 0.390 ± 0.060 Å, which is shown by the red dotted line in each panel of Fig. 5. It is evident from Fig. 5 that the CGM of cluster galaxies that are closer or reside in massive clusters is considerably deficient in gas as compared to non-cluster galaxies. However, the difference in gas content decreases as one moves further away from the clusters or towards low-mass clusters. These observations indicate that the CGM is strongly influenced by the large-scale environments of galaxies.

DISCUSSION
Here we report on the first statistical detection of cool, neutral gas traced by H i 1215 at more than 99% CL along with a marginal detection (> 95% CL) of C iv in the median stacked spectra of 2146 and 688 quasar−cluster pairs, respectively, at a median 9 The velocity difference between the pixels where optical depth is 5% and 95% of the total optical depth. 10 With the null probability that the two distributions are drawn from the same parent population being ≲ 5%. cl / 500 ≈ 7. We obtain an upper limit of ∼ 0.006Å (3 ) for O vi based on the non-detection in the high-SNR composite spectrum of 1765 quasar−cluster pairs. Furthermore, we determine the CF of 0.21 +0.01 −0.01 for th = 0.1 Å for Ly , 0.10 +0.03 −0.04 for th = 0.05 Å for C iv, and 0.10 +0.01 −0.02 for th = 0.05 Å for O vi in the cluster outskirts (median cl / 500 ≈ 7). We also find that the CGM of cluster galaxies exhibits a significant environmental dependence. In this section, we discuss the distribution and origin of the cool, neutral gas and metals detected in the cluster outskirts, and the role of environment in determining the gas reservoirs in the CGM of galaxies.

Distribution of neutral gas in cluster outskirts: Ly REW-profile
The distribution of cool, neutral gas as a function of impact parameter for low-mass haloes ( vir ≲ 10 13 M ⊙ ) has been studied extensively in the literature (e.g., Tumlinson et al. 2013;Wilde et al. 2021;Dutta et al. 2023), with several studies indicating a decline in the Ly REW as the impact parameter from the hosting haloes increases. At high mass scales (> 10 14 M ⊙ ) of galaxy clusters, on the other hand, observations and simulations of the hot ICM show that the temperature, density, and metal content of hot gas vary significantly from the core to the outskirts (Nagai & Lau 2011;Simionescu et al. 2011;Urban et al. 2017), which may influence the distribution of cool, neutral gas surrounding clusters. For instance there are studies suggesting that the cool gas traced by the Ly absorbers is shock heated to high temperatures in the ICM, avoiding the inner region of the clusters (see Yoon et al. 2012). However, there are observations of cool gas in the hot ICM. For example, the detection of extended H emission surrounding the galaxy NGC 1275, located at the center of the Perseus cluster (see Conselice et al. 2001). The H emitting gas may have originated from the cooling of the ICM, as suggested by Fabian (1994). Additionally, simulations of galaxy clusters revealed that fast, radiatively cooling outflows can also cause condensation of cool gas (e.g., see Qiu et al. 2021). Therefore, the REW-profile of Ly as a function of cluster impact parameter, which has not been explored thoroughly using statistically significant samples of clusters, is crucial. Using 23 background quasar sightlines, Yoon et al. (2012) found that Ly absorbers avoid the central X-ray emitting regions of the Virgo cluster. Moreover, relatively stronger Ly absorbers (REW > 0.3 Å) are found to be associated with Virgo substructures. Excluding the Ly absorbers associated with the substructures, they reported a weakly increasing Ly REW with an increasing clustocentric impact parameter, which is in contrast with CGM observations. A similar analysis for the Coma, using 8 background quasars, also revealed no significant trend between Ly REW and impact parameter (see Yoon & Putman 2017). Similarly, Muzahid et al. (2017) found no significant trend between the H i column density and normalized cluster impact parameter up to 5 500 (see their figure 3) by combining their results with those of Yoon et al. (2012), Yoon & Putman (2017), and Burchett et al. (2018).
In this study, the Ly REW-profile, constructed using a statistically significant sample of 2146 quasar-cluster pairs, shows a significant decline with the impact parameter (and normalized impact parameter) as reported in the literature for galactic halos. A singlecomponent power-law can adequately explain the REW-profiles (see Fig. 4). 11 The best-fitting relations are: 1215 = (0.069 ± 0.009 Å) × ( cl 2 Mpc Recently, Dutta et al. (2023) showed that the Ly REW-profile for low-galaxies can be decomposed into two components: (i) a Gaussian/log-linear component, prominent only at small scales (≲ 200 ), and (ii) a power-law component, reminiscent of galaxyabsorber clustering, describing the large-scale features. Owing to the smaller sample size, particularly at cl < 500 ( 200 ), and much weaker overall absorption signal compared to the CGM of galaxies (see next section for details), we could not investigate the possible existence of similar characteristics in the Ly REW-profile around low-clusters. Fig. 6 shows the Ly REW-profile plotted against the cluster impact parameter normalized by the 200 12 for our full sample of 2146 quasar−cluster pairs. The solid black line represents the best-fit single-component power-law to our data points. We compare this with the log-linear + power-law (Model 1, dotted magenta line) and Gaussian + power-law (Model 2, dashed blue line) models describing the Ly REW profile for the CGM of ≈ 0.15 galaxies from Dutta et al. (2023). Fig. 6 indicates that the Ly REW profile for the clusters differs from that of the CGM of galaxies within 200 . However, beyond 200 the REW profiles of clusters and galaxies are consistent with each other within the 1 allowed uncertainties, albeit the clusters' measurements are systematically lower. Moreover, we notice that the power-law slopes are (≈ −0.62) also consistent outside the virial radius where the signal is expected to be dominated by large-scale clustering (see Dutta et al. 2023).

Distribution of neutral gas in cluster outskirts: Ly covering fraction
We determine the covering fraction of Ly using a subsample of quasar-cluster pairs with simultaneous coverage of the Ly and Ly lines. We obtained CF = 0.21 +0.01 −0.01 (0.16 +0.01 −0.01 ) for a threshold REW of 0.1 Å (0.3 Å) for the full sample. The Ly CF-profiles obtained for the full sample, shown on the left panels of Fig. 4, show a declining trend with increasing cl and cl / 500 (see also Fig. A1). Overall, the Ly CF remains in the range of 10-25% which is broadly consistent with the predictions from hydrodynamical simulation (see figure 3 of Butsky et al. 2019) over a length scale of 3 Mpc from clusters. The Ly covering fraction increases rapidly to ≈ 1 in the simulation for cl < 500 kpc (≈ 0.7 200 ). We are yet to probe this impact parameter range observationally.
When a LOS velocity of ≈ 1000 kms −1 is used to associate Ly absorbers with Coma, the CF is measured to be ≈ 33% for th > 0.1 Å (≈ 25% for th > 0.3 Å) both within and outside vir (see Yoon & Putman 2017). For the less massive Virgo cluster, the authors reported a marginally (≈ 1 ) higher Ly CF outside the vir (0.38 +0.14 −0.12 for th > 0.1 Å). The covering fraction measurements for Virgo and Coma presented by Yoon & Putman (2017) have large uncertainties, and are broadly consistent with our measurements. Here we caution that our CF measurements could be systematically 11 Suggested by an F-test with > 95% confidence. 12 We use the Python package COLOSSUS (Diemer 2018) to compute the 200 for a given and 500 of a cluster. smaller as we required the presence of Ly (or higher order Lyman series lines) for detection, which will miss the stand-alone, weak Ly absorbers. Next we compare the CF of Ly absorbers in the outskirts of galaxy clusters with the CF measured in the CGM of galaxies. Several studies in the literature have constrained the Ly CF in the CGM of galaxies across a wide range of redshift, stellar mass, SFR, and impact parameter (see Tumlinson et al. 2013;Liang & Chen 2014;Keeney et al. 2018;Wilde et al. 2021). For instance, Tumlinson et al. (2013) found CF of 100% and 75% for star-forming and passive galaxies, respectively, for Ly absorbers with (H i) > 10 13 cm −2 within 0.5 vir . Meanwhile, the CF of Ly absorbers with th > 0.05 Å in Liang & Chen (2014) remained above 60% out to 8 vir . Keeney et al. (2018) reported 50% CF for (H i) > 10 13 cm −2 Ly absorbers in both < * and ≥ * galaxies up to 4 vir . Recently, Wilde et al. (2021) showed that the CF of Ly absorption with (H i) = 10 13−15 cm −2 is nearly 80% over 4 vir from galaxies. Overall, most CGM studies have reported a high CF of >50% for Ly absorbers out to 4 vir . In contrast, our study finds a low CF of Ly absorption systems with W th > 0.1Å (i.e., (H i) > 10 13.3 cm −2 ) ranging from ≈ 25% at < 500 to ≈ 15% at around 6 − 10 500 (see Fig. 4). We argue that the lower CF of the Ly absorbers in the outskirts of clusters as compared to the CGM of galaxies may account for the lower REWs measured for clusters, particularly within 200 , in contrast to the CGM of galaxies, as shown in Fig. 6.

Distribution of metals in the cluster outskirts
We produce stacks of all the prominent metal lines accessible to COS. However, except for a marginal C iv, none of the metal lines are detected (see Table 1). By examining the relevant parts of the COS spectra, we determine the covering fraction of the C iv and O vi, the two most commonly studied metal lines. For C iv, we obtain CF = 0.10 +0.03 −0.04 for th = 0.05 Å. The CF for O vi is found to be similar (i.e., 0.10 +0.01 −0.02 for th = 0.05 Å). The covering fraction of C iv and/or O vi are not well constrained in the earlier studies for the cluster outskirts. Using a sample of only 7 cluster pairs, Tejos et al. (2016) obtained a CF of 0.14 +0.3 −0.1 for O vi with th = 0.06 Å. The CF value is consistent with our measurement. Burchett et al. (2018) found no O vi absorption within ±1500 kms −1 of the 5 clusters in their study. This is also consistent with the low O vi CF measured here.
The distribution of metals traced by C iv and O vi absorption have been explored in the simulation of Butsky et al. (2019). They found CF of C iv within ≈ 2 200 to be < 0.1, whereas for O vi, it is > 0.15. The authors attributed this difference to the lower ionization potential of C iv (48 eV) compared to the O vi (114 eV), which makes it challenging to sustain the triply-ionized state for carbon at temperatures of > 10 5 K. However, for both the ions, the covering fraction increases significantly in the inner 500 kpc regions. The C iv covering fraction we estimate in the impact parameter range of 500-3000 kpc is somewhat higher compared to the simulation prediction (i.e., only a few percent; see figure 3 of Butsky et al. 2019). The O vi CF of ≈ 10%, however, is consistent with our measurements (see Fig. 4).
The non-detection of O vi in the stack of 1765 quasar-cluster pairs is intriguing. It is generally believed that the gas temperature in cluster outskirts is more suitable for the tracers of the warmhot phase, such as the O vi. However, the non-detection of O vi in the high SNR stacked spectra puts a stringent upper limit on the (O vi) of 10 12.7 cm −2 . Combined with the marginal detection of C iv, we obtain an upper limit on the gas temperature of 10 5.3 K, under collisional ionization equilibrium (CIE; Gnat & Sternberg 2007) assuming single-phase gas with solar relative abundances. It indicates that cluster outskirts are not potential reservoirs of warmhot gas.
We point out that the marginal detection of C iv in the stack of 688 quasar-cluster pairs while the non-detection of O vi, albeit with a ≈ 2.5 times larger sample size and a covering fraction similar to C iv, could be due to two factors: (i) The REW of the detected C iv absorption systems tends to be somewhat higher compared to O vi. (ii) The LOS velocity distribution of the detected absorbers shows a somewhat larger scatter for O vi.
It is noteworthy that among the 49 C iv absorption systems, 44 systems show associated Ly absorption (median REW = 0.8 Å), while the remaining 5 C iv systems do not have spectral coverage for Ly . Likewise, all 114 O vi absorbers in our study are associated with Ly (median REW = 0.7 Å) and/or Ly absorption. The association of the strong Ly absorption with the metal lines suggests that a fraction of these systems may be arising from the CGM of cluster galaxies close to the sightline near a background quasars. In the next section we will estimate the contribution of the circumgalactic gas to the stacked signals.

Contribution of the CGM to the stacked signals
In order to understand whether the Ly absorption signals in the stacks are dominated by the CGM of cluster member galaxies, we adopt an approach similar to Mishra & Muzahid (2022). Similar to Section 3.4, for each cluster, we first identify galaxies within a radius of 300 kpc surrounding the background quasar with photometric and/or spectroscopic redshifts consistent with that of the cluster. Since the majority of our clusters are from DESI, we use the photometric redshift catalog of galaxies from the DESI Legacy Imaging Surveys (Zou et al. 2019) to identify galaxies around the clusters. The typical uncertainty of the photo-estimates of the DESI galaxies is 0.017; we therefore use a LOS velocity window of ±4500 kms −1 (≡ × Δ /1 + ) for linking galaxies to a cluster. Next, we eliminate those clusters from our analysis when one or more galaxies are found. For 1140 of the 2146 quasar-cluster pairs, we identify a total of 2251 cluster galaxies with a median ≈ 19. 13 We exclude those 1140 quasar-cluster pairs, for which the background quasars may probe the CGM of cluster galaxies, and call the remaining sample of 1006 (2146 − 1140) quasar-cluster pairs the ExCGM4500 sample. Furthermore, since the LOS velocity of ±4500 kms −1 corresponds to a large cosmological distance of ≈ ±60 Mpc at the median redshift of 0.14, we use one more conservative window of 1000 kms −1 and label the corresponding sample as ExCGM1000. We identify a total of 694 galaxies associated with 476 quasar−cluster pairs within the velocity window of ±1000 kms −1 in ExCGM1000 sample.
The Ly REWs measured for the ExCGM1000 and Ex-CGM4500 subsamples are 0.028 ± 0.006 Å and 0.027 ± 0.008 Å, respectively. Both are significant at the 99% CL, and are consistent with the REW measured for the full sample within the ≈ 1.5 allowed uncertainties. As can be seen from Fig. 7, the Ly REWs measured for the sub-samples after excluding the possible CGM contributions are consistent with the full sample within 1 , indicating that the CGM of bright cluster galaxies in the outskirts does not dominate the Ly absorption signal.

Cosmic filaments as the origin of the stacked signals
To investigate the origin of the observed cool, neutral gas and metals around clusters, we re-examine the samples of 1160 and 649 quasar−cluster pairs used for the covering fraction analysis of Ly and C iv (see Section 3.3). As mentioned earlier, we identify 241 and 49 Ly and C iv absorption systems, respectively, within ±500 kms −1 around 218 and 45 quasar−cluster pairs. We isolate the subsamples of 942 (= 1160 -218) and 604 (= 649 -45) quasar−cluster pairs without any individual absorption detection respectively for Ly and C iv and produce stacks. The top panels of Fig. 8 show the stacked Ly and C iv profiles for the full samples of 1160 and 649 quasar-cluster pairs respectively. Notably, we detect Ly and C iv signals at >99% and 95% CL, respectively. However, no significant absorption is seen for either Ly or C iv when we exclude the pairs with visually identified Ly and C iv absorbers (see the bottom panels of Fig. 8). This indicates that the signals in the full stacks (top panels of Fig. 8) are primarily dominated by the individual absorbers used in the covering fraction analysis. These findings indicate that the distribution of cool gas and metals in the cluster outskirts is patchy, and the detected absorption signals in the stacked spectra may arise from cosmic filaments feeding the clusters with cool gas rather than from a volume-filling gas phase. The observed low covering fractions of gas and metals, particularly in comparison with the CGM of galaxies, are also consistent with this picture. In line with our argument, Wakker et al. (2015), in a study to probe nearby galaxy filaments using Ly absorption, found a strong anti-correlation between Ly REW and filament impact parameter, with absorbers with (H i) > 10 13 cm −2 primarily detected within ≈ 500 kpc from the filament axis (see also Bouma et al. 2021).
We also search for the galaxies around the 241 Ly absorption systems as identified in Section 3.3. We use the similar initial DESI galaxy catalog used for creating samples in Section 4.4. Around each Ly system, we look for galaxies with photometric and/or spectroscopic redshifts within |Δ | < 500, 1000, 2000, and 4000 kms −1 of the absorbers and a projected separation of < 300 kpc. Out of 241 Ly systems, we find 13% (32 / 241), 28% (68 / 241), 39% (94 / 241), and 56% (136 / 241) systems have associated bright galaxies (median ≈ 19) within |Δ | < 500, 1000, 2000, and 4000 kms −1 , respectively. This suggests that most of the detected Ly absorbers are likely unrelated to the halos of bright galaxies, and may instead originate from the cosmic filaments. However, conducting deeper galaxy surveys around the cluster outskirts is crucial to search for the presence of faint, low-mass galaxies around these absorbers. Alternatively, the observed absorption signal could arise from stripped off gas due to multiple environmental processes faraway from any galaxies as discussed below. Fig. 5 shows that the CGM of galaxies located in the proximity of clusters (i.e., cl / 500 < 4) is more gas-deficient compared to their counterparts located farther away and compared to non-cluster galaxies with matched SFR, ★ , gal , and redshift. Similarly, at comparable distances from clusters, the CGM of galaxies residing in more massive halos (i.e., 500 > 1.2 × 10 14 M ⊙ ) is more gasdeprived compared to galaxies in the outskirts of low-mass clusters and compared to non-cluster galaxies. Since the galaxies in the comparable bins are matched in galaxy properties (e.g., SFR, ★ ), we argued that this difference in the gas content is not governed by internal galactic processes but by environmental effects. The dearth of cool, neutral gas in the CGM of galaxies residing close to cluster centres and/or residing in the outskirts of massive clusters can be understood in terms of efficient ram pressure stripping of circumgalactic matter that is loosely bound to the galaxies. Both the density of the ambient medium ( amb ) and the relative velocity of infalling galaxies with respect to the ambient medium ( rel ) are expected to be higher for massive clusters and/or in the inner regions of clusters. This naturally facilitates higher ram pressure ( ram = amb × 2 rel ) stripping for galaxies residing in massive clusters and/or closer to cluster centers. Indeed, the recent simulation by Rohr et al. (2023) revealed that the majority of the jellyfish galaxies lose the bulk of their cold gas (< 10 4.5 K) via ram pressure stripping within 1-2 Gyr since infall when they are at a distance as far as 3 − 4 200 . The authors further showed that satellite galaxies of all types deposit more than 10 10 M ⊙ of cold gas via ram pressure stripping in hosts of mass > 10 13 M ⊙ over a time-scale of 5 Gyr. Such stripped gas faraway from galaxies can explain the observed absorption signals if the gas remains cold.

Stripping as the origin of the stacked signals
The other important environmental process that can cause gas loss closer to cluster centres (< 2 − 3 200 ) is "overshooting" (e.g., Pimbblet 2011;McGee et al. 2009;Bahé et al. 2013;Borrow et al. 2023). It is found in simulation that a substantial fraction of galaxies in cluster outskirts are actually backsplash galaxies, which have undergone a significant gas loss due to severe tidal disruption and stripping caused by ram pressure during their first pericentre passage. The current position of such backsplash galaxies and the locations at which their gas has been stripped in the outskirts can be considerably different. Such "overshot" galaxies can give rise to the signals detected here. This process, nevertheless, is not efficient beyond a few virial radii (see e.g., Bahé et al. 2013).
The results in the right panel of Fig. 5 are intriguing as they indicate that galaxies residing in more massive clusters exhibit a significant deficiency of circumgalactic gas compared to those in less massive clusters, even at distances as far as ≈ 5 500 . The Ly absorption strength for the galaxies residing in low-mass clusters is somewhat lower, albeit consistent with non-cluster galaxies within 1 . The lack of cool gas in the CGM of galaxies residing in high-mass clusters is likely due to "pre-processing". Gas removal through "pre-processing" happens when a galaxy becomes a satellites of a group (substructure) in the outskirts of a massive cluster prior to infall (e.g., Wetzel et al. 2013;Jaffé et al. 2016;Hough et al. 2023). According to the hierarchical structure formation model, more substructures are expected around massive halos, which facilitates efficient "pre-processing" of the CGM of galaxies residing in the outskirts of massive clusters.
It is expected that all of these processes discussed above (ram pressure stripping, overshooting, and pre-processing) will contribute to the cool gas cross-section in the cluster outskirts to a varied degree, depending on the host mass and distance from the core. However, the overall low covering fraction of the neutral gas in cluster outskirts suggests that the stripped circumgalactic gas is either mixing with the ambient hotter medium or is heated up by AGN feedback or that the galaxies supplying cool gas via stripping are predominantly falling onto clusters through large-scale cosmic filaments.

CONCLUSION
By cross-matching the HSLA quasar catalog with seven cluster catalogs from the literature, we have built a sample of 3191 quasarcluster pairs with a median cluster redshift and mass of 0.19 and 1.3 × 10 14 M ⊙ , respectively. The quasars are probing the clusters with impact parameters in the range 0.7-10 Mpc with a median of 4.8 Mpc (median cl / 500 ≈ 7). Using spectral stacking, we present the first statistical detection of cool, neutral gas (traced by Ly ) and metals (traced by C iv) in the outskirts of clusters. In addition, we construct the Ly REW-profile out to 10 500 . We visually identify 241 Ly , 49 C iv, and 114 O vi absorption systems within ±500 kms −1 of the cluster redshifts, which allow us to determine the covering fractions of Ly and the other metal ions. We leverage the combined results of the stacking analysis and the distribution of individual absorption systems to investigate the nature and origin of the cool, neutral gas and metals in the outskirts of the clusters.
The key results of this paper are summarized as follows: • In the median stacked spectra of 2146 and 688 quasar−cluster pairs, we detect significant Ly and marginal C iv absorption with REWs of 0.043 ± 0.006 Å (at > 99% CL) and 0.020 ± 0.007 Å (at > 95% CL) respectively. We do not detect O vi in the high-SNR stacked spectrum of 1765 quasar−cluster pairs giving a 3 upper limit on the REW of < 0.006 Å (Fig. 2).
• The REW-profile of Ly exhibits a decreasing trend with increasing cl ( cl / 500 ) which is well described by a power-law with a slope of −0.74 ± 0.17 (−0.60 ± 0.15; see Fig. 3). Similarly, the covering fractions of Ly and O vi also gradually decrease with increasing cl and cl / 500 , while no such trend is observed for the C iv covering fraction profile (Fig. 4).
• We compare the Ly REW-profile obtained for the cluster outskirts to the CGM of field galaxies (Fig. 6) and observe an appreciable difference, particularly in the inner regions (≲ 200 ). We attribute this difference to the observed lower covering fraction for Ly in the cluster outskirts.
• We show that the Ly REW-profiles for the 1000 and 4500 sub-samples are consistent with the REW-profile of the full sample within 1 indicating that the signal is not dominated by the CGM of bright galaxies (Fig. 7).
• We do not find any detectable absorption in the stacked spectra when we remove the visually identified absorbers from our analysis (Fig. 8). Only 13% of the visually identified Ly absorbers are related to bright ( ≈ 19) galaxies with |Δ | < 500 kms −1 and gal < 300 kpc. Based on these, we argued that the distribution of the detected cool gas in the cluster outskirts is filamentary, which is also consistent with the observed low covering fractions.
• We find that the CGM of galaxies residing in cluster outskirts is significantly gas poor compared to their field counterparts. This effect is more pronounced for galaxies that are closer to clusters or that are in massive clusters (Fig. 5).
We discuss the roles of environmental processes such as ram pressure stripping, overshooting, and pre-processing to understand the observed gas deficiencies in cluster galaxies. Our findings suggest that at smaller clustocentric distances, ram pressure stripping and overshooting are likely the dominant processes, whereas at large distances (≳ 4 − 5 500 ) pre-processing dominates the stripping of circumgalactic gas. We argue that the cool gas detected in cluster outskirts arises in part from stripped-off CGM and from large-scale filaments feeding the clusters with cool gas. In the future, we plan to perform detailed photoionization modelling of the individual absorbers detected in the outskirts of clusters to constrain the physical conditions and metal enrichment and to determine the relative contributions of the two channels to the observed absorption signals.
the Bok telescope, Steward Observatory, University of Arizona; and the Mayall telescope, Kitt Peak National Observatory, NOIRLab. NOIRLab is operated by the Association of Universities for Research in Astronomy (AURA) under a cooperative agreement with the National Science Foundation. Pipeline processing and analyses of the data were supported by NOIRLab and the Lawrence Berkeley National Laboratory (LBNL). Legacy Surveys also uses data products from the Near-Earth Object Wide-field Infrared Survey Explorer (NEOWISE), a project of the Jet Propulsion Laboratory/California Institute of Technology, funded by the National Aeronautics and Space Administration. Legacy Surveys was supported by: the Director, Office of Science, Office of High Energy Physics of the U.S. Department of Energy; the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility; the U.S. National Science Foundation, Division of Astronomical Sciences; the National Astronomical Observatories of China, the Chinese Academy of Sciences and the Chinese National Natural Science Foundation. LBNL is managed by the Regents of the University of California under contract to the U.S. Department of Energy.

DATA AVAILABILITY
The spectral data underlying this paper are available in the HSLA Archive at https://archive.stsci.edu/ missions-and-data/hsla. The DESI photo-z galaxy catalog is publicly available at http://cdsarc.u-strasbg.fr/viz-bin/ cat/J/ApJS/242/8.