Searching for cold gas traced by MgII quasar absorbers in massive X-ray-selected galaxy clusters

,


Introduction
Galaxy clusters are an important laboratory for testing models of gravitational structure formation, constraining the parameters of cosmological models, measuring the mean matter density of the Universe, and for investigating galaxy evolution as well as plasma physics (Kravtsov & Borgani 2012).Despite the advances of decades of observations and simulations, the interactions between the gas phases that make up the baryon content of clusters remain largely unknown.
The baryon budget in clusters within R 5001 , although not well constrained, is made up of approximately 70% hot intracluster medium (ICM), 13% cold gas from stars and galaxies, and 17% warm hidden baryons that are yet to be observed (Ettori 2003;Fukugita & Peebles 2004;Gonzalez et al. 2007a,b;Kravtsov et al. 2005).The ICM in galaxy clusters has been thoroughly studied in X-rays with instruments such as ROSAT, XMM-Newton, Chandra (see e.g.Sarazin 1986;Rosati et al. 2002;Böhringer & Werner 2010, and references therein), and now eROSITA (Merloni 2012;Predehl et al. 2021;Liu et al. 2022).In the meantime, optical and near-infrared observations from the Dark Energy Survey and the Hyper Supreme-Cam survey provide a different perspective on the growth and evolution of baryons in clusters (Aihara et al. 2018;DESI Collaboration et al. 2016).Galaxy cluster outskirts are found to be multi-phase structures, where new infalling material flows inwards along filaments and cooler clumps become an important component (Reiprich et al. 2013).The gas density distributions in this multiphase structure depend on several physical properties, including temperature and dynamical state.Indeed, there exists a steep relation between gas mass and temperature that implies a decrease in the total gas mass content of cooler clusters relative to highermass systems (Croston et al. 2008).
Alongside these developments, structure formation studies have progressively focused on the cycle of baryons.The physical processes by which gas is accreted onto galaxies, transformed into stars, and then expelled from galaxies into the circumgalactic medium (CGM) are of paramount importance for galaxy formation and evolution (Tumlinson et al. 2017).The vast majority of the CGM studies so far have focused on field galaxies, while almost ∼50% of galaxies in the local Universe are in clusters or groups (Eke et al. 2004).In addition, simulations predict that a major fraction of the baryons at lower redshift are found in gas 10 5 -10 6 K, the so-called warm-hot intergalactic medium (WHIM; Cen & Ostriker 1999;Davé et al. 2001).Focusing on the Virgo cluster as a close-by laboratory, early results from Yoon & Putman (2013) showed that the CGM surveys must consider the role of the environment.Clusters provide rather complex astrophysical and dynamical systems, where many different physical processes take place, such as galactic winds, active galactic nucleus (AGN) feedback, and gas stripping, leading to mixing and redistribution of metals in the ICM (Simionescu et al. 2009;Kirkpatrick & McNamara 2015).Nevertheless, as mentioned above, the baryon budget in clusters is not well constrained, which is mainly because it is unclear whether or not all baryonic constituents have been identified and quantified.
The ICM is an essential component for the assessment of the cosmic baryon and metal budgets, given that the ICM contributes even more to the overall baryon budget than stars (e.g.Gonzalez et al. 2013;Péroux & Howk 2020).Within these large structures, the ICM material is composed and mixed with the CGM material of many group members.Galaxies are infalling and orbiting within halos, while their properties are likely affected by encounters with the warm and hot gas sitting in the gravitational well (e.g.Popesso et al. 2015).The large amounts of hot gas trapped within their deep potential wells make clusters of galaxies shine in the X-ray waveband (e.g.Jones & Forman 1999).However, determining the abundance of elements other than iron probed by X-ray observations is currently challenging and uncertain (Mernier et al. 2018;Frebel 2018).
From a theoretical standpoint, multi-dimensional numerical plasma simulations in the ICM of clusters, groups, and massive galaxies show how cold-gas filaments can condense out of a hot halo.These filaments fall through the hot gas where cool blobs can survive the passage through the hot medium due to Kevin-Helmholtz instability (Sharma et al. 2012;McCourt et al. 2012).Modern cosmological hydrodynamic simulations show that cool-warm gas from the CGM of galaxies is stripped from the galactic potentials through ram-pressure forces of the hot ICM (Ayromlou et al. 2019;Yun et al. 2019).On the other hand, cosmic web filaments interact with galaxy clusters permeating gas streams of relatively metal-poor cool-warm gas.These and other mechanisms create an inhomogeneous and turbulent multi-phase ICM with a range of physical and chemical properties (Bahé et al. 2013;Simionescu et al. 2019;Kravtsov & Borgani 2012;Kunz et al. 2022).Models of interactions between the supermassive black holes and the large-scale atmospheres from these massive structures suggest that the cooling of the CGM is rapid and inhomogeneous (Donahue & Voit 2022).The conditions of the ambient gas in these haloes quickly develop a multiphase structure, where the colder and denser gas sinks into the central galaxy and is accreted by the supermassive black hole (Tremblay et al. 2016).This cold gas is then again pushed outwards through AGN outflows, but can rain back down after the hydrodynamic phenomena that lifted it subsides (McCourt et al. 2012;Sharma et al. 2012).This cycle has been observed and was interpreted as being driven by AGN feedback (Revaz et al. 2008).
The hotter component (T ∼ 10 6−7 K) of this multi-phase gas emits at X-ray wavelengths.As a reference, the hot phase of the Milky Way's CGM has been observed both in emission and absorption, and has been characterised in great detail and found to be made up of three distinct components: the hot, warm, and warm-hot components (Das et al. 2021;Mathur et al. 2021;Locatelli et al. 2023).Beyond the Milky Way observations, the gas is diffuse and challenging to detect in emission due to low sur-face brightness (Cantalupo et al. 2019); it is best probed using quasar absorption line spectroscopy, which has proven to be a powerful probe of these environments (Péroux & Howk 2020).Absorption lines detected against bright background quasars offer the most compelling way to study the distribution, chemical properties, and kinematics of CGM gas (Hamanowicz et al. 2020;Szakacs et al. 2021).In these quasar absorbers, the minimum gas density that can be detected is set by the brightness of the background source and thus the detection efficiency is independent of redshift (Tripp et al. 1998).While individual absorption measurements are limited to a pencil-beam along the line of sight, a sample of sightlines allows us to statistically measure the mean properties of galaxy clusters.These techniques have been extensively used to investigate the gaseous halos of isolated galaxies (Szakacs et al. 2021).However, without a quantitative description of the most massive structures in the low-redshift Universe, a full census of baryons and metals in the Universe cannot be obtained.
Among quasar absorbers, the low-ionisation MgII doublet is known to trace cold 10 4 K gas.Because of its distinct doublet feature, MgII has been used extensively in a large number of spectroscopic surveys.In the last two decades, MgII absorption system surveys have been used to study the physical properties of large samples of galaxies over a wide range of luminosities and morphologies (Lanzetta & Bowen 1990;Nestor et al. 2005;Narayanan 2007;Lopez et al. 2008;Seyffert et al. 2013;Anand et al. 2021).Absorption by MgII in cosmological galaxy formation simulations has been analysed (Nelson et al. 2020;Augustin et al. 2021), and is shown to closely relate to star forming regions, galactic outflow, and galactic discs (Bowen & Chelouche 2011).These are known to be associated with galaxies and their CGM (Bouché et al. 2006;Zhu et al. 2015), and individual absorption lines allow us to characterise the spatial distribution and physical properties of the cold-gas clouds located in the vicinity of galaxies (Lan & Mo 2018;Zhu et al. 2015).
Recently, large spectroscopic surveys have focused on the incidence of MgII in and around galaxy clusters, with different methods and samples, showing a wide range of detections (Anand et al. 2022;Mishra & Muzahid 2022).In particular, Mishra et al. (2023) probed cold neutral gas in the outskirts of low-redshift galaxy clusters.The present study extends previous works to higher cluster masses.Indeed, we build on these works to analyse the cold gas traced by MgII absorbers in a large sample of massive X-ray-selected clusters.Aiming at this higher mass range provides fresh clues as to the evolution of the physical properties of galaxy clusters, including gas mass fractions as a function of galaxy cluster masses.
The present paper is organised as follows: Section 2 presents observational data used in this study.Section 3 details the analysis performed, while Section 4 focuses on the results.In section 5, we discuss our findings in the context of other works and state-of-the-art simulations.We summarise our findings and present our conclusions in Section 6.Here, we adopt an H 0 = 67.74km s −1 Mpc −1 , Ω M = 0.3089, and Ω Λ = 0.6911 cosmology.The equivalent width (EW) always refers to the MgII 2796Å line unless stated otherwise.

