The eROSITA Final Equatorial-Depth Survey (eFEDS):

For all the clusters in our sample we apply the procedure presented in Ghirardini et al. (2021a) and Liu A. et al. (submitted) in order to recover surface brightness and density profiles. While we direct the reader to Ghirardini et al. (2021a) for details, we here provide a short summary of the major steps of the analysis. We start from the clean event files, we extract images and exposure maps in the 0.5–2.0 keV energy band using eSASS tasks evtool and expmap respectively with size as large as 4 times $R_{500}$ around each cluster position. We directly fit the 2D distribution of the X-ray photons in the produced images of these clusters, allowing for multiple cluster fits at the same time (when more than one cluster is present in the image), for fitting of cluster centers, and for modeling of point sources detected. As in our previous works we model only bright point sources (with more than 0.1 count per second), because modelling all of them would be computationally too expensive. These point sources are modelled as delta functions convolved with the PSF, thus taking into account effects caused by PSF wings. Faint point sources are masked by removing a circular region around them whose radius is large enough to assure that point source signal outside is consistent with the background level. The clusters are modelled in 2D using the Vikhlinin et al. (2006) model:

where the priors on our parameters are $\epsilon<5$ (as suggested by Vikhlinin et al., 2006), $\beta>1/3$, and $\alpha>0$, and we freeze $r_{s}=r_{c}$. If more clusters are present in the image, then this model above is also used for these other clusters. The center of the clusters is not fixed, but allowed to vary using a gaussian prior centered on the detection location and $\sigma$ of 20 arcsec. The resulting model cluster images are convolved with the PSF of eROSITA. The instrumental background, see Freyberg et al. (2020), (particle-induced background and camera noise) model folded with the un-vignetted exposure map, and the sky background (including the cosmic X-ray background and the soft background component from the galactic halo and local bubble) folded with the vignetted exposure map are added to the total model.

We then fit the image obtained from the eROSITA observations with the model image in 2D using the Monte Carlo Markov Chain (MCMC) code emcee (Foreman-Mackey et al., 2013) to find the best-fit model parameters. We assume a Poisson log-likelihood function $\sum N_{i}-\mu_{i}\log N_{i}$, where $N_{i}$ are the model predicted counts, and $\mu_{i}$ are the observed counts in each pixel of the image.

The fitting of the images can be interpreted as a true density once we consider the emissivity of the gas, which is determined, as in our previous works, by making use of the spectral information. Similarly as in Ghirardini et al. (2021a) and Liu A. et al. (submitted) we fit cluster spectra within $R_{500}$ to obtain the conversion factor from count rate to emissivity, thus determining the physical properties of the clusters. In this work, we make use of only the products of imaging data, that are only density profiles, surface brightness profiles, and luminosity within $R_{500}$. The final luminosities used in this work are provided in Bahar et al. in preparation.

The central value of the gas density, $n_{0}$, is another indicator of the relaxation state of clusters, since relaxed systems are expected to have higher central densities (Hudson et al., 2010). We use the value of the electron density computed at 0.02 $R_{500}$. This value is quite close to the center, and for the vast majority of eFEDS clusters this is well within the field-of-view (FoV) averaged PSF value of eROSITA. However, as we also specified for the concentration parameter, the electron density profile is deconvolved by eROSITA’s PSF, see Sect. 2.1. Therefore, the $n_{0}$ we present in this work is PSF-independent, and can be used to compare with previous results. We do not consider the self-similar evolution of the density profile (McDonald et al., 2017; Ghirardini et al., 2021b) as a redshift correction on this morphological parameter because we will model it separately in order to identify the actual evolution of the central density with redshift.

The concentration parameter, $c_{\rm SB}$, has been introduced by Santos et al. (2008) as an indicator of the presence of a peaked X-ray surface brightness, which has been shown to correlate with the relaxation state of clusters. It is defined as the ratio between the integrated surface brightness in two different circular apertures. In this work we use two definitions: the original one introduced by Santos et al. (2008),

and one scaled with $R_{500}$ (Maughan et al., 2012)

