X-ray surface brightness and gas density profiles of galaxy clusters up to 3 × R_500c with SRG/eROSITA

ABSTRACT

Using the data of the SRG/eROSITA all-sky survey, we stacked a sample of ∼40 galaxy cluster images in the 0.3–2.3 keV band, covering the radial range up to 10 × R_500c. The excess emission on top of the Galactic and extragalactic X-ray backgrounds and foregrounds is detected up to ∼3 × R_500c. At these distances, the surface brightness of the stacked image drops below ∼1 per cent of the background. The density profile reconstructed from the X-ray surface brightness profile agrees well (within ∼30 per cent) with the mean gas profile found in numerical simulations, which predict the local gas overdensity of ∼ 20–30 at 3 × R_500c and the gas fraction close to the universal value of |$\frac{\Omega _b}{\Omega _m}\approx 0.15$| in the standard Λ cold dark matter model. Taking at face value, this agreement suggests that up to ∼3 × R_500c the X-ray signal is not strongly boosted by the gas clumpiness, although a scenario with a moderately inhomogeneous gas cannot be excluded. A comparison of the derived gas density profile with the electron pressure profile based on the Sunyaev-Zeldovich (SZ) effect measurements suggests that by r ∼ 3 × R_500c the gas temperature drops by a factor of ∼ 4–5 below the characteristic temperature of a typical cluster in the sample within R_500c, while the entropy keeps growing up to this distance. Better constraints on the gas properties just beyond 3 × R_500c should be possible with a sample larger than used for this pilot study.

1 INTRODUCTION

Clusters of galaxies represent the high-mass end of the virialized halos in the present-day Universe (see Kravtsov & Borgani 2012, for a review). To the first approximation, their properties depend on their mass and redshift. In particular, their radial profiles of total and gas densities and gas temperature or entropy profiles, once properly scaled, are expected to be approximately universal (e.g. Markevitch et al. 1998; Voit 2005; Vikhlinin et al. 2006; Arnaud et al. 2010; Walker et al. 2012; Eckert et al. 2012, 2013; Ghirardini et al. 2019, among others) although the details of the formation history of individual objects lead to inevitable deviations. From this perspective, averaged (stacked) profiles are useful since these deviations might be suppressed. The X-ray emission from massive clusters has been routinely detected up to R_500c¹ or even up to ≈1.5–2R_500c in individual cases (see e.g. Reiprich et al. 2013; Walker et al. 2019; Walker & Lau 2022 for reviews). The same is true for SZ signal from clusters (e.g. Plagge et al. 2010; Planck Collaboration V 2013). Combining the data from many clusters helps to extend the radial range probed. To this end, the data of the SRG/eROSITA all-sky survey (Predehl et al. 2020, 2021; Sunyaev et al. 2021) are especially valuable since they provide uniform coverage of the regions around clusters with the same sensitivity and spatial resolution. Here we perform a pilot study of stacking a small sample of clusters using eROSITA data.

Throughout the paper, we adopted the Λ cold dark matter cosmology with Ω_M = 0.3, Ω_Λ = 0.7, and H₀ = 70 km s⁻¹ Mpc⁻¹.

2 SAMPLE SELECTION

For the purpose of this study, we need a sample of clusters with known redshifts and reliable mass estimates. We further would like to suppress spatial variations of the X-ray sky backgrounds and foregrounds to get an accurate estimate of the excess emission associated with clusters, i.e. the stacked image should not be dominated by a few very bright objects. Simultaneously, we want the clusters to be spatially well-resolved and sufficiently bright. A combination of these requirements calls for a sample of massive nearby clusters, from which a few nearest objects are excluded.

The galaxy clusters from the CHEX-MATE catalogue (namely, Tier-1; CHEX-MATE Collaboration 2021; Campitiello et al. 2022) appear to meet all our requirements. The CHEX-MATE is an unbiased, signal-to-noise limited sample of 118 galaxy clusters detected by Planck via Sunyaev–Zel’dovich effect (Sunyaev & Zeldovich 1972). It is composed of clusters at 0.05 < |$z$| < 0.2 with masses |$2\times 10^{14} \, \mathrm{M}_{\odot } \lt M_{500c} \lt 9\times 10^{14} \, \mathrm{M}_{\odot }$| from the PSZ2 catalogue (Planck Collaboration XXVII 2016). We restrict the sample to objects with a Galactic longitude 0° < l < 180° (available to us according to the data-sharing agreement) and a Galactic latitude |b| > 15° to avoid background variations due to proximity to the Galactic plane and the Cygnus X star formation region. For the same reason, we further exclude a few clusters located near the North Polar Spur. The final list of 38 clusters is presented in Table A1. The sample median mass is equal to |$M_{500c} = 4.13 \cdot 10^{14} \, \mathrm{M}_{\odot }$|⁠, the median redshift is |$z$| = 0.113, and the hydrogen column density is N_H = 2.17 × 10²⁰ cm⁻².