The foreground cluster sample
We start our analysis from the largest spectroscopically confirmed sample of massive clusters compiled from the SPectroscopic IDentification of ERosita Sources (SPIDERS) program of the Sloan Digital Sky Survey (SDSS; Clerc et al. 2016Clerc et al. , 2020;;Kirkpatrick et al. 2021;Ider Chitham et al. 2020).We drew the cluster candidates for SPIDERS from a subset of CODEX (Finoguenov et al. 2020), an X-ray-selected catalogue of clusters from the Rosat All Sky Survey (RASS), specifically in the 5 350 square degree BOSS imaging footprint; 2 740 CODEX clusters are included in the SPIDERS sample (Kirkpatrick et al. 2021).These are complemented with the X-CLASS cluster sample (based on serendipitous XMM-Newton observations) with 124 uniquely identified and confirmed clusters (Clerc et al. 2014).Visual inspections of individual spectra are carried out by trained astronomers to verify the existence of a cluster.After a first run with the redMaPPer (Rykoff et al. 2016) algorithm for automatic membership assignment, a minimum of three visual inspecting members are required to converge on a final redshift determination, the final condition being that these members lie close to each other in a velocity-distance diagram.Due to the uncertainty associated with the position of the detections, several measures are taken to ensure the closest possible estimation of redshift.The likelihood that a candidate will be validated is highly dependent on redshift, cluster richness, and the number of spectra available.For example, in the spectroscopic redshift determination, there are 622 instances where the mean spread in the measurements is δ z = 0.00049 (147 km s −1 ), while the maximum is δ z = 0.0055 (1 650 km s −1 ).We conservatively assume that all the objects with velocity offsets from the cluster centre greater than 5000 km s −1 are rejected as members of the cluster.
Figure 1 shows the distribution of foreground cluster masses.Scaling from their X-ray luminosities, the average mass of these systems amounts to M 500c = 2.7×10 14 M ⊙ (converted from M 200c = 3.9×10 14 M ⊙ ).Here, we mainly express the cluster mass in units of M 200c , that is, the mass within a radius R 200c .This corresponds to the radius enclosing an average density that is 200 times the critical density of the Universe at the cluster's redshift and ranges within 0.048 < R 200c < 0.160 deg or within 1Mpc< R 200c < 3Mpc depending on redshift.This radius is expressed as: In order to compare these results to the latest related published research, we use the public open-source package COLOS- SUS to convert between M 200c and M 500c (Diemer 2018), and continue to do so throughout the paper for consistency.
The observed galaxy cluster distribution sample covers a redshift span of 0.03 < z < 0.677 (a subset of which -systems with z>0.3-is shown in Figure 2), with a typical statistical uncertainty on each redshift of ∆ z /(1 + z) = 6 × 10 −4 (Kirkpatrick et al. 2021).The number of spectroscopic members per system ranges between 3 and 75, with a mean of 12 members.

The background quasar sample
We selected spectra of background quasars from SDSS.The spectra were obtained with the 2.5m Sloan Telescope as part of the BOSS and eBOSS surveys (Smee et al. 2013), covering a wavelength range of 3600 Å to 10 400 Å at a spectral resolution of λ/∆λ ∼2000, or with the SDSS-I instrument, covering 3800 Å to 9100 Å, with the same spectral resolution.The Data Release 16 of SDSS comprises a complete selection of spectroscopically confirmed quasars (Lyke et al. 2020;Ahumada et al. 2020), and includes 750 414 confirmed quasars with 0.8 < z < 2.2 (Table 1).

