<figure><img src="image_1.png" /></figure>

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

A.Y Fresco, C Peroux, A Merloni, J Comparat, R Szakacs, S Weng

► To cite this version:

A.Y Fresco, C Peroux, A Merloni, J Comparat, R Szakacs, et al. Searching for cold gas traced by MgII quasar absorbers in massive X-ray-selected galaxy clusters. Astron.Astrophys., 2024, 684, pp.A136. 10.1051/0004-6361/202346866. hal-04559381

HAL Id: hal-04559381

https://hal.science/hal-04559381

Submitted on 5 May 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_1.png" /></figure>

HAL
open science

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

A.Y Fresco, C Peroux, A Merloni, J Comparat, R Szakacs, S Weng

► To cite this version:

A.Y Fresco, C Peroux, A Merloni, J Comparat, R Szakacs, et al.. Searching for cold gas traced by MgII quasar absorbers in massive X-ray-selected galaxy clusters. Astronomy & Astrophysics - A&A, 2024, 684, pp.A136. 10.1051/0004-6361/202346866. hal-04559381

HAL Id: hal-04559381

https://hal.science/hal-04559381v1

Submitted on 5 May 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Searching for cold gas traced by MgII quasar absorbers in massive X-ray-selected galaxy clusters*

A. Fresco<fnref n="1" />, C. Péroux<fnref n="2" />,<fnref n="3" />, A. Merloni<fnref n="1" />, J. Comparat<fnref n="1" />, R. Szakacs<fnref n="2" />, and S. Weng<fnref n="2" />,<fnref n="4" />,<fnref n="5" />,<fnref n="6" />

Received 10 May 2023 / Accepted 29 January 2024

## ABSTRACT

*Context.* Almost 50% of galaxies in the local Universe are in clusters or groups coexisting with both hot and cold gas components. In the present study, we observationally probed the cold-gas content of X-ray-selected massive galaxy clusters with spectroscopic redshift measured from the SDSS/SPIDERS survey. This paper focuses on the most massive structures: galaxy clusters with a mean mass of $M_{500c} = 2.7 \times 10^{14} M_{\odot}$.

*Aims.* We used a large number of background quasar optical spectra from SDSS DR16 to probe the diffuse $T = 10^4$ K gas in their intracluster medium.

*Methods.* We first analysed a sample of spectra with known MgII absorbers, and then blindly stacked about 16 000 archival spectra at the redshifts of the foreground galaxy clusters.

*Results.* We tentatively (3.7σ significance) detect MgII in the clusters with an equivalent width EW(MgII λ2796) of 0.056 ± 0.015 Å, corresponding to a column density of log [N(MgII)/cm⁻<fnref n="2" />] = 12.12 ± 0.1. We tested our methodology by generating 22 000 mock SDSS spectra with MgII absorbers from Illustris-TNG50 cosmological magnetohydrodynamical simulations, combining photo-ionisation modelling and ray tracing. We also performed bootstrapping stacking at different cluster redshifts and stacked quasar spectra with no intervening clusters in the line of sight to measure the significance of our detection.

*Conclusions.* These results are in line with the findings of recent, similar observational studies but challenge predictions from Illustris-TNG simulations. Together, our findings indicate that large amounts of cold gas may be found in the most massive structures of the Universe.

**Key words.** galaxies: abundances – galaxies: clusters: intracluster medium – galaxies: evolution – galaxies: halos – quasars: absorption lines

## 1. 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_{500}^1$, 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. 2007; 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 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 multi-phase 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 higher-mass 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

<footnote n="1">¹ Max-Planck-Institut für Extraterrestrische Physik (MPE), Giessenbachstrasse 1, 85748 Garching, Germany
e-mail: yryfresquito@gmail.com</footnote>

<footnote n="2">² European Southern Observatory (ESO), Karl-Schwarzschildstrasse 2, 85748 Garching bei München, Germany</footnote>

<footnote n="3">³ Aix-Marseille Université, CNRS, LAM (Laboratoire d'Astrophysique de Marseille) UMR 7326, 13388 Marseille, France</footnote>

<footnote n="4">⁴ Sydney Institute for Astronomy, School of Physics A28, University of Sydney, NSW 2006, Australia</footnote>

<footnote n="5">⁵ ARC Centre of Excellence for All Sky Astrophysics in 3 Dimensions (ASTRO 3D), Mount Stromlo Rd, Australian Capital Territory 2611, Australia</footnote>

<footnote n="6">⁶ ATNF, CSIRO Space and Astronomy, PO Box 76, Epping, NSW 1710, Australia</footnote>

<footnote>* All the data used in this work are publicly available.</footnote>

<footnote n="1">¹ We define $R_{500}$ as the radius where the density is 500 times the critical density of the Universe at a given redshift.</footnote>
# Searching for cold gas traced by MgII quasar absorbers in massive X-ray-selected galaxy clusters*

A. Fresco<fnref n="1" />, C. Péroux<fnref n="2" />,<fnref n="3" />, A. Merloni<fnref n="1" />, J. Comparat<fnref n="1" />, R. Szakacs<fnref n="2" />, and S. Weng<fnref n="2" />,<fnref n="4" />,<fnref n="5" />,<fnref n="6" />

Received 10 May 2023 / Accepted 29 January 2024

## ABSTRACT

*Context.* Almost 50% of galaxies in the local Universe are in clusters or groups coexisting with both hot and cold gas components. In the present study, we observationally probed the cold-gas content of X-ray-selected massive galaxy clusters with spectroscopic redshift measured from the SDSS/SPIDERS survey. This paper focuses on the most massive structures: galaxy clusters with a mean mass of $M_{500c} = 2.7 \times 10^{14} M_{\odot}$.

*Aims.* We used a large number of background quasar optical spectra from SDSS DR16 to probe the diffuse $T = 10^4$ K gas in their intracluster medium.

*Methods.* We first analysed a sample of spectra with known MgII absorbers, and then blindly stacked about 16 000 archival spectra at the redshifts of the foreground galaxy clusters.

*Results.* We tentatively (3.7σ significance) detect MgII in the clusters with an equivalent width EW(MgII λ2796) of 0.056 ± 0.015 Å, corresponding to a column density of log [N(MgII)/cm⁻<fnref n="2" />] = 12.12 ± 0.1. We tested our methodology by generating 22 000 mock SDSS spectra with MgII absorbers from Illustris-TNG50 cosmological magnetohydrodynamical simulations, combining photo-ionisation modelling and ray tracing. We also performed bootstrapping stacking at different cluster redshifts and stacked quasar spectra with no intervening clusters in the line of sight to measure the significance of our detection.

*Conclusions.* These results are in line with the findings of recent, similar observational studies but challenge predictions from Illustris-TNG simulations. Together, our findings indicate that large amounts of cold gas may be found in the most massive structures of the Universe.

**Key words.** galaxies: abundances – galaxies: clusters: intracluster medium – galaxies: evolution – galaxies: halos – quasars: absorption lines

## 1. 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_{500}^1$, 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. 2007; 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 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 multi-phase 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 higher-mass 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

<footnote n="1">¹ Max-Planck-Institut für Extraterrestrische Physik (MPE), Giessenbachstrasse 1, 85748 Garching, Germany
e-mail: yryfresquito@gmail.com</footnote>

<footnote n="2">² European Southern Observatory (ESO), Karl-Schwarzschildstrasse 2, 85748 Garching bei München, Germany</footnote>

<footnote n="3">³ Aix-Marseille Université, CNRS, LAM (Laboratoire d'Astrophysique de Marseille) UMR 7326, 13388 Marseille, France</footnote>

<footnote n="4">⁴ Sydney Institute for Astronomy, School of Physics A28, University of Sydney, NSW 2006, Australia</footnote>

<footnote n="5">⁵ ARC Centre of Excellence for All Sky Astrophysics in 3 Dimensions (ASTRO 3D), Mount Stromlo Rd, Australian Capital Territory 2611, Australia</footnote>

<footnote n="6">⁶ ATNF, CSIRO Space and Astronomy, PO Box 76, Epping, NSW 1710, Australia</footnote>

<footnote>* All the data used in this work are publicly available.</footnote>

<footnote n="1">¹ We define $R_{500}$ as the radius where the density is 500 times the critical density of the Universe at a given redshift.</footnote>
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 multi-phase 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 \sim 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. 2024). Beyond the Milky Way observations, the gas is diffuse and challenging to detect in emission due to low surface 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. (2024) 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 Sect. 4 focuses on the results. In Sect. 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 Sect. 6. Here, we adopt an $H_0 = 67.74$ km s$^{-1}$ Mpc$^{-1}$, $\Omega_M = 0.3089$, and $\Omega_\Lambda = 0.6911$ cosmology. The equivalent width (EW) always refers to the MgII 2796 Å line unless stated otherwise.

## 2. Observational data

### 2.1. 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. 2016, 2020;
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 multi-phase 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 \sim 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. 2024). Beyond the Milky Way observations, the gas is diffuse and challenging to detect in emission due to low surface 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. (2024) 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 Sect. 4 focuses on the results. In Sect. 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 Sect. 6. Here, we adopt an $H_0 = 67.74$ km s$^{-1}$ Mpc$^{-1}$, $\Omega_M = 0.3089$, and $\Omega_\Lambda = 0.6911$ cosmology. The equivalent width (EW) always refers to the MgII 2796 Å line unless stated otherwise.