3 DATA AND DATA PROCESSING

Orbital observatory SRG (Sunyaev et al. 2021), featuring two focusing X-ray telescopes, eROSITA (Predehl et al. 2021, 0.3–8.0 keV) and ART-XC (Pavlinsky et al. 2021, 4–30 keV) was launched in 2019 July and started to perform the all-sky survey mission in 2019 December. By now, four complete all-sky surveys have been completed, resulting in unprecedentedly deep X-ray maps of the whole sky being obtained.

The eROSITA telescope consists of seven nearly identical independent modules, each having its own mirror system, a position-sensitive detector, and electronics (Predehl et al. 2021), and for the imaging analysis we use the data of all seven telescope modules (TMs). Here we use the data accumulated over four consecutive scans, with the total effective exposure amounting to 1225 s (i.e. ≈8600 s in equivalent exposure for one TM). Initial reduction and processing of the data were performed using standard routines of the esass software (Brunner et al. 2018; Predehl et al. 2021), while the imaging and spectral analysis were carried out with the background models, vignetting, point spread function (PSF), and spectral response function calibrations built upon the standard ones via slight modifications motivated by results of calibration and performance verification observations (e.g. Churazov et al. 2021; Khabibullin et al. 2023).

For each galaxy cluster from Table A1, we made an X-ray image in the 0.3–2.3 keV band of a region 10R_500c × 10R_500c centred on a cluster, where the values of R_500c and |$z$| taken from the original sample were used to set the angular size of the image. As outlined below, a few operations have been performed before adding the images to the stack.

3.1 Removal of sources

In the 0.3–2.3 keV band, the eROSITA detector’s intrinsic background makes up about ∼33 per cent of the total count rate, while distant active galactic nuclei, and the Galactic diffuse emission contribute ∼42 per cent and ∼25 per cent, respectively, in the areas with small N_H (namely, for the sample median N_H = 2.17 × 10²⁰ cm⁻²). Of course, these numbers fluctuate from field to field and are also subject to temporal variations of the detector background.

To reduce fluctuations in surface brightness due to bright sources, we removed compact sources with the flux exceeding |$3\times 10^{-14}\, {\rm erg\, cm^{-2}\, s^{-1} }$| in the standard 0.5–2 keV energy band. The procedure for detection and modelling of the point sources is identical to the one implemented and described in Churazov et al. (2021). Fluxes are converted from the detection band to the standard 0.5–2 keV energy band via a constant conversion factor, calculated for a typical cosmic X-ray background (CXB) spectrum. Although some sources might have different spectra, e.g. stars, they are not numerous at the flux limit of the constructed source catalogues. The shape and normalization of the logN–logS distributions of the obtained catalogues are fully consistent with the earlier measurements and models.

Similar procedure is applied to mildly extended sources [with the characteristic radius ≲ 3 arcmin, i.e. ≲10-times the characteristic radius for the β-profile model of the core PSF; see e.g. Churazov et al. (2021) for the detailed description], which are detected and modelled with the β-profile (Cavaliere & Fusco-Femiano 1976) of the surface brightness together with the point sources at the first stage. The limiting flux for the mildly extended sources depends on their size and the background level for the locations of the clusters in our sample, and a rather conservative significance limit corresponding to the flux of ≈10⁻¹³ erg s⁻¹ cm⁻² has been chosen to ensure its uniformity across the sample while keeping suppression of the shot noise from the individual sources at the required level.

After automatic detection, modelling and removal of the point and mildly extended sources, more extended sources, typically having sizes bigger than a few arcmin and irregular shapes, are visually identified on the smoothed residual images (as exemplified by the supernova remnant candidates found blindly in the X-ray survey data; see Khabibullin et al. 2023). Those are not numerous and were masked manually after a few iterations with varying mask sizes.

Fig. 1 illustrates the source removal procedure. Fig. 2 compares the stacked images with and without the removal of sources. The right-hand panel of Fig. 2 shows the major strength of the survey data when the intrinsic size of the telescope FoV is not imposing any constraints on the size of the studied region.

3.2 Stray light