Pairing background quasars with foreground clusters
A key component of the analysis is to determine which quasar sightlines pass close to a foreground galaxy cluster in projection on the plane of the sky.We determined the number of quasarcluster pairs by cross-matching the foreground cluster sample with the background quasar catalogue.We made use of the cluster optical centre from the DR16 catalogue.We selected all SDSS spectra from the parent sample with a distance in right ascension and declination within three times R 200 (apparent R 200c , ∼ 3 Mpc) of the centre of a cluster, and a redshift of z quasar > z cluster + 0.01 for a total of 434,736 spectra.Among these clusterquasar pairs, we select the ones for which the foreground cluster will have a redshift such that MgII is covered by the SDSS spectra (z>0.3), and a quasar redshift of z < 5.This leads to a total of 18 694 cluster-quasar pairs.Table 1 presents the numbers of objects resulting from the various selection cuts: (i) the quasar redshift must be higher than the redshift of the intervening cluster; (ii) the corresponding MgII absorption line should fall within the SDSS wavelength coverage and the quasar redshift must be z < 5; and (iii) the final spectrum must be free of wavelength gaps 0. 5 1.0 1.5 2.0 2.5 3.0 3.5 4 and other errors.Indeed, some of the SDSS spectra have gaps in their wavelength coverage due to bad detector columns or high noise peaks related to cosmic rays for example.These objects were removed from the sample.Figure 3 shows the distribution of the redshift of the High-quality SDSS spectra, the redshift of the known MgII absorbers in these spectra, and the redshift of the foreground clusters.We note that the cluster centre is not always well constrained because of the uncertainty in the redshift of some of the members.
The signal-to-noise ratio (S/N) is computed following Mas-Ribas et al. (2017) for each spectrum as S/N= mean(flux)/spread, where the spread is the standard deviation defined as the dispersion of the flux relative to the mean, and the mean calculated over the flux of the spectra.The distribution of the S/N of the spectra is presented in Figure 4.
For each of the pairs, we calculate the angular separation between the foreground cluster and the background quasar projected on the sky plane.  1 cuts: (i) the quasar redshift must be higher than the redshift of the intervening cluster; (ii) the corresponding MgII absorption line should fall within the SDSS wavelength coverage and the quasar redshift must be z < 5; and (iii) the final spectrum must be free of wavelength gaps and other errors.
the R 200 radius of the cluster.We stress that while we select in units of R 200 , all the figures presented in the present paper are in units of R 500 in order to ease comparison with results from the literature.

Sample of known MgII quasar absorbers
Our analysis also makes use of the MgII quasar absorption catalogue of Anand et al. (2021).The authors developed an automated pipeline to detect intervening metal absorption line systems with a matched kernel convolution technique and adaptive S/N criteria.By processing one million quasars from the SDSS DR16, these authors compiled a sample of about 160 000 MgII absorbers in the redshift range 0.3 < z abs < 2.3.After crossmatching our sample of cluster-quasar pairs with this absorption catalogue, we get a match for 4 150 spectra.
Figure 3 shows the redshift distribution of these MgII quasar absorber matches in orange while the quasar spectra with foreground clusters are shown in blue, and the background quasars in green.We use the publicly available catalogue to estimate the velocity difference, ∆v, between the absorption redshift (z abs ) and the cluster redshift (z cluster ).We express ∆v as follows: where c is the speed of light in km/s and z cluster is the redshift of the clusters for the individual cluster-quasar pairs.These velocity differences can be very large, because the absorption redshifts are not necessarily close to the clusters.Figure 5 displays the full range of ∆v within [-5000:5000] km s −1 .Due to the high velocity dispersion expected in galaxy clusters (Girardi et al. 1993), we chose to focus on the range of [-2000:2000] km s −1 around the clusters.Figure 5 indicates that a lower cut, namely [-1000:1000] km s −1 , would make little difference to the sample.In total, we found 32 absorbers associated with the clusters in our sample, resulting in an overall incidence of 32/16 224 ≈ 0.2% where 16 224 refers to the total number of quasar spectra studied (see Table 1).

Normalising the background quasar spectra
In order to retrieve the quasar absorption systems, we removed the intrinsic spectral signature of the background quasars.To this end, we modelled each of the SDSS quasar spectra with the highly flexible Python QSO fitting code, PyQSOfit (Guo et al. 2018).This algorithm fits the spectral features, broad emission lines and continuum slope of the quasars.Figure 6  the ratio between the original data and the fitted continuum, resulting in a normalised spectrum.For some objects, the emission lines cannot be perfectly modelled and after normalisation, some residuals remain.We note in particular the noise of the SDSS spectra increases at each end of the wavelength window, resulting in noisy edges in the normalised spectra.Nonetheless, the method proved very efficient to remove the features of quasar emission.

Stacking background quasar spectra at the position of foreground clusters
In order to detect the weakest tracers of the cold gas, we take advantage of the large number of spectra available, and stack the background quasar spectra at the redshift of the foreground clusters.We first shift the quasar spectra to the cluster rest frame.We use the fact that all the SDSS spectra have a constant log step of 0.0001 (approximately 69 km s −1 ) to avoid interpolating between pixels.We note here that for this specific cluster sample Clerc et al. (2020); Kirkpatrick et al. (2021) indicate that the typical statistical uncertainty on cluster redshifts is ∆z/(1 + z) = 6 × 10 −4 , which is about 180 km/s.This mean value does not reflect the distribution of the cluster redshift uncertainties, which may depart from a Gaussian profile.Following the method of Fresco et al. (2020), we performed a median stack in the velocity region that can be attributed to the foreground cluster within a window of [-2000:2000] km s −1 .After this process, a Gaussian smoothing of 1 pixel was applied using gaussian smooth from specutils.
In this analysis, we stacked two independent sets of quasar spectra.First, we stacked quasar spectra with known MgII absorbers that happen to be close in velocity space to a foreground cluster.We refer to this sample as the MgII-selected sample; this sample is used as a test for our methodology and to check whether or not the stacked spectra recover the features even with the large velocity dispersions inherent to massive galaxy clusters, which will likely induce velocity offsets at the impact parameter where the quasar's line-of-sight pierces the cold gas traced by MgII.Table 2 summarises the (small) number of avail-able spectra in the MgII-selected sample and the corresponding S/N of the resulting stacked spectrum.
Second, we 'blindly' (i.e.without prior knowledge of the presence of MgII absorbers) stack all the quasar spectra that have a foreground cluster; these correspond to the so-called blind sample and are described in Table 2.