## 2. Observational data

### 2.1. 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. 2016, 2020;
<figure><img src="image_3.png" /><figcaption>**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.</figcaption></figure>

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 $\delta_z = 0.00049$ (147 km s$^{-1}$), while the maximum is $\delta_z = 0.0055$ (1650 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 \times 10^{14} M_{\odot}$ (converted from $M_{200c} = 3.9 \times 10^{14} M_{\odot}$). 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 $1 \text{ Mpc} < R_{200c} < 3 \text{ Mpc}$ depending on redshift. This radius is expressed as:

$$R_{200} = \left( \frac{3M_{200c}}{4\pi 200 \rho_{cr}} \left( \frac{H_z}{100} \right)^{-2} \right)^{1/3} \quad (1)$$

In order to compare these results to the latest related published research, we use the public open-source package COLOSSUS to convert between $M_{200c}$ and $M_{500c}$ (Diemer 2018), and continue to do so throughout the paper for consistency.

<figure><img src="image_4.png" /><figcaption>**Fig. 2.** Distribution of cluster redshifts from the parent SPIDERS sample (Clerc et al. 2020) used in this study. The clusters constitute a sub-sample 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$.</figcaption></figure>

<table><caption>**Table 1.** Number of spectra resulting from the various selection cuts.</caption><thead><tr><td>Sample</td><td>Number of objects</td></tr></thead><tbody><tr><td>Parent sample</td><td>750 414</td></tr><tr><td>$z_{\text{quasar}} > z_{\text{cluster}}$ and up to $3 \times R_{200c}$</td><td>434 736</td></tr><tr><td>MgII range and $z_{\text{quasar}} < 5$</td><td>18 694</td></tr><tr><td>High-quality SDSS spectra</td><td>16 224</td></tr></tbody></table>

**Notes.** Table 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 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 Fig. 2), with a typical statistical uncertainty on each redshift of $\Delta z/(1+z) = 6 \times 10^{-4}$ (Kirkpatrick et al. 2021). The number of spectroscopic members per system ranges between 3 and 75, with a mean of 12 members.

## 2.2. 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 $\lambda/\Delta\lambda \sim 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).

## 2.3. 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 quasar-cluster 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}$.
<figure><img src="image_3.png" /><figcaption>**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.</figcaption></figure>

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 $\delta_z = 0.00049$ (147 km s$^{-1}$), while the maximum is $\delta_z = 0.0055$ (1650 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 \times 10^{14} M_{\odot}$ (converted from $M_{200c} = 3.9 \times 10^{14} M_{\odot}$). 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 $1 \text{ Mpc} < R_{200c} < 3 \text{ Mpc}$ depending on redshift. This radius is expressed as:

$$R_{200} = \left( \frac{3M_{200c}}{4\pi 200 \rho_{cr}} \left( \frac{H_z}{100} \right)^{-2} \right)^{1/3} \quad (1)$$

In order to compare these results to the latest related published research, we use the public open-source package COLOSSUS to convert between $M_{200c}$ and $M_{500c}$ (Diemer 2018), and continue to do so throughout the paper for consistency.

<figure><img src="image_4.png" /><figcaption>**Fig. 2.** Distribution of cluster redshifts from the parent SPIDERS sample (Clerc et al. 2020) used in this study. The clusters constitute a sub-sample 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$.</figcaption></figure>

<table><caption>**Table 1.** Number of spectra resulting from the various selection cuts.</caption><thead><tr><td>Sample</td><td>Number of objects</td></tr></thead><tbody><tr><td>Parent sample</td><td>750 414</td></tr><tr><td>$z_{\text{quasar}} > z_{\text{cluster}}$ and up to $3 \times R_{200c}$</td><td>434 736</td></tr><tr><td>MgII range and $z_{\text{quasar}} < 5$</td><td>18 694</td></tr><tr><td>High-quality SDSS spectra</td><td>16 224</td></tr></tbody></table>

**Notes.** Table 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 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 Fig. 2), with a typical statistical uncertainty on each redshift of $\Delta z/(1+z) = 6 \times 10^{-4}$ (Kirkpatrick et al. 2021). The number of spectroscopic members per system ranges between 3 and 75, with a mean of 12 members.

## 2.2. 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 $\lambda/\Delta\lambda \sim 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).

## 2.3. 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 quasar-cluster 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}$.
<figure><img src="image_5.png" /><figcaption>**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 2881 MgII known absorbers redshift from Anand et al. (2021) found in these quasars, and the blue as the redshift of the foreground clusters (1066 clusters).</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Fig. 4.** Distribution of the S/N of the spectrum of 16 224 background quasars (see text) from the SDSS DR16 sample with foreground clusters. The median S/N of spectra is S/N ~ 3.79 as indicated by the vertical dotted line. The 10th percentile is a dashed-dotted line, and the 90th percentile is the solid line.</figcaption></figure>

(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 cluster-quasar 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 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 Fig. 4.

For each of the pairs, we calculate the angular separation between the foreground cluster and the background quasar projected on the sky plane. The angular separation is normalised by 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.

## 2.4. 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 cross-matching our sample of cluster-quasar pairs with this absorption catalogue, we get a match for 4150 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:

$$ \Delta v = \frac{\Delta z \times c}{1 + z_{\text{cluster}}} \qquad (2) $$

where c is the speed of light in km s$^{-1}$ 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).

# 3. Analysis

## 3.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 displays an example spectrum showing both the quasar continuum fit and 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
<figure><img src="image_5.png" /><figcaption>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 2881 MgII known absorbers redshift from Anand et al. (2021) found in these quasars, and the blue as the redshift of the foreground clusters (1066 clusters).</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 4. Distribution of the S/N of the spectrum of 16 224 background quasars (see text) from the SDSS DR16 sample with foreground clusters. The median S/N of spectra is S/N ~ 3.79 as indicated by the vertical dotted line. The 10th percentile is a dashed-dotted line, and the 90th percentile is the solid line.</figcaption></figure>

(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 cluster-quasar 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 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 Fig. 4.

For each of the pairs, we calculate the angular separation between the foreground cluster and the background quasar projected on the sky plane. The angular separation is normalised by 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.

## 2.4. 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 cross-matching our sample of cluster-quasar pairs with this absorption catalogue, we get a match for 4150 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:

$$ \Delta v = \frac{\Delta z \times c}{1 + z_{\text{cluster}}} \qquad (2) $$

where c is the speed of light in km s$^{-1}$ 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).

# 3. Analysis

## 3.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 displays an example spectrum showing both the quasar continuum fit and 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
<figure><img src="image_7.png" /><figcaption>Fig. 5. Velocity difference, $\Delta v$, between the known MgII absorber from Anand et al. (2021) within $3 \times R_{200}$ from the cluster centre and the cluster redshift, as defined in Eq. (2). The figure displays the full range of $\Delta 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.</figcaption></figure>

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.

