© The Authors 2024

1 Introduction

In the current understanding of the “bottom–up” structure formation of the Universe in the standard ΛCDM cosmology, small overdensities collapse first, overcoming the cosmological expansion and merging to form larger halos (Springel 2005). In this scenario, the gas encapsulated in dark matter halos forms the first stars and galaxies as it cools and condenses. The effects of tidal forces, mergers, and interactions in their surroundings regulate the galaxy formation and evolution process. The majority of galaxies in the Universe are found in dense environments as galaxy groups and include a large fraction of the universal baryon budget (Mulchaey 2000). The interaction between galaxies and the intra-group medium (IGrM; the gas encapsulated within the galaxy groups) plays a crucial role in shaping the properties and evolution of galaxies. For instance, the feedback from supernovae, star formation, and central supermassive black holes impacts the evolution of galaxies within the galaxy groups.

Although there is no clear definition of galaxy groups, dark matter halos with fewer than 50 galaxies and/or masses between 5 × 10¹²–10¹⁴ M_⊙ are classified as galaxy groups (Crain et al. 2009) in the literature. In addition to their member galaxies, galaxy groups contain diffuse baryonic matter in the form of hot plasma with temperatures (T¹) ranging from 0.1–2 keV, all encapsulated by the potential well provided by dark matter. Groups are further categorized as loose (or poor) groups, compact groups, and fossil groups, depending on their optical properties (Hickson 1997; Mulchaey & Zabludoff 1998; Mulchaey 2000; Voevodkin et al. 2010).

Despite the abundance of galaxy groups and the essential role they play in the assembly process of dark matter halos in the Universe (Crain et al. 2009), their detection has been challenging due to their low richness, faint X-ray signal, and shallow potential wells. A variety of methods have been employed to search for galaxy groups. In the optical domain, clustering and friends-of-friends (FoF) algorithms have been used to catalog groups in spectroscopic or photometric galaxy redshift surveys (e.g., Hickson 1982; Robotham et al. 2011; Tempel et al. 2017; Gozaliasl et al. 2022). However, due to their relatively low richnesses, group catalogs compiled using optical observations may suffer from large contamination fractions from the random superposition of galaxies along the line-of-sight, otherwise known as projection effects (see Costanzi et al. 2019; Grandis et al. 2021; Myles et al. 2021). On the other hand, in the X-ray domain, the emission from IGrM makes them appear as extended sources in the X-ray sky, where they typically exhibit rapidly increasing X-ray emission profiles from the outskirts of the system to the center. Because of their characteristic surface brightness profiles, if their emission is above the background level, they can be easily identified and do not suffer significantly from projection effects. Aside from being a reliable tool for detecting groups, X-ray observations also enable the measurement of the physical properties of the hot ionized IGrM through imaging and spectral analysis.

The effects of the non-thermal astrophysical phenomena governing galaxy formation are easier to study using galaxy groups compared to clusters, as the input energy associated with these phenomena is comparable to the binding energies of groups (Balogh et al. 2001). For investigating the non-thermal phenomena, entropy ($S=T/n_{\text{e}}^{2/3}$) measurements of IGrM are often used, which retain a historical record of the thermodynamic state of the gas and reflect the changes in the cooling and heating processes in galaxy groups (Voit et al. 2005; Pratt et al. 2010). For example, OverWhelmingly Large Simulations (OWLS, Schaye et al. 2010; McCarthy et al. 2010) showed that the outflows from the central active galactic nuclei (AGN) elevate the entropy of the IGrM by mechanically expelling low-entropy gas instead of directly heating it. This process mitigates rapid cooling and prevents excessive star formation (Bryan 2000; Balogh et al. 2001). On the other hand, the feedback generated by supernova-driven winds associated with galaxies can also increase the entropy of the IGrM through direct heating (see Eckert et al. 2021, for a review).

Another important observational constraint on the nonthermal astrophysical phenomenon is the shape of the stellar mass function. Observational studies have demonstrated that the stellar mass exhibits a cutoff at M_* ~ 10¹¹ M_⊙ (Davidzon et al. 2017). Early studies trying to reproduce the cutoff with stellar feedback (e.g., through the energy and momentum released by supernovae explosions and stellar winds) were unsuccessful, as the injected energy was proven to not be enough to prevent cooling and regulate the star formation efficiency (Benson et al. 2003). Therefore, it is commonly agreed that feedback from another source, such as AGN, is needed to reproduce the observed cutoff in the stellar mass function (Harrison 2017). Moreover, there are further observational constraints on the energy released by the non-gravitational feedback, such as the Si and Fe abundance profiles of the intracluster and intragroup mediums. These measurements cannot be reproduced even with the assumption of 100 percent efficient stellar feedback heating (Kravtsov & Yepes 2000), a scenario that is rejected by the measurements of the galactic outflows (Martin 1999). Therefore, the total amount of energy that can be injected through stellar feedback on the IGrM is constrained relatively well by the abundant observational data. Consequently, by measuring the thermodynamic properties of galaxy groups, one effectively constrains the energetics of the remaining source of energy, AGN. For the higher mass groups (M_500c > 10^13.5 M_⊙), the impact of stellar feedback is at a negligible level such that entropy measurements put direct constraints on the impact of AGN (Le Brun et al. 2014). For the low-mass groups (M_500c < 10¹³ M_⊙), constraints from multiple observables should be combined to isolate the impact of AGN on its surroundings (e.g., see Altamura et al. 2023).

The AGN heating in galaxy clusters and groups is observationally confirmed by shocks, ripples, and cavities detected in X-ray wavelengths (e.g., Fabian et al. 2006; Randall et al. 2011) as well as the detection of radio-loud AGN in a significant proportion of the brightest cluster and group galaxies of the cool core galaxy clusters and groups (e.g., Burns 1990; Best et al. 2007; Smolcic et al. 2011). Furthermore, deeper radio observations have revealed that nearly every central galaxy in X-ray bright groups hosts radio emission (Kolokythas et al. 2019). In fact, radio observations of galaxy groups are highly complementary to the X-ray view of groups for investigating the impact of AGN on IGrM (Eckert et al. 2021). Simultaneously studying their X-ray and radio properties allows for putting constraints on the radio mode feedback from the central engine (e.g., Pasini et al. 2022; Böckmann et al. 2023). Nevertheless, combining multi-wavelength datasets comes with challenges. For instance, crossmatching X-ray and radio catalogs makes it challenging to have a good handle on the selection effects, which is crucial for achieving unbiased conclusions about galaxy groups at the population level. Given the challenges and caveats, in this work, we only focus on putting constraints on the impact of non- gravitational feedback mechanisms through X-ray observations and leave the investigation of the multi-wavelength picture of the eRASS1 galaxy groups sample to future work.

Entropy of IGrM can be measured using X-ray observations, where the electron density and temperature measurements can be made using the imaging and spectroscopic capabilities of X-ray telescopes, such as SRG/eROSITA, XMM-Newton, and Chandra. Ponman et al. (1999) measured the entropy of 25 bright galaxy clusters and groups at a radius of 0.1 r_virial² using ROSAT and GINGA observations and reported that the entropy measurements at the core lie above the expected power-law relation with temperature for the first time. Subsequently, Lloyd-Davies et al. (2000), Finoguenov et al. (2002), and Ponman et al. (2003) measured the entropy profiles of galaxy clusters and groups using ROSAT and ASCA observations that provided the first hint that in galaxy groups, the excess entropy is not limited only to the core but can also be prominent at larger radii. Voit et al. (2005) formulated a baseline entropy profile that can be used for evaluating the impact of non-gravitational processes for galaxy clusters and groups using four sets of simulations that only include gravitational processes. Using X-ray instruments with a higher spatial and spectral resolution, such as XMM-Newton and Chandra, significantly improved our understanding of the excess entropy in galaxy groups by accurately measuring their entropy profiles. Johnson et al. (2009) investigated entropy profiles of galaxy groups by analyzing XMM-Newton observations of 28 nearby galaxy groups from the Two-Dimensional XMM-Newton Group Survey. They divided their sample into two subsamples (cool core and non-cool core) based on the temperature gradient at the core of their groups and found that the entropy profiles of the groups in their non-cool core sample exhibit less scatter compared to entropy profiles of their cool core sample. Around the same time, Sun et al. (2009, S09 hereafter) conducted a comprehensive study on the thermodynamic gas properties of 43 nearby galaxy groups using the archival Chandra observations, where they constrained the temperature, electron density, and metallic- ity profiles of 23 groups out to r_500c³ accurately thanks to the outstanding imaging capabilities of Chandra and the relatively deep archival observations of some systems in their sample. Furthermore, they compared their entropy profiles with the baseline entropy profile of Voit et al. (2005) and found that even though the entropy excess reduces as a function of the radius, it remains significant out to r_500c . Subsequently, the detailed analysis of the outskirts of RX J1159+5531, UGC 03957, and Virgo using Suzaku observations revealed that the entropy excess can go beyond r_500c (Humphrey et al. 2012; Thölken et al. 2016; Simionescu et al. 2017). More recently, Panagoulia et al. (2014) analyzed 66 galaxy groups from the NORAS and REFLEX samples to investigate the properties of IGrM at the core and found that entropy profiles of galaxy groups at the core follow a power-law relation and do not exhibit any entropy floor.

Previous studies in the literature on the thermodynamic properties of galaxy groups have been conducted using relatively small (fewer than 100) and highly incomplete samples that lack well-defined selection functions. Notably, eROSITA opens a new window for galaxy group studies by providing the largest pure X-ray selected sample with a well-defined selection function, which is crucial for achieving robust conclusions that reflect the physical properties of the galaxy group well at the population level. Furthermore, the superb soft X-ray band sensitivity and the scanning observing strategy of eROSITA make it an excellent instrument for investigating the physical properties of the hot gas in galaxy groups, as the emission of IGrM peaks at the soft X-ray band and the brightest galaxy groups above the detection capabilities of the current instruments are at low redshift and well extended.

In this work, we examine the effect of the feedback on the thermodynamics of galaxy groups detected by eROSITA in its first All-Sky Survey. We accomplish this by performing joint imaging and spectral analysis on the eRASS:4 (the four consecutive eROSITA All-Sky Surveys stacked together) observations of the galaxy groups in our sample. The extended ROentgen Survey with an Imaging Telescope Array (eROSITA), the soft X- ray telescope on board the Spectrum-Roentgen-Gamma (SRG) mission (Sunyaev et al. 2021), was launched on July 13, 2019 (Predehl et al. 2021). The first All-Sky Survey with eROSITA was successfully executed on June 11, 2020, after 184 days of operation. In this first All-Sky Survey (eRASS1), eROSITA detected a total of 12247 optically confirmed galaxy groups and clusters spanning the redshift range 0.003 < z < 1.32 with a sample purity level of 86% in the Western Galactic half of the survey (359.9442 deg > l > 179.9442 deg), where the data rights belong to the German eROSITA consortium (Merloni et al. 2024; Bulbul et al. 2024; Kluge et al. 2024).

In this paper, we combine the imaging and spectroscopic information of 1178 eROSITA-detected galaxy groups and obtain average entropy measurements at 0.15r_500c, r_2500c and r_500c to investigate the effects of AGN feedback and compare our findings with the state-of-the-art numerical simulations from MillenniumTNG (Hernández-Aguayo et al. 2023; Pakmor et al. 2023), Magneticum (Hirschmann et al. 2014), and OverWhelmingly Large Simulations (Schaye et al. 2010; McCarthy et al. 2010). The findings represent the first study of a comprehensive group sample with a well-defined selection function. This paper is organized as follows: In Sect. 2, we describe the construction of the galaxy groups sample from the primary eRASS1 sample, and in Sect. 3, we describe the X-ray data reduction and the analysis of the groups. In Sect. 4, we provide a discussion of major systematics and the details of the quantification and incorporation of them in our results, and in Sect. 5, we provide our final results on the entropy measurements of the sample and comparisons with the previous measurements. In Sect. 6, we provide a comparison between our measurements and the predictions of the state-of-the-art simulations. Lastly, we provide a summary of our findings and list our conclusions in Sect. 7. Throughout this paper, we adopt a flat ΛCDM cosmology using the Planck Collaboration XIII (2016) results, namely Ω_m = 0.3089, Ω_b = 0.0486, σ₈ = 0.8147, and H₀ = 67.74 km s⁻¹ Mpc⁻¹. Quoted error bars correspond to a 1-σ confidence level unless noted otherwise.

2 Sample of galaxy groups

This work utilizes a subsample of the X-ray selected, optically identified primary eRASS1 galaxy cluster and group sample detected in the Western Galactic hemisphere of the first eROSITA All-Sky Survey (Bulbul et al. 2024; Kluge et al. 2024). Below, we briefly describe the detection of galaxy clusters and groups in eRASS1 observations and a brief summary of the optical and X-ray cleaning performed in Bulbul et al. (2024). Subsequently, we provide the details of the additional selection and cleaning applied to the eRASS1 galaxy clusters and groups catalog.

The X-ray emitting celestial objects in the eRASS1 master X- ray catalog (Merloni et al. 2024) are detected using the eROSITA source detection pipeline, which is part of the eROSITA Science Analysis Software System (eSASS, Brunner et al. 2022). The pipeline locates detection candidates and calculates detection and extent likelihood (ℒ_det and ℒ_ext ) parameters by comparing the spatial distribution and the abundance of the photons around the candidate with the local background. To construct the primary galaxy groups and clusters sample (Bulbul et al. 2024), a ℒ_ext > 3 cut is applied to increase the completeness of the galaxy groups and clusters sample (see Bulbul et al. 2022, for the motivation). The DESI Legacy Survey DR9 and DR10 datasets are used in the optical identification processes by the eROMaPPer pipeline, which is based on the matched-filter red-sequence algorithm from redMaPPer (Rykoff et al. 2014, 2016) tailored and optimized for the identification of eROSITA extended sources (Ider Chitham et al. 2020; Kluge et al. 2024). If available, spectroscopic redshifts (z_spec) are prioritized over photometric redshift (z_λ) by the eROMaPPer pipeline (see Kluge et al. 2024, for further details). In this primary sample, 12 705 extended sources in the redshift range of 0.01 to 1.35 are identified as galaxy clusters or groups with a contamination fraction of 14% (Bulbul et al. 2024; Kluge et al. 2024).