Measuring equivalent widths and column densities
We first search blindly -that is without prior knowledge of the presence of any MgII in the quasar spectrum-for the minimum flux within the full velocity range [-2000:2000] km/s around the MgII 2796 Å line.The wavelength of the minimum flux is set to the MgII 2796 Å line.We then use the precise 750 km s −1 velocity separation between the MgII lines of the doublet to materialise the position of MgII 2803 Å.To quantify the amount of cold gas, we measure the EW of the MgII lines.By focusing on the MgII 2796 Å line, we obtain a conservative measurement of the strength of the MgII doublet.The observed spectra can be seen in Figures 7 and 8 for the MgII-sample and blind sample, respectively.The EW 2796 is measured by integrating the area of the spectrum over the wavelength range (i.e.summing the absorbed flux values per pixel) defined by vertical lines around our first line of the doublet, which we fixed at the MgII 2796 Å line -500km/s blueward and +375km/s redward.The continuum level is defined as the observed spectrum divided by the fitted continuum as previously described in section 3.1.We then measured the second line of the MgII doublet, the MgII 2803 Å line, which is located at +750km/s redward.The same method is applied to the MgII 2803 Å line, which is integrated from -375km/s (blueward) to +500km/s (redward) of the rest-frame wavelength.These measures result in estimates of the strength of the detected MgII absorbers.We made use of the Linetools software package (Prochaska et al. 2017), which calculates the EW 2796 in observed wavelength space.We then converted the observed-frame EW 2796 to rest-frame EW, following the usual relation: where the redshift ⟨z cluster ⟩ is the average cluster redshift of the MgII selected sample as listed in Table 2.The estimates of the uncertainties are calculated from the mock spectra, as described in section 4.4.
We then computed the MgII column density according to the linear relation between the EW and the column density using the following equation from Zhu et al. (2014): where λ rest is the rest wavelength of each of the MgII lines of the doublet, λ rest = 2796 Å and λ rest = 2803 Å.The corresponding oscillator strengths are f osc = 0.615 and f osc = 0.306.For saturated absorbers with EW 2796 > 0.15 Å, we derive a lower limit on the MgII column density (see Table 2).

The MgII-selected sample
We used the method described in the previous section to compute the EW of the detected MgII feature in the MgII-selected stack., is located at -100km/s from the cluster redshift, but we note that this slight offset is eight times smaller than the typical velocity dispersion expected in galaxy clusters (∼800 km/s, Kirkpatrick et al. 2021).This means that the expected shift from the velocity dispersion alone likely explains the small offset observed here.
To remove the effects inherent to the a priori knowledge of MgII, we also compute the frequency-weighted EWs and column densities as follows: EW = EW rest f c , where the covering fraction, f c , is defined as the number of spectra with known ab-sorbers over the total number of spectra available for each bin of mass and angular separation.These values are listed in Table 2.

The blind sample
Figure 8 displays the stack of all the spectra, with a total of 16 224 cluster-quasar pairs.This stack shows an absorption of the MgII doublet at the redshift of the cluster centre, with a rest EW for MgII λ2796Å of EW 2796 =0.056 ± 0.015Å (3.7σ significance), and a column density of log [N(MgII)/cm −2 ]=12.12±0.1.The mean mass of clusters contributing to the stacked spectra is 9.7 × 10 14 M⊙ (M 500 ), and a mean angular separation of the cluster-quasar pairs is about 2 × R 200 .We note that the absorption is redshifted from the cluster centre by 350 km s −1 from the absorption line centre.This shift might be related to the presence of inflows and/or outflows, although it is challenging to determine this from absorption line studies alone because the orientation of the clusters and sightlines might complicate the interpretation.The uncertainty in the measurements was calculated by randomising the redshift of the clusters and repeating the stacking 500 times, as described in section 4.4.
The detection is broad, ranging over 1 000 km/s in width for each line.This is likely the result of the offset in velocities in the different clusters along the background sightlines.tion, we computed the stack of the blind sample after removing the known MgII absorbers.We find that the absorption features remain largely unchanged, which further supports the idea that the detected absorption feature is not solely the result of a few known strong absorption systems.

Mock Sloan spectra with simulated MgII absorbers
In order to quantify the systematic uncertainty of the method and derive an error on the EW measurements, we use 22 000 S/N=200 mock Sloan stacks representative of the S/N achieved with our observed stack.The synthetic spectra are convolved with the line spread function (LSF) of the the SDSS instrument to achieve the same spectral resolution as the observations.Similarly, noise properties typical of Sloan spectra are included.We used Gaussian noise with a spread similar to Sloan's typical error, including increased noise towards each end of the spectrum.We additionally inserted an MgII absorber doublet from the TNG50 cosmological magnetohydrodynamical simulations.
To compute MgII, we take the total magnesium mass per cell as tracked during the simulation, and use CLOUDY (Ferland et al. 2017) to calculate the ionisation state assuming both collisional and photo-ionisation following the modelling approach of Nelson et al. (2020).We note that the CLOUDY modelling indicates that at densities of < 10 −2.5 cm −3 , the photo-ionisation dominates, but overall both photo-ionisation and collisional ionisation processes are expected to play a role at densities typical of the ICM (Peterson et al. 2001).We then ray-trace through the simulated gas distribution to create synthetic absorption spectra akin to those in real observations (Szakacs et al. 2023, Nelson et al. in prep.).This is similar in its rationale to several other techniques for creating absorption spectra from hydrodynamical simulations, such as specwizard (Theuns et al. 1998;Schaye et al. 2003), Trident (Hummels et al. 2017), and pygad (Gad 2021).The absorbers are inserted at a given wavelength position but with different EWs, ranging from 0.05 to 5 Å.We used these mock spectra to run the same search and measurement algorithm as that used for analysing the observations.Our results show that in 22 000 mock MgII stacks, only 3% of the detected absorption as found by the minimum-flux pixel approach is not associated with a simulated input MgII line.This means that the code finds the minimum flux outside of the [-2000:2000]km/s area where the MgII doublet should be located at a given redshift.
We then measured the EW 2796 of the simulated spectra with the method used to analyse the observations.We find that the EW measurements from our analysis are consistent with the input EW from the simulations.The analysis also indicates that the mean difference between the measured values and the input EW is 0.013 Å for EWs typical of the one we measure in the blind stack.This is similar to our estimated error in the blind stack of 0.015Å.