## 3.2. 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 $\Delta z/(1+z) = 6 \times 10^{-4}$, which is about 180 km s$^{-1}$. 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 available 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.

## 3.3. 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$^{-1}$ 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 Figs. 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 −500 km s$^{-1}$ blueward and +375 km s$^{-1}$ redward. The continuum level is defined as the observed spectrum divided by the fitted continuum as previously described in Sect. 3.1. We then measured the second line of the MgII doublet, the MgII 2803 Å line, which is located at +750 km s$^{-1}$ redward. The same method is applied to the MgII 2803 Å line, which is integrated from −375 km s$^{-1}$ (blueward) to +500 km s$^{-1}$ (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:

$$ EW_{\text{rest}} = \frac{EW_{\text{obs}}}{(1 + \langle z_{\text{cluster}} \rangle)}, \qquad (3) $$

where the redshift $\langle z_{\text{cluster}} \rangle$ 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 Sect. 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):

$$ N[\text{MgII}] = 1.13 \times 10^{20} \left( \frac{EW_{\text{rest}}}{f_{\text{osc}} \lambda_{\text{rest}}^2} \right) \text{cm}^{-2}, \qquad (4) $$

where $\lambda_{\text{rest}}$ is the rest wavelength of each of the MgII lines of the doublet, $\lambda_{\text{rest}} = 2796$ Å and $\lambda_{\text{rest}} = 2803$ Å. The corresponding oscillator strengths are $f_{\text{osc}} = 0.615$ and $f_{\text{osc}} = 0.306$. For saturated absorbers with EW$_{2796} > 0.15$ Å, we derive a lower limit on the MgII column density (see Table 2).

# 4. Results

## 4.1. 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. The corresponding measurements are listed in Table 2 and the stacked spectrum is displayed in Fig. 7. The detection of the MgII doublet in the stack validates the methodology by showing the very strong absorption of the MgII doublet present at the redshift of the clusters. Also, despite the high-velocity dispersions within the clusters, this detection demonstrates that the cold $10^4$ K gas can be traced by MgII absorbers observed in the spectrum of background quasars. The detected
<figure><img src="image_7.png" /><figcaption>Fig. 5. Velocity difference, Δv, between the known MgII absorber from Anand et al. (2021) within 3 × R$_{200}$ from the cluster centre and the cluster redshift, as defined in Eq. (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.</figcaption></figure>

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.

## 3.2. 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$^{-1}$. 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 available 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.

## 3.3. 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$^{-1}$ 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 Figs. 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 −500 km s$^{-1}$ blueward and +375 km s$^{-1}$ redward. The continuum level is defined as the observed spectrum divided by the fitted continuum as previously described in Sect. 3.1. We then measured the second line of the MgII doublet, the MgII 2803 Å line, which is located at +750 km s$^{-1}$ redward. The same method is applied to the MgII 2803 Å line, which is integrated from −375 km s$^{-1}$ (blueward) to +500 km s$^{-1}$ (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:

$$ EW_{\text{rest}} = \frac{EW_{\text{obs}}}{(1 + \langle z_{\text{cluster}} \rangle)}, \qquad (3) $$

where the redshift $\langle z_{\text{cluster}} \rangle$ 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 Sect. 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):

$$ N[\text{MgII}] = 1.13 \times 10^{20} \left( \frac{EW_{\text{rest}}}{f_{\text{osc}} \lambda_{\text{rest}}^2} \right) \text{cm}^{-2}, \qquad (4) $$

where $\lambda_{\text{rest}}$ is the rest wavelength of each of the MgII lines of the doublet, $\lambda_{\text{rest}} = 2796$ Å and $\lambda_{\text{rest}} = 2803$ Å. The corresponding oscillator strengths are $f_{\text{osc}} = 0.615$ and $f_{\text{osc}} = 0.306$. For saturated absorbers with EW$_{2796}$ > 0.15 Å, we derive a lower limit on the MgII column density (see Table 2).

# 4. Results

## 4.1. 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. The corresponding measurements are listed in Table 2 and the stacked spectrum is displayed in Fig. 7. The detection of the MgII doublet in the stack validates the methodology by showing the very strong absorption of the MgII doublet present at the redshift of the clusters. Also, despite the high-velocity dispersions within the clusters, this detection demonstrates that the cold $10^4$ K gas can be traced by MgII absorbers observed in the spectrum of background quasars. The detected
<figure><img src="image_8.png" /><figcaption>**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.</figcaption></figure>

<table><caption>**Table 2.** Properties of the two stacks performed.</caption><thead><tr><th>Sample</th><th># Quasar Spectra</th><th>Mean MgII Redshift</th><th>S/N MgII sample</th><th>EW MgII 2796 Å</th><th>EW MgII 2803 Å</th><th>MgII Doublet Ratio</th><th>Col Den 2796 log N(MgII) cm<sup>-2</sup></th><th>Col Den 2803 log N(MgII) cm<sup>-2</sup></th><th>Frequency Weighted 2796 log N(MgII) cm<sup>-2</sup></th></tr></thead><tbody><tr><td>MgII selected</td><td>32/16224 = 0.2%</td><td>0.51</td><td>73</td><td>0.35 ± 0.015</td><td>0.19 ± 0.015</td><td>1.8</td><td>≥12.92</td><td>12.95 ± 0.1</td><td>≥10.22</td></tr><tr><td>Blind</td><td>16224 = 100%</td><td>0.45</td><td>195</td><td>0.056 ± 0.015</td><td>0.024 ± 0.015</td><td>2.3</td><td>12.12 ± 0.1</td><td>12.05 ± 0.1</td><td>12.12 ± 0.1</td></tr></tbody></table>

**Notes.** Table 2 lists 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.

lines are strong and show indications of saturation. Given the saturated EWs in the range $0.15 < EW_{2796} < 2$ Å (Churchill et al. 2000), we calculate a lower limit for the column density using Eq. (4). The resulting $EW_{2796}$ is $EW_{2796} = 0.35 \pm 0.015$ Å for the MgII 2796Å line, corresponding to a column density of $log [N(MgII)/cm^{-2}] = 12.92$. The bluest line of the MgII doublet, 2796Å, is located at $-100$ km s$^{-1}$ 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$^{-1}$, 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 col-

umn densities as follows: $EW = EW_{rest}f_c$, where the covering fraction, $f_c$, is defined as the number of spectra with known absorbers over the total number of spectra available for each bin of mass and angular separation. These values are listed in Table 2.

## 4.2. The blind sample

Figure 8 displays the stack of all the spectra, with a total of 16224 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 $\lambda2796$ Å of $EW_{2796} = 0.056 \pm 0.015$ Å ($3.7\sigma$ significance), and a column density of $log [N(MgII)/cm^{-2}] = 12.12 \pm 0.1$. The mean mass of clusters contributing to the stacked
<figure><img src="image_8.png" /><figcaption>**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.</figcaption></figure>

<table><caption>**Table 2.** Properties of the two stacks performed.</caption><thead><tr><th>Sample</th><th># Quasar Spectra</th><th>Mean MgII Redshift</th><th>S/N MgII sample</th><th>EW MgII 2796 Å</th><th>EW MgII 2803 Å</th><th>MgII Doublet Ratio</th><th>Col Den 2796 log N(MgII) cm<sup>-2</sup></th><th>Col Den 2803 log N(MgII) cm<sup>-2</sup></th><th>Frequency Weighted 2796 log N(MgII) cm<sup>-2</sup></th></tr></thead><tbody><tr><td>MgII selected</td><td>32/16224 = 0.2%</td><td>0.51</td><td>73</td><td>0.35 ± 0.015</td><td>0.19 ± 0.015</td><td>1.8</td><td>≥12.92</td><td>12.95 ± 0.1</td><td>≥10.22</td></tr><tr><td>Blind</td><td>16224 = 100%</td><td>0.45</td><td>195</td><td>0.056 ± 0.015</td><td>0.024 ± 0.015</td><td>2.3</td><td>12.12 ± 0.1</td><td>12.05 ± 0.1</td><td>12.12 ± 0.1</td></tr></tbody></table>