A small fraction of incoming photons can reach eROSITA detectors after one scattering by the telescope mirrors instead of nominal two scatterings. These photons (usually called ‘stray light’) produce spurious extended halos (up to ∼3°) around each X-ray source (see e.g. fig. 15 in Sunyaev et al. 2021). To estimate the contribution of the stray light, we followed the approach described in appendix A of Churazov et al. (2023). Namely, every image has been convolved with a kernel that approximately characterizes the spatial distribution of stray light photons in the 0.4–2.3 keV band. The resulting convolved images are stacked similarly to cluster images and subtracted from them. This way a first-order correction for the stray light contribution is achieved.

3.3 Stacking

To do actual stacking, we re-mapped every X-ray image as if the cluster was located at D_ref = 343 Mpc away from us and had |$R_{\rm 500c} = 1\, {\rm Mpc}$| corresponding to the angular size of ≈10 arcmin. The choice of these parameters is rather arbitrary and motivated largely by the desire to have a simple conversion of radii into the units of R_500c or Mpc. The X-ray surface brightness was also adjusted to reflect the differences in distance, mass, redshift (see Appendix  B) and the area distortions due to the tangential (gnomonic) projection of individual images. Given that the redshift and mass dependencies are accounted for during the stacking procedure, the above reference distance D_ref is simply a coefficient that relates the physical and angular sizes in Euclidean space.

When doing stacking, we keep track of the photon counting noise and the estimated surface density of unresolved sources so that their contributions to the noise in the final images can be estimated. In parallel, the (properly scaled) exposure maps are also stacked in the same manner as the X-ray images.

4 RESULTS

Fig. 2 shows the final stacked image with a few characteristic radii indicated by white lines. Namely, we adopt the following relation between these radii R_ta:R_200m:R_200c:R_500c = 8.1:2.7:1.6:1 according to Nelson, Lau & Nagai (2014) and Diemer & Kravtsov (2014). By design, R_500c in the stacked image corresponds to 1 Mpc (or 10 arcmin at the adopted reference distance). Given that all images were originally made in equatorial orientation and then co-added, the final 2D image does not bear much information except for a demonstration that no prominent peculiar structures have been left in the image and the resulting image appears symmetric. Instead, most of the important information is contained in the radial surface brightness profile, which is related to the 3D emissivity profile. The latter profile is assumed to be spherically symmetric. However, 2D stacked images could be used to compare the expected and actual levels of noise by calculating the radial profile in several wedges as we do below.

4.1 The X-ray surface brightness and gas density profiles

To describe the observed surface brightness profile of the stacked image, one could use a modified β-model (e.g. Vikhlinin et al. 2006, 2009) to parametrize the radial dependence of the emission measure, i.e. the product of the electron and proton number densities

In cosmological simulations, dark matter and gas density profiles in the outer cluster regions show varying logarithmic slope (e.g. Diemer & Kravtsov 2014; O’Neil et al. 2021). This behaviour is not fully captured by formula (1), so we further modify it by introducing an additional multiplier

so that the resulting profile can change its slope near r_s and r_s2. In total, the model has 10 free parameters. While the values of these parameters are correlated, the equation (2) provides us with sufficient flexibility to describe the radial behaviour of the emission measure found in simulations. For an X-ray emissivity that weakly depends on temperature (see below), the expected surface brightness profile is described by the line-of-sight integral of the above expression.

We extract the surface brightness profile from the stacked X-ray image (Fig. 2) and then fit it with the projected model given by equation (2) plus a spatially constant background component. Fig. 3 shows the observed surface brightness profile together with the best-fitting model. The parameters of the model are given in Appendix  C. The uncertainty in the constant background level is estimated to be of the level of ≈1 per cent. This value was obtained by dividing the stacked image into four 90° wedges and calculating the standard deviation between best-fitting constants in those wedges.

The best-fitting model to the surface brightness in units of |${\rm counts\, s^{-1}\, arcmin^{-2}}$| has been converted into the proton density² profile as follows:

where θ = 2.9 × 10⁻⁴ corresponds to 1 arcmin and n_p/n_e = 0.83. The value |$\epsilon _c=5\times 10^{-14} \, {\rm cm^5 counts~s^{-1}}$| used in the above expression corresponds to a typical combination of cluster parameters in the sample (see Appendix  B and Fig. B1). The spread in the absorption column densities, mean temperatures, and redshifts among the clusters in the sample translates into variations of ϵ_c between 4 and 5. This translated into ∼11 per cent uncertainty in the recovered typical gas density. A larger spread could be present when the entire radial range is considered (including the cool core regions or the very outskirts). For this pilot study, we keep the assumption of a constant ϵ_c and plan to do a first-order correction using temperature and abundance profiles from numerical simulations in the subsequent work. A brief discussion on the expected level of uncertainties is given in Appendix  B.