To construct a final clean and secure galaxy group sample, we apply further cuts based on the X-ray and optical properties of the primary sample. While the literature lacks a precise definition for galaxy groups, we classify an object as a group if its mass ranges between 5 × 10¹² < M_500c < 10¹⁴ M_⊙. The upper end of our group definition (10¹⁴M_⊙) corresponds to a plasma temperature of T ~ 2 keV and is commonly used in previous X-ray studies for distinguishing galaxy clusters from galaxy groups (Lovisari et al. 2021). For incorporating this mass criterion, we use the M_500c estimates obtained in Sect. 3.3 using a Bayesian X-ray observable estimation framework that jointly estimates the soft-band (0.5–2 keV) X-ray luminosity (L_X), temperature (T) and the mass (M_500c) of galaxy clusters and groups from their count-rate profiles (see Sect. 3.3 for the details of our L_X − T − M_500c estimation). After applying a mass cut of M_500c < 10¹⁴M_⊙, we select 2526 galaxy group candidates with a median redshift of 0.11.

To further reduce the contamination, we apply other cleaning methods using the deeper eRASS:4 data. Contaminants in the ℒ_ext > 3 sample of the eRASSl clusters and groups catalog can be classified into two categories: misclassified sources (mostly AGN) and spurious sources. Given that our preliminary sample has a median redshift of 0.11, the “real” galaxy groups in our sample are expected to be relatively well extended in the sky, whereas misclassified point sources, by definition, should have a low extent. We make use of this fact and conservatively remove 841 objects that have EXT < 20 arcsec and ℒ_ext < 5.5. These cuts remove most point sources, leaving 1685 group candidates in the sample.

Once the misclassified point sources are removed, spurious sources are left to be cleaned from our group sample. We use the count measurements in the 0.3–1.8 keV band (see Appendix A for the details on the choice of the energy band) obtained from eRASS:4 observations as described in Sect. 3.2 to clean the spurious sources. We first apply an X-ray count cut of 10. This cut removes 423 objects from our group candidates. Furthermore, we remove 10 more sources from the remaining sample with count measurements 1σ consistent with the background level. This procedure removes most of the spurious sources since one would expect the galaxy groups to be more prominent and bright as the survey gets deeper. On the other hand, the spurious sources are expected to have low counts and be consistent with the background level since they are mostly due to background fluctuations or superpositions of undetected AGN in the eRASS1 observations. Applying these cuts, we remove a large fraction of contaminants in our sample and obtain a highly pure sample with 1252 galaxy groups.

One of the major benefits of the strict cleaning procedure described above is its applicability to simulations. Our cleaning procedure relies on the detection pipeline outputs (e.g., ℒ_ext and EXT), and therefore the selection process is fully reproducible in the simulations of the eRASS1 digital twin (Comparat et al. 2020; Seppi et al. 2022). This allows us to construct a robust selection function for our sample using the eROSITA’s digital twin simulations. We further note that the cleaning applied in this work to remove spurious sources has no impact on the selection function⁴.

Following the cleaning procedure, we visually inspect the groups that are centrally peaked (r_c < 32 arcsec)⁵ and have low extent likelihoods (ℒ_ext < 10) using the eRASS:4 and the Legacy Survey data. As a result of the visual inspection, we further flagged and removed 74 falsely classified point sources from our sample. These objects make a small fraction of our clean sample (6%); therefore, their removal has a negligible impact on the selection function, especially compared to the systematic uncertainties of X-ray simulations at group scales used to construct the selection function. The final sample, consisting of 1178 galaxy groups, has a median redshift of 0.11 and a median mass (M_500c) of 6.3 × 10¹³ M_⊙. The Redshift and mass distributions of the final sample can be seen in the right panel of Fig. 1. Moreover, the 2D projected distribution of these groups in the eROSITA sky is shown in Fig. 2. We note that some of the “cleaning” procedures described above (e.g., the EXT and ℒ_ext cuts) not only clean spurious sources but also unavoidably remove some of the faint real groups from the sample according to the expected purity of the eRASS1 cluster sample in the cost of achieving a more secure groups sample. Nevertheless, the resulting extra selection is taken into account in our analysis by incorporating a selection function built for our final sample. After the cleaning procedure described above, our sample ended up having three objects with mass estimates slightly below the lower bound of our group mass definition (5 × 10¹² M_⊙ < M _500c < 10¹⁴ M_⊙). We eventually decided to keep them in our group sample since the removal of three objects has little to no impact on our final results, and their “true” masses can well be within our group mass definition due to the intrinsic scatter of the L_X − M_500c relation.

The final galaxy group sample described above is obtained to construct a well-defined selection function using eROSITA’s twin simulations. A good handle on selection effects is key for achieving universal conclusions about the properties and the governing physics of studies of groups. The deeper eRASS:4 observations of an unprecedented number of galaxy groups we use in this work are particularly well-suited for studying the baryonic physics in galaxy groups because of the higher statistics allowed by the deeper survey data and large field of view necessary to measure the X-ray properties out to large radii. In the next section, we present our eRASS:4 analysis of the galaxy groups in the sample.

3 Data analysis

3.1 X-ray data reduction and analysis

Taking advantage of the higher signal-to-noise, deeper survey observation, we used the eRASS:4 observations of the eRASS1 selected galaxy groups with the processing version 020 (briefly described in Appendix C of Merloni et al. 2024), which is an updated version of 010 processing used for the first data release (DR1). The main updates on the 020 version (internal catalog version 221031) are the improved boresight correction, low- energy detector noise suppression, improved subpixel resolution, and updated pattern and energy tasks. We further reduce the calibrate event files using the using the eROSITA Science Analysis Software System (eSASS, Brunner et al. 2022)⁶ with the version id eSASSusers_211214 that is the same version used for Bulbul et al. (2024) and Merloni et al. (2024) for DR1. Time variable (solar incident angle dependent) optical light contamination (light leak) is observed in the data from the telescope modules (TMs) 5 and 7, which has a large impact on the calibration of the low-energy band of the spectrum (see Predehl et al. 2021; Coutinho et al. 2022; Merloni et al. 2024, for further details). In this work, we analyze the hot gas properties of galaxy groups whose emission peaks at the soft X-ray band that suffers from the contamination; therefore, we followed a conservative approach and only used the data from TMs 1, 2, 3, 4, and 6, removing the data from TMs 5 and 7 that suffer from contamination due the optical light leak. Furthermore, we obtain a clean event list by applying the standard flag 0xE000F000 to select all the possible patterns (singles, doubles, triples, and quads) and run the flaregti eSASS task to have flare filtered good-time-intervals.

Fig. 1 Left: mass and redshift distributions of the galaxy group sample used in this work consisting of 1178 objects where colors of the data points represent the eROSITA counts of the groups in soft X-ray band (0.3-1.8 keV) Right: mass and redshift histograms of the group sample where the median redshift (z) is 0.11 and median mass (M500c) is 6.5 × 1013 M⊙.

Fig. 2 Projected locations of the 1178 galaxy groups in the primary catalog in the eROSITA and Legacy Survey DR9N and DR10 13,116 deg2 common footprints. The redshift confirmed by the follow-up algorithm eROMaPPer is color coded (Kluge et al. 2024), while the sizes of the detections are scaled with the angular sizes (r500c) of the groups (see Sect. 3.3 for the r500c estimation procedure). The inhomogeneity of the source density in this figure is due to the exposure variation across the eROSITA-DE X-ray sky (see Fig. 2 in Bulbul et al. 2024).

3.2 Imaging analysis

For the imaging analysis, we use an energy band of 0.3–1.8 keV to maximize the signal-to-noise ratio for a soft X-ray emitting source such as groups (see Appendix A for the details of this optimization scheme). We extract images, vignetted and nonvignetted exposure maps in this band centered around each group in the catalog using the evtool and expmap tasks in eSASS with a standard eROSITA pixel size of 4 arcsec and the FLAREGTI option. For the extraction region, we used an image size of ~8r_500C,eSASS × 8r_500C,eSASS that covers well the region from the source center beyond the Virial radius for the local background measurements. The radius, r_500c,_eSASS, was estimated using the flux reported in column ML_FLUX_1 of the eRASS1 X-ray catalog (Merloni et al. 2024) and an L_X − M_500c relation of the eROSITA Depth Final Equatorial Survey (eFEDS) clusters and groups (Bahar et al. 2022; Chiu et al. 2022). These estimates were only used to determine the image size that has a negligible impact on the results. After generating X-ray images and exposure maps, we used the eRASS:4 point source catalog in the 0.2–2.3 keV band to mask or co-fit the point sources in the field of view in the rest of the imaging analysis. Following the same procedure in our eFEDS analysis (Ghirardini et al. 2021; Liu et al. 2022; Bahar et al. 2022), we masked the faint point sources with ML_RATE_l < 0.1 cnts/s) out to the radii where their emission becomes consistent with the background. On the other hand, we co-fit bright point sources (ML_RATE_l > 0.1 cnts/s) in the surface brightness analysis. In addition to the bright point sources, we also modeled and co-fit the closest extended sources to the central galaxy group to clean the image from contaminating X-ray emission. The remaining ones in the field are conservatively masked out to their 2r_500c,eSASS. Example eROSITA images of a bright nearby group (1eRASS J024933.9- 311126) and a group at the median redshift of our group sample (1eRASS J045547.0-572404) are shown in Fig. 3.

Following a forward modeling approach, we fit the X-ray images using a Bayesian fitting pipeline to deproject the surface brightness emission. We assume a Poisson likelihood for the X-ray counts and sample the likelihood using the emcee package (Foreman-Mackey et al. 2013) that employs the Goodman & Weare (2010) Affine Invariant Markov chain Monte Carlo (MCMC) technique. The fit is performed to account for the cross-talk between the emission from the nearby co-fitted extended and point sources. For extended sources, we model emissivity using a modified Vikhlinin et al. (2006) profile: $${\text{Λ}}_{\text{ep}}(T,Z)n_{\text{e}}^2(r)=N_{\text{Λ},n_{\text{e}}^2}{\left(\frac{r}{r_c}\right)}^{-\alpha }{\left(1+{\left(\frac{r}{r_c}\right)}^2\right)}^{-3\beta +\alpha /2}{\left(1+{\left(\frac{r}{r_s}\right)}^3\right)}^{-ϵ/3},$$(1)

where Λ_ep(T, Z) is the band-averaged cooling function, Z is the metallicity, and $N_{\text{Λ},n_e^2},r_c,r_s,\alpha ,\beta ,ϵ$ are the free parameters of the emissivity profile, normalization, core radius, scale radius, and the power law exponents, respectively. These parameters are allowed to vary in the fits. At each MCMC step, the profile is projected along the line of sight following the equation $$S_X=\frac{1}{4\pi {(1+z)}^4}{\displaystyle \int n_{\text{e}}}{n}_{\text{p}}{\text{Λ}}_{\text{ep}}(T,Z)\text{d}l,$$(2)

where the number density of protons (n_p) are related to the number density of electrons (n_e) via n_p = n_e/1.2 (Bulbul et al. 2010). Then, the projected count rate is convolved with the eROSITA point spread function (PSF), multiplied with the exposure map, and compared with the masked X-ray image. In addition to the electron density profile parameters, two additional parameters are left free for the central position of the groups, which adds up to eight free parameters for every extended source modeled in the image. For the bright point sources, we only allowed the normalization of their profile to vary while keeping the centroids fixed due to the high positional accuracy of eROSITA (Brunner et al. 2022; Merloni et al. 2024). Lastly, for each image, the background count rates are assumed to be constant across the image, and two more parameters are allowed to be free for the vignetted and unvignetted backgrounds.

As an output of this fitting procedure, we obtain the best- fit de-projected emissivity profiles (${\text{Λ}}_{\text{ep}}(T,Z)n_{\text{e}}^2(r)$), count-rate profiles, and the associated uncertainties that include the crosstalk between the co-fitted nearby extended and point sources. We show the performance of our pipeline in Fig. 4 for a bright group (1eRASS J045547.0-572404, the second group in Fig. 3) with the images before and after subtracting the emission from the modeled extended and point sources on the left. It is clear from the bottom left figure that after the removal of the modeled profiles, the image is free from any X-ray source, and only noise remains. The surface brightness profile of the observed field (in blue) and the best-fit model (in red) are shown on the right panel of the same figure. The PSF convolved best-fit surface brightness model represents the eROSITA data well. The peaked emission at large radii (at 500 and 1000 arcsec) shows the contribution of the modeled point sources to the overall emission, which are successfully modeled and removed from the total source model.

Given the PSF of eROSITA being relatively large, we also ran tests on the robustness of our fitting procedure around the core region (0.15r_500C) of groups by simulating and fitting synthetic galaxy group observations. The synthetic observations are obtained by first generating galaxy group profiles in a nonparametric way employing a covariance matrix obtained from XXL observations following Comparat et al. (2020)⁷. These profiles are then convolved with the eROSITA PSF, and the X-ray observations are obtained by creating the Poisson realizations of the PSF convolved surface brightness distributions. Through this procedure, we have fitted 30 simulated groups at a redshift of z = 0.11 (the median redshift of our sample) and 30 groups at a redshift of z = 0.2 (85% of the groups in our sample are at z < 0.2). As a result of these tests, we found that our fitting procedure is capable of robustly deconvolving the profiles with PSF and recovering the input surface brightness profiles around 0.15r_500c. We also found that because of the PSF smoothing, the recovered profiles of the objects with intrinsically larger surface brightness fluctuations may deviate more from the input simulated profiles; however, at the sample level, these fluctuations cancel out such that our measurements, on average, are unbiased. Furthermore, we have also investigated the possible impact of an undetected central compact source on the surface brightness measurements of the groups in our sample at 0.15r_500c by comparing the fitted surface brightness profiles of groups with the PSF profile. From this investigation, we find that an undetected point source at the center of a group can only change n_e(0.15r_500c) a few percent, which is within the total error budget of our n_e measurements that includes statistical and systematic uncertainties (see Sect. 4 for details on the systematic uncertainties taken into account in this study). Therefore, we conclude that the galaxy groups we use in this work are well extended in the sky such that an undetected point source at the center of the group has little to no impact on the electron density measurements at 0.15r_500c.

Fig. 3 Left: eRASS:4 soft band 0.3–1.8 keV images of two bright groups (1eRASS J024933.9-311126 and 1eRASS J045547.0-572404) at redshifts 0.023 and 0.105. Right: legacy Survey DR10 images of the same groups with the eRASS:4 X-ray contours overlayed.

3.3 Estimation of X-ray observables