**Notes.** Table 2 lists 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.

lines are strong and show indications of saturation. Given the saturated EWs in the range $0.15 < EW_{2796} < 2$ Å (Churchill et al. 2000), we calculate a lower limit for the column density using Eq. (4). The resulting $EW_{2796}$ is $EW_{2796} = 0.35 \pm 0.015$ Å for the MgII 2796Å line, corresponding to a column density of $log [N(MgII)/cm^{-2}] = 12.92$. The bluest line of the MgII doublet, 2796Å, is located at $-100$ km s$^{-1}$ 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$^{-1}$, 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 col-

umn densities as follows: $EW = EW_{rest}f_c$, where the covering fraction, $f_c$, is defined as the number of spectra with known absorbers over the total number of spectra available for each bin of mass and angular separation. These values are listed in Table 2.

## 4.2. The blind sample

Figure 8 displays the stack of all the spectra, with a total of 16224 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 $\lambda2796$ Å of $EW_{2796} = 0.056 \pm 0.015$ Å ($3.7\sigma$ significance), and a column density of $log [N(MgII)/cm^{-2}] = 12.12 \pm 0.1$. The mean mass of clusters contributing to the stacked
<figure><img src="image_9.png" /></figure>

**Fig. 7.** Stack of quasar spectra with a sightline within $3 \times R_{200}$ of the centre of a cluster, and with known MgII absorbers from the MgII absorption catalogue of Anand et al. (2021) satisfying $|d\nu| < 2000 \text{ km s}^{-1}$. The figures display the normalised flux in velocity space. The first dashed line at $-100 \text{ km s}^{-1}$ is fixed at the location of the minimum flux within $[-2000:2000] \text{ km s}^{-1}$ of the redshift of the known MgII absorber. The second dashed line is placed at $+750 \text{ km s}^{-1}$ from the first dashed line, which corresponds to the separation of the two MgII lines in the doublet. The solid vertical lines are fixed at $-500 \text{ km s}^{-1}$ from 2976 and $+500 \text{ km s}^{-1}$ from 2803 line. The green line shows the fitted continuum used for the EW measurement.

<figure><img src="image_10.png" /></figure>

**Fig. 8.** Blind sample stack of all the background quasar spectra available at the cluster redshift, showing the normalised flux in log in velocity space. This sample comprises a total of 16 224 spectra and leads to an S/N of the stack of 195. The detection has a rest EW of $EW_{2796} = 0.056 \pm 0.015$ Å (3.7$\sigma$ significance), corresponding to a column density of log [$N(\text{MgII})/\text{cm}^{-2}$] = 12.12$\pm$0.1. The green line shows the fitted continuum used for the EW measurement.

spectra is $9.7 \times 10^{14} M_{\odot}$ ($M_{500}$), and a mean angular separation of the cluster-quasar pairs is about $2 \times 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 Sect. 4.4.

The detection is broad, ranging over $1000 \text{ km s}^{-1}$ in width for each line. This is likely the result of the offset in velocities in the different clusters along the background sightlines. In addition, 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.

### 4.3. 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 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 Illustris-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} \text{ 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] \text{ km s}^{-1}$ 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 Å.

### 4.4. 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 Sect. 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
<figure><img src="image_9.png" /></figure>

**Fig. 7.** Stack of quasar spectra with a sightline within $3 \times R_{200}$ of the centre of a cluster, and with known MgII absorbers from the MgII absorption catalogue of Anand et al. (2021) satisfying $|d\nu| < 2000 \text{ km s}^{-1}$. The figures display the normalised flux in velocity space. The first dashed line at $-100 \text{ km s}^{-1}$ is fixed at the location of the minimum flux within $[-2000:2000] \text{ km s}^{-1}$ of the redshift of the known MgII absorber. The second dashed line is placed at $+750 \text{ km s}^{-1}$ from the first dashed line, which corresponds to the separation of the two MgII lines in the doublet. The solid vertical lines are fixed at $-500 \text{ km s}^{-1}$ from 2976 and $+500 \text{ km s}^{-1}$ from 2803 line. The green line shows the fitted continuum used for the EW measurement.

<figure><img src="image_10.png" /></figure>

**Fig. 8.** Blind sample stack of all the background quasar spectra available at the cluster redshift, showing the normalised flux in log in velocity space. This sample comprises a total of 16 224 spectra and leads to an S/N of the stack of 195. The detection has a rest EW of $EW_{2796} = 0.056 \pm 0.015$ Å (3.7$\sigma$ significance), corresponding to a column density of log [$N(\text{MgII})/\text{cm}^{-2}$] = 12.12$\pm$0.1. The green line shows the fitted continuum used for the EW measurement.

spectra is $9.7 \times 10^{14} M_{\odot}$ ($M_{500}$), and a mean angular separation of the cluster-quasar pairs is about $2 \times 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 Sect. 4.4.

The detection is broad, ranging over $1000 \text{ km s}^{-1}$ in width for each line. This is likely the result of the offset in velocities in the different clusters along the background sightlines. In addition, 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.

### 4.3. 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 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 Illustris-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} \text{ 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] \text{ km s}^{-1}$ 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 Å.

### 4.4. 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 Sect. 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
<figure><img src="image_11.png" /></figure>

**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 Fig. 8 to illustrate the difference.

<figure><img src="image_12.png" /></figure>

**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<sub>2796</sub> corresponding to 0.056±0.015 Å, and with a vertical dashed line the mean of EW<sub>2796</sub> of 0.003 Å. We measure a standard deviation from the distribution of EW<sub>2796</sub> of 0.015 Å, which, considering our tentative detection, gives us a significance of 3.7σ.

one final median array shown as the black dotted line in Fig. 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 over-plot in red the blind stack from Fig. 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 Fig. 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.

# 5. Discussion

## 5.1. 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⁴ 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⁴ 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⁴ 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 Fig. 11 the EW (left y-axis) and column density (right y-axis) of MgII-absorbing 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.

## 5.2. Cold gas in galaxy clusters

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_\odot$) have an overdensity of absorbers with EW<sub>2796</sub> > 1 Å that is twice higher than for moderate mass clusters ($M \sim 2 \times 10^{13} M_\odot$). The effect is even more pronounced at smaller distances ($d < 1$ Mpc) 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
<figure><img src="image_11.png" /></figure>

**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 Fig. 8 to illustrate the difference.

<figure><img src="image_12.png" /></figure>

**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<sub>2796</sub> corresponding to 0.056±0.015 Å, and with a vertical dashed line the mean of EW<sub>2796</sub> of 0.003 Å. We measure a standard deviation from the distribution of EW<sub>2796</sub> of 0.015 Å, which, considering our tentative detection, gives us a significance of 3.7σ.

one final median array shown as the black dotted line in Fig. 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 over-plot in red the blind stack from Fig. 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 Fig. 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.

# 5. Discussion

## 5.1. 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⁴ 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⁴ 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⁴ 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 Fig. 11 the EW (left y-axis) and column density (right y-axis) of MgII-absorbing 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.

## 5.2. Cold gas in galaxy clusters

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_\odot$) have an overdensity of absorbers with EW<sub>2796</sub> > 1 Å that is twice higher than for moderate mass clusters ($M \sim 2 \times 10^{13} M_\odot$). The effect is even more pronounced at smaller distances ($d < 1$ Mpc) 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
<figure><img src="image_13.png" /><figcaption>**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.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Fig. 12.** Comparison between our results and measurements of column density and EW of MgII absorption in LRGs from Lan & Mo (2018), LRGs and cluster samples from Anand et al. (2022), and our tentative detection from the blind sample, as a function of projected distance.</figcaption></figure>

quasar spectra in this latter work, the authors report that the MgII absorber detection rate per quasar is $2.70 \pm 0.66$ times higher inside the clusters than outside them. This shows that Mg II absorbers are abundant in clusters compared to the field.

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_{\odot}$. 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.