We also note here that in the outer regions of a cluster, several other effects can impact the X-ray signal. These include non-equilibrium ionization (e.g. Yoshikawa & Sasaki 2006), photoionization (e.g. Churazov et al. 2001; Khabibullin & Churazov 2019; Štofanová et al. 2022), and resonant scattering of the line (e.g. Gilfanov, Syunyaev & Churazov 1987) and/or continuum photons. For the latter process, the difference in temperatures between the cluster core and outskirts implies that different ions are present in these two regions. Therefore, the photons from the continuum will be scattered in the cluster’s distant cluster outskirts. In this case, the expected signal is more than an order of magnitude below the surface brightness levels relevant to this study.

4.2 Unresolved CXB sources

Apart from pure photon-counting noise, additional variations in the derived surface brightness profiles could come from fluctuations of various backgrounds and foregrounds, and, also, from ‘peculiarities’ of individual clusters. Some of these contributions can be straightforwardly included in the error budget. For instance, this is true for Poisson fluctuations in the number of unresolved sources that compose CXB. The surface density of unresolved sources was estimated and the expected noise level was corrected to include these fluctuations (single-halo approximation). Note that clustering of X-ray sources (e.g. Vikhlinin & Forman 1995) is not accounted for here, and may increase noise associated with unresolved X-ray sources. To this end, the number counts of the Chandra deep fields were used (Luo et al. 2017). Given that we use the same threshold for compact sources in all fields, this type of noise scales with the solid angle of the studied region similarly to the photon counting noise, albeit its amplitude does not depend on the exposure time. Effectively, for this sample, the noise increases by a factor ∼1.9 compared to the pure photon counting noise.

The comparison of wedge-to-wedge variations in the radial profile extracted from the final stacked image (four 90° wedges were used) shows that the scatter is within a factor of less than 2 from expectations. Part of this extra noise comes from the remaining azimuthal asymmetries in the cluster core regions, where statistical errors are very small. In the most important radial range between ∼0.6 and 3 × R_500c, the scatter between wedges agrees well with expectations. We therefore concluded that other types of uncertainties, e.g. large-scale sky background variations, do not dominate the error budget.

5 DISCUSSION

Given our choice³ of |$R_{\rm 500c}=1\, {\rm Mpc}$| as a fiducial value, the corresponding mass is |$M_{\rm 500c}=500\rho _c\frac{4}{3}\pi R_{\rm 500c}^3\approx 2.85\times 10^{14}\, {\rm M_\odot }$|⁠. Below, we will use this mass to scale models for comparison with observational results.

The radial gas density profile⁴ that extends up to ∼3 × R_500c is the most principal result of this study. For estimates, we assume that the uncertainties in the recovered density profile at a given radius can be characterized by the scatter between best-fitting model profiles at the same radius extracted in four 90° wedges (scaled down by a factor of 2). On top of this, comes the uncertainty in the conversion coefficient ϵ_c discussed above.

5.1 Density

The density profile (number density of protons normalized by the mean density |$n_{p,z=0}=n_{{\rm H},z=0}=1.92\times 10^{-7}\, {\rm cm^{-3}}$| at |$z$| = 0) corresponding to the best-fitting model is shown in Fig. 4 with the blue line (a set of best-fitting parameters is given in Appendix  C). The light-blue shaded region shows the r.m.s. scatter between the best-fitting models in four 90° wedges. This r.m.s. value was divided by 2, to show the expected amplitude of variations in the profile under the assumption that the variations of the surface brightness in the four wedges are statistically independent. The dashed vertical line marks the position of the last significant point in the surface brightness profile (see Fig. 3). At this radius, the r.m.s. is ∼0.3 of the best-fitting model value.

The gas density profiles around massive halos have been calculated in numerical simulations (e.g. Lau et al. 2015; O’Neil et al. 2021) and, also, derived from X-ray data (e.g. Vikhlinin et al. 2006; Croston et al. 2008; Eckert et al. 2012; Pratt et al. 2022). For example, the red line in Fig. 4 shows the gas overdensity profile of the stacked halos between 10¹⁴ and |$10^{14.5}\, {\rm M_\odot }$| at |$z$| = 0 in TNG300-1 from O’Neil et al. (2021; see their fig. 5). The ratio of the two profiles is shown in the bottom panel with the purple line. Within the radial range (∼ 0.2–3) × R_500c, the two profiles agree within ±30 per cent. The brown curve shows the ratio of the gas density (again normalized by n_{p, |$z$| = 0}) from the X-ray data to the total matter overdensity from O’Neil et al. (2021). As expected, this ratio (known as f_gas) is close to the mean baryon fraction Ω_b/Ω_m ≈ 0.15 beyond R_500c.