The shallow nature of the eROSITA survey only allows the measurement of the physical properties of a few nearby bright galaxy groups. The eRASS1 group sample should be binned into smaller samples to achieve sufficient signal-to-noise for reliably constraining the physical properties of the faint galaxy groups through joint spectral analysis (see Sect. 3.4). For an optimal binning scheme, a low-scatter temperature estimator should be used such that groups with similar temperatures can be binned together. Moreover, mass estimates of the galaxy groups are needed to extract spectra within a physically motivated scale radii (r_500c). For these purposes, we use L_X − M_500c and L_X − T scaling relations and calculate the soft-band (0.5–2 keV) X-ray luminosity (L_X), temperature (T), and mass ( M_500c) estimates of the galaxy clusters and groups from the count-rate profiles measured in Sect. 3.2. For a self-consistent treatment of the eROSITA groups, we employ the L_X − M_500c and L_X − T relations calibrated using the eFEDS observations (Chiu et al. 2022; Bahar et al. 2022).

The selection effects and the mass function need to be accounted for to obtain unbiased estimations of the physical properties of an underlying population from intrinsically scattered scaling relations. For this purpose, we built a Bayesian framework that simultaneously estimates the L_X and T, M_500c observables from the observed count-rate profiles, Ĉ_R (r) (see Sect. 3.2 for the details of the count-rate profile measurement procedure). The formulation of the Bayesian estimation framework is as follows. To simultaneously estimate L_X, T and M_500c observables, the joint probability density function, P(L_X, T, M_500c|D, θ_all), given the data, D, and a set of model parameters, θ_all is needed to be computed. This can be expanded as $$P\left(L_{\text{X}},T,M_{500c}\mid D,{\theta }_{\text{all\hspace{0.17em}}}\right)\equiv P\left(L_{\text{X}},T,M_{500c}\mid {\widehat{C}}_{R,500c},I,z,\mathscr{H},{\theta }_{\text{LT}},{\theta }_{\text{LM}}\right),$$(3)

where the count-rate within r_500c (Ĉ_R,500c), detection information of the galaxy group (I), redshift (z), and sky position (H), represent the data (D); and the scaling relations parameters (θ_LT and θ_LM) represent the model parameters (θ_all). We note that by definition, the Ĉ_R,500c term has an intrinsic dependence on M_500c such that Ĉ_R,500c is different for every M_500c in the L_X − T − M_500c parameter space. Taking this into account, our framework allows all the information of the measured count-rate profiles to be included in our analysis rather than count-rate measurements within fixed radii.

Using the Bayes rule, Eq. (3) can be rewritten as $$\begin{array}{l}P\left(L_{\text{X}},T,M_{500c}\mid {\widehat{C}}_{R,500c},I,z,\mathscr{H},{\theta }_{\text{LT}},{\theta }_{\text{LM}}\right)\\ =\frac{P\left(L_{\text{X}},T,M_{500c},{\widehat{C}}_{R,500c},I\mid z,\mathscr{H},{\theta }_{\text{LT}},{\theta }_{\text{LM}}\right)}{{\displaystyle \iiint P}\left(L_{\text{X}},T,M_{500c},{\widehat{C}}_{R,500c},I\mid z,\mathscr{H},{\theta }_{\text{LT}},{\theta }_{\text{LM}}\right)d{L}_{\text{X}}dTd{M}_{500c}}.\end{array}$$(4)

Furthermore, one can expand the common term in the numerator and the denominator as $$\begin{array}{l}P\left(L_{\text{X}},T,M_{500c},{\widehat{C}}_{R,500c},I\mid z,\mathscr{H},{\theta }_{\text{LT}},{\theta }_{\text{LM}}\right)\\ \approx P\left({\widehat{C}}_{R,500c}\mid L_{\text{X}},T,M_{500c},z\right)P\left(I\mid L_{\text{X}},z,\mathscr{H}\right)\\ \times P\left(L_{\text{X}},T\mid M_{500c},z,{\theta }_{\text{LT}},{\theta }_{\text{LM}}\right)P\left(M_{500c}\mid z\right),\end{array}$$(5)

where the first term, P(Ĉ_R,500c|L_X, T, M_500c, z), stands for the measurement uncertainty of the count-rate. The conditional probability distribution for this term is obtained by first calculating the true count rate (Ĉ_R,500c) for every point in the L_X − T parameter space by assuming the source emitting an unabsorbed APEC (Smith et al. 2001) spectrum in Xspec (Arnaud 1996) at a redshift z with an abundance of 0.3Z_⊙, a temperature of T and a luminosity of L_X. Then the true count-rates (C_R,500c) are compared with the observed count-rates calculated (Ĉ_R,500c) at every r_500c value in the mass parameter space (M_500c) and the value of the conditional probability is obtained. The P(I|L_X, z, H) term in Eq. (5) is the selection function term that is a function of soft-band X-ray luminosity, redshift, and sky position where the sky position includes the local background surface brightness, exposure, and the neutral hydrogen column density information.

The selection function is obtained by simulating the eROSITA X-ray All-Sky observations using the baryon painting method (Comparat et al. 2019, 2020) and applying the same routines of the eSASS source detection pipeline to construct one- to-one correspondence of the catalogs and selection (Seppi et al. 2022; Clerc et al. 2024). The P(M_500c|z) term is the mass function term for which the analytical formulation of Tinker et al. (2008) is used in this work. Lastly, the P(L_X, T|M_500c, z, θ_LT, θ_LM) term is the intrinsically scattered scaling relation term that gives the L_X and T distributions at a given mass and redshift.

Ideally, one would use a jointly fit, intrinsically scattered L_X − T − M_500c scaling relation for the P(L_X, T|M_500c, z, θ_LT, θ_LM); however, there is no such relation in the literature yet that is calibrated by taking into account the selection effects and covers a similar mass range with eROSITA. For this reason, we expanded this term as $$P\left(L_{\text{X}},T\mid M_{500c},z,{\theta }_{\text{LT}},{\theta }_{\text{LM}}\right)\approx P\left(T\mid L_{\text{X}},z,{\theta }_{\text{LT}}\right)P\left(L_{\text{X}}\mid M_{500c},z,{\theta }_{\text{LM}}\right)$$(6)

and used the Bahar et al. (2022) L_X − T and Chiu et al. (2022) L_X − M_500c relations that are calibrated for eROSITA by taking into account the selection effects. The P(T |L_X, z, θ_LT) term was obtained from P(L_X |T, z, θ_LT) using Bayes theorem: $$P\left(T\mid L_{\text{X}},z,{\theta }_{\text{LT}}\right)=\frac{P\left(L_{\text{X}}\mid T,z,{\theta }_{\text{LT}}\right)P(T\mid z)}{{\displaystyle \int P}\left(L_{\text{X}}\mid T,z,{\theta }_{\text{LT}}\right)P(T\mid z)d{L}_{\text{X}}},$$(7)

where self-consistently, the same temperature function is used for the P(T |z) term as in Eq. (6) in Bahar et al. (2022).

As a final step, we substitute the terms in Eqs. (5)–(4) and calculate the joint probability density function, P(L_X, T, M_500c|Ĉ_R,500c, I, z, 𝓗,θ_LT,θ_LM), for each galaxy group. Subsequently, we marginalize over the nuisance parameters and obtain L_X, T, and M_500c estimates given the data and the scaling relations parameters⁸. We provide the distributions of mass, temperature, soft-band (0.5–2 keV) X-ray luminosity estimates obtained through this Bayesian framework along with the distributions of redshift and count (0.3–1.8 keV) of the final galaxy group sample in Figs. 1, 5 and 6.

The estimated M_500c are then converted to r_500c and r_2500c by assuming an average dark matter concentration of c_500c = r_500c/r_s,NFW = 4.2⁹ and scaling the r_500c estimates accordingly (S09). These characteristic radii are employed to determine the spectral extraction region (r < r_500c) and serve as characteristic radii (0.15r_500c, r_2500c, r_500c) for the entropy measurements.

The net effect of accounting for selection effects and the mass function when estimating L_X, T, and M_500c from scaling relations depends on two factors: the scatter in the scaling relations and the interplay between the selection and mass functions. With zero scatter, there’s a one-to-one relationship between observables, allowing straightforward conversions. As the scatter of the relation increases, the scaling relation estimates that ignore the selection effects will be more vulnerable to being biased. Furthermore, the net effect also depends on the shapes of the selection and mass functions. The interplay between the selection function and the mass function across the L_X − T −M_500c parameter space is often not trivial; however, to the first order, if we consider X-ray selection as a redshift dependent L_X cut, not accounting for selection and mass functions would lead to both T and M_500c being overestimated.

Fig. 4 Example of the eROSITA imaging analysis for the group 1eRASS J045547.0-572404. Left: residual image of the group before and after subtracting the co-fit extended sources and point sources in the field. Nearby co-fit clusters and groups are shown with red circles, co-fit nearby bright AGN are shown with green circles, and the fitted galaxy group is shown with a blue rectangle. The smooth noise level indicates that the contaminant emission is modeled properly in the analysis. Right: surface brightness profile of the same group. The observed surface brightness profile of the image is plotted in blue, the best-fit model of the image is plotted in red, the PSF profile over the measured background is plotted in green, the PSF deconvolved surface brightness profile of the galaxy group is plotted in orange, the measured background level is shown with a horizontal dashed red line, and three characteristic radii of the group (0.15r 500c, r2500c and r500c from the core to the outskirt respectively) are shown with dashed gray lines.

3.4 Grouping and the spectral analysis

Measuring temperature through X-ray spectroscopy requires considerably more photons than measuring surface brightness properties with imaging analysis. Given the shallow nature of the eROSITA All-Sky Survey, the photon counts of most of the galaxy groups in our sample are insufficient for temperature measurements, even though the flux or luminosity of these objects can be reliably measured from X-ray images. For instance, more than half of the groups in the sample, shown on the right panel of Fig. 5, have fewer than 100 counts within the 0.3–1.8 keV band, which is not sufficient for measuring their temperature reliably.

The two canonical ways to overcome the problem of insufficient photon counts for spectral analysis are co-fitting or stacking. A plethora of examples of both techniques exist in the literature. For example, McDonald et al. (2014) co-fit radially extracted spectra of 80 South Pole Telescope (SPT) selected massive clusters, while Bulbul et al. (2014) and Zhang et al. (2024) stacked megaseconds of XMM-Newton and eROSITA spectra respectively to achieve a high signal-to-noise level and reveal faint spectral features. In this work, we employ the co-fitting technique to maintain the spectral information of individual groups that would be averaged out when stacked. This method is the most suitable for the primary goal of this work. Compared to stacking, this approach is computationally expensive; however, with the recently developed high-performance Central Processing Units (CPUs) and the improvements in parallel computing, we are able to employ the co-fitting technique in this work.

We grouped the sample such that the IGrM temperatures of the galaxy groups are similar in each bin. Moreover, we required the statistical constraining power (photon counts) of the groups in the same bin to be similar to each other to avoid the source with the highest count from biasing the measurements. In other words, our aim was to minimize the temperature and photon count variation − Δ T_500c ~ std(T_500c)¹⁰ and ΔC_500c ~ std(C_500c) – in each bin while trying to achieve a sufficient signal-to-noise ratio. To achieve this, we grouped the sample using the Voronoi binning technique (Cappellari & Copin 2003). To apply the tessellation technique, we pixelated our temperature proxy T_500c,sc (surface brightness inferred temperature estimate; see Sect. 3.3 for the details) and count measurements C_500c such that each pixel was occupied by only one galaxy group. The axes are then re-scaled, and the resulting image is given to the Vorbin package, the Python implementation of the Voronoi binning technique. The free parameters, namely, the axes scaling factors and the target S/N, are fine-tuned until the temperature variation (Δ T_500c) and the photon count variation (ΔC_500c) in the Voronoi bins are sufficiently small. The final binning scheme, shown in Fig. 6, is achieved using an S/N target, S /N = 22.36 (equivalent to 500 counts). We further present the distributions of the total counts and the number of galaxy groups in the Voronoi bins in Fig. 6. Using this binning scheme, we obtained 271 bins with a median count of 905, sufficient for obtaining reliable spectroscopic temperature measurements at T < 2 keV for each bin. We note that the binning scheme can be slightly different if a different target S/N or axis scaling factors are chosen; however, the impact of the chosen binning scheme is negligible on the final results as long as the resulting Δ T_500c and ΔC_500c are similar.

After the grouping, the source and local background spectra of the galaxy groups in our sample are extracted using the eSASS task, srctool. The source spectra are extracted from the circular regions centered around the galaxy group and have a radius of r_500c (see Sect. 3.3 for the details of the M_500c/r_500c estimation procedure). Similarly, the local background spectra are extracted from annuli that are centered around the galaxy group and have a radial range of 4r_500c < r < 6r_500c. The best-fit count-rate profiles are used to determine the masking radius of the bright point sources and nearby extended objects co-fitted during the imaging analysis. The remaining point sources and nearby extended sources within the extraction region are masked as described in Sect. 3.2.

We extract ancillary response files (ARFs) and redistribution matrix files (RMFs) using the srctool task for the background and the source region in different settings to be assigned to various components of the source and background models. The ARF assigned to the source component is extracted with the exttype=BETA and psftype=2D_PSF settings to consider the energy-dependent PSF and vignetting corrections. Over the extraction regions, the flux distribution of the vignetted X-ray background is assumed to be flat, and the exttype=TOPHAT and psftype=NONE settings are used for extracting the ARFs assigned to the vignetted X-ray background components of the source and background regions.

Similar to our eFEDS analysis (Ghirardini et al. 2021; Liu et al. 2022; Bahar et al. 2022; Bulbul et al. 2022), the local background model consists of two major components: particle- induced instrumental background (see Bulbul et al. 2020; Freyberg et al. 2020, for further details) and X-ray background including the Galactic foreground, and unresolved point sources in the sky. The total model includes a spectral model component with an absorbed thermal component. The spectral analysis was performed using PyXspec, the Python interface of the standard X-ray spectral analysis package Xspec (version 12.12.1, Arnaud 1996), which employs the AtomDB atomic database (version 3.0.9, Foster et al. 2012). The Xspec model of the X-ray foreground consists of an unabsorbed APEC (Smith et al. 2001) for the local hot bubble (Yeung et al. 2023, T ~ 0.084 keV), two absorbed APECs for the hot and cold components of the galactic halo (Ponti et al. 2023; Bulbul et al. 2012, T ~ 0.49 and 0.157 keV respectively). To model the cosmic X-ray background, we use an absorbed power-law for the unresolved AGN (Cappelluti et al. 2017, Γ = 1.45). For the shape of the instrumental background, we use the best-fitting model of Yeung et al. (2023) obtained by calibrating the filter wheel closed (FWC) data (see Appendix A.1. and A.2. of Yeung et al. 2023, for the details of the modeling of FWC data). This instrumental background component is folded with unvignetted ARF (Freyberg et al. 2020), while the cosmic X-ray background and Galactic foreground are folded with the respective vignetted ARF in the fits.