Figure 11 displays our results together with the EWs from the above surveys. The left panel of Fig. 11 displays the EW from the blind sample in our work compared to the values reported by Anand et al. (2022) and Mishra & Muzahid (2022) at different values of projected distance. The right panel of Fig. 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 \times 10^{14}$ to $8.4 \times 10^{15}$, with an EW$_{2796}$ of $0.056 \pm 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 Fig. 11.

Our higher EW measurements and limits could be associated with the higher mass in our sample, as shown in Fig. 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 \sim 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 \sim 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} \ge 0.5$ Å) in the lower-redshift ($z \sim 1$) studies (i.e. at higher mass) from the MusE GAs FLOw and Wind survey or MEGAFLOW (Schroetter et al. 2016, 2019) than at $z \sim 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
<figure><img src="image_13.png" /><figcaption>**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.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Fig. 12.** Comparison between our results and measurements of column density and EW of MgII absorption in LRGs from Lan & Mo (2018), LRGs and cluster samples from Anand et al. (2022), and our tentative detection from the blind sample, as a function of projected distance.</figcaption></figure>

quasar spectra in this latter work, the authors report that the MgII absorber detection rate per quasar is 2.70 ± 0.66 times higher inside the clusters than outside them. This shows that Mg II absorbers are abundant in clusters compared to the field.

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_{\odot}$. 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.

Figure 11 displays our results together with the EWs from the above surveys. The left panel of Fig. 11 displays the EW from the blind sample in our work compared to the values reported by Anand et al. (2022) and Mishra & Muzahid (2022) at different values of projected distance. The right panel of Fig. 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 \times 10^{14}$ to $8.4 \times 10^{15}$, with an EW$_{2796}$ of $0.056 \pm 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 Fig. 11.

Our higher EW measurements and limits could be associated with the higher mass in our sample, as shown in Fig. 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 \sim 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 \sim 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} \ge 0.5$ Å) in the lower-redshift ($z \sim 1$) studies (i.e. at higher mass) from the MusE GAs FLOw and Wind survey or MEGAFLOW (Schroetter et al. 2016, 2019) than at $z \sim 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¹¹.² $M_\odot$ 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 Fig. 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 Fig. 12. Similarly, Anand et al. (2022) perform a 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}$.

## 5.3. 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_\odot$. 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 Fig. 13, we compare our results with the TNG50 simulation (Pillepich et al. 2019; Nelson et al. 2019) from the Illustris-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 Illustris-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³ resolution elements in a ~50 Mpc (comoving) box. Nelson et al. (2020) make predictions of the physi-

cal 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 2024). 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 (2024) 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.

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 VI (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 Eq. (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_\odot$) displayed in the figure is smaller than our smallest mass cluster ($3.28 \times 10^{14} M_\odot$). 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_\odot) > 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.
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¹¹.² $M_\odot$ 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 Fig. 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 Fig. 12. Similarly, Anand et al. (2022) perform a 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}$.

## 5.3. 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_\odot$. 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 Fig. 13, we compare our results with the TNG50 simulation (Pillepich et al. 2019; Nelson et al. 2019) from the Illustris-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 Illustris-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 \times 2160^3$ resolution elements in a ~50 Mpc (comoving) box. Nelson et al. (2020) make predictions of the physi-

cal 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 2024). 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 (2024) 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.

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 VI (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 Eq. (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 ($\sim 10^{13.5} M_\odot$) displayed in the figure is smaller than our smallest mass cluster ($3.28 \times 10^{14} M_\odot$). 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_\odot) > 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.
<figure><img src="image_15.png" /><figcaption>**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_\odot$) displayed in the figure is smaller than our smallest mass cluster ($3.28 \times 10^{14} M_\odot$). Considering this mass difference, we note that there is likely an excess of MgII-absorbing cold gas in the TNG50 simulations.</figcaption></figure>

## 5.4. 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 low- and medium-resolution quasar spectra in extremely large numbers, notably the DESI experiment (DESI Collaboration 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.

## 6. 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 can be summarised as follows:

- From the sample of known MgII absorbers within [-2000:2000] km s⁻¹ of the clusters, we detect strong absorption rest EW for MgII λ2796 of EW<sub>2796</sub> = 0.35 ± 0.015 Å. The uncertainty is calculated through bootstrapping. As the absorption lines are saturated (at EW<sub>2796</sub> > 0.15 Å), we derive a lower limit for the column density of

log [N(MgII)/cm⁻²] ≥ 12.92. These results validate the technique and demonstrate that despite the high velocity dispersions within the clusters, the cold 10⁴K gas can be traced by MgII absorbers observed in the spectra of background quasars.

- We then stacked quasar spectra without prior knowledge of the presence of absorption systems (so-called blind sample). We tentatively (≈3.7σ significance) detect MgII absorption signal by stacking a total of 16 224 quasar absorber-cluster pairs within [-2000:2000] km s⁻¹ from the clusters, with EW<sub>2796</sub> = 0.056 ± 0.015 Å, corresponding to a column density of log [N(MgII)/cm⁻²] = 12.12 ± 0.1 cm⁻².

- 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 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.

*Acknowledgements.* This research made use of Astropy (http://www.astropy.org), a community-developed core Python package for Astronomy (Astropy Collaboration 2018). We thank Dylan Nelson for generating the TNG50 MgII spectra. A. Fresco gratefully acknowledges the support provided by the fellowship from Becas Carlos Antonio Lopez (BECAL), Paraguay.

## References

Ahumada, R., Prieto, C. A., Almeida, A., et al. 2020, ApJS, 249, 3

Aihara, H., Arimoto, N., Armstrong, R., et al. 2018, PASJ, 70, S4

Anand, A., Nelson, D., & Kauffmann, G. 2021, MNRAS, 504, 65

Anand, A., Kauffmann, G., & Nelson, D. 2022, MNRAS, 513, 3210

Astropy Collaboration (Price-Whelan, A. M., et al.) 2018, AJ, 156, 123

Augustin, R., Péroux, C., Hamanowicz, A., et al. 2021, MNRAS, 505, 6195

Ayromlou, M., Nelson, D., Yates, R. M., Kauffmann, G., & White, S. D. M. 2019, MNRAS, 487, 4313

Bahé, Y. M., McCarthy, I. G., Balogh, M. L., & Font, A. S. 2013, MNRAS, 430, 3017

Bird, S., Garnett, R., & Ho, S. 2017, MNRAS, 466, 2111

Böhringer, H., & Werner, N. 2010, A&A Rev., 18, 127

Bouché, N., Murphy, M. T., Péroux, C., Csabai, I., & Wild, V. 2006, MNRAS, 371, 495

Bowen, D. V., & Chelouche, D. 2011, ApJ, 727, 47

Butsky, I. S., Burchett, J. N., Nagai, D., et al. 2019, MNRAS, 490, 4292

Cantalupo, S., Pezzulli, G., Lilly, S. J., et al. 2019, MNRAS, 483, 5188

Cen, R., & Ostriker, J. P. 1999, ApJ, 519, L109

Churchill, C. W., Mellon, R. R., Charlton, J. C., et al. 2000, ApJ, 543, 577

Clerc, N., Adami, C., Lieu, M., et al. 2014, MNRAS, 444, 2723

Clerc, N., Merloni, A., Zhang, Y. Y., et al. 2016, MNRAS, 463, 4490

Clerc, N., Kirkpatrick, C. C., Finoguenov, A., et al. 2020, MNRAS, 497, 3976

Croston, J. H., Pratt, G. W., Böhringer, H., et al. 2008, A&A, 487, 431

Das, S., Mathur, S., Gupta, A., & Krongold, Y. 2021, ApJ, 918, 83

Davé, R., Cen, R., Ostriker, J. P., et al. 2001, ApJ, 552, 473

DESI Collaboration (Aghamousa, A., et al.) 2016, ArXiv e-prints [arXiv:1611.00036]
<figure><img src="image_15.png" /><figcaption>**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_\odot$) displayed in the figure is smaller than our smallest mass cluster ($3.28 \times 10^{14} M_\odot$). Considering this mass difference, we note that there is likely an excess of MgII-absorbing cold gas in the TNG50 simulations.</figcaption></figure>

## 5.4. 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 low- and medium-resolution quasar spectra in extremely large numbers, notably the DESI experiment (DESI Collaboration 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.