In the same figure, we show the density profile model obtained by Pratt et al. (2022, 2023) from the XMM-Newton observations of 118 clusters. It shows an excellent agreement with the eROSITA-based profile in the range (0.03–1) × R_500c. At large radii (R_500c < R < 3R_500c), the models start to deviate from each other. In this regime, the extrapolation of Pratt et al. (2022) model predicts higher density than derived from the eROSITA data.

5.2 Temperature and entropy

Measuring stacked spectra to determine the gas temperature beyond R_500c is much more difficult than measuring the X-ray surface brightness in a broad energy band. However, one can use the ratio of pressure profiles derived from SZ data to the density profile from X-ray data to measure the temperature. Various approximations of the pressure profiles (usually referred to as Universal Pressure Profile or UPP) are available, starting from the work by Arnaud et al. (2010), which was based on X-ray data and used a parametric model (generalized NFW) motivated by numerical simulations (Nagai, Kravtsov & Vikhlinin 2007). Subsequent studies (e.g. Planck Collaboration V 2013; Sayers et al. 2013; Pointecouteau et al. 2021; Melin & Pratt 2023) used SZ data that can probe larger radii than typically mapped with X-rays. In the top panel of Fig. 5, we plot the universal pressure profile from Planck Collaboration V (2013; the red line) scaled to our reference value of the mass |$M_{\rm 500c}=2.85\times 10^{14}\, {\rm M_\odot }$|⁠. This pressure profile was converted from the electron pressure profile P_e to the ‘proton’ pressure, i.e. |$P_{\rm p}=n_{\rm p} kT=\frac{\mu _{\rm e}}{\mu _{\rm p}}P_{\rm e}$|⁠, where μ_e and μ_p are the mean atomic weights per electron and proton, respectively. For comparison, we plot the density profile n_p(r) from Fig. 4 with the blue line. The ratio of these two curves is the effective temperature that is shown in the bottom panel. For comparison, the horizontal dashed line shows the expected value of temperature from the scaling relation of Vikhlinin et al. (2006) for |$M_{\rm 500c}=2.85\times 10^{14}\, {\rm M_\odot }$|⁠. This value corresponds to the mass-weighted temperature T_mg that on average is some 20 per cent lower than the peak temperature in the cluster radial profile. In the derived temperature profile, the broad peak of |$T\sim 4.8\, {\rm keV}$| is near ∼0.25R_500c. Within R_500c the temperature profile agrees with typical trends found in observations of nearby clusters and simulations (e.g. Markevitch et al. 1998; De Grandi & Molendi 2002; Loken et al. 2002). At 3 × R_500c, the temperature drops to |$\sim 1\, {\rm keV}$|⁠, i.e. a factor of ∼5 lower than the peak value. Other variants of the UPP yield qualitatively similar results, although the scatter is substantial. In particular, the wiggles in the pressure profile found by Anbajagane et al. (2022) are just beyond 3 × R_500c. This scatter will be discussed in more detail in the subsequent work.

The gas entropy profile is yet another important quantity that characterizes the history of gas heating via shocks and other dissipative processes (e.g. Tozzi & Norman 2001; Voit, Kay & Bryan 2005; Cavagnolo et al. 2009; Pratt et al. 2010; Walker et al. 2012; Ghirardini et al. 2019). The same approach of combining the X-ray and SZ data can be used to estimate the entropy of the gas. We follow the standard convention and show the value of |$K_{\rm e}=P_{\rm e}/n_{\rm e}^{5/3}$| in units of |${\rm keV\, cm^2}$| in the middle panel of Fig. 5. The radial dependence of entropy across the radial range shown in Fig. 5 can be approximated with a power law |$K_{\rm e}=1150 \left({r/R_{\rm 500c}} \right)^{0.8}\, {\rm keV\, cm^2}$| (the dashed red line). While this dependence is shallower than expected in non-radiative models (Tozzi & Norman 2001; Voit et al. 2005), deviations from the steeper variants of the entropy profile approximation (outside the innermost region) are within a factor of ∼2. For example, the dashed orange line corresponds to |$K_{\rm e}=73+1260 \left({r/R_{\rm 500c}} \right)^{1.21}\, {\rm keV\, cm^2}$|⁠, which is based on the full ACCEPT sample (see table 4 in Cavagnolo et al. 2009).

5.3 Clumpiness