Uncertainty assessment from bootstrapping
In addition to using the mock spectra with simulated absorbers to quantify the possible systematic errors of our method, we also performed a non-parametric bootstrapping experiment by stacking the full sample of cluster-quasar pairs -at cluster redshifts that have been randomised-500 times.We then repeated the measurement process mentioned in section 3, where we first find the minimum in the flux, and then measure the EW of the possible absorption feature.We do a binning of the 500 stacks for each flux point, and calculate the median of each bin, obtaining one final median array shown as the black dotted line in Figure 9, where the grey lines in the background represent each random stack, and the overlayed black dotted line represents the medium flux per velocity bin of the 500 random stacks.We further overplot in red the blind stack from Figure 8.The results show that the stacks at random redshifts do not display absorption features at the level of the blind stack shown in red.This demonstrates that the detection is significant.Indeed, the measurement of the EW of the blind stack is higher than the distribution of the measured EW of the absorption features of the random stacks as seen in Figure 10.As expected, the distribution peaks around zero, because we do not expect any absorption feature at the random redshift we are considering.Assuming a normal distribution, we get a standard deviation from the EW distribution of 0.015Å.Taking this as our measurement uncertainty gives us a 3.7σ significance for the MgII 2796Å absorption feature.

Evidence of cold gas in X-ray-selected clusters
The detection of strong MgII absorption in the MgII-selected sample validates our approach of stacking many medium-S/N background quasar spectra to increase the sensitivity to cold gas in clusters.This detection further demonstrates that, despite the high velocity dispersions within the clusters, the 10 4 K gas can be traced by the MgII absorbers observed in the spectrum of background quasars.Our findings based on the blind stack of the full sample also indicate the presence of some cold 10 4 K gas traced by MgII in the intracluster environment.We note that the CGM has a rather loose definition, and so whether these MgII absorbers are associated with the circumgalactic gas of individual galaxy members or with the intracluster gas is somewhat subjective.More important in this work is the total amount of 10 4 K cold gas that is being detected in the dense cluster regions.These observational results are in line with expectation from various numerical (Sharma et al. 2012;McCourt et al. 2012) and hydrodynamical simulations (Pillepich et al. 2018;Nelson et al. 2018;Tremmel et al. 2019;Butsky et al. 2019), which show that the interactions of supermassive black holes and large-scale atmospheres of massive clusters cause the CGM to cool down rapidly and homogeneously (Donahue & Voit 2022).The low entropy in groups of galaxies and galaxy clusters explains how the cold gas filaments can condense inside a hot halo, where local thermal instability makes it possible for cool gas blobs to survive the hot medium (Sharma et al. 2012;McCourt et al. 2012).
To put our results into perspective, we plot in Figure 11 the EW (left y-axis) and column density (right y-axis) of MgIIabsorbing gas as a function of the foreground cluster mass (right panel) and the projected distance expressed in kiloparsecs (kpc) (left panel) in comparison with similar studies, which we discuss further below.Lopez et al. (2008) were the first to look for cold gas in foreground galaxy clusters using background quasar spectra.They based their study on a sample of 442 cluster-quasar pairs.To this end, the authors made use of the third data release of the SDSS with high-redshift cluster or group candidates from the Red-Sequence Cluster Survey (Gladders & Yee 2005).Lopez et al. (2008) found that there is proportionally less cold gas in more massive clusters than in low-mass systems when using models of galaxy counts.This refers to the relation between overdensities of MgII absorbers in clusters with much denser galaxy environments.When considering the stellar baryon fraction, studies still disagree on the exact slope of the stellar and total baryon fraction as a function of cluster halo mass (Gonzalez et al. 2013).Nevertheless, the Lopez et al. (2008) results are in agreement with our findings: the most massive structures of their sample (M > 10 14 M ⊙ ) have an overdensity of absorbers with EW 2796 > 1Å that is twice higher than for moderate mass clusters (M ∼ 2 × 10 13 M ⊙ ).The effect is even more pronounced at smaller distances (d< 1Mpc) from the cluster centre.At the same time, Lopez et al. (2008) point out that the more massive clusters contain five times more galaxies than the less massive ones.While searching for this cold gas through MgII absorption lines, these latter authors find that a subsample of their massive clusters yield a stronger and more significant signal.Recently, Lee et al. (2021) performed a cluster-quasar cross-correlation with SDSS DR14 quasars and redMaPPer clusters, with a total of 82 000 cluster-quasar pairs.Although there was no stacking of quasar spectra in this latter work, the authors report that the MgII absorber detection rate per quasar is 2.70±0.66times higher inside the clusters than outside them.This shows that Mg II absorbers are abundant in clusters compared to the field.