## 6. 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 can be summarised as follows:

- From the sample of known MgII absorbers within [-2000:2000] km s⁻¹ of the clusters, we detect strong absorption rest EW for MgII λ2796 of EW<sub>2796</sub> = 0.35 ± 0.015 Å. The uncertainty is calculated through bootstrapping. As the absorption lines are saturated (at EW<sub>2796</sub> > 0.15 Å), we derive a lower limit for the column density of

log [N(MgII)/cm⁻²] ≥ 12.92. These results validate the technique and demonstrate that despite the high velocity dispersions within the clusters, the cold 10⁴K gas can be traced by MgII absorbers observed in the spectra of background quasars.

- We then stacked quasar spectra without prior knowledge of the presence of absorption systems (so-called blind sample). We tentatively (≈3.7σ significance) detect MgII absorption signal by stacking a total of 16 224 quasar absorber-cluster pairs within [-2000:2000] km s⁻¹ from the clusters, with EW<sub>2796</sub> = 0.056 ± 0.015 Å, corresponding to a column density of log [N(MgII)/cm⁻²] = 12.12 ± 0.1 cm⁻².

- 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 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.

*Acknowledgements.* This research made use of Astropy (http://www.astropy.org), a community-developed core Python package for Astronomy (Astropy Collaboration 2018). We thank Dylan Nelson for generating the TNG50 MgII spectra. A. Fresco gratefully acknowledges the support provided by the fellowship from Becas Carlos Antonio Lopez (BECAL), Paraguay.

## References

Ahumada, R., Prieto, C. A., Almeida, A., et al. 2020, ApJS, 249, 3

Aihara, H., Arimoto, N., Armstrong, R., et al. 2018, PASJ, 70, S4

Anand, A., Nelson, D., & Kauffmann, G. 2021, MNRAS, 504, 65

Anand, A., Kauffmann, G., & Nelson, D. 2022, MNRAS, 513, 3210

Astropy Collaboration (Price-Whelan, A. M., et al.) 2018, AJ, 156, 123

Augustin, R., Péroux, C., Hamanowicz, A., et al. 2021, MNRAS, 505, 6195

Ayromlou, M., Nelson, D., Yates, R. M., Kauffmann, G., & White, S. D. M. 2019, MNRAS, 487, 4313

Bahé, Y. M., McCarthy, I. G., Balogh, M. L., & Font, A. S. 2013, MNRAS, 430, 3017

Bird, S., Garnett, R., & Ho, S. 2017, MNRAS, 466, 2111

Böhringer, H., & Werner, N. 2010, A&A Rev., 18, 127

Bouché, N., Murphy, M. T., Péroux, C., Csabai, I., & Wild, V. 2006, MNRAS, 371, 495

Bowen, D. V., & Chelouche, D. 2011, ApJ, 727, 47

Butsky, I. S., Burchett, J. N., Nagai, D., et al. 2019, MNRAS, 490, 4292

Cantalupo, S., Pezzulli, G., Lilly, S. J., et al. 2019, MNRAS, 483, 5188

Cen, R., & Ostriker, J. P. 1999, ApJ, 519, L109

Churchill, C. W., Mellon, R. R., Charlton, J. C., et al. 2000, ApJ, 543, 577

Clerc, N., Adami, C., Lieu, M., et al. 2014, MNRAS, 444, 2723

Clerc, N., Merloni, A., Zhang, Y. Y., et al. 2016, MNRAS, 463, 4490

Clerc, N., Kirkpatrick, C. C., Finoguenov, A., et al. 2020, MNRAS, 497, 3976

Croston, J. H., Pratt, G. W., Böhringer, H., et al. 2008, A&A, 487, 431

Das, S., Mathur, S., Gupta, A., & Krongold, Y. 2021, ApJ, 918, 83

Davé, R., Cen, R., Ostriker, J. P., et al. 2001, ApJ, 552, 473

DESI Collaboration (Aghamousa, A., et al.) 2016, ArXiv e-prints [arXiv:1611.00036]
Diemer, B. 2018, ApJS, 239, 35

Donahue, M., & Voit, G. M. 2022, Phys. Rep., 973, 1

Dutta, R., Fumagalli, M., Fossati, M., et al. 2020, MNRAS, 499, 5022

Eke, V. R., Baugh, C. M., Cole, S., et al. 2004, MNRAS, 348, 866

Ettori, S. 2003, MNRAS, 344, L13

Ferland, G. J., Chatzikos, M., Guzmán, F., et al. 2017, Rev. Mex. A&A, 53, 385

Finoguenov, A., Rykoff, E., Clerc, N., et al. 2020, A&A, 638, A114

Frebel, A. 2018, Ann. Rev. Nucl. Part. Sci., 68, 237

Fresco, A. Y., Péroux, C., Merloni, A., Hamanowicz, A., & Szakacs, R. 2020, MNRAS, 499, 5230

Fukugita, M., & Peebles, P. J. E. 2004, ApJ, 616, 643

Gad, A. F. 2021, ArXiv e-prints [arXiv:2106.06158]

Girardi, M., Biviano, A., Giuricin, G., Mardirossian, F., & Mezzetti, M. 1993, ApJ, 404, 38

Gladders, M. D., & Yee, H. K. C. 2005, ApJS, 157, 1

Gonzalez, A. H., Zaritsky, D., & Zabludoff, A. I. 2007, ApJ, 666, 147

Gonzalez, A. H., Sivanandam, S., Zabludoff, A. I., & Zaritsky, D. 2013, ApJ, 778, 14

Guo, H., Shen, Y., & Wang, S. 2018, Astrophysics Source Code Library[record ascl:1809.008]

Hamanowicz, A., Péroux, C., Zwaan, M. A., et al. 2020, MNRAS, 492, 2347

Hummels, C. B., Smith, B. D., & Silvia, D. W. 2017, ApJ, 847, 59

Ider Chitham, J., Comparat, J., Finoguenov, A., et al. 2020, MNRAS, 499, 4768

Jones, C., & Forman, W. 1999, ApJ, 511, 65

Kirkpatrick, C. C., & McNamara, B. R. 2015, MNRAS, 452, 4361

Kirkpatrick, C. C., Clerc, N., Finoguenov, A., et al. 2021, MNRAS, 503, 5763

Kraljic, K., Laigle, C., Pichon, C., et al. 2022, MNRAS, 514, 1359

Kravtsov, A. V., & Borgani, S. 2012, ARA&A, 50, 353

Kravtsov, A. V., Nagai, D., & Vikhlinin, A. A. 2005, ApJ, 625, 588

Kunz, M. W., Jones, T. W., & Zhuravleva, I. 2022, in Plasma Physics of the Intracluster Medium, eds. C. Bambi, A. Santangelo, Handbook of X-ray and Gamma-ray Astrophysics (Springer, Singapore)

Laganá, T. F., Martinet, N., Durret, F., et al. 2013, A&A, 555, A66

Lan, T.-W., & Mo, H. 2018, ApJ, 866, 36

Lanzetta, K. M., & Bowen, D. 1990, ApJ, 357, 321

Lee, J. C., Hwang, H. S., & Song, H. 2021, MNRAS, 503, 4309

Liu, A., Bulbul, E., Ghirardini, V., et al. 2022, A&A, 661, A2

Locatelli, N., Ponti, G., Zheng, X., et al. 2024, A&A, 681, A78

Lopez, S., Barrientos, L. F., Lira, P., et al. 2008, ApJ, 679, 1144

Lyke, B. W., Higley, A. N., McLane, J. N., et al. 2020, ApJS, 250, 8

Marinacci, F., Vogelsberger, M., Pakmor, R., et al. 2018, MNRAS, 480, 5113

Mas-Ribas, L., Miralda-Escudé, J., Pérez-Ràfols, I., et al. 2017, ApJ, 846, 4

Mathur, S., Gupta, A., Das, S., Krongold, Y., & Nicastro, F. 2021, ApJ, 908, 69

McCourt, M., Sharma, P., Quataert, E., & Parrish, I. J. 2012, MNRAS, 419, 3319

Merloni, A. 2012, Science from the Next Generation Imaging and Spectroscopic Surveys, 43

Merloni, A., Alexander, D. A., Banerji, M., et al. 2019, The Messenger, 175, 42

Mernier, F., Biffi, V., Yamaguchi, H., et al. 2018, Space Sci. Rev., 214, 129

Mishra, S., & Muzahid, S. 2022, ApJ, 933, 229

Mishra, S., Muzahid, S., Dutta, S., Srianand, R., & Charlton, J. 2024, MNRAS, 527, 3858

Naiman, J. P., Pillepich, A., Springel, V., et al. 2018, MNRAS, 477, 1206

Nandra, K., Barret, D., Barcons, X., et al. 2013, ArXiv e-prints [arXiv:1306.2307]

Narayanan, A. 2007, Am. Astron. Soc. Meeting Abstr., 211, 72.04

Nelson, D., Pillepich, A., Springel, V., et al. 2018, MNRAS, 475, 624

Nelson, D., Pillepich, A., Springel, V., et al. 2019, MNRAS, 490, 3234

Nelson, D., Sharma, P., Pillepich, A., et al. 2020, MNRAS, 498, 2391

Nelson, D., Pillepich, A., Ayromlou, M., et al. 2023, A&A, submitted [arXiv:2311.06338]

Nestor, D. B., Turnshek, D. A., & Rao, S. M. 2005, ApJ, 628, 637

Nielsen, N. M., Kacprzak, G. G., Pointon, S. K., et al. 2020, ApJ, 904, 164

Noterdaeme, P., Petitjean, P., Carithers, W. C., et al. 2012, ViZieR Online Data Catalog, J/A+A/547/L1

Parks, D., Prochaska, J. X., Dong, S., & Cai, Z. 2018, MNRAS, 476, 1151

Péroux, C., & Howk, J. C. 2020, ARA&A, 58, 363

Peroux, C., Merloni, A., Liske, J., et al. 2023, The Messenger, 190, 42

Peterson, J. R., Paerels, F. B. S., Kaastra, J. S., et al. 2001, A&A, 365, L104

Pillepich, A., Nelson, D., Hernquist, L., et al. 2018, MNRAS, 475, 648

Pillepich, A., Nelson, D., Springel, V., et al. 2019, MNRAS, 490, 3196

Planck Collaboration VI. 2020, A&A, 641, A6

Popesso, P., Biviano, A., Finoguenov, A., et al. 2015, A&A, 574, A105

Predehl, P., Andritschke, R., Arefiev, V., et al. 2021, A&A, 647, A1

Prochaska, J. X., Tejos, N., Crighton, N., et al. 2017, https://doi.org/10.5281/zenodo.1036773

Ramesh, R., & Nelson, D. 2024, MNRAS, 528, 3320

Reiprich, T. H., Basu, K., Ettori, S., et al. 2013, Space Sci. Rev., 177, 195

Revaz, Y., Combes, F., & Salomé, P. 2008, A&A, 477, L33

Rosati, P., Borgani, S., & Norman, C. 2002, ARA&A, 40, 539

Rykoff, E. S., Rozo, E., Hollowood, D., et al. 2016, ApJS, 224, 1

Saeedzadeh, V., Jung, S. L., Rennehan, D., et al. 2023, ArXiv e-prints [arXiv:2304.03798]

Sarazin, C. L. 1986, Rev. Mod. Phys., 58, 1

Schaye, J., Aguirre, A., Kim, T.-S., et al. 2003, ApJ, 596, 768

Schroetter, I., Bouché, N., Wendt, M., et al. 2016, ApJ, 833, 39

Schroetter, I., Bouché, N. F., Zabl, J., et al. 2019, MNRAS, 490, 4368

Seyffert, E. N., Cooksey, K. L., Simcoe, R. A., et al. 2013, ApJ, 779, 161

Sharma, P., McCourt, M., Quataert, E., & Parrish, I. J. 2012, MNRAS, 420, 3174

Simionescu, A., Werner, N., Böhringer, H., et al. 2009, A&A, 493, 409

Simionescu, A., ZuHone, J., Zhuravleva, I., et al. 2019, Space Sci. Rev., 215, 24

Slee, S. A., Gunn, J. E., Uomoto, A., et al. 2013, AJ, 146, 32

Springel, V., Pakmor, R., Pillepich, A., et al. 2018, MNRAS, 475, 676

Szakacs, R., Péroux, C., Zwaan, M., et al. 2021, MNRAS, 505, 4746

Szakacs, R., Péroux, C., Nelson, D., et al. 2023, MNRAS, 526, 3744

Theuns, T., Leonard, A., Efstathiou, G., Pearce, F. R., & Thomas, P. A. 1998, MNRAS, 301, 478

Tremblay, G. R., Oonk, J. B. R., Combes, F., et al. 2016, Nature, 534, 218

Tremmel, M., Quinn, T. R., Ricarte, A., et al. 2019, MNRAS, 483, 3336

Tripp, T. M., Lu, L., & Savage, B. D. 1998, ApJ, 508, 200

Tumlinson, J., Peebles, M. S., & Werk, J. K. 2017, ARA&A, 55, 389

van de Voort, F., Springel, V., Mandelker, N., van den Bosch, F. C., & Pakmor, R. 2019, MNRAS, 482, L85

Weinberger, R., Springel, V., Pakmor, R., et al. 2018, MNRAS, 479, 4056

Wen, Z. L., Han, J. L., & Liu, F. S. 2012, ApJS, 199, 34

XRISM Science Team 2022, ArXiv e-prints [arXiv:2202.05399]

Yoon, J. H., & Putman, M. E. 2013, ApJ, 772, L29

Yun, K., Pillepich, A., Zinger, E., et al. 2019, MNRAS, 483, 1042

Zhu, G., Ménard, B., Bizyaev, D., et al. 2014, MNRAS, 439, 3139

Zhu, G. B., Comparat, J., Kneib, J.-P., et al. 2015, ApJ, 815, 48

Zou, H., Gao, J., Xu, X., et al. 2021, ApJS, 253, 56
Diemer, B. 2018, ApJS, 239, 35

Donahue, M., & Voit, G. M. 2022, Phys. Rep., 973, 1

Dutta, R., Fumagalli, M., Fossati, M., et al. 2020, MNRAS, 499, 5022

Eke, V. R., Baugh, C. M., Cole, S., et al. 2004, MNRAS, 348, 866

Ettori, S. 2003, MNRAS, 344, L13

Ferland, G. J., Chatzikos, M., Guzmán, F., et al. 2017, Rev. Mex. A&A, 53, 385

Finoguenov, A., Rykoff, E., Clerc, N., et al. 2020, A&A, 638, A114

Frebel, A. 2018, Ann. Rev. Nucl. Part. Sci., 68, 237

Fresco, A. Y., Péroux, C., Merloni, A., Hamanowicz, A., & Szakacs, R. 2020, MNRAS, 499, 5230

Fukugita, M., & Peebles, P. J. E. 2004, ApJ, 616, 643

Gad, A. F. 2021, ArXiv e-prints [arXiv:2106.06158]

Girardi, M., Biviano, A., Giuricin, G., Mardirossian, F., & Mezzetti, M. 1993, ApJ, 404, 38

Gladders, M. D., & Yee, H. K. C. 2005, ApJS, 157, 1

Gonzalez, A. H., Zaritsky, D., & Zabludoff, A. I. 2007, ApJ, 666, 147

Gonzalez, A. H., Sivanandam, S., Zabludoff, A. I., & Zaritsky, D. 2013, ApJ, 778, 14

Guo, H., Shen, Y., & Wang, S. 2018, Astrophysics Source Code Library[record ascl:1809.008]

Hamanowicz, A., Péroux, C., Zwaan, M. A., et al. 2020, MNRAS, 492, 2347

Hummels, C. B., Smith, B. D., & Silvia, D. W. 2017, ApJ, 847, 59

Ider Chitham, J., Comparat, J., Finoguenov, A., et al. 2020, MNRAS, 499, 4768

Jones, C., & Forman, W. 1999, ApJ, 511, 65

Kirkpatrick, C. C., & McNamara, B. R. 2015, MNRAS, 452, 4361

Kirkpatrick, C. C., Clerc, N., Finoguenov, A., et al. 2021, MNRAS, 503, 5763

Kraljic, K., Laigle, C., Pichon, C., et al. 2022, MNRAS, 514, 1359

Kravtsov, A. V., & Borgani, S. 2012, ARA&A, 50, 353

Kravtsov, A. V., Nagai, D., & Vikhlinin, A. A. 2005, ApJ, 625, 588

Kunz, M. W., Jones, T. W., & Zhuravleva, I. 2022, in Plasma Physics of the Intracluster Medium, eds. C. Bambi, A. Santangelo, Handbook of X-ray and Gamma-ray Astrophysics (Springer, Singapore)

Laganá, T. F., Martinet, N., Durret, F., et al. 2013, A&A, 555, A66

Lan, T.-W., & Mo, H. 2018, ApJ, 866, 36

Lanzetta, K. M., & Bowen, D. 1990, ApJ, 357, 321

Lee, J. C., Hwang, H. S., & Song, H. 2021, MNRAS, 503, 4309

Liu, A., Bulbul, E., Ghirardini, V., et al. 2022, A&A, 661, A2

Locatelli, N., Ponti, G., Zheng, X., et al. 2024, A&A, 681, A78

Lopez, S., Barrientos, L. F., Lira, P., et al. 2008, ApJ, 679, 1144

Lyke, B. W., Higley, A. N., McLane, J. N., et al. 2020, ApJS, 250, 8

Marinacci, F., Vogelsberger, M., Pakmor, R., et al. 2018, MNRAS, 480, 5113

Mas-Ribas, L., Miralda-Escudé, J., Pérez-Ràfols, I., et al. 2017, ApJ, 846, 4

Mathur, S., Gupta, A., Das, S., Krongold, Y., & Nicastro, F. 2021, ApJ, 908, 69

McCourt, M., Sharma, P., Quataert, E., & Parrish, I. J. 2012, MNRAS, 419, 3319

Merloni, A. 2012, Science from the Next Generation Imaging and Spectroscopic Surveys, 43

Merloni, A., Alexander, D. A., Banerji, M., et al. 2019, The Messenger, 175, 42

Mernier, F., Biffi, V., Yamaguchi, H., et al. 2018, Space Sci. Rev., 214, 129

Mishra, S., & Muzahid, S. 2022, ApJ, 933, 229

Mishra, S., Muzahid, S., Dutta, S., Srianand, R., & Charlton, J. 2024, MNRAS, 527, 3858

Naiman, J. P., Pillepich, A., Springel, V., et al. 2018, MNRAS, 477, 1206

Nandra, K., Barret, D., Barcons, X., et al. 2013, ArXiv e-prints [arXiv:1306.2307]

Narayanan, A. 2007, Am. Astron. Soc. Meeting Abstr., 211, 72.04

Nelson, D., Pillepich, A., Springel, V., et al. 2018, MNRAS, 475, 624

Nelson, D., Pillepich, A., Springel, V., et al. 2019, MNRAS, 490, 3234

Nelson, D., Sharma, P., Pillepich, A., et al. 2020, MNRAS, 498, 2391

Nelson, D., Pillepich, A., Ayromlou, M., et al. 2023, A&A, submitted [arXiv:2311.06338]

Nestor, D. B., Turnshek, D. A., & Rao, S. M. 2005, ApJ, 628, 637

Nielsen, N. M., Kacprzak, G. G., Pointon, S. K., et al. 2020, ApJ, 904, 164

Noterdaeme, P., Petitjean, P., Carithers, W. C., et al. 2012, ViZieR Online Data Catalog, J/A+A/547/L1

Parks, D., Prochaska, J. X., Dong, S., & Cai, Z. 2018, MNRAS, 476, 1151

Péroux, C., & Howk, J. C. 2020, ARA&A, 58, 363

Peroux, C., Merloni, A., Liske, J., et al. 2023, The Messenger, 190, 42

Peterson, J. R., Paerels, F. B. S., Kaastra, J. S., et al. 2001, A&A, 365, L104

Pillepich, A., Nelson, D., Hernquist, L., et al. 2018, MNRAS, 475, 648

Pillepich, A., Nelson, D., Springel, V., et al. 2019, MNRAS, 490, 3196

Planck Collaboration VI. 2020, A&A, 641, A6

Popesso, P., Biviano, A., Finoguenov, A., et al. 2015, A&A, 574, A105

Predehl, P., Andritschke, R., Arefiev, V., et al. 2021, A&A, 647, A1

Prochaska, J. X., Tejos, N., Crighton, N., et al. 2017, https://doi.org/10.5281/zenodo.1036773

Ramesh, R., & Nelson, D. 2024, MNRAS, 528, 3320

Reiprich, T. H., Basu, K., Ettori, S., et al. 2013, Space Sci. Rev., 177, 195

Revaz, Y., Combes, F., & Salomé, P. 2008, A&A, 477, L33

Rosati, P., Borgani, S., & Norman, C. 2002, ARA&A, 40, 539

Rykoff, E. S., Rozo, E., Hollowood, D., et al. 2016, ApJS, 224, 1

Saeedzadeh, V., Jung, S. L., Rennehan, D., et al. 2023, ArXiv e-prints [arXiv:2304.03798]

Sarazin, C. L. 1986, Rev. Mod. Phys., 58, 1

Schaye, J., Aguirre, A., Kim, T.-S., et al. 2003, ApJ, 596, 768

Schroetter, I., Bouché, N., Wendt, M., et al. 2016, ApJ, 833, 39

Schroetter, I., Bouché, N. F., Zabl, J., et al. 2019, MNRAS, 490, 4368

Seyffert, E. N., Cooksey, K. L., Simcoe, R. A., et al. 2013, ApJ, 779, 161

Sharma, P., McCourt, M., Quataert, E., & Parrish, I. J. 2012, MNRAS, 420, 3174

Simionescu, A., Werner, N., Böhringer, H., et al. 2009, A&A, 493, 409

Simionescu, A., ZuHone, J., Zhuravleva, I., et al. 2019, Space Sci. Rev., 215, 24

Smee, S. A., Gunn, J. E., Uomoto, A., et al. 2013, AJ, 146, 32

Springel, V., Pakmor, R., Pillepich, A., et al. 2018, MNRAS, 475, 676

Szakacs, R., Péroux, C., Zwaan, M., et al. 2021, MNRAS, 505, 4746

Szakacs, R., Péroux, C., Nelson, D., et al. 2023, MNRAS, 526, 3744

Theuns, T., Leonard, A., Efstathiou, G., Pearce, F. R., & Thomas, P. A. 1998, MNRAS, 301, 478

Tremblay, G. R., Oonk, J. B. R., Combes, F., et al. 2016, Nature, 534, 218

Tremmel, M., Quinn, T. R., Ricarte, A., et al. 2019, MNRAS, 483, 3336

Tripp, T. M., Lu, L., & Savage, B. D. 1998, ApJ, 508, 200

Tumlinson, J., Peebles, M. S., & Werk, J. K. 2017, ARA&A, 55, 389

van de Voort, F., Springel, V., Mandelker, N., van den Bosch, F. C., & Pakmor, R. 2019, MNRAS, 482, L85

Weinberger, R., Springel, V., Pakmor, R., et al. 2018, MNRAS, 479, 4056

Wen, Z. L., Han, J. L., & Liu, F. S. 2012, ApJS, 199, 34

XRISM Science Team 2022, ArXiv e-prints [arXiv:2202.05399]

Yoon, J. H., & Putman, M. E. 2013, ApJ, 772, L29

Yun, K., Pillepich, A., Zinger, E., et al. 2019, MNRAS, 483, 1042

Zhu, G., Ménard, B., Bizyaev, D., et al. 2014, MNRAS, 439, 3139

Zhu, G. B., Comparat, J., Kneib, J.-P., et al. 2015, ApJ, 815, 48

Zou, H., Gao, J., Xu, X., et al. 2021, ApJS, 253, 56