Gas clumpiness⁵ is one of the issues that complicate the conversion of the X-ray surface brightness to the gas density. In the simplest case, when the emissivity does not depend on temperature, the X-ray surface brightness depends on 〈n²〉 that is larger than 〈n〉², which makes the density estimates ∝ 〈n²〉^1/2 biased high. Some level of clumpiness is inevitably present both in the simulations and real objects and can affect all thermodynamic properties, although the magnitude of the effect is still a question of debate (e.g. Simionescu et al. 2011; Zhuravleva et al. 2013; Ettori 2015). In particular, numerical simulations suggest that the probability density distribution of the gas density in radial shells can be well described by a lognormal distribution, once a small number of high-density clumps is removed (e.g. Kawahara et al. 2007; Khedekar et al. 2013; Zhuravleva et al. 2013), i.e.

where n is the gas density divided by the median density value. In this case, the clumpiness factor

The bias in the value of the gas density derived from X-ray flux is equal to |$C^{1/2}=e^{\sigma ^2/2}$|⁠. The known uncertainties shown by shaded regions in the bottom panel of Fig. 4 allow for a factor of 2 larger densities derived from X-ray data compared to that found in numerical simulations. This corresponds to the value of |$\sigma =\sqrt{2\ln 2} \approx 1.2$|⁠. Typical values of σ found in numerical simulations at (2 − 2.5) × R_500c are in the range of 0.4–0.6 [see fig. 4 in Zhuravleva et al. (2013) or fig. 4 in Khedekar et al. (2013)]. Therefore, the derived value can be considered as a confirmation that much more extreme levels of X-ray-emitting gas clumpiness are not required by the data. The curve |$\rho^{stacked}_{gas, X}/\rho_{gas, Sim}$| shown in Fig. 4 remains below ∼2 at all radii, implying that a rough upper limits of σ ≲ 1.2 applies to the entire radial range probed. We note here that the calculated ratio |$\rho^{stacked}_{gas, X}/\rho_{gas, Sim}$|_, is better suited to demonstrate overall agreement between observations and simulations (and set crude upper limits on the gas clumpiness), than to derive the actual level of clumpiness at a given radius.

A more elaborate version of the same estimate could include a possible correlation of the density and temperature variations (e.g. Khedekar et al. 2013). However, this is important only if there are strong variations of the emissivity with temperature. Given the effective (i.e. calculated taking into account eROSITA’s response function) emissivity curves shown in the Appendix  B, it is clear that the departures towards higher temperatures are unlikely to affect the above conclusions. On the other hand, for deviations towards a much smaller temperature, say below 0.1 keV, the effective emissivity drops dramatically, effectively reducing the contribution of very cold clumps to the chosen X-ray band. Therefore, the above conclusion that the X-ray-emitting gas is not extremely clumpy stays unchanged. A much larger sample combined with a matching sample of simulated clusters would be needed to (i) bring down the uncertainties at large radii, (ii) verify the consistency of directly observed quantities, e.g. surface brightness in observations and simulations, and (iii) verify that the level of clumpiness found in simulations can be robustly recovered over the range or radii.

5.4 Can we go beyond 3 × R_500c?

Figs 4 and 5 suggest that X-rays remain a good proxy of the gas at least up to ∼3 × R_500c. At some larger distances from the cluster, the X-ray signal will plausibly be dominated by a highly inhomogeneous combination of individual halos (e.g. groups of galaxies) and the warm-hot intergalactic medium emission. The signal there will also reflect a complicated structure of the accretion/merger shocks (Vazza et al. 2009; Aung, Nagai & Lau 2021; Zhang et al. 2021) and/or a non-equilibrium ionization (e.g. Yoshikawa & Sasaki 2006). At even larger distances, e.g. comparable to the cluster turn-around radius, the excess X-ray signal may simply reflect the overall overdensity of all ‘objects’ compared to volume-averaged quantities for a representative chunk of the Universe. In other words, the X-ray signal might switch from the quadratic dependence of the baryon density |$\propto \rho _{\rm gas}^2$| in the inner part to a more linear dependence ∝ ρ_gas in the far outskirts. On top of this, come the changes in gas emissivity.