To account for the X-ray absorption, we use the TBAB S (Wilms et al. 2000) interstellar medium (ISM) absorption model in Xspec (Arnaud 1996). We use the HI4PI survey (HI4PI Collaboration 2016) for calculations of the hydrogen column density (n_H). The n_H values at the positions of eRASS1 galaxy groups are relatively low because of their locations at higher Galactic latitudes; therefore, using the total hydrogen column density (n_H,tot) rather than the neutral hydrogen column density (n_H,I) has a negligible impact on our results at the sample level¹¹ as also noted in Bulbul et al. (2024). We use the solar abundances of Asplund et al. (2009) when measuring the metallicity of the groups. We use C-statistic (Cash 1979) for the statistical interpretation of our spectra that provides unbiased estimates of the model parameters at the low and high count regimes (Kaastra 2017). We employ the co-fitting technique for the spectral analysis. This required us to explore likelihoods with relatively high dimensional parameter space. For this purpose, we chose to employ the MCMC fitting technique. Xspec has a built- in MCMC sampler; however, the amount of control it allows the user over the priors is limited. For this reason, we employ the widely used MCMC sampler emcee (Foreman-Mackey et al. 2013) rather than the built-in Xspec sampler to explore high dimensional likelihoods. We achieved this by developing an interface that allows cross-talk between PyXspec and emcee and updates the model parameters at every MCMC step accordingly.

Galaxy groups are low-mass objects with relatively low plasma temperatures (T < 2 keV) due to their shallower potential wells. They share this low-temperature parameter space with other background/foreground components, such as the cold (T ~ 0.157 keV) and hot (T ~ 0.49 keV) components of the galactic halo (Ponti et al. 2023). This results in degeneracies between the source and background/foreground components at the low count regime. At a given energy band, it is relatively easy to separate the source and background components using the 2D distribution of photons through imaging analysis since we expect the local background rate to be relatively flat, whereas the source emission roughly follows a projected Vikhlinin profile (Vikhlinin et al. 2006). In this work, we make use of this fact and combine the spatial and spectral information of photons following a novel approach with the aim of lifting the aforementioned degeneracies the best we can. We achieve this by using the observed countrates (0.3–1.8 keV) measured through imaging analysis as priors in the spectral analysis in our pipeline that combines PyXspec and emcee.

To obtain the average temperatures of the binned galaxy groups (see the first paragraph of this section for the details of the grouping), we only link the temperature parameter of the (APEC) model and co-fit all the source and background spectra of the galaxy groups in each Voronoi bin. In total, 2 × N_i,gr spectra are co-fit (one spectrum for the source region and one for the background region) where N_i,gr is the number of galaxy groups in the i’th Voronoi bin. During the fitting, the temperatures, the normalizations of the X-ray background components, and the normalizations of the unvignetted particle background components are allowed to be free with flat priors in the logarithmic parameter space.

Many studies in the literature show strong degeneracy between the temperature and metallicity measurements for galaxy clusters and groups hosting multi-phase gas. This manifests as the so-called ‘Fe bias’ (Buote & Fabian 1998; Buote 2000; Gastaldello et al. 2021; Mernier et al. 2022) where temperature and abundance measurements change depending on the number of gas components fitted. In our work, we take this effect into account by allowing the metallicities of each galaxy group to be free with a Gaussian prior centered around 0.3 Z_⊙ with a standard deviation of 0.025 while measuring average temperatures. The normalization of the source emission of each galaxy group and the co-fit temperature are left free with log-uniform priors. Consequently, following the spectral co-fitting analysis procedure described above, we obtained 271 average temperature measurements within r_500c for each Voronoi bin shown in Fig. 6.

Fig. 5 Left: soft-band (0.5–2 keV) X-ray luminosity (LX,500c) and redshift (z) distributions of the galaxy group sample used in this work consisting of 1178 objects. The luminosity of the groups cover a range of 3.9 ×1040–1.4 × 1043 ergs s−1, and the redshift span of the sample is 0.003–0.48. Right: count (C500c) and mass (M500c) distributions of the 1178 galaxy group used in this work. Measured counts of galaxy groups range between 10–14 380, and their masses span a range of 2.1 × 1012–1014 M⊙. Median values of the LX,500c, z, C500c and M500c observables are 6.3 × 1042 ergs s–1, 0.11, 95, 6.5 × 1013 M⊙ respectively.

Fig. 6 Left: Voronoi binning scheme used for grouping the sample obtained from the distribution of count (within r500c) measurements (C500c) and scaling relation based temperature estimates (T500c,sc). Top right: histogram of total counts in 271 Voronoi bins with a median of 905 counts. Bottom right: Histogram of galaxy groups in Voronoi bins.

3.5 Electron density, temperature, and entropy profiles

From the imaging analysis described in Sect. 3.2, we obtain deprojected emissivity profiles of all the 1178 galaxy groups. The spectral analysis, described in Sect. 3.4, yields average temperature measurements of the 271 galaxy group bin within r_500c. The electron number density profile measurements of the IGrM have a non-negligible dependence on the assumptions on the temperature and metallicity profiles. Furthermore, temperature profiles of the galaxy groups are needed for obtaining entropy profiles. For most binned groups, only a single average temperature and metallicity measurement can be achieved within r_500c due to low S/N eROSITA data. We overcome this limitation and incorporate the impact of temperature variation on the thermodynamic properties as a function of radius, as described below.

To account for the radial temperature variation, we first determine the average shape of the 3D temperature profile of groups (T(r)/T(r < r_500c)) using the temperature profile measurements of the tier 1 and 2 groups presented in S09 (see Sect. 4.1 for the details). The average and individual shapes of the T (r)/T (r < r_500c) profiles are shown in Fig. 7. We then rescale the average shape with the integrated temperature measurements and obtain the average temperature profiles for each binned group. The temperature profile measurements presented in S09 for a sample of 43 groups are obtained by analyzing deep Chandra observations; therefore, the overall shape of the profiles is relatively well-constrained. Following this procedure, we obtain the average temperature profiles of the 271 galaxy group bins. We note that our approach of obtaining temperature profiles of groups is equivalent to fixing the shape of an assumed temperature profile and fitting spectra by allowing the normalization of the profile to be free. We further note that the observed temperature measurement discrepancy between telescopes (e.g., Liu et al. 2023) does not affect our work given that the temperature measurements of eROSITA for galaxy groups agree very well with the Chandra and XMM-Newton temperatures (Migkas et al. 2024).

For the metallicity profile, we consider the following studies in the literature that have reasonably large galaxy group samples: Sun (2012), Mernier et al. (2017), and Lovisari & Reiprich (2019). We find that their measurements agree relatively well within the scatter of the metallicity profile reported in Mernier et al. (2017) (see Fig. 8). For this reason, we adopt the median Mernier et al. (2017) metallicity profile (Z_M17) for our measurements and consider the scatter of their profile as our systematic uncertainty. We quantify the impact of this uncertainty on the thermodynamic profile measurements (n_e(r) and S (r)) and consider the resulting difference as part of the overall error budget (see Sect. 4.1 for the details of the assumed metallicity profile and quantification of the impact of this choice).

We then calculate the electron density profiles from our deprojected emissivity measurements by constructing APEC models in Xspec (similarly done in Liu et al. 2022) having the temperatures equal to the temperature profiles of the binned groups and the metallicities equal to the assumed metallicity profile at all radii out to 2r_500c. During the n_e profile calculation, in addition to the errors of the imaging analysis, the uncertainties of the average metallicity and temperature measurements are propagated as well using the MCMC chains of the spectral analysis, taking into account the covariances. Then, the electron density profiles of the objects within each Voronoi bin are averaged to get the average electron density profiles of the binned groups.

Lastly, the entropy profiles of the binned groups are obtained by combining the average electron density and the temperature profiles using the equation below: $$S(r)=T(r)/n_{\text{e}}{(r)}^{2/3}.$$(8)

The entropy profiles are then sampled at three characteristic radii, and the final entropy measurements of 271 galaxy group bins are obtained. The full shape of the entropy profiles of binned groups, along with the other thermodynamic profiles such as electron density and pressure (P = n_e T), will be presented in Bahar et al. (in prep.).

Besides the systematics resulting from the metallicity profile, we also consider other major systematics that have a non- negligible impact on the thermodynamic properties measured in this work. We discuss and quantify the impact of these systematics on our measurements in Sect. 4 and take them into account as part of the total error budget when we draw conclusions in the next section.

Fig. 7 Normalized temperature profiles, T(r)/T(r < r500c), of 23 groups (gray) presented in S09 and their average (blue). The blue dashed line provides an average conversion ratio between the temperature profile, T(r), and the characteristic temperature measurements, T(r < r500c).

Fig. 8 Average metalicity profiles of galaxy groups reported by Mernier et al. (2017) and Lovisari et al. (2015) and the stacked metallicity profiles of Sun (2012) for galaxy groups in two temperature bins (T = 0.75–1.3 keV and 0.75–1.3 keV).

4 Assumptions, corrections, and systematics

Having fair comparisons between the thermodynamic properties of groups observed with different X-ray observatories and simulations is a challenging task as various systematics should be taken into account in the measurements, such as the systematic uncertainties on the metallicity and temperature profiles, systematic uncertainties on the group masses, the flux calibration mismatches between instruments and the systematics resulting from the use of different atomic database versions. In this section, we provide a list of assumptions, corrections, and systematics taken into account in this work, along with our approach to account for them. We also list a summary of the description and implementation of the assumptions, corrections, and systematics in Table 1.

4.1 Assumptions on temperature and metallicity profiles

Accounting for the temperature and metallicity radial variation is key to having reliable thermodynamic profiles. Shallow survey observations and low signal-to-noise data of most groups in our sample are insufficient to measure temperature profiles reliably. There are various studies in the literature on the average temperature profile of the hot gas in clusters (e.g., McDonald et al. 2014; Ghirardini et al. 2019), while the studies focusing on the shape of the average temperature profile of groups with a large enough sample are limited. An in-depth study of 43 nearby galaxy groups with deep Chandra observations by S09 (among which 23 of them have good temperature constraints out to r_500c ) is one of the few studies we compared within this work. In this work, we used the temperature profile measurements of these 23 groups to get the average shape of the temperature profile of groups. To get the average shape, we first calculated the characteristic temperatures of the groups, T(r < r_500c), by projecting and integrating all the temperature profiles within a cylindrical volume of radius r_500c . We achieved this by following the Mazzotta et al. (2004) weighting and projection formulas $$w=n_{\text{e}}{n}_{\text{p}}{T}^{\alpha }$$(9)

and $$T_{\text{sl}}=\frac{{\displaystyle \int w}TdV}{{\displaystyle \int w}dV},$$(10)

where for α we used 0.76, which is calibrated for eROSITA (ZuHone et al. 2023). Furthermore, we used our average group electron density profile for n_e (Bahar et al., in prep.). We then divided the temperature profiles with the characteristic temperatures and obtained the normalized temperature profiles, T(r)/T(r < r_500c). Lastly, we took the average of these profiles and renormalized them to get the average shape of the temperature profiles of groups. The average and individual normalized profiles of 23 groups are shown in Fig. 7.

Unlike clusters, the band-averaged cooling function of groups has a strong metallicity dependence because of the significant contribution from the line emission at temperatures, T < 2 keV. Therefore, the radial change in metallicity from the center to the outskirts should be accounted for to calculate electron density profiles accurately.

During the last decade, the shape and the strength of emission lines have significantly changed in most commonly used plasma emission codes (see Sect. 4.4) that had a strong influence on the metallicity measurements of groups (e.g., Mernier et al. 2018). For this reason, it is important to use the most recent publications and account for uncertainties in metallicity profiles in the systematics error budget. Among the metallicity profile measurements in the literature, the results presented in Sun (2012), Mernier et al. (2017) and Lovisari & Reiprich (2019) stand out as the most recent studies with moderately large galaxy group samples with sufficiently deep observations. Fig. 8 presents the stacked metallicity profiles of Sun (2012) and the average metal- licity profiles of Mernier et al. (2017) and Lovisari & Reiprich (2019) that are renormalized based on the iron abundance ratio of Asplund et al. (2009). In Sun (2012), the author reports stacked abundance profiles of 39 galaxy groups in three temperature bins (0.75–1.3, 1.3–1.9 and 1.9–2.7 keV). In this work, we consider only the results of the first two temperature bins (0.75–1.3 keV and 1.3–1.9 keV), which are relevant to our sample that has a median temperature of T(r < r_500c) = 1.45 keV.

Overall, the average metallicity profile reported in Mernier et al. (2017) lies between the Sun (2012) and Lovisari & Reiprich (2019) measurements and the Sun (2012) measurements in the 1.3–1.9 keV temperature bin lies above the Mernier et al. (2017) profile and the average measurements of Lovisari & Reiprich (2019) lie below the Mernier et al. (2017) profile. When calculating the thermodynamic properties, the differences in metallicity measurements must be accounted for as systematics because of the strong dependence of emissivity on metallicity at group scales. Given the large spread of metallicity measurements, we take the average profile of Mernier et al. (2017) as our default profile and conservatively consider the shaded area as the systematics of the average profile measurements. To account for the impact of the choice of average metallic- ity profile, we construct APEC spectra in Xspec and obtain deprojected electron density profiles of the 1178 galaxy groups in our sample (similarly done in Liu et al. 2022) using the scaled temperature profiles and three metallicity profiles (low- scattered, median, and up-scattered Z_M17 profiles) shown with red in Fig. 8. We present the ratios of the electron densities in Fig. 9 that are obtained by using the aforementioned three metallicity profiles (low-scattered Z_M17: n_e,L, median Z_M17: n_e, and up-scattered Z_M17: n_e,U) for 1178 groups at the three characteristic radii (0.15r_500c, r_2500c and r_500c). The ratios (n_e,U/n_e and n_e,L/n_e) ranging between 0.84-1.24 in Fig. 9 indicating a non-negligible difference between the electron density measurements. The ratios deviate more from unity as the characteristic temperature decreases, and the measurement radius increases. This is due to line emission, coupled with metalicity, which plays a more important role as the temperature decreases. The procedure described above is followed for the final results, and the radius/temperature dependent systematics due to the choice of metallicity profile are quantified and propagated to our final entropy measurements presented in Sect. 5.

Fig. 9 Ratio of the electron densities obtained by assuming low/up-scattered Mernier et al. (2017) metallicity profiles to the electron densities obtained by the median Mernier et al. (2017) metallicity profile at three radii, 0.15r500c, r2500c, and r500c as a function of characteristic temperature, T(r < r500c). Green, blue, and purple data points represent the ratio between the electron densities obtained by assuming the lower envelope (ne,L) of the red shared area and the dark red median line (ne) in Fig. 8. Orange, yellow, and red data points represent the ratio between the electron densities obtained by assuming the upper envelope of the red shared area (ne,U) and the dark red median line (ne) in Fig. 8.