Cold gas in galaxy clusters
Using a photometric sample of galaxy clusters identified in the Legacy Imaging Survey of Dark Energy Survey Instrument (DESI) by Zou et al. (2021), a follow-up study led by Anand et al. (2022) cross-correlates the MgII absorption catalogue (Anand et al. 2021) with cluster and luminous red galaxies (from SDSS BOSS survey) with halo masses in the range between 10 13.8 and 10 14.8 M ⊙ .Furthermore, Mishra & Muzahid (2022) cross-match the SDSS cluster catalogue of Wen et al. (2012) and the SDSS DR16 quasar catalogue of Lyke et al. (2020).The median impact parameter of the clusters from the quasar sightlines is 2.4 Mpc (median 3.6 R 500 ).Mishra & Muzahid (2022) were the first authors to stack the background quasars to increase the sensitivity of the experiment; in their work, they measured the total EW and assumed that the observed MgII line falls on the linear part of the curve of growth, meaning that the EW 2796 is two-thirds of the total EW.2022) and Mishra & Muzahid (2022) at different values of projected distance.The right panel of Figure 11 shows the EWs as a function of cluster mass.Clearly, the works of Anand et al. (2022) and Mishra & Muzahid (2022) probe a lower mass range than our study.We compare the median of the two closest numbers in projected distance from Anand et al. (2022), and the closest point from Mishra & Muzahid (2022) with our tentative detection (red circle).Both Mishra & Muzahid (2022) and Anand et al. (2022) also report a clear trend of decreasing EW 2796 with increased projected distance.Despite the differences between the approaches, the reported EW is in line with results from Anand et al. (2022).The cluster masses in our study range from M 500 = 1 × 10 14 to 8.4 × 10 15 , with an EW 2796 of 0.056±0.015Å of MgII 2796 line.Using three times the number of quasar cluster pairs, the rest EW 2796 from Mishra & Muzahid (2022) remain one order of magnitude smaller than our measured rest EWs from both our MgII-selected sample and our blind sample, as shown in the lower panel of Figure 11.
Our higher EW measurements and limits could be associated with the higher mass in our sample, as shown in Figure 11.Indeed, Sharma et al. (2012);McCourt et al. (2012); Donahue & Voit (2022) propose that although cooling is negligible near the virial radius, it becomes more important at higher densities and smaller radii, which is in line with these observations.Additionally, complementary observational studies have analysed the incidence of MgII absorbers in the context of isolated and group environments.Specifically, MgII metal line emission studies from galaxies at z∼1 find extended MgII emission in a large blind galaxy survey in the Muse Analysis of Gas around Galaxies or MAGG (Dutta et al. 2020).Their results show that the MgII is associated with multiple galaxies, and their measurements of the MgII emission flux in these groups of galaxies are on average five times stronger than in isolated galaxies.Their findings therefore favour the scenario where the hydrodynamic interactions among group members are the primary reason for the increased strength of MgII lines.Similarly, observational studies of the cold gas in the CGM at higher redshifts (z∼2) by Nielsen et al. (2020) reveal a contrast, in that these galaxies are actively forming at cosmic noon, where the CGM is less massive than at lower redshifts.The measurements of MgII EW are found to be larger (EW 2796 ≥ 0.5Å) in the lower-redshift (z∼1) studies (i.e. at higher mass) from the MusE GAs FLOw and Wind survey or MEGAFLOW (Schroetter et al. 2016(Schroetter et al. , 2019) ) than at z∼2. Figure 12 contrasts our results for EW 2796 with observations of lower-mass structures, including photometric redshift clusters from the legacy imaging survey of the Dark Energy Survey (Anand et al. 2022), luminous red galaxies (LRGs) from SDSS DR16 (Anand et al. 2022), and SDSS DR14 (Lan & Mo 2018) as a function of projected distance.The LRGs with typical masses of ∼ 10 11.2 M ⊙ from Anand et al. (2022) show a smaller average MgII EW than the clusters.Although these values seem to converge at large impact parameters, they report this difference as being due to the larger halo masses and denser environments in the cluster sample.On the other hand, the LRGs from Lan & Mo (2018) show larger EW 2796 at lower impact parameter, and a steeper decline at higher projected distance compared to clusters from the legacy imaging survey of the Dark Energy Survey and the LRGs from Anand et al. (2022).We may again attribute the difference in EW 2796 we report in this study to the mass differences, because our sample of clusters has a higher mean mass than the samples used in these previous works (see lower panel of Figure 11).
While there still exists disagreement in the trend of total baryon fraction as a function of cluster mass (Gonzalez et al. 2013;Laganá et al. 2013), the difference in mass between all the previously mentioned studies and our work, in connection with the higher measurements from EW from our higher-cluster-mass sample, reveals a trend: increased baryon fraction at higher cluster masses.These differences in gas measurements for halos of different mass can also be seen when comparing galaxies versus galaxy clusters as shown in Figure 12.Similarly, Anand et al. (2022) perform direct comparison of LRGs with legacy imaging survey clusters, finding consistently higher MgII EWs in the cluster sample than in the LRGs, where the clusters are between two and three times more massive than LRGs within R 500 .