When calculating the concentration parameter for our sample, the effects caused by eROSITA’s PSF are fully taken into account. In fact, the concentration is not simply the count ratio in the two apertures, but comes from the two dimensional image fitting we have described in Sect. 2.1 which produces PSF-deconvolved surface brightness profiles. This also means that the concentration we measure here, being corrected for instrumental effects, can be directly compared with previous results with different instruments.

The ellipticity parameter, $\epsilon$, is defined as the ratio between the minor and major axes. Contrary to several past studies, we do not compute the two axes using the second-order moments (variance) of the flux distribution in the cluster image because statistically the variance is the “width” of the distribution only when the distribution is gaussian. Rather, it is known that the photon distribution in cluster images is not gaussian but is similar to a $\beta$-model. Therefore we measure the ellipticity by extending the analysis described in Sect. 2.1, allowing the density profile to be elliptical and rotated on the plane-of-the sky. This is accomplished by introducing a rotation angle and the ellipticity in the model construction of the image. In short, the 2 dimensional density profile along the x-axis is squeezed by changing the core radius as $r_{c,x}=\epsilon\times r_{c,y}$, and then it is rotated by an angle $\theta$. The best fitting parameter $\epsilon$ is the ellipticity we adopt in this work. We note that on these two parameters we put uniform prior between 0 and 1 on the ellipticity, and uniform prior between 0 and $\pi$ on the rotation angle, not between 0 and 2$\pi$ to avoid problems caused by the $\pi$ rotational symmetry of the ellipse.

The power ratios $P_{m}$ consist of a 2-dimensional decomposition of the surface brightness distribution within a specific aperture, of $R_{500}$, and they account for radial fluctuations since the higher the order of the power ratio the higher the sensitivity to smaller fluctuations. They are introduced by Buote & Tsai (1995), and are defined as $P_{m0}=P_{m}/P_{0}$, where

and

The values of $a_{m}$ and $b_{m}$ are calculated as

where the X-ray surface brightness $S_{B}(x)$ is calculated at the position expressed in polar coordinates $x=(r,\phi)$. In this work, we include only the power ratios up to order 4, that is $P_{10}$, $P_{20}$, $P_{30}$, and $P_{40}$. We remind the reader that power ratios can only be directly computed from the images, therefore it is not possible to be corrected for eROSITA’s PSF. Since the results will be affected by the PSF, comparison with the literature will be challenging due to instrumental differences. For this reason we do not compare the obtained values with previous studies.

In Fig. 3 we show the parameter-parameter planes and in Fig. 4 we show the correlation matrix between the parameters. As observed in previous works (e.g. Lovisari et al., 2017) there is strong correlation in-between core sensitive parameters, like central density and concentration, and in-between core insensitive parameters, like power ratios and photon asymmetry.

In order to investigate the redshift dependence of the morphological parameters, we first divide our sample into four homogeneous sub-samples, selecting the cluster population according to their quartile in the redshift distribution, $z<0.23$, $0.23<z<0.35$, $0.35<z<0.46$, and $z>0.46$, therefore with same number of clusters in each bin.

It must be noted that the detection (or non-detection) of a significant redshift evolution of samples covering a large redshift and luminosity span can be caused by the significant correlation between luminosity and redshift. Therefore, the redshift and luminosity evolution should be simultaneously modelled. In fact, as shown in Fig. 2, we observe a clear trend where at higher redshifts we observe a larger number of luminous clusters due to eROSITA’s selection function. Therefore, we further divide the sample based on their luminosity as in the following: $L_{500}<1.1\times 10^{43}$ erg/s, $1.1\times 10^{43}$ erg/s ¡ $L_{500}<2.6\times 10^{43}$ erg/s, $2.6\times 10^{43}$ erg/s ¡ $L_{500}<5.9\times 10^{43}$ erg/s, and $L_{500}>5.9\times 10^{43}$ erg/s, where also in this case we have used the quartile of the distribution to have approximately the same number of clusters in each bin. The cumulative distributions of the morphological parameters for the sample of eFEDS clusters, split based on the above criteria, are shown in Fig. 5 and Fig. 6.