4.2 Correction for the flux discrepancy

A discrepancy of 15% is reported in the luminosity measurements in the soft-band (0.5–2 keV) of a subsample of massive galaxy clusters observed with both eROSITA and Chandra (Bulbul et al. 2024). The observed flux difference is constant with no flux or luminosity dependence. Some of this difference can be explained by the photon loss in the latest processing due to the higher CCD thresholds (Merloni et al. 2024); however, further investigation is required to understand the observed flux discrepancy, which could be due to various calibration effects. We account for the flux discrepancy while comparing our entropy measurements with those reported in the literature. Among our two main X-ray observables n_e and T, only electron density is impacted by the flux discrepancy since the spectroscopic T measurements are not sensitive to the overall flux normalization. To roughly estimate the impact, we assumed the shape of the measured electron density profile to be the same for different instruments, used the fact that $L\propto n_{\text{e}}^2$, and obtained a ratio of n_e,eRO/n_e,Cha = 0.85^0.5 ~ 0.92 between the electron density measurements of Chandra and eROSITA. The 8% underestimation of n_e corresponds to a 5% overestimation of entropy. This fraction is factored in the Chandra measurements in S09 when comparing with the eROSITA results in Fig. 10.

4.3 Systematics related to mass measurements

We note that entropy measurements at overdensity radii, 0.15r_500c, r_2500c, and r_500c are sensitive to the assumed masses of the galaxy groups. This dependence is due to the entropy profile of galaxy groups being a strong function of the radial distance and the overdensity radius, which is a mass-dependent quantity (e.g. Bulbul et al. 2010; Ghirardini et al. 2019). Therefore, any disagreement in radius and mass may lead to a bias in the measured thermodynamic profiles and their comparisons between different methods. Masses of galaxy groups can be estimated in different ways, such as by assuming hydrostatic equilibrium, using the shear information of the lensed galaxies, or using scaling relations; however, these methods have advantages and disadvantages along with introduced biases. Comparison of the mass estimation techniques for galaxy groups is beyond the scope of this paper; therefore, in this paper, we account for the bias introduced while comparing our results with the literature. We find that our scaling relations based r_500c estimates are ∼11% higher than the hydrostatic equilibrium based r_500c estimates in S09 for 9 crossmatched groups. Comparison between our characteristic radii estimates (r_500c) and those in the literature (r_500c,S09) is provided in Fig. 11. A discrepancy of ∼11% approximately corresponds to a bias of ~37% (1.11³ = 1.37) on M_500c. This result is close but slightly below the 45% hydrostatic mass bias observed in galaxy groups (Nagai et al. 2007), see also Sect. 6.2 of S09. We note that the masses of S09 are obtained with an outdated version of ATOMDB. Lovisari et al. (2015) and Sun (2012) independently confirmed that using a more recent ATOMDB version (v2.0.1) increases the temperatures by ∼15%. Such an increase would reduce the mass mismatch to a ∼22% level and the radius mismatch to ∼7% level.

One should note that constraining hydrostatic mass bias, especially in galaxy groups, is challenging due to our limited knowledge of the magnitude of non-thermal pressure support, magnified in galaxy groups due to powerful AGN feedback compared to galaxy clusters. Moreover, the fraction of discrepancy may be due to other systematics in the measurements, such as the representation of the L_X – M relation with a single powerlaw spanning a large mass range. We provide the r_500c estimates of the 9 crossmatched groups in Table 2 along with the slopes of the average entropy profiles of the binned groups at the three characteristic radii in Table 3.

The mass estimate dependence is also responsible for the observed scatter in the entropy of the sample. Fig. 12 shows a strong correlation between the scatter of the average entropy measurements at a given temperature and the average concentration ($c_{\text{SB},r_{500c}}=SB\left(r<0.1r_{500c}\right)/SB\left(r<r_{500c}\right)$) of the sample (see Sanders et al., in prep. for the concentration measurements). This is expected since the M_500c (or r_500c) estimates used in this work are obtained using an LX – M scaling relation. At a given mass, the cool-core galaxy groups with higher luminosity would have higher entropy measurements while the others scatter around the sample’s median. The average entropy and temperature measurements of the sample are presented in Fig. 12, where the colors of the data points indicating average concentration obtained by averaging the measurements provided in Sanders et al. (in prep.). A few extreme cases with a larger concentration, electron number density, and entropy are easily noticeable in Fig. 12. By construction, the r_500c (or M_500c) estimates are, on average, unbiased at the sample scale. Therefore, a small scatter does not significantly impact the conclusions in this work. On average, these effects cancel out such that at all three radii, the error-weighted average entropy plotted in magenta coincides with the intermediate concentration values ${\log }_{10}\left(c_{\text{SB},r_{500c}}\right)~-0.7$.

Fig. 10 Comparison between the redshift evolution scaled average eROSITA entropy measurements (white diamonds) with the Chandra measurements of 43 galaxy groups (black error bars) presented in S09 at three radii (0.15r500c, r2500c, r500c) as a function of characteristic temperature, T(r < r500c). Blue error bars represent the statistical uncertainties of the average measurements, and the red error bars represent the overall error budget resulting from the statistical and systematic uncertainties (see Sect. 4 for the details of the accounted systematics). For consistency, coreexcised temperatures presented in S09 are converted to core-included temperatures, and the entropy measurements presented in S09 are normalized such that the new data points are measured at the same angular radii with eROSITA flux calibration (see Sect. 5 for the details of these corrections).

Fig. 11 Comparison between the scaling relation based r500c estimates of the crossmatched galaxy groups obtained in this work (y-axis: r500c) with the estimates of S09 (x-axis: r500c,S09) obtained by assuming hydrostatic equilibrium.

4.4 Systematics related to the atomic databases

Over time, our knowledge of the strength of atomic transitions has changed significantly for T < 2 keV plasma. The change in our knowledge of the strength and the width of the Fe-L complex is the most relevant to this work, given its key role in determining the temperature and metallicity of the IGrM. The change in our knowledge resulted in significant updates on the predictions of the spectral models, APEC (Smith et al. 2001) and CIE (an updated version of the MEKAL model, Mewe et al. 1985, 1986; Liedahl et al. 1995), that retrieve data from atomic databases, AtomDB (Foster et al. 2012) and SPEXACT (Kaastra et al. 1996, 2020), around the Fe-L complex (e.g., see Fig. 1 in Gastaldello et al. 2021).

These updates may cause measurement mismatches between the studies that employed different versions of the atomic databases. For this reason, the possible impact of using different versions of the atomic databases should be taken into account while comparing results with the literature. In fact, Sun (2012) reports that for the T < 2 keV plasma, the temperature measurements presented in S09 increase by 10–20% when a more recent AtomDB version, v2.0.1, is used instead of the one employed in S09, v1.3.1. The increase in the temperature is also independently confirmed by Lovisari et al. (2015), who found their temperature measurements to be 13% higher than the temperatures of S09 for the seven crossmatched objects. In this work, we take this change into account and apply a 15% correction to the entropy and temperature measurements of S09 to their T < 2 keV groups while comparing their results with the eROSITA measurements (e.g., in Fig. 10). With this correction, we aim to have the fairest comparison of our results.

Lovisari et al. (2015) further report that in the case of such an update in the atomic database, the normalization of the APEC spectrum reduces ∼10% for a group (NGC3402) residing at the low-T parameter space (T ∼ 1 keV), where the impact of the change is expected to be the most prominent. The density of the plasma scales as the square root of the APEC normalization, which results in a ∼3% increase in entropy. Given that the 3% increase in entropy is an upper limit (as explained in Section A.3 of Lovisari et al. 2015) and is well within the entropy error bars presented in S09, we did not propagate the impact of this change to the S09 entropy measurements.

Lastly, Sun (2012) reported that abundance measurements drop ∼20% after updating the AtomDB version (from 1.3.1 to 2.0.1). This reduces the stacked metallicity data points of Sun (2012) by ∼20%; however, given that we are employing Mernier et al. (2017) profile in our work, our results are not affected by this. We further note that even if we apply such a correction to Sun (2012) measurements, they will be well within our conservative systematic error bars for the metallicity profile (shaded area in Fig. 8); therefore, it does not pose any challenge to our measurements.

It can be seen from Fig. 1 of Gastaldello et al. (2021) that the width and the normalization of the Fe-L complex seemed to be converging. However, it is hard to know how far we are from the absolute calibration. Therefore, future work should keep systematics related to the atomic databases in mind while comparing results from the literature.

5 Entropy and characteristic temperature measurements

In this work, we constrain the entropy of the IGrM, utilizing the deep eRASS:4 observations of the galaxy groups detected in the eRASS1 survey. In this section, we present our measurements of the characteristic temperature, T(r < r_500c), and entropy at three overdensity radii, 0.15r_500c, r_2500c, and r_500c, from the combined analysis of 1178 galaxy groups. We then compare our findings with the previously reported results in the literature.

Following the procedures described in Sects. 2 and 3, we obtain average entropy and temperature measurements of a sample of 1178 eROSITA selected galaxy groups in 271 bins. Our measurements for the average temperature and average entropy at three characteristic radii scaled by the self-similar redshift evolution (E(z)^4/3) can be seen in Fig. 13, where the error- weighted average entropy measurements are shown in white circles, the statistical uncertainties are shown with blue error bars and the overall error budget resulting from the statistical and systematic (see Sect. 4) uncertainties are shown with red error bars. Error-weighted average of the redshift evolution scaled entropy measurements, their statistical uncertainties, and the impact of systematics are also provided in Table 3 as a function of IGrM temperature. Besides the average entropy measurements, we provide slopes of the average entropy profiles in Table 3 at the three characteristic radii for future work to account for the mass measurement systematic while comparing with our results (see Sect. 4.3 for the details of accounting the systematics in mass measurements). We note that the few outlying entropy measurements in Fig. 13 are due to the extreme cool-core objects and have negligible impact on our error-weighted entropy measurements (see Sect. 4.3 for a more detailed discussion on the outliers and their impact on our final results).

We find that the characteristic temperature measurements of our binned group sample span a range of ~0.73–2.68 keV. Furthermore, we find that the characteristic temperature measurements of the binned groups correlate well with the entropy measurements at the three characteristic radii such that, an increase in the characteristic temperature from 0.73 keV to 2.68 keV corresponds to an increase from the redshift scaled entropy levels of 121, 271, 341 keV cm² to 404, 722, 1135 keV cm² respectively. This trend can also be clearly seen from the error-weighted average measurements (white circles) in Fig. 13 such that the error-weighted averages at these three radii increase from the levels of 145, 335, 455 keV cm² to 316, 635, 905 keV cm² respectively. It can further be noticed that as the characteristic temperature of the binned groups decreases, the statistical uncertainty of the error-weighted average profiles (blue error bars) decreases. Conversely, the decrease in the temperature corresponds to an increase in the overall error budget (red error bars). This is due to the fact that the line emission from the Fe-L complex at the low-temperature parameter space becomes significant, which provides an additional spectral feature for measuring the temperature and reduces the statistical uncertainty. Meanwhile, at the same parameter space, the systematic uncertainty of the metallicity profile has a very large impact on the electron density measurements, which results in the entropy measurements having large systematic uncertainties. At higher temperatures, the systematic uncertainties of the electron density measurements become smaller due to the reduced line emission, and the less significant Fe-L complex increases the statistical uncertainties.

Furthermore, we find that the redshift scaled average entropy and the average temperature measurements of our sample at the three radii follow power-law relations within the uncertainties. To quantify the normalizations and the slopes of the relations, we fit a power-law model to the measurements of the form S (r)E(z)^4/3 = A(T/ T_piv)^B where the S (r) term stands for entropy measured at three characteristic radii (0.15r_500c, r_2500c, and r_500c), the T term stands for the characteristic temperature (T(r < r_500c)), T_piv term stands for the pivot temperature, and the A and B terms stand for the normalization and the slope of the relation respectively. For our galaxy group sample, we took a pivot value of Tpiv = 1.44 keV and obtained the best- fit power-law models to the average measurements at the three radii as shown in Fig. 14. Furthermore, we noticed that the slope of the S (r) – T relations seem to be changing around T = 1.44 keV such that the warm/higher-mass groups in our sample (T = 1.44–2.68 keV) are steeper than the S (r) – T relations of the cool/lower-mass groups (T = 0.73–1.44 keV). To quantify the difference, we separately fit the average measurements of the warm and cool groups using the same relation, assuming pivot values of T_piv = 1.79 and 1.08 keV, respectively. Best-fit parameters of the S (r) – T relation at three radii for all galaxy groups in our sample, along with the best-fit parameters for the warm, and cool groups, are listed in Table 4.

As a result of the fitting procedure, we find that for all galaxy groups (T = 0.73–2.68 keV), the slopes of the S (r_2500c) – T(r < r_500c) and S (r_500c) – T(r < r_500c) relations ($B={0.70}_{-0.07}^{+0.07}\text{\hspace{0.17em}and\hspace{0.17em}}{0.69}_{-0.09}^{+0.09}$ respectively) are in very good agreement with each other whereas we find that the slope of the S (0.15r_500c) - T(r < r_500c) relation ($B={0.86}_{-0.06}^{+0.06}$) being slightly steeper. We also compare our relations with the self-similar predictions presented in Voit et al. (2005). We find that our best-fit normalizations of the S (r)–T relations are >5σ larger than the self-similar predictions at all radii. This agrees well with the previous findings on the entropy excess in galaxy groups and indicates that non-gravitational processes such as AGN feedback and radiative cooling are non-negligible and play an important role in shaping the entropy profiles of galaxy groups. Moreover, we compared the best-fit slopes of the $S\left(0.15r_{r500c}\right)-T\left(r<r_{500c}\right),S\left(r_{2500c}\right)-T\left(r<r_{500c}\right),S\left(r_{500c}\right)-T\left(r<r_{500c}\right)$ relations for all galaxy groups in our sample with the self-similar prediction (B_self = 1 at all radii) and found that our slopes are 2.3, 4.3, and 3.4σ shallower than the self-similar slope.