It is nevertheless interesting to ask, what are the prospects of detecting the gas beyond 3 × R_500c if the density follows predictions of numerical simulations, and the emissivity of the gas does not change dramatically, i.e. the X-ray signal is set by the line-of-sight integral |$\epsilon _c \rho _{\rm gas}^2$|⁠, where ϵ_c ≈ const is the same as below 3 × R_500c. At 3 × R_500c the local overdensity is ∼25–30 and it declines approximately as R⁻³. Assuming that the noise scales as ∝ N^−1/2A^−1/2, where N is the size of the sample, and |$A=2\pi \frac{\delta R}{R}R^2$| is the area of annulus with the radius R and a constant relative thickness |$\frac{\delta R}{R}$|⁠. Since the surface brightness declines as ρ²R ∝ R⁻⁵, the size of the sample should increase as N ∝ R⁸. Therefore, in order to go from ∼3 × R_500c to ∼4 × R_500c, the size of the sample has to be increased by an order of magnitude. While still feasible, it is clear that at these distances the detection of the proper X-ray signal will be very difficult against various systematic effects. Nevertheless, doing a larger sample makes sense in order to reduce the noise near 3 × R_500c. This will be a subject of the subsequent publication. The plan is to supplement the above analysis with stacking that uses R_200m instead of R_500c, which might be more appropriate for the outskirts.

As a caveat, we mention that the full calibration of the stray light profiles is yet to be completed and therefore some adjustments in the outer radial bins are possible in the future. Here, the main conclusion is that with eROSITA it is possible to reach levels ≲ 1 per cent of the sky background in the 0.3–2.3 keV that are needed to place constraints on the X-ray emission near 3 × R_500c. With this conclusion in mind, it is possible to consider larger samples and a fully ‘forward’ approach to the fitting procedure for a sample of cluster profiles rather than the stacked image approach adopted here.

6 CONCLUSIONS

Exploiting the ‘unlimited’ field-of-view of SRG/eROSITA in the all-sky survey, we stacked X-ray images of 38 nearby galaxy clusters drawn from the CHEX-MATE sample. The images cover the radial range up to the turn-around radii of these clusters. The excess X-ray emission on top of the sky background is detected up to ∼3 × R_500c. The recovered gas density profile agrees well with the results of cosmological simulations in the standard Λ cold dark matter model and does not show evidence of strong gas clumpiness (for the gas with temperatures in the range relevant to the X-ray band). Combining the density profile with the universal pressure profile from SZ data yields the temperature and entropy profiles that do not feature any strong change in their radial trends up to the maximum radius probed.

Acknowledgement

We are grateful to the referee for careful reading of the manuscript and valuable suggestions and comments.

This work is based on observations with the eROSITA telescope onboard SRG space observatory. The SRG observatory was built by Roskosmos in the interests of the Russian Academy of Sciences represented by its Space Research Institute (IKI) in the framework of the Russian Federal Space Program, with the participation of the Deutsches Zentrum für Luft-und Raumfahrt (DLR). The eROSITA X-ray telescope was built by a consortium of German Institutes led by MPE, and supported by DLR. The SRG spacecraft was designed, built, launched, and operated by the Lavochkin Association and its subcontractors. The science data are downlinked via the Deep Space Network Antennae in Bear Lakes, Ussurijsk, and Baikonur, funded by Roskosmos. The eROSITA data used in this work were converted to calibrated event lists using the esass software system developed by the German eROSITA Consortium and analysed using proprietary data reduction software developed by the Russian eROSITA Consortium. AAS was partly supported by the project number 0033-2019-0005 of the Russian Ministry of Science and Higher Education. IK acknowledges support by the COMPLEX project from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program grant agreement ERC-2019-AdG 882679.

DATA AVAILABILITY

X-ray data analysed in this article were used with the permission of the Russian SRG/eROSITA consortium. The data will become publicly available as a part of the corresponding SRG/eROSITA data release along with the appropriate calibration information. All other data are publicly available and can be accessed at the corresponding public archive servers.

Footnotes

References

APPENDIX A: LIST OF CLUSTERS

Table A1 lists 38 clusters drawn from the CHEX-MATE sample (CHEX-MATE Collaboration 2021) used for stacking. Selection criteria are discussed in Section 2.

APPENDIX B: RE-SCALING IMAGES FOR STACKING

The observed energy flux from a volume element dV of a cluster in the energy range between E₁ and E₂ is

where |$\Lambda (E^{^{\prime }},T)$| is the gas emissivity of the cluster at the rest-frame energy E′, n_e, and n_p are the electron and proton densities, respectively, D_a is the angular diameter distance, and |$z$| is the redshift. Therefore,

since |$E^{^{\prime }}=E(1+z)$|⁠. In terms of the surface brightness (in counts per second per unit solid angle) recorded by a telescope with effective area A(E),

where l is the distance along the line of sight and ϵ_c(T, |$z$|⁠, E₁, E₂) for eROSITA is given in (Churazov et al. 2021, see their appendix C).⁶