We perform the Kolmogorov-Smirnov statistical test (KS test) to check whether the morphological parameters in each redshift (and luminosity) bin are drawn from the same parent distribution as the unbinned cluster sample and whether there are clear indications for evolution with redshift. We compute the $p$-value from the KS test, which gives the probability for clusters in each redshift quartile to be indistinguishable from the unbinned cluster sample. Small value indicates that the morphological parameters in a specific redshift bin are clearly different than the unbinned distribution, while large values indicate they are indistinguishable. We find that the central density and the Gini coefficient have very small $p$-values, indicating that these two parameters might significantly evolve with redshift, while for the other parameters, with large $p$-values, such evidence is not present. On the other hand, all the morphological parameters, except for the cuspiness, ellipticity, and concentration, show a significant luminosity dependence. However, we remind that such dependence might be mimicked by or canceled out by the degeneracy between luminosity and redshift, therefore conclusive analysis on the redshift evolution or luminosity dependence can be derived only when we model at the same time both dependencies, considering eROSITA selection function and luminosity function. In the next section, we provide the detailed analysis of this modeling.

eROSITA’s great survey capabilities allow to reach down to luminosities as faint as $10^{41}$ erg/s (corresponding approximately to a group-scale mass of $10^{13}M_{\odot}$), and allows to study high redshift clusters up to $z$ = 1.2. For the first time, we are able to constrain both luminosity and redshift evolution of the morphological parameters simultaneously. To understand the underlying evolution in both redshift and luminosity simultaneously in morphological parameters, we model a power-law relation between each parameter $\mathcal{M}$, luminosity $L_{500}$, and redshift $z$, taking into account both the selection function (probability of detecting a cluster with given luminosity and redshift) and the luminosity function (the probability of detecting a cluster with a given luminosity and redshift at a fixed cosmology), thus taking into account the observed distribution in luminosity and redshift space of our clusters. We note that the scaling relation used to connect intrinsic morphological parameters $\mathcal{M}$ with luminosity and redshift is written as:

where $L_{\rm piv}=2.6\times 10^{43}$ erg/s and $z_{\rm piv}=0.35$ are the median luminosity and redshift respectively, of the eFEDS sample. $\gamma$ represents luminosity dependent slope, and $\beta$ represents the slope of the redshift evolution. The likelihood derivation is provided in detail in Appendix A. The final likelihood for each cluster is written in Eq. (31). The best-fit values of all morphological parameters are given in Table 1. We show the morphological parameters as a function of luminosity and redshift, while the best-fit models are shown in gray shaded regions in Fig. 7. We highlight a few important results obtained from our joint modeling:

The observed increase in central density values with luminosity can be easily explained: first we notice that gas density does not scale with mass according to the self-similar model, however literature results (Pratt et al., 2009; Lovisari et al., 2015) show that gas mass fraction increases with cluster mass, and the only way to obtain a higher gas mass fraction at a fixed scaled radius is by increasing the gas density profile, and therefore increasing the central density as well. On the other hand, McDonald et al. (2017) find lack of evolution in the distribution of the central density combining the following three cluster samples: 49 low-$z$ X-ray selected clusters (Vikhlinin et al., 2009), 90 SPT selected clusters from $z=0.2$ to $z=1.2$ (McDonald et al., 2013), and the 8 most massive clusters in the SPT cluster survey (Bleem et al., 2015) above $z=1.2$. They interpret this as lack of evolution in the core properties of clusters in the last $\sim$10 Gyr. A similar result was later reported in Ghirardini et al. (2021b), when comparing the core properties of the same 8 most massive clusters in SPT (Bleem et al., 2015) above $z=1.2$, and a Planck selected sample of 12 low redshift ($z<0.1$) clusters (X-COP, Eckert et al., 2017). On the contrary, Sanders et al. (2018) analyzed a sample of 83 SPT selected clusters from $z=0.2$ to $z=1.2$ and found no significance difference with respect to the self-similar evolution model. One important remark is that McDonald et al. (2017) selected clusters in order to satisfy evolutionary requirements, meaning that clusters at low redshift where chosen in a mass range expected from the evolutionary scenario (Fakhouri et al., 2010), and therefore they might not be representative of the cluster population, but of how single clusters are expected to evolve.