Furthermore, we compared the best-fit S (r) – T relations of the warm and cool groups and found that the slopes of the warm/higher-mass groups, $B={0.98}_{-0.10}^{+0.10},{0.87}_{-0.11}^{+0.11}\text{\hspace{0.17em}and\hspace{0.17em}}{0.72}_{-0.17}^{+0.17}$ are much steeper than the slopes of the cool/lower-mass groups, $B={0.59}_{-0.16}^{+0.16},{0.32}_{-0.22}^{+0.22}\text{\hspace{0.17em}and\hspace{0.17em}}{0.56}_{-0.28}^{+0.28}$respectively. We then compared our best-fit normalizations of the warm and cool groups with the self-similar predictions presented in Voit et al. (2005) and found that the normalizations for warm and cool groups are >5σ larger than the self-similar predictions at all radii. Moreover, we also compared the slopes of the warm and cool groups with the self-similar prediction and found that the best-fit slopes of the warm/higher-mass groups, $B={0.98}_{-0.10}^{+0.10},{0.87}_{-0.11}^{+0.11}\text{\hspace{0.17em}and\hspace{0.17em}}{0.72}_{-0.17}^{+0.17}$, agree much better with the self-similar prediction (B_self = 1) compared to the slopes of the cool/lower-mass groups, $B={0.59}_{-0,16}^{+0.16},{0.32}_{-0.22}^{+0.22}\text{\hspace{0.17em}and\hspace{0.17em}}{0.56}_{-0.28}^{+0.28}$ 0.15r_500c, r_2500c, and r_500c respectively. This comparison suggests that even though the non-gravitational processes significantly increase the overall entropy levels of galaxy groups at all temper- ature/mass scales, they result in the slope of the S(r) – T relation for cool/lower-mass groups deviating more from the self-similar prediction compared to the slope of the warm/higher-mass groups that live at a temperature/mass parameter space closer to clusters.

We note that the observed flattening for the cool groups could also be due to other reasons, such as the selection effects or relatively large systematic uncertainties. The X-ray selection probability of galaxy clusters and groups has a dependency on the emission profile of the cluster and group (Clerc et al. 2024). Therefore, for the most robust calibration of the S (r) – T scaling relation, a selection function that takes entropy at given radii as input should be used, which is currently not implemented in our framework for galaxy groups. In the future, a better understanding of the temperature and metallicity profiles of galaxy groups, along with the use of a profile-dependent selection function for fitting the S (r) – T relation, would make the picture clearer and help us understand the origin of the observed flattening for the cool/lower-mass groups.

Comparing the entropy measurements at the three radii, with the previous measurements in the literature as a function of temperature, T (r < r_500c), is challenging given the limited number of studies reporting these quantities and the information provided in these studies to account for the discrepancies in the measurement radii being limited across the literature. In this section, we present our findings from the comparison performed between our measurements and the results reported by S09 and Johnson et al. (2009).

To have a fair comparison, we quantify and account for the systematic differences between our results and those presented in S09. The average temperature of galaxy groups in S09 are measured using core-excised apertures (0.15r_500c < r < r_500c), whereas in this work, we measure temperatures including the core out to r500c to maximize the signal-to-noise in our measurements. We convert the core-excised characteristic temperature measurements, T(0.15r_500c < r < r_500c), reported in S09 to the core-included temperatures by applying a conversion factor of T(r < r_500c)/T(0.15r_500c < r < r_500c) = 1.07. This conversion factor is obtained by projecting the average group temperature profile (see Sect. 4.1) using the temperature weighting and projection formulation of ZuHone et al. (2023) calibrated for eROSITA and the average group electron density profile (Bahar et al., in prep.). Additionally, we apply a correction factor to account for the difference in the characteristic radii estimation to the S09 entropies as described in Sect. 4.3. The correction factors are obtained by comparing our scaling relations based on r500c estimates with the masses obtained assuming hydrostatic equilibrium (see Sect. 4.1) and quantifying the impact of the mismatch on the average entropy measurements at three characteristic radii. Furthermore, we apply a 5% correction to the entropy measurements of S09 to account for the 15% flux discrepancy between Chandra and eROSITA (see Sect. 4.2) that is reported in Bulbul et al. (2024). Lastly, we applied a 15% correction to the entropy and temperature measurements of the T < 2 keV groups in the S09 sample due to the use of an old version of the atomic database AtomDB (see Sect. 4.4 for the details on the correction).

The comparison between our temperature and entropy measurements with the findings of S09 are shown in Fig. 10. After the systematic differences are accounted for, we find our measurements to agree well with the S09 results within ~1σ confidence at all three radii. Fig. 10 shows that our results only deviate from the measurements of S09 at the cooler temperatures with T(r < r_500c) < 1.15 keV. We argue that this may be due to the completeness of our sample being much higher than the completeness of the S09 sample at these temperatures. The S09 sample is based on the archival Chandra follow-up observations of 43 ROSAT-detected bright groups. For this reason, by construction, the limited number of measurements reported in S09 at the low-temperature parameter space is for bright groups that have relatively high n_e and low S . This effect is visible in Fig. 10 such that S09 has only three measurements for S (r_500c) and 11 for S (r_2500c) and S (0.15r_500c) at the cool temperatures end, whereas our sample includes temperature measurements of a sample of 323 galaxy groups binned into 61 groups. Furthermore, we note that the sample studies analyzing individually followed-up systems, such as S09, exclude morphologically disturbed systems. Nevertheless, in our work, with the aim of achieving comprehensive conclusions about galaxy groups, we did not make such a distinction and have analyzed all the galaxy groups in our sample by allowing the centroid of the surface brightness profile to be free. Parallel to having a preferential selection toward bright objects, excluding morphologically disturbed systems can also potentially decrease the average entropy of a sample at a given mass/temperature. That is in agreement with the direction of the slight discrepancy we see in Fig. 10 and can possibly explain the mismatch. A more quantitative statement on the relationship between the morphology and gas properties of galaxy clusters and groups, as well as the relationship between the eROSITA selection and the morphology of the extended objects, will be explored in Sanders et al. (in prep.).

Furthermore, Johnson et al. (2009) analyzed XMM-Newton observations of 28 nearby galaxy groups and found that the entropy measurements of their sample at the core, S (0.1r_500c), follow a power-law relation with their core-excluded temperature measurements, T(0.1r_500c < r < 0.3r_500c). They further report that the best-fit slope of their S (0.1r_500c)−T(0.1r_500c < r < 0.3r_500c) relation (0.79 ± 0.06) agrees with the slope of a similar relation, S (0.15r_500c)−T(0.15r_500c < r < r_2500c), presented in S09, 0.78 ± 0.12. In our work, we measure entropy at the core at a slightly different radius (0.15r_500c) than the one used in Johnson et al. (2009), and we use core-included temperature measurements rather than the core-excised temperatures. Nevertheless, we compare our most similar relation, S (0.15r_500c)−T(r < r_500c), to the relations fitted in S09 and Johnson et al. (2009) and find that our slope 0.86 ± 0.08 is in good statistical agreement with the slopes of S09 and Johnson et al. (2009). Moreover, we further compare our best-fit slope for the S (r_2500c)−T(r < r_500c) relation with the best-fit slope of another similar relation, S (r_2500c)−T(0.15r_500c < r < r_500c), presented in S09 and find that our slope (0.7 ± 0.07) agrees well with the slope reported in S09 (0.76 ± 0.06). Lastly, we compare the best-fit slope of our S (r_500c)−T(r < r_500c) relation with the best-fit slope of S09 for the S (r_500c)−T(0.15r_500c < r < r_500c) relation and find that our slope (0.69 ± 0.09) is in good statistical agreement with the slope reported in S09 (0.8 ± 0.2).

Fig. 12 Average entropy measurements of binned groups (diamonds) at three radii (0.15 r500c, r2500c, r500c) as a function of characteristic temperature, T(r < r500c). The colors of the diamond data points represent the average concentration ($c_{\text{SB},r_{500c}}=SB\left(r<0.1r_{500c}\right)/SB\left(r<r_{500c}\right)$) of the groups within the corresponding Voronoi bin. Error-weighted averages of the diamond data points are shown in magenta.

Fig. 13 Redshift evolution scaled average entropy measurements of binned groups (gray diamonds) at three radii (0.15r500c, r2500c, r500c) as a function of characteristic temperature, T(r < r500c). The colors of the gray data points represent the average masses (M500c) of the groups within the corresponding Voronoi bin. Error-weighted averages of the data points are plotted as white circles, statistical uncertainties of the averages are shown with blue error bars, and the overall error budget of the average measurements resulting from the statistical and systematic uncertainties are shown with red error bars (see Sect. 4 for the details of the accounted systematics).

Fig. 14 S (r) – T relations at 0.15r500c, r2500c and r500c. The white circles represent redshift evolution scaled average entropy measurements, blue error bars represent statistical uncertainties of the averages, and the red error bars represent the overall error budget of the average measurements resulting from the statistical and systematic uncertainties (see Sect. 4 for the details of the accounted systematics). The best-fit power-law models to the data points are plotted as light green lines, and the uncertainties of the best-fit lines are shown with dark green shaded regions.

6 Comparison with the numerical simulations

In this section, we compare our results with the state-of-the-art simulations to place constraints on the physics of the AGN feedback. We quantify the agreement or disagreement between our measurements and the predictions of the various AGN feedback models implemented in these cosmological hydrodynamical simulations, including MillenniumTNG¹² (Hernández-Aguayo et al. 2023; Pakmor et al. 2023), Magneticum¹³ (Hirschmann et al. 2014) and the Overwhelmingly Large Simulations (OWL simulations, Schaye et al. 2010).

We compare our results with the Magneticum, MillenniumTNG, and OWL simulations that include different implementations for AGN feedback. The entropy measurements of the OWL simulations are plotted in Fig. 2 of McCarthy et al. (2010); therefore, we extracted their measurements from their paper and directly used them in this work. In contrast, the entropy and the characteristic temperature profile measurements for MillenniumTNG and Magneticum simulations are not publicly available. For this reason, we extract thermodynamic profiles of the gaseous halos from the MillenniumTNG and Magneticum simulations. A brief description of the MillenniumTNG and Magneticum simulations and the extraction process is as follows.

6.1 MillenniumTNG simulations

The MTNG740 flagship full physics run of the MillenniumTNG project (Pakmor et al. 2023) is used in this work, and it simulates a 500 Mpc/h cosmological box with the Planck Collaboration XIII (2016) cosmology. Its galaxy formation model is close to the IllustrisTNG model (Weinberger et al. 2017; Pillepich et al. 2018) and includes primordial and metal line cooling, a sub-grid model for star formation and the interstellar medium, mass return from stars via AGB stars and supernovae, an effective model for galactic winds, as well as a model for the formation, growth, and feedback from supermassive black holes. At a baryonic mass resolution of 3 × 10⁷ M_⊙ MTNG740 reproduces well many properties of observed galaxies and galaxy clusters for halos with M_500c > 2.3 × 10¹² M_⊙.

The hydrodynamical profiles of groups shown in this work are computed the same way as the galaxy cluster profiles in Pakmor et al. (2023). MTNG outputs a number of useful quantities: gas cell mass (m), electron abundance (x), and internal energy (ϵ). Assuming a primordial hydrogen mass fraction of X_H = 0.76, and an adiabatic index of γ = 5/3, we compute the volume-weighted electron number density, n_e, and temperature T, for each gas particle i as $$V_i{n}_{\text{e},i}=x_i{m}_i\frac{X_H}{m_p}$$(11) $$T_i=(\gamma -1){ϵ}_i/k\frac{4m_p}{1+3X_H+4X_H{x}_i},$$(12)

where m_p is the proton mass and k is the Boltzmann constant. We then computed the averages of the two quantities in radial bins b̂_j: $$n_{\text{e},j}=V^{-1}\left({\widehat{b}}_j\right){\displaystyle \sum _{i\in {\widehat{b}}_j}{V}_i}{n}_{\text{e},i}.$$(13) $$T_j=m^{-1}\left({\widehat{b}}_j\right){\displaystyle \sum _{i\in {\widehat{b}}_j}{m}_i}{T}_i.$$(14)

The radial bins were chosen as 23 logarithmically spaced intervals between 0.001 and 10 r_200c, where r_200c is the radius from the halo center encompassing 200 times the critical density of the Universe, ρ_c. As a derived quantity, we also computed the entropy as $$S\left({\widehat{b}}_j\right)=\frac{T\left({\widehat{b}}_j\right)}{n_{\text{e}}{\left({\widehat{b}}_j\right)}^{2/3}}.$$(15)

6.2 Magneticum simulations

The Magneticum simulations used in this work were performed with the TreePM/SPH code P-Gadget3, an extended version of P-Gadget2 (Springel 2005). In addition to hydrodynamics and gravity, the simulations also account for a variety of physical baryonic processes, including radiative cooling and heating from a time-dependent UV background (Haardt & Madau 2001), star formation and feedback (Springel & Hernquist 2003), metal enrichment from stellar evolution (Tornatore et al. 2004, 2007), and black hole growth and gas accretion powering energy feedback from AGN (Springel et al. 2005; Di Matteo et al. 2005; Fabjan et al. 2010). In this study, we consider a sample of 22,254 galaxy groups and clusters with mass M_500c > 5 × 10¹² M_⊙ identified using the SubFind algorithm (Springel et al. 2001; Dolag et al. 2009) in the "Box2_hr" simulation box. This covers a comoving volume of (352h⁻¹cMpc)³ and is resolved with 2 × 1584³ particles (which correspond to mass resolutions of m_DM = 6.9 × 10⁸ h⁻¹ M_⊙ and m_gas = 1.4 × 10⁸ h⁻¹ M_⊙, for DM and gas respectively). A ΛCDM cosmology with h = 0.704, Ω_b = 0.0451, Ω_m = 0.272, Ω_Λ = 0.728 and σ₈ = 0.809 is employed (from the 7-year results of the Wilkinson Microwave Anisotropy Probe; Komatsu et al. 2011). See also Biffi et al. (2022) for further details.

We compute gas thermodynamical profiles by selecting the hot, diffuse gas component in each halo of the sample, characterized by temperature higher than 5 × 10⁵ K and non-star-forming, representing the X-ray emitting intra cluster and group mediums. The gas mass-weighted temperature and electron density radial profiles are then generated by considering three-dimensional radial bins. The profiles of massive halos are resolved with 50 linear radial bins up to 1.5 r_500c. For smaller systems resolved with fewer gas particles, we instead adopt an equal-particle binning of the radial profiles, to ensure statistically reliable estimates. In all cases, we ensure a minimum of 150 selected gas particles in each radial shell.

6.3 Overwhelmingly Large Simulations