Comparison with simulations
Butsky et al. (2019) used the RomulusC simulations (Tremmel et al. 2019) to probe the nature of the multi-phase cool-warm (10 4 < T < 10 6 K) gas in and around a galaxy cluster of mass 10 14 M ⊙ .Their study makes predictions for the covering fractions of key absorption-line tracers, both in the ICM and CGM of cluster galaxies using synthetic spectra.The authors find there is a significant quantity of multi-phase gas in the cool (10 4−5 K) gas at all clustocentric radii.The results from Butsky et al. (2019) indicate that the column density of all ions declines from the cluster centre out to 1 Mpc, but remains relatively flat towards the edge of the halo.In a more recent set of CGM simulations around massive galaxies and groups of galaxies dubbed Romulus, Saeedzadeh et al. (2023) expand their previous simulations, showing that the presence of cold gas in the CGM can be described as filaments of inflowing cooling gas as well as gaseous tails from possible satellites.Furthermore, condensation patches can also originate from density perturbations cooling rapidly, which provides further theoretical support to the observations of cold gas in clusters.
In Figure 13, we compare our results with the TNG50 simulation (Pillepich et al. 2019;Nelson et al. 2019) from the TNG suite (Marinacci et al. 2018;Springel et al. 2018;Naiman et al. 2018;Nelson et al. 2018;Pillepich et al. 2018).The aim of TNG is to study the physical processes that drive galaxy formation and to investigate how galaxies evolve within large-scale structures.TNG50 is a gravomagnetohydrodynamics (MHD) cosmological simulation including a comprehensive model for galaxy formation physics (Weinberger et al. 2018;Pillepich et al. 2018) at the highest resolution.TNG50 includes 2×2160 3 resolution elements in a ∼ 50 Mpc (comoving) box.Nelson et al. (2020) make predictions of the physical properties of the cold gas traced by MgII in the CGM of galaxies and groups at z = 0.5 based on these TNG-50 simulations.By identifying discrete structures of cool, MgII-rich gas, the authors found that cold gas in these massive halos is made up of thousands to tens of thousands of small (kpc), discrete clouds.The results also indicate a tendency whereby the most massive halos have the highest covering fraction, while the highest column densities are located near the halo centre, where the gas densities are greater.When specifically addressing the physics of small scales in the CGM, simulations are facing the complication of resolving these lower-density structures.This convergence issue arises from limitations in computational resources that limit the maximum resolution, and the priorities on the coding side that are usually set on denser structures.To counter this issue, van de Voort et al. ( 2019) used standard mass refinement and additional uniform spatial refinement to instead increase the resolution in the CGM of a Milky-Way-mass galaxy.Their findings demonstrate a drastic change in the radial profile of neutral hydrogen column density compared with previous simulations.The authors found both the HI covering fraction and column density to increase.Recently, major enhanced CGM resolution cosmological simulations with the zoom-in approach using the TNG galaxy formation model have become available (GIBLE, Ramesh & Nelson 2023).The results indicate that by improving the mass resolution, the cold gas in CGM regions is better resolved.The increase in resolution in the work of Ramesh & Nelson (2023) leads to a larger number of small clouds.While it is worth mentioning that the galactic scales simulated in these works are not directly comparable to our study, the convergence issue is relevant.We stress that there are large mass differences between the models and the data: the TNG50 highest mass bin (∼ 10 13.5 M ⊙ ) displayed in the figure is smaller than our smallest mass cluster (3.28 × 10 14 M ⊙ ).Considering this mass difference, we note that there is likely an excess of MgII-absorbing cold gas in the TNG50 simulations.
Figure 13 shows our results in terms of the column density of MgII as a function of projected distance from the cluster.The conversion of angular size from R 500 in degrees to linear size in kpc at a given redshift was performed using cosmological parameters from Planck Collaboration et al. (2020).We note that the TNG50 simulations probe significantly smaller impact parameters than our observations.Interestingly, the EW in TNG50 changes by more than three orders of magnitude from small (10 kpc) to large (1 Mpc) scales.However, it should be noted that Nelson et al. (2020) provide predictions in terms of column density derived from integrating the gas density of a halo along the wavelength axis to calculate its surface density.The conversion to EWs is done in the present study using an empirically based conversion described in equation 4, which might be less reliable at these extreme column-density values.Equally, we stress that there are large mass differences between the models and the data: the TNG50 highest mass bin (∼ 10 13.5 M ⊙ ) displayed in the figure is smaller than our smallest mass cluster (3.28 × 10 14 M ⊙ ).Considering this mass difference, and the observed increase in the cold column density as a function of cluster mass, we note here a possible disagreement and a likely excess of MgII-absorbing cold gas in the TNG50 simulations compared to observations.To address these findings, a more appropriate comparison would make use of the new TNG-Cluster simulations introduced in Nelson et al. (2023).Indeed, these simulations significantly increase the statistical sampling of the most massive and rarest objects in the Universe, specifically galaxy clusters evolved in cosmic time to z=0 masses of log(M 200c / M ⊙ ) > 14.3 -15.4.These objects would be most relevant to the observational results presented here.Nevertheless, the analysis of the MgII gas properties in these newly released simulations is beyond the scope of the present paper.

Looking forward
The analysis presented here stems from the ambitious endeavours of the SDSS.This series of multi-object surveys, based on extremely large surveys of thousands of quasar absorbers, brought quasar studies to a new era (e.g.Noterdaeme et al. 2012;Bird et al. 2017;Parks et al. 2018).Such surveys advanced the field significantly because they produced homogeneous data products for well over one million low-resolution quasar spectra.In the near future, dedicated spectroscopic surveys on 4m class telescopes will provide a wealth of new lowand medium-resolution quasar spectra in extremely large numbers, notably the DESI experiment (DESI Collaboration et al. 2016), the WEAVE-QSO survey (Kraljic et al. 2022), and a surveys with the 4MOST experiment (including Merloni et al. 2019;Peroux et al. 2023).In particular, DESI is projected to obtain 3 million quasar spectra and 800 000 MgII absorbers.In collaboration with new generations of X-ray missions, including the full eROSITA survey (Merloni 2012), X-ray imaging and Spectroscopy mission (XRISM; XRISM Science Team 2022), and ESA/Athena (Nandra et al. 2013), the study presented here can be expanded greatly by not only increasing the number of background quasar spectra, but also increasing the spectral resolution of these data.We note that the latter is key as the EW limit of detectable absorbers scales linearly with spectral resolution, and a 1 dex increase in EW provides ten times more absorbers.

Conclusions
In this work, we used background quasar spectra from the SDSS/SPIDERS survey to explore the cold-gas content of foreground massive X-ray-selected galaxy clusters.To this end, we stacked approximately 16 000 quasar spectra in the quasar sightline.Our main results are can be summarised as follows: -From the sample of known MgII absorbers within [-2000:2000] km s −1 of the clusters, we detect strong absorption rest EW for MgII λ2796 of EW 2796 =0.35±0.015Å.The uncertainty is calculated through bootstrapping.
As the absorption lines are saturated (at EW 2796 >0.15 Å), we derive a lower limit for the column density of log [N(MgII)/cm −2 ]≥12.92.These results validate the technique and demonstrate that despite the high velocity dispersions within the clusters, the cold 10 4 K gas can be traced by MgII absorbers observed in the spectra of background quasars.
-While we are probing a different cluster mass range, our results are in line with observational findings previously published in the literature (Lee et al. 2021;Anand et al. 2022;Mishra & Muzahid 2022).We also report an excess of MgII gas in the predictions from TNG50 hydrodynamical cosmological simulations (Nelson et al. 2020) compared to our work and others with lower-mass objects in their sample (Anand et al. 2022;Lan & Mo 2018).
Using a similar approach of stacking in upcoming surveys with increased spectral resolution and number of background spectra (DESI, DESI Collaboration et al. 2016), WEAVE (Kraljic et al. 2022) and 4MOST (Merloni et al. 2019), combined with a new generation of X-ray facilities (most notably eROSITA (Merloni 2012) on board SRG), will likely provide additional insight into the cold-gas content of galaxy clusters.In terms of simulations, TNG50 profiles could be extended out to an angular separation of 4000 kpc to overlap with the observations.Further, using the new TNG-Cluster simulations (Nelson et al. 2023) would provide a sample of halo masses that overlap with our observations, which would allow us to assess the importance of the reported halo mass mismatch.

Fig. 1 .
Fig. 1.Distribution of the number of background quasar spectra of the full sample used in the analysis as a function of the foreground cluster mass within R 500 .The dotted line represents the mean mass of the clusters.