Nurgaliev et al. (2017) compared a sample of 36 ROSAT selected clusters (400d, Burenin et al., 2007) in redshift range 0.35 $<z<$ 0.9 with a sample of 90 clusters selected from the SPT survey in redshift range 0.25 $<z<$ 1.2. Measuring photon asymmetry, they find no statistical difference between the morphological properties of their X-ray and SZ selected cluster samples. Furthermore, they find absence of mass or redshift evolution in photon asymmetry. They also studied a sample of 85 simulated clusters, applying X-ray and SZ selection, and since the measured X-ray morphology is indistinguishable, they conclude that X-ray and SZ surveys are probing the same cluster population. They interpret this lack of evolution in the morphological properties as a lack of direct correlation between dynamical state of the clusters and these properties. Their results are justified by theoretical studies that find that substructures statistics can vary significantly on very short time scales during cluster mergers. The authors claim that this absence of difference between the X-ray and SZ cluster samples indicates that high resolution (1 arcmin) SZ surveys are not biased toward selecting preferentially merging clusters, while other SZ instruments with lower resolution might still have a bias since it is more likely that multiple clusters along the line of sight contribute to the integrated signal.

In summary, it is difficult to properly quantify the effect of PSF on the measured evolution with redshift for some parameters (Gini coefficient, photon asymmetry, and power ratios). The reduction in angular size of clusters at high redshift should imply that cluster images become smoother at higher redshift, and smoother images will result in larger Gini coefficients and smaller photon asymmetry values. The fact that we observe the opposite trend implies that the underlying evolution of these parameters should be even stronger than the measured values. The same logic then can be applied to the power ratios as their evolution should, also in that case, be stronger than the measured ones. The only parameter that would indicate the opposite trend, with more cool cores at high redshift, is the central density. However in this case the expected evolution is not very significant, and at 3$\sigma$ marginally consistent with no evolution.

The fitting procedure we have applied on our data to recover redshift evolution and luminosity dependence of our clusters, see Sect. 4.2, can be also exploited to define new morphological parameters which are constructed from the original distributions after the luminosity and redshift dependence are factored out:

using $\gamma$ and $\beta$ from the corresponding row of Table 1. By construction, these newly defined parameters lack any redshift evolution or luminosity dependence, and therefore can be used to properly and correctly investigate the presence or absence of bimodality in our parameter distribution. In Fig. 8 we show the changes in the distribution of the corrected morphological parameters once this correction is applied. Some correlations become tighter, while others become more loose. In particular the correlation between core properties, central density, both concentrations, and cuspiness, become tighter, indicating that all these four parameters measure just the cluster core properties, and therefore the connection among them is almost one to one. On the other hand, correlation between parameters sensitive to large scale fluctuations, power ratios, Gini coefficient, and photon asymmetry in particular, become much less significant, indicating that these parameters are measuring complementary properties of the large scale emission of clusters.

In the literature, it is common to use morphological parameters to characterize the physical state of the ICM, determining which clusters are relaxed and which one are disturbed. In Lovisari et al. (2017), the authors use a morphological parameter, which was previously introduced in Rasia et al. (2013). This parameter combines the concentration and centroid shift values and then is used to distinguish between the relaxed and disturbed systems. The authors measure the deviation from the mean value of concentration and centroid shift, adding these together but changing the sign of the centroid shift deviation to take into account the anti-correlation between these two parameters. However, one important issue is that this parameter does not take into account the strength of the anti-correlation and, since the correlation between these parameters is very strong, this parameter is likely affected by double counting issues.

Here in this work, we introduce a new parameter, the relaxation score, or $R_{\rm score}$, that contains the information stored in all the morphological parameters by taking into account both the sign (positive or negative) of the correlations of the morphological parameters with respect to concentration as well as the strength of these correlations. The definition of $R_{\rm score}$ is:

where $\mathcal{MN}(\mathbf{\mu},\mathbf{\Sigma})$ is the multivariate normal PDF in $n$-dimensional space of our parameter distribution, where $\mathbf{\mu}$ is an array containing the mean values for all the morphological parameters $\mathcal{M}_{1}\dots\mathcal{M}_{n}$, and $\mathbf{\Sigma}$ is the covariance matrix computed from our data. In other words, we compute the cumulative distribution function (CDF) in $n$-dimensional space. For the quantities that anti-correlate with concentration, such as power ratios and photon asymmetry, we consider their reciprocal in the calculation of this quantity, Eq. (15), in order to build this parameters from quantities that correlate, not from a mix of correlating and anti-correlating quantities. As noted before, to avoid being biased towards redshift and luminosity dependence in our parameters, we consider the corrected morphological parameters as introduced in Eq. (14) both for computation of means and covariance matrix, and for computation of the $R_{\rm score}$ itself. By construction, the $R_{\rm score}$ should be higher for objects with large concentration, central density, ellipticity, cuspiness, or Gini coefficient, and lower for objects with high photon asymmetry, or high power ratios. This is reflected in Fig. 8.

We show the distribution of this newly introduced parameter in Fig. 9. At a first glance, it looks unimodal, however we refer to the next section for the appropriate investigation on the best fitting distribution for this new parameter. It is worth to point out that this parameter has been constructed with redshift and luminosity corrected morphological parameters, therefore, by construction, it does not depend on redshift or luminosity.

Fig. 10 shows the correlation between the concentration, computed using the (Maughan et al., 2012)’s definition with $R_{500}$, and the relaxation score. We observe a clear correlation between the two parameters, where relaxation score increases with concentration. The correlation becomes insignificant above $c_{\rm SB,R_{500}}=0.27$, the threshold adopted by Lovisari et al. (2017) to identify relaxed clusters. We interpret this as an indication that concentration is a good indicator to probe the very central state of clusters, but it does not correlate well with the relaxation state. In fact, a relaxed cluster will generally have a cool core in its center, while the opposite may not be true: a merger along the line of sight (e.g., Dupke et al., 2007), or a merger in its initial stage, might have a centrally peaked density profile. The mergers whose axis is along the line-of-sight might still not be visible from the cluster center prospective, while it is visible on a larger scale. Therefore we claim that clusters with large concentration can still be disturbed clusters, as supported by recent theoretical papers (e.g., Rasia et al., 2015; Biffi et al., 2016). On the other hand, by construction, the relaxation score is able to capture the ICM dynamics from small to large scales, and therefore can be used to identify the cases where a cool core is still existing in an disturbed cluster. Furthermore, in Fig. 10 we observe that the thresholds for disturbed clusters adopted by Lovisari et al. (2017), are in agreement with the median of the log-normal best fit shown with the horizontal black line, meaning that almost all the clusters to the left of the dashed line are also below the horizontal line. Therefore we make use of this line to establish a threshold in terms of relaxation score: we define clusters to be relaxed if $R_{\rm score}>0.0137$, while they are defined to be disturbed if the opposite happens.

This definition allows us to compute the fraction of relaxed clusters as a function of redshift, by dividing the cluster sample in 5 redshift bins, and in each bin counting clusters based on their relaxation score. We have to add a further complication to this analysis: the relaxation score distribution for most clusters is wide and quite skewed, and therefore computing the fraction using just the median values will bias our results. Therefore we have to count the number of clusters that satisfy our threshold by first randomizing the relaxation score values by using the MCMC chain obtained when estimating it. In other words we perform 10000 bootstraps iteration of the relaxation score for each cluster, each time calculating the fraction of clusters that are relaxed.

The obtained evolution of the fraction of objects that we classify as relaxed is shown in Fig. 11. We observe a fraction of relaxed clusters that is about 50% in the lowest redshift bin, and at the higher redshift bins it becomes about 30–35%, consistent with cool-core fraction estimated in previous works (e.g., Rossetti et al., 2017; Sanders et al., 2018). We therefore find a hint (slightly more than 1$\sigma$) of possible evolution with redshift of the fraction of relaxed clusters. However we might argue that the difference in the first bin can also be not a consequence of redshift evolution in the fraction of relaxed clusters, but may be caused by the incompleteness due to our source finding and confirmation technique in the detection of very extended, low-surface-brightness galaxy groups or clusters at low redshift, similarly as it was shown in a focused reanalysis of ROSAT observations (Xu et al., 2018).