The Overwhelmingly Large Simulations (Schaye et al. 2010) used in this work is a cosmological hydrodynamical simulation that was performed using an extended version of the ThreePM/SPH code Gadget3 (Springel 2005). The simulations include radiative cooling (Wiersma et al. 2009a), star formation (Schaye & Dalla Vecchia 2008), stellar evolution and chemical enrichment (Wiersma et al. 2009b), kinetic supernovae feedback (Dalla Vecchia & Schaye 2008) and feedback from supermassive black holes (Springel et al. 2005; Booth & Schaye 2009). The simulations employ a flat ACDM cosmology with h = 0.73, Ω_b = 0.0418, Ω_m = 0.238, Ω_Λ = 0.762 and σ₈ = 0.74 that are taken from the analysis of 3-year results of WMAP (Spergel et al. 2007). The entropy and temperature measurements of galaxy groups in OWL simulations are taken from McCarthy et al. (2010) where ~200 galaxy groups are selected by applying a mass criterion of M_200c > 10¹³ M_⊙. The temperature and electron density profiles are obtained by emission-weighting the gas properties, and the entropy measurements are obtained by combining the temperature and electron density measurements.

They ran their simulations two times. In their first run (AGN run), they included all the sub-grid physics listed above, and in their second run (REF run), they turned off the AGN feedback. Both runs are indistinguishable in all aspects, with the only difference being that the former includes feedback from super-massive black holes as prescribed in (Booth & Schaye 2009). We refer the reader to Schaye et al. (2010) and McCarthy et al. (2010) for further details on the OWL simulations and the measurements of the thermodynamic properties of the galaxy groups in OWLS.

6.4 Comparisons of observations with OWL, MillenniumTNG, and Magneticum simulations

A fair comparison between simulations and observations is only possible when the compared samples are subject to a similar selection and the simulated temperatures are weighted with a well-calibrated temperature weighting scheme that would reproduce the observed spectroscopic-like temperatures. We achieve this by using the eROSITA selection function that is obtained by applying the same cleaning procedure described in Sect. 2 to the mock catalogs from the eROSITA’s digital twin (Seppi et al. 2022). The eRASS1 selection function, encapsulating the selection and cleaning information, provides detection probabilities of the group- and cluster-scale dark matter halos as a function of M_500c and ɀ (see Clerc et al. 2024, for details). The selection function obtained through this procedure is applied to the MillenniumTNG and Magneticum clusters and groups, and the entropy profiles of the simulated samples are measured at the three characteristic radii. Furthermore, the spectroscopic-like characteristic temperatures, T(r < r_500c), are obtained by weighting and projecting the simulated 3D temperature profiles of MillenniumTNG and Magneticum halos, using the temperature weighting scheme calibrated for eROSITA (ZuHone et al. 2023). Lastly, we compare the entropy and characteristic temperature measurements of simulated groups with the observations.

Before comparing our results with the simulations employing different AGN feedback implementations, we first test the reference (REF) run of the OWL simulations of McCarthy et al. (2010) with the observations. The predictions of the REF run serve as a baseline simulation and can be used to test the overall impact of the AGN feedback on the entropy profiles and characteristic temperatures of groups (Schaye et al. 2010). Their REF run is performed with their full model, identical to their final run, but does not include AGN feedback. To have a fair comparison, we convert the core-excised characteristic temperature measurements in simulations (T(0.15r_500c < r < r_500c)) to the core-included temperatures, T(r < r500c), following the same approach presented in Sect. 5. After the correction, we compare our measurements and the predictions of the REF run of the OWL simulations at the three characteristic radii, as shown in Fig. 15. We find that the main visible effect of AGN feedback is to lift the entropy values of the groups outside the core (e.g., at r_2500c or r_500c) for all groups, with a more marked effect at the lowest temperatures. The excess entropy induced by the AGN feedback in observations leads to significant disagreement between our measurements and the predictions of their REF run outside the core. On the other hand, the consistency with the AGN run (the twin simulations of the REF run, but the AGN feedback is turned on) suggests that the observations significantly favor the presence of strong AGN feedback in galaxy groups. It can further be seen from Fig. 15 that at the core 0.15r_500c, the data points of McCarthy et al. (2010) for the AGN feedback run are indistinguishable from data points of the reference run, and their results agree well with ours at this radius. Contrarily, outside the core (at r_2500c or r_500c), the reference and the feedback runs are significantly different at the low-temperature parameter space, indicating the strong AGN feedback imprint. This suggests that the entropy measurements at the core probe the AGN feedback much less efficiently than those outside the core at the group scale halos, which is also shown in McCarthy et al. (2010) for core-excised quantities. We note that the agreement at the core of groups between the reference and AGN feedback runs are obtained from the OWL simulations, where a thermal AGN feedback model is implemented. It is unclear whether such a conclusion would hold true for simulations implementing a feedback model that includes both kinetic and thermal feedback that leads to significant turbulent stirring in the core, such as MillenniumTNG. Therefore, to verify the universality of the conclusion, further tests are required in the future on simulations implementing kinetic and thermal feedback.

After comparing our measurements with the (REF) run of the OWL simulations, we compare our results with the three simulations employing different AGN feedback implementations: Magneticum, MillenniumTNG, and OWL (AGN run). The comparison between our measurements and the predictions of the simulations at the three characteristic radii is shown in Fig. 16.

At the cores of groups (0.15r_500c), we find the eROSITA observations agreeing relatively well with the OWL simulations between the IGrM temperatures T = 0.73–1.79 keV whereas at the warmer IGrM temperatures T = 1.79–2.68 keV the observations fall above their predictions. We further find that the Magneticum and MillenniumTNG simulations overpredict the average entropy for the cool/lower-mass groups (T = 0.73–1.44 keV), whereas the agreement becomes better in the MillenniumTNG at the warmer IGrM temperatures, close to the cluster ICM temperature range. Even though our measurements disagree with the average entropy predictions of Magneticum and MillenniumTNG, our measurements still lie within the 1σ scatter of the simulated profiles.

At the mid-region of groups (r_2500c), the entropy measured by eROSITA agrees well with the Magneticum and OWL simulation for cool/lower-mass groups, whereas the MillenniumTNG simulations underpredict the entropy for the groups in this region at a ~2.5σ level. We find that our measurements also agree well with the Magneticum simulations for warm/higher-mass groups. Furthermore, we find that as the temperatures/masses of the groups increase (T > 1.44 keV) and approach the cluster temperatures/masses, the offset between observations and the predictions of MillenniumTNG and OWL simulations becomes more significant, starting from a statistical disagreement level of 3.5σ at T = 1.44 keV and going up to a level of 8.5σ at T = 2.68 keV.

At the group outskirts (r_500c), we observe an overall agreement between the eROSITA observations and the Magneticum simulations at all temperatures. While the entropy measurements of galaxy groups in MillenniumTNG simulations are slightly below observations, we find that they are consistent at a 2σ confidence level at all temperatures. Although our measurements are in 1σ agreement with the OWL simulations for cooler/lower-mass groups T < 1.44 keV, the departure from observations becomes more significant (>2σ) at warmer temperatures, close to the ICM regime.

The AGN feedback implementations in most large-scale simulations are variations of the (Springel et al. 2005; Di Matteo et al. 2005) model. Even though similar in spirit, different adaptions and extensions to the original model have been made for the different simulations, leading to a strongly varying impact on the hot gas among the three simulations covered in this work. While the main feature behind the model for the OWL simulation is to achieve an effective AGN feedback by accumulating the injected energy until a ∆T_heat is reached, Magneticum and MillenniumTNG introduce a transition from quasar-mode to a stronger radio mode feedback (e.g., see Sijacki et al. 2007; Fabjan et al. 2010) instead. However, while in Magneticum the radio mode feedback is still modeled as a isotropic, thermal feedback with an increased efficiency (e.g., see Hirschmann et al. 2014), in MillenniumTNG it is modeled as a kinetic feedback (e.g., see Pillepich et al. 2018; Weinberger et al. 2017). In addition to these differences, the detailed setting within the individual models cannot only significantly change the distribution and properties of gas within galaxies, but also impact the hot gas and especially the entropy within and around clusters and groups (e.g., see figure 3 in Fabjan et al. 2010). Subsequently, in an updated version of the OWL simulations (Le Brun et al. 2014), three different ∆T_heat parameters (10⁸, 10^8.5 and 10^8.7 keV) of the Booth & Schaye (2009) model, coupled with the burstiness and the energeticness of the feedback, are tested. The authors report that ∆T_heat can be used to tune the impact of the feedback in the IGrM gas. They further show that for groups, the normalization and slope of the S − T relation (or the S − M relation) does not depend much on ∆T_heat at the core (0.15r_500c). In comparison, an increase in the ∆T_heat parameter corresponds to a larger normalization and possibly a different slope at the outskirts (r_25ooc and r_5ooc) due to the low-entropy gas getting ejected (see Fig. 6 in Le Brun et al. 2014). This can also be seen in Fig. 7 of Le Brun et al. (2014) where the gas gets ejected in the case of more bursty and energetic feedback and results in lower densities at all radii out to 1.6r_500c.

In light of these indications found within simulations, the discrepancy between the observed entropy and the predictions of MillenniumTNG and OWL simulations at the outskirts can be justified by the effectiveness and energetics of the individual AGN feedback models being not sufficient to match our observations at the group scale. The weaker AGN feedback found in MillenniumTNG, leading to lower entropy in the outskirts, may also be related to the elimination of the magnetic fields from the TNG physics model, as noted in Pakmor et al. (2023). The agreement between the observations and the Magneticum simulations in the group outskirts suggest that the underlying model is able to effectively treat a broad range in mass, extending the previously reported good agreement of entropy profiles for galaxy clusters (Planelles et al. 2014; Rasia et al. 2015) to group scales.

The agreement between the entropy measurements in eROSITA observations and the numerical simulations at the cores is better than the outskirts. Slightly higher entropy in simulations is harder to explain using the Booth & Schaye (2009) model since the normalization and the shape at the core (0.15r_500c) do not seem to be affected by the change of the ∆T_heat parameter in the runs presented in Le Brun et al. (2014).

We note that the interpretations above for the OWL simulation assume that the Booth & Schaye (2009) model can reproduce the gas properties by tuning the ∆T_heat parameter relatively accurately, which may or may not be the case. Furthermore, as also discussed above, the AGN feedback implementation in MillenniumTNG and Magneticum simulations are even different (among them as well as compared to OWL). For these reasons, a more detailed analysis would be needed to better understand which aspects of the underlying AGN feedback models are driving the differences between the simulations and the eROSITA observations, which is beyond the scope of this paper. Furthermore, we note that the selection function we used in this work is a function of mass and redshift, and the use of it significantly improves the robustness of our comparison. In future studies, a profile-dependent selection function can be used that would make the selection procedure applied to the simulations even more realistic. Ultimately, the most accurate comparison would be achieved by producing synthetic eROSITA observations through fully forward modeling cosmo-logical hydrodynamic simulations to which the same detection and data analysis pipeline used for real observations can be applied.

Fig. 15 Comparison between the redshiſt evolution scaled average eROSITA entropy measurements (white circles) with the predictions of the REF (AGN feedback off) and AGN (AGN feedback on) runs of OWL simulations (yellow triangles and green squares) presented in McCarthy et al. (2010) at three radii (0.15r500c, r2500c, r500c) as a function of characteristic temperature, T(r < r500c). Blue error bars represent the statistical uncertainties of the average measurements, and the red error bars represent the overall error budget resulting from the statistical and systematic uncertainties (see Sect. 4 for the details of the accounted systematics). For consistency, core-excised temperatures presented in McCarthy et al. (2010) are converted to core-included temperatures (see Sect. 6 for the details of the correction).

Fig. 16 Comparison of the redshift evolution scaled average eROSITA entropy measurements (yellow circles) with the predictions of three simulations (Magneticum, MillenniumTNG, and OWL simulations) at 0.15r5ooc (top left), r25ooc (top right), and r5ooc (bottom) as a function of characteristic temperature, T(r < r500c). Gray crosses represent the average entropy measurements of the binned groups, and the blue error bars represent the overall error budget of the average measurements resulting from the statistical and systematic uncertainties (see Sect. 4 for the details of the accounted systematics). The solid lines represent the predictions of simulations at the group regime (5 × 1012 < M500c < 1014 M⊙), the dashed lines represent the predictions at the low-mass cluster regime (1014 < M50oc < 3 × 1014 M⊙) and the shaded regions represent the scatter of the measurements.

7 Conclusions

Our work places the tightest constraints on the impact of AGN feedback on the average thermodynamic properties by populating the low-mass galaxy groups down to cool IGrM temperatures of 0.7 keV. eROSITA’s superb sensitivity in the soft X-ray band led to the detection of a large number of galaxy groups with a well-understood selection function. When stacked, the eROSITA data provides unprecedented statistical power for the measurements of X-ray properties of these sources. We used a sample of 1178 galaxy groups to place constraints on the impact of AGN feedback on the thermodynamic properties of the IGrM. We selected the galaxy group sample based on the primary sample of the eRASS1 clusters and groups and applied a rigorous selection and cleaning. The cleaning procedure was designed to provide a pure sample with a well-defined selection function while maximizing the sample size. A Bayesian imaging analysis was carried out for all the 1178 galaxy groups in the sample, where the nearby clusters and bright point sources were co-fit. The galaxy groups with similar statistical and physical properties, such as count and temperature, were then grouped together into 271 bins. A joint Bayesian spectral fitting was performed on the groups in the same bin to increase the statistical power in each bin with the sources with similar properties.

Constraining baryonic physics at group scales is a highly challenging task. Systematic effects must be considered to achieve a reliable and robust conclusion. We have quantified and discussed three major systematics for thermodynamic profile studies and conservatively took them into account for our final measurements and conclusions. Apart from the robust thermodynamic property measurements, we have also provided a detailed comparison of our findings with the state-of-the-art simulations by accounting for the selection effects. We assessed the agreement between our measurements and the simulations, employing various AGN feedback implementations to pave the way for more realistic AGN feedback modeling in numerical simulations.

This paper is the first in-depth study of galaxy groups with eROSITA focusing on the impact of AGN feedback on the entropy and temperature measurements within an overdensity radius. The main conclusions of this work are as follows:

• With a sample of 1178 galaxy groups and average thermodynamic properties of 271 binned groups, our work stands as the most comprehensive study of the hot gas in galaxy groups in terms of sample size, diversity, and statistics. The selection effects have been considered for the first time while comparing the measured thermodynamic profiles of galaxy groups with the numerical simulations;

• Overall, the entropy measurements at three characteristic radii, 0.15r_500c, r_2500c, and r_500c, and the characteristic temperature of the galaxy groups detected in the eROSITA first All-Sky Survey observations are in good agreement with the previously reported results in S09 and Johnson et al. (2009), within the uncertainties (see Sect. 4). The largest mismatch between the eROSITA measurements presented in this work and the Chandra measurements presented in S09 is at the low-temperature parameter space (T(r < r_500c) < 1.15 keV). We argue that this is because the completeness and the selection of the two samples is different in this temperature range. The archival sample used in S09 consists mostly of bright and morphologically undisturbed galaxy groups detected in ROSAT observations and was followed up by Chandra. In contrast, the group sample used in this work has a more uniform and well-defined selection;