Fig. 2 .
Fig.2.Distribution of cluster redshifts from the parent SPIDERS sample(Clerc et al. 2020) used in this study.The clusters constitute a subsample of 1066 objects from the parent sample of 2740 for which MgII absorption lines would fall within the observed SDSS spectral range.The dotted line represents the mean redshift of the clusters at z=0.41.

Fig. 3 .Fig. 4 .
Fig. 3. Histogram of the number of quasars as a function of their redshift.Here, the green colour represents the quasar redshift (16,224 quasars), the orange colour represents the 2,881 MgII known absorbers redshift from Anand et al. (2021) found in these quasars, and the blue as the redshift of the foreground clusters (1,066 clusters).

Fig. 5 .
Fig.5.Velocity difference, ∆v, between the known MgII absorber fromAnand et al. (2021) within 3 × R 200 from the cluster centre and the cluster redshift, as defined in equation 2. The figure displays the full range of ∆v within [-5000:5000] km s −1 .In the subsequent analysis, we chose to focus on the range of[-2000:2000] km s −1 around the cluster as indicated by the dashed lines.

Fig. 6 .
Fig. 6.Example observed-frame SDSS quasar spectrum and data/continuum model ratio spectra after removing the features of quasar emission and converting to unit average flux density with the package PyQSOfit (Guo et al. 2018).On top, we show the original SDSS spectrum, with the multi-component fit performed with PyQSOfit in colour.The bottom panel displays the resulting quasar spectrum normalised to an arbitrary value of 1 by dividing the observed spectrum by the continuum.

Fig. 7 .Fig. 8 .
Figure8displays the stack of all the spectra, with a total of 16 224 cluster-quasar pairs.This stack shows an absorption of the MgII doublet at the redshift of the cluster centre, with a rest EW for MgII λ2796Å of EW 2796 =0.056 ± 0.015Å (3.7σ significance), and a column density of log [N(MgII)/cm −2 ]=12.12±0.1.The mean mass of clusters contributing to the stacked spectra is 9.7 × 10 14 M⊙ (M 500 ), and a mean angular separation of the cluster-quasar pairs is about 2 × R 200 .We note that the absorption is redshifted from the cluster centre by 350 km s −1 from the absorption line centre.This shift might be related to the presence of inflows and/or outflows, although it is challenging to determine this from absorption line studies alone because the orientation of the clusters and sightlines might complicate the interpretation.The uncertainty in the measurements was calculated by randomising the redshift of the clusters and repeating the stacking 500 times, as described in section 4.4.The detection is broad, ranging over 1 000 km/s in width for each line.This is likely the result of the offset in velocities in the different clusters along the background sightlines.In addi-

Fig. 9 .
Fig.9.Plot of 500 random stacks from the bootstrapping in grey.The stacks are from the same sample of 16 224 cluster-quasar pairs stacked at randomised redshifts from the full redshift distribution.The overlayed black dotted line represents the median flux per velocity bin of the 500 random stacks.On top, we over-plot in red the original blind stack from Figure8to illustrate the difference.

Fig. 10 .
Fig.10.Histogram of the measurements of the EWs of the 500 random stacks from bootstrapping for the MgII 2793 and 2803 Å lines stacked at random redshifts from the cluster redshift distribution.We mark with a vertical solid line the position of our measurement of EW 2796 corresponding to 0.056±0.015Å, and with a vertical dashed line the mean of EW 2796 of 0.003 Å.We measure a standard deviation from the distribution of EW 2796 of 0.015Å, which, considering our tentative detection, gives us a significance of 3.7σ.

Figure 11
Figure11displays our results together with the EWs from the above surveys.The left panel of Figure11displays the EW from the blind sample in our work compared to the values reported byAnand et al. (2022) andMishra & Muzahid (2022) at different values of projected distance.The right panel of Figure11shows the EWs as a function of cluster mass.Clearly, the works ofAnand et al. (2022) andMishra & Muzahid (2022) probe a lower mass range than our study.We compare the median of the two closest numbers in projected distance fromAnand et al. (2022), and the closest point fromMishra & Muzahid (2022) with our tentative detection (red circle).BothMishra & Muzahid (2022) andAnand et al. (2022) also report a clear trend of decreasing EW 2796 with increased projected distance.Despite the differences between the approaches, the reported EW is in line with results fromAnand et al. (2022).The cluster masses in our study range from M 500 = 1 × 10 14

Fig. 11 .Fig. 12 .
Fig. 11.Equivalent width and column density from the blind sample compared with results from the literature (Anand et al. 2021; Mishra & Muzahid 2022).The left panel displays the detection in the blind sample as a function of projected distance (in units of R 500 ), with the detection in our blind stack depicted as a red circle.In the right panel, we plot the mean of the two points of Anand et al. (2021) and the point of Mishra & Muzahid (2022) closest to our results in projected distance space (see left panel).This second figure highlights the difference in mass range between the samples.

Fig. 13 .
Fig. 13.Comparison between our blind stack with hydrodynamical cosmological simulations from TNG50(Nelson et al. 2020).The shape and colour of the symbols are as in previous figures.The purple symbols display predictions of the MgII column density in TNG50 as a function of projected distance in kpc.We stress that there are large mass differences between the models and the data: the TNG50 highest mass bin (∼ 10 13.5 M ⊙ ) displayed in the figure is smaller than our smallest mass cluster (3.28 × 10 14 M ⊙ ).Considering this mass difference, we note that there is likely an excess of MgII-absorbing cold gas in the TNG50 simulations.

Table 1 .
The angular separation is normalised by Number of spectra resulting from the various selection cuts.

Table 2 .
Properties of the two stacks performed.Table2lists the number of background quasar spectra stacked.The column of the number of quasar spectra corresponds to the fraction of quasars over the entire number of spectra with its representative percentage.The column density for the MgII-selected sample is frequency-weighted, considering the number of spectra used.Possible saturation (at EW 2796 > 0.15 Å) leads to a lower limit on the column density of the MgII-selected sample.MgII absorption is tentatively detected in the blind stack.This tentative detection provides a measure of the column density, which relates to the amount of cold gas in galaxy clusters.