When stacking cluster images with different masses and at different redshifts, the aim is to rescale the angular size and the observed surface brightness so that if these clusters obey the most simple form of scaling relations they will appear similar to each other. In particular, one can relate the mass to some characteristic radius and adjust the angular size and density accordingly. While for clusters’ outskirts scaling with radii at fixed matter overdensity is more appropriate, e.g. R_200m (see e.g. O’Neil et al. 2021), we used here the radii at fixed overdensity relative to the critical density of the Universe at a given redshift, namely R_500c. Since the mass and redshift ranges used here are small 0.05 < |$z$| < 0.2, the impact of choosing a different scaling radius is expected to be small. Accordingly, the observed surface brightness I_obs for every cluster was renormalized assuming that the typical gas density at R_500c is ∝ E(⁠|$z$|⁠)² and the physical size is ∝ R_500c. Thus,

From equation (3) it follows that the recovered density scales as |$\epsilon _c^{-1/2}$|⁠. The sensitivity of ϵ_c to the Galactic hydrogen column density, hot gas metallicity, and temperature are illustrated in Fig. B2. The black curve (Z/Z_⊙ = 0.2) can be considered typical (e.g. Ghizzardi et al. 2021; Angelinelli et al. 2022) for the outskirts of clusters located far above the Galactic plane (hence low N_H). Corresponding variations of |$\epsilon _c^{-1/2}$| with gas temperature are modest down to |$kT \sim 0.3\, {\rm keV}$|⁠. For instance, the uncertainties in temperature at 3 × R_500c (Fig. 5) derived under the assumption of ϵ_c(T) = const, imply |$kT\approx 1\pm 0.3\, {\rm keV}$|⁠. From Fig. B2 we can infer that corresponding variations of |$\epsilon _c^{1/2}$| are ∼10 per cent, i.e. a factor of ∼3 smaller. The same is true for the uncertainty arising from the redshift dependence of ϵ_c (see Fig. B1). These arguments show that there is no internal inconsistency in the calculation of gas density and temperature arising from the assumption of ϵ_c(T) = const. However, they do not provide direct proof that the X-ray-emitting gas is homogeneous in temperature and/or density. A more accurate answer would be possible with future microcalorimeter missions, especially those with full imaging capabilities like LEM (Kraft et al. 2022) that can measure both the mean temperature in the cluster outskirts as well as a distribution over gas temperatures using Fe-L complex and lines of lighter elements like Ne and O below 1 keV. In any case, X-ray observations place constraints only on the X-ray-emitting gas. Cooler gas, invisible in X-rays, might be present too. However, the broad agreement of the density profile derived in Section 4 with simulations suggests that the hot phase dominates at least up to 3 × R_500c.

APPENDIX C: FIDUCIAL MODEL PARAMETERS

The multiparameter model used to fit the X-ray surface brightness (see equation 2) is flexible enough to describe the data equally well with a slightly different combination of parameters. In this sense, only the overall shape of the recovered density distribution makes sense rather than the values of individual parameters. With this clause in mind, the following set of parameters yields the resulting density profile in the form

where r, r_c, r_s, r_s2 are units of R_500c. The best-fitting values of the parameters are: |$N=1.4\times 10^{-2}\, {\rm cm^{-3}}$|⁠, r_c = 0.187, α = −0.7, β = 0.193, r_s = 20.2, g = 3.29, e = 5.73, r_{s, 2} = 20.8, g₂ = 0.698, e₂ = 7.99.

The relevant radial range, after exclusion of the outer regions where the signal is below 1 per cent of the background is ∼(10⁻² − 3) × R_500c.

For the entropy radial profile in the range ∼(3 × 10⁻² − 3) × R_500c, the following power-law fit provides reasonable approximation:

APPENDIX D: WEDGE-TO-WEDGE VARIATIONS

In Section 4, an azimuthally averaged X-ray surface brightness profile is discussed and approximated with an analytic function (see also Appendix  C). Here we show (Fig. D1) deviations from this function in individual 90° wedges. In this plot, I_{X, m} stands for the full analytic model that includes both the cluster emission and a flat background. Inside the central region (10 arcmin corresponds to R_500c), where the cluster emission dominates, deviations at the level of ∼10 per cent are present, reflecting the residual variations of the surface brightness in the stacked image. These deviations might be caused by the uncertainties in the definition of cluster centres, intrinsic ellipticities of individual objects, etc. At large radii, the sky background (CXB and Milky Way foreground) dominates and the plot illustrates that the final scatter of the measured data points in wedges amounts to ∼2 per cent of the total background.