• We compared our entropy measurements with the reference (REF) run of OWL simulations that include various non-gravitational processes but not the AGN feedback. From this comparison, we conclude that the impact of AGN feedback on the entropy profiles of groups is significant at the out-skirts (r_2500c and r_500c) and less pronounced near the core (0.15r_500c). This result suggests that for groups, AGN feedback has a larger impact on the outskirts than the core, contrary to its impact on the observed gas properties of clusters of galaxies (Le Brun et al. 2014). Due to their shallower potential wells, the feedback from the central black hole is able to move the low-entropy gas to much larger radii (e.g., to the outskirts) of galaxy groups;

• Our measurements have significant constraining power on the impact of AGN feedback on the thermodynamic properties of the IGrM gas and serve as a reference for the feedback implementations in numerical simulations and theoretical models. We compared our results with three state-of-the-art cosmological hydrodynamic simulations, including MillenniumTNG, Magneticum, and OWL, employing various AGN feedback implementations by accounting for the selection effects. In the cores of the galaxy groups, the entropy agrees well with OWL simulations out to T = 1.79 keV. Although the sample-averaged entropy from the MillenniumTNG and Magneticum simulations are higher, the measurements are within the scatter of the respective simulations. At the out-skirts, (r_2500c and r_500c), the observed entropy agrees well with the Magneticum simulations for groups with IGrM temperatures down to 0.79 keV. In OWL simulations, although we observed a similar entropy-flattening trend for the cooler groups, the departure from observations becomes significant toward the galaxy cluster regime at higher temperatures. A similar trend has been observed in MillenniumTNG simulations; however, in this case, the entropy offset is relatively significant for all temperature regimes. Overall, the AGN feedback implementation in Magneticum simulations reproduces our observations the best at the three characteristic radii.

This work demonstrates the potential of eROSITA for exploring the baryonic physics at the galaxy groups out to large radii. Deeper data with the eROSITA All-sky Survey will allow for the detection of a larger sample of galaxy groups pushing down in mass and IGrM temperature floors. Employing these groups in similar future studies with larger statistical power will enable the testing of hydrodynamical simulations in very cool temperatures below <0.5 keV, which are currently unreachable with the eRASS1 sample.

Acknowledgement

This work is based on data from eROSITA, the soft X-ray instrument aboard SRG, a joint Russian-German science mission supported by the Russian Space Agency (Roskosmos), in the interests of the Russian Academy of Sciences represented by its Space Research Institute (IKI), and the Deutsches Zentrum für Luft und Raumfahrt (DLR). The SRG spacecraft was built by Lavochkin Association (NPOL) and its subcontractors and is operated by NPOL with support from the Max Planck Institute for Extraterrestrial Physics (MPE). The development and construction of the eROSITA X-ray instrument were led by MPE, with contributions from the Dr. Karl Remeis Observatory Bamberg & ECAP (FAU Erlangen-Nuernberg), the University of Hamburg Observatory, the Leibniz Institute for Astrophysics Potsdam (AIP), and the Institute for Astronomy and Astrophysics of the University of Tübingen, with the support of DLR and the Max Planck Society. The Argelander Institute for Astronomy of the University of Bonn and the Ludwig Maximilians Universität Munich also participated in the science preparation for eROSITA. The eROSITA data shown here were processed using the eSASS/NRTA software system developed by the German eROSITA consortium. E. Bulbul, A. Liu, V. Ghirardini, C. Garrel, S. Zelmer, and X. Zhang acknowledge financial support from the European Research Council (ERC) Consolidator Grant under the European Union’s Horizon 2020 research and innovation program (grant agreement CoG Dark-Quest No 101002585). N. Clerc was financially supported by CNES. S. Bose is supported by the UK Research and Innovation (UKRI) Future Leaders Fellowship [grant number MR/V023381/1]. C. Hernández-Aguayo acknowledges support from the Excellence Cluster ORIGINS which is funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy – EXC-2094 – 390783311. This work made use of the following Python software packages: Astropy¹⁴ (Astropy Collaboration 2022), Colossus¹⁵ (Diemer 2018), emcee¹⁶ (Foreman-Mackey et al. 2013), Matplotlib¹⁷ (Hunter 2007), NumPy¹⁸ (Harris et al. 2020), pyfftw¹⁹ (Frigo & Johnson 2005) PyXspec²⁰ (Arnaud 1996), SciPy²¹ (Virtanen et al. 2020), Vorbin²² (Cappellari & Copin 2003),

Appendix A Energy band selection for the imaging analysis

We chose the energy band for the imaging analysis by optimizing the signal-to-noise ratio given multiple X-ray foreground and background components. The models we used during this optimization scheme include 1) a 0.1 keV unabsorbed APEC component for the Local Hot Bubble, whose flux is scaled to 2.3 × 10⁻¹³ erg s⁻¹ cm⁻² deg⁻² in the 0.3–0.7 keV band (Yeung et al. 2023); 2) a 0.18 keV absorbed APEC component for the Galactic Halo, whose absorbed flux is scaled to 1.1 × 10⁻¹² erg s⁻¹ cm⁻² deg⁻² in the 0.5–2.0 keV band (Henley & Shelton 2013); 3) a Γ = 1.4 power law component for the Cosmic X-ray Background component, whose absorbed flux is scaled to 6 × 10⁻¹³ erg s⁻¹ cm⁻² deg⁻² in the 0.5–2.0 keV band, which corresponds to a 10⁻¹³ erg s⁻¹ cm⁻² point source flux cut; and 4) the instrumental background from eROSITA filter-wheel-closed observations. The foreground HI column density is fixed to 2.7 × 10²⁰ cm⁻² for this study, which is the averaged value of the sample. With these foreground and background configuration, we find that the signal-to-noise ratio of a T = 1 keV APEC source component at the redshift of 0.18 reaches the maximum in the 0.3–1.8 keV band at a source-to-background ratio of 1. We note that the upper boundary of the optimal band could increase if we adopt a higher source-to-background ratio since the log N − log S curve of the selected sample is a power low. Nevertheless, we select a source-to-background ratio of 1 and we adopt the 0.3–1.8 keV band for our imaging analysis to maintain high signal-to-noise for the faint sources.

References

All Tables

All Figures

	Fig. 1 Left: mass and redshift distributions of the galaxy group sample used in this work consisting of 1178 objects where colors of the data points represent the eROSITA counts of the groups in soft X-ray band (0.3-1.8 keV) Right: mass and redshift histograms of the group sample where the median redshift (z) is 0.11 and median mass (M500c) is 6.5 × 1013 M⊙.
In the text

Fig. 2 Projected locations of the 1178 galaxy groups in the primary catalog in the eROSITA and Legacy Survey DR9N and DR10 13,116 deg2 common footprints. The redshift confirmed by the follow-up algorithm eROMaPPer is color coded (Kluge et al. 2024), while the sizes of the detections are scaled with the angular sizes (r500c) of the groups (see Sect. 3.3 for the r500c estimation procedure). The inhomogeneity of the source density in this figure is due to the exposure variation across the eROSITA-DE X-ray sky (see Fig. 2 in Bulbul et al. 2024).

In the text

	Fig. 3 Left: eRASS:4 soft band 0.3–1.8 keV images of two bright groups (1eRASS J024933.9-311126 and 1eRASS J045547.0-572404) at redshifts 0.023 and 0.105. Right: legacy Survey DR10 images of the same groups with the eRASS:4 X-ray contours overlayed.
In the text

Fig. 4 Example of the eROSITA imaging analysis for the group 1eRASS J045547.0-572404. Left: residual image of the group before and after subtracting the co-fit extended sources and point sources in the field. Nearby co-fit clusters and groups are shown with red circles, co-fit nearby bright AGN are shown with green circles, and the fitted galaxy group is shown with a blue rectangle. The smooth noise level indicates that the contaminant emission is modeled properly in the analysis. Right: surface brightness profile of the same group. The observed surface brightness profile of the image is plotted in blue, the best-fit model of the image is plotted in red, the PSF profile over the measured background is plotted in green, the PSF deconvolved surface brightness profile of the galaxy group is plotted in orange, the measured background level is shown with a horizontal dashed red line, and three characteristic radii of the group (0.15r 500c, r2500c and r500c from the core to the outskirt respectively) are shown with dashed gray lines.

In the text

Fig. 5 Left: soft-band (0.5–2 keV) X-ray luminosity (LX,500c) and redshift (z) distributions of the galaxy group sample used in this work consisting of 1178 objects. The luminosity of the groups cover a range of 3.9 ×1040–1.4 × 1043 ergs s−1, and the redshift span of the sample is 0.003–0.48. Right: count (C500c) and mass (M500c) distributions of the 1178 galaxy group used in this work. Measured counts of galaxy groups range between 10–14 380, and their masses span a range of 2.1 × 1012–1014 M⊙. Median values of the LX,500c, z, C500c and M500c observables are 6.3 × 1042 ergs s–1, 0.11, 95, 6.5 × 1013 M⊙ respectively.

In the text

	Fig. 6 Left: Voronoi binning scheme used for grouping the sample obtained from the distribution of count (within r500c) measurements (C500c) and scaling relation based temperature estimates (T500c,sc). Top right: histogram of total counts in 271 Voronoi bins with a median of 905 counts. Bottom right: Histogram of galaxy groups in Voronoi bins.
In the text

	Fig. 7 Normalized temperature profiles, T(r)/T(r < r500c), of 23 groups (gray) presented in S09 and their average (blue). The blue dashed line provides an average conversion ratio between the temperature profile, T(r), and the characteristic temperature measurements, T(r < r500c).
In the text

	Fig. 8 Average metalicity profiles of galaxy groups reported by Mernier et al. (2017) and Lovisari et al. (2015) and the stacked metallicity profiles of Sun (2012) for galaxy groups in two temperature bins (T = 0.75–1.3 keV and 0.75–1.3 keV).
In the text

Fig. 9 Ratio of the electron densities obtained by assuming low/up-scattered Mernier et al. (2017) metallicity profiles to the electron densities obtained by the median Mernier et al. (2017) metallicity profile at three radii, 0.15r500c, r2500c, and r500c as a function of characteristic temperature, T(r < r500c). Green, blue, and purple data points represent the ratio between the electron densities obtained by assuming the lower envelope (ne,L) of the red shared area and the dark red median line (ne) in Fig. 8. Orange, yellow, and red data points represent the ratio between the electron densities obtained by assuming the upper envelope of the red shared area (ne,U) and the dark red median line (ne) in Fig. 8.

In the text

Fig. 10 Comparison between the redshift evolution scaled average eROSITA entropy measurements (white diamonds) with the Chandra measurements of 43 galaxy groups (black error bars) presented in S09 at three radii (0.15r500c, r2500c, r500c) as a function of characteristic temperature, T(r < r500c). Blue error bars represent the statistical uncertainties of the average measurements, and the red error bars represent the overall error budget resulting from the statistical and systematic uncertainties (see Sect. 4 for the details of the accounted systematics). For consistency, coreexcised temperatures presented in S09 are converted to core-included temperatures, and the entropy measurements presented in S09 are normalized such that the new data points are measured at the same angular radii with eROSITA flux calibration (see Sect. 5 for the details of these corrections).

In the text

	Fig. 11 Comparison between the scaling relation based r500c estimates of the crossmatched galaxy groups obtained in this work (y-axis: r500c) with the estimates of S09 (x-axis: r500c,S09) obtained by assuming hydrostatic equilibrium.
In the text

Fig. 12 Average entropy measurements of binned groups (diamonds) at three radii (0.15 r500c, r2500c, r500c) as a function of characteristic temperature, T(r < r500c). The colors of the diamond data points represent the average concentration ($c_{\text{SB},r_{500c}}=SB\left(r<0.1r_{500c}\right)/SB\left(r<r_{500c}\right)$) of the groups within the corresponding Voronoi bin. Error-weighted averages of the diamond data points are shown in magenta.

In the text

Fig. 13 Redshift evolution scaled average entropy measurements of binned groups (gray diamonds) at three radii (0.15r500c, r2500c, r500c) as a function of characteristic temperature, T(r < r500c). The colors of the gray data points represent the average masses (M500c) of the groups within the corresponding Voronoi bin. Error-weighted averages of the data points are plotted as white circles, statistical uncertainties of the averages are shown with blue error bars, and the overall error budget of the average measurements resulting from the statistical and systematic uncertainties are shown with red error bars (see Sect. 4 for the details of the accounted systematics).

In the text

Fig. 14 S (r) – T relations at 0.15r500c, r2500c and r500c. The white circles represent redshift evolution scaled average entropy measurements, blue error bars represent statistical uncertainties of the averages, and the red error bars represent the overall error budget of the average measurements resulting from the statistical and systematic uncertainties (see Sect. 4 for the details of the accounted systematics). The best-fit power-law models to the data points are plotted as light green lines, and the uncertainties of the best-fit lines are shown with dark green shaded regions.

In the text

Fig. 15 Comparison between the redshiſt evolution scaled average eROSITA entropy measurements (white circles) with the predictions of the REF (AGN feedback off) and AGN (AGN feedback on) runs of OWL simulations (yellow triangles and green squares) presented in McCarthy et al. (2010) at three radii (0.15r500c, r2500c, r500c) as a function of characteristic temperature, T(r < r500c). Blue error bars represent the statistical uncertainties of the average measurements, and the red error bars represent the overall error budget resulting from the statistical and systematic uncertainties (see Sect. 4 for the details of the accounted systematics). For consistency, core-excised temperatures presented in McCarthy et al. (2010) are converted to core-included temperatures (see Sect. 6 for the details of the correction).

In the text

Fig. 16 Comparison of the redshift evolution scaled average eROSITA entropy measurements (yellow circles) with the predictions of three simulations (Magneticum, MillenniumTNG, and OWL simulations) at 0.15r5ooc (top left), r25ooc (top right), and r5ooc (bottom) as a function of characteristic temperature, T(r < r500c). Gray crosses represent the average entropy measurements of the binned groups, and the blue error bars represent the overall error budget of the average measurements resulting from the statistical and systematic uncertainties (see Sect. 4 for the details of the accounted systematics). The solid lines represent the predictions of simulations at the group regime (5 × 1012 < M500c < 1014 M⊙), the dashed lines represent the predictions at the low-mass cluster regime (1014 < M50oc < 3 × 1014 M⊙) and the shaded regions represent the scatter of the measurements.

In the text