© The Authors 2025

1. Introduction

Clusters of galaxies are permeated by hot (∼10⁷ − 10⁸ K), rarefied (∼10⁻² − 10⁻⁴ ions per cm³), optically thin gas, known as the intracluster medium (ICM), which is believed to be approximately in equilibrium in the total gravitational potential of the system, dominated by the dark-matter halo. The study of the ICM is key to measuring the parameters that define the cosmological background (e.g., Pratt et al. 2019) and to understanding non-gravitational physical processes, such as the mechanisms that prevent gas from cooling and forming stars in the cluster core. A major limitation to the determination of the cosmological parameters and to the understanding of the ongoing heating mechanisms in the ICM is the difficulty in gathering information on the velocity field (that is, the turbulent and bulk motions) of the gas.

The ICM emits in X-rays via thermal bremsstrahlung and the emission lines of the heavy metals (e.g., Sanders 2023, for a review), and distorts the cosmic microwave background (CMB) through the inverse Compton scattering of ICM electrons, known as the Sunyaev-Zeldovich (SZ; Sunyaev & Zeldovich 1972) effect (see Mroczkowski et al. 2019, for a review). Our current understanding of the ICM thus comes from spectral analyses of the broadband X-ray thermal continuum and assessment of the main CMB distortions. In recent years, several observational campaigns have aimed to advance our understanding of the thermodynamic properties of the ICM from the center out to the outskirts of clusters. These efforts have been focused on nearby, massive clusters, leveraging the joint power of X-ray and SZ effect data from the XMM-Newton and Planck observatories, respectively.

The limited sensitivity of microwave detectors to the weak distortions in the CMB spectrum caused by the bulk and turbulent motions of electrons in the ICM, known as the kinetic SZ effect (see Baldi et al. 2018; Mroczkowski et al. 2019; Altamura et al. 2023), and the degeneracy of this signal with the intrinsic properties of the ICM (like the gas density and geometry) has prevented the characterization of the velocity field in a large sample of objects. Though the poor spectral resolution of X-ray CCD-like spectrometers (such as XMM-Newton) has hindered the detection of the broadening and shifting of the X-ray emitting lines, several attempts have constrained the motions of the ICM to a few hundred km s⁻¹ in the inner regions, except for a few merging systems with velocities above 500 km s⁻¹ (e.g., Tamura et al. 2014; Liu & Tozzi 2019; Sanders et al. 2020; Gatuzz et al. 2022a,b, 2024). Overall, due to the limitations of the available instruments, we have relatively little information about the bulk and turbulent motions in the ICM both in the inner region, where, according to numerical simulations, galaxy-scale processes routinely inject significant kinetic energy into the ICM, and near the virial radius, where substantial kinetic energy is expected to be injected by large-scale structure processes (e.g., Vazza et al. 2012; Angelinelli et al. 2020).

Precious information on nonthermal motions in the ICM came from the spectra obtained by the X-ray Multi-Mirror Reflecting Grating Spectrometer (RGS) on board XMM-Newton in the 0.3−2.5 keV energy range. However, due to the slitless nature of the detector, the extent of the source contributes to the line width. By assuming that the objects under examination (typically the cores of galaxy clusters) are point sources, upper limits of ≈500 km s⁻¹ on the velocity broadening have been obtained (Sanders et al. 2011; Pinto et al. 2015; Bambic et al. 2018; Sanders 2023, for a review).

With the advent of the high-spectral-resolution X-ray spectrometer Resolve on board the X-ray Imaging and Spectroscopy Mission (XRISM; Tashiro et al. 2018), measurements of the bulk and turbulent motions in the ICM are now possible (see, e.g., Ota et al. 2018; Bartalesi et al. 2024, hereafter B24). However, due to the relatively small collective area and field of view (FoV) of XRISM, they are limited to the cores of nearby, massive clusters, where the signal is higher. The published measurements of the turbulent velocity dispersion of 164 ± 10 km s⁻¹ (obtained with Hitomi, the predecessor of XRISM; Hitomi Collaboration 2016) and 169 ± 10 km s⁻¹ (Miller et al. 2025) in the cores of the massive clusters Perseus and A2029, respectively, placed an upper limit of 0.04 on the nonthermal-to-total energy ratio in the central regions of these objects, which is an indication of dynamically unimportant turbulence in the inner regions of massive galaxy clusters.

Our current knowledge on the mass distribution of galaxy clusters primarily comes from gravitational lensing and, specifically in the outskirts, from statistical measurements of the weak lensing (WL) signal, that is, the mild distortion of the shape of background galaxies (for a review see Umetsu 2020). However, WL mass measurements for nearby massive objects can be complicated due to their large angular size in the sky. Given that the median value of the baryon fraction and the scatter around it agree far beyond the core between samples of clusters formed in cosmological hydrodynamical simulations with different solvers and sub-grid models, an independent, albeit indirect, estimate of the cluster mass can be obtained by assuming a “universal” value of the gas mass fraction as inferred from these simulations (Ghirardini et al. 2018; Eckert et al. 2019; Ettori & Eckert 2022). Since the first measurements of the mass of galaxy clusters were obtained, solid evidence has been found that the masses obtained under the assumption of hydrostatic equilibrium are biased low compared to the WL masses (e.g., Pratt et al. 2019; Ettori et al. 2019). This mass discrepancy, the so-called hydrostatic mass bias, has been attributed to our inability to account for some unresolved, residual kinetic energy in the ICM. This discrepancy can potentially explain much of the current tension in determining cosmological parameters: between those obtained from cluster halo abundances calibrated on hydrostatic masses and those inferred from modeling the primary anisotropies of the CMB (Planck Collaboration XXIX 2014).

The results of hydrodynamic cosmological simulations indicate that both rotation and turbulence are expected to contribute to the velocity field of the ICM (e.g., Lau et al. 2009; Suto et al. 2013; Baldi et al. 2017; Braspenning et al. 2025). However, there is no general consensus on how the nonthermal support is expected to be split between rotation and turbulence, because of the dependence on the specific implementation of some baryonic physics processes in hydrodynamic codes. For instance, as far as rotation is concerned, this uncertainty in the predictions of cosmological simulations is illustrated by Fig. 3 of B24, where the average ICM rotation velocity profiles obtained in two different cosmological simulations are shown (see also Suto et al. 2013; Braspenning et al. 2025).

The possibility that the ICM rotates significantly is interesting also for reasons that go beyond gauging the nonthermal support in clusters. In particular, given that the ICM is weakly magnetized (e.g., Bruggen 2013), the rotation of the ICM could be relevant to the energy balance of the gas in galaxy clusters because the magnetorotational instability (Balbus & Hawley 1991) could be at work in a magnetized rotating ICM (Nipoti & Posti 2014; Nipoti et al. 2015). The nonlinear evolution of the magnetorotational instability is expected to lead to turbulent heating, which could contribute to regulating the thermal evolution of the ICM, cooperating with feedback from active galactic nuclei in offsetting the radiative cooling of the ICM in the central regions of cool-core clusters (e.g., McNamara & Nulsen 2012; Hlavacek-Larrondo et al. 2022).

In an attempt to improve our understanding of the velocity field of the ICM, in this paper we present cluster models that allow for the presence of rotation and turbulence in the ICM. Using a Bayesian model-data comparison, these models are tuned to reproduce the radial profiles of the thermodynamic properties of the ICM as resolved in the XMM-Newton and Planck maps for the clusters of the XMM-Newton Cluster Outskirts Project (X-COP) sample (Eckert et al. 2017), and to be consistent with the available virial mass estimates. Among the studied clusters, A2255 turns out to be the most promising object in which ICM rotation could be detectable with currently available X-ray spectrometers: we present a mock spectrum of an XRISM pointing in the center of this cluster based on our model with turbulent and rotating ICM.

The paper is organized as follows. In Sects. 2 and 3 we describe, respectively, the observational data and the models; Sect. 4 details the statistical method used in our analysis; in Sect. 5, we present the results of our analysis, we compare them with those of previous works, and we show and discuss the mock spectrum of A2255; Sect. 6 summarizes our main findings. Throughout this work, we assume a flat Λ cold dark matter (CDM) cosmological model¹, with present-day matter density parameter Ω_m, 0 = 0.3 and Hubble constant H₀ = 70 km s⁻¹ Mpc⁻¹. We indicate with M_Δ the mass enclosed within a sphere and centered in the cluster center with a radius r_Δ such that the average density in the sphere is Δ times the critical density of the Universe, ρ_crit ≡ 3H²(z)/(8πG), where H(z) is the Hubble parameter. In the Bayesian analysis, we adopt the 68% highest posterior density interval (HPDI²) as the 1σ credible interval of the marginal posterior (e.g., Sect. 2.3 of Gelman et al. 2013).

2. Published data from X-ray, SZ, and WL analyses of X-COP clusters

We analyzed the clusters belonging to the sample of X-COP, a large observational campaign designed to advance our understanding of the virialization region of galaxy clusters (Eckert et al. 2017). The X-COP sample consists of 12 nearby (0.04 < z < 0.09), massive (2 × 10¹⁴ M_⊙ ≲ M₅₀₀ ≲ 3 × 10¹⁵ M_⊙) galaxy clusters. The clusters were selected to have the highest signal-to-noise ratio of the SZ effect as resolved in the Planck maps (Planck Collaboration XXIX 2014). The cluster sample was further refined by excluding the objects that appeared significantly disturbed. The remaining clusters were followed up with XMM-Newton to obtain detailed spectral data out to ∼r₅₀₀. As discussed in Eckert et al. (2022a), in clusters with a very peaked density profile near the cluster center, the point spread function (PSF) of XMM-Newton contaminates the spectra extracted from a given radial bin with some emission coming from neighboring bins. Given that in the model-data comparison we neglect this effect, in this work we excluded A2029, which has a very peaked central density profile. Thus, our sample of observed clusters consists of the 11 objects listed in Table 1.

2.1. ICM thermodynamic properties from X-rays

We present here the properties of the X-COP clusters derived from the X-ray spectral analysis of Ghirardini et al. (2019, hereafter G19) and available at the X-COP website³.

For each cluster, N_bins,X X-ray spectra were extracted from the relative XMM-Newton mosaic: the i-th spectrum is taken in a circular annulus with inner and outer radii ${\widehat{R}}_{\text{i}}$ and ${\widehat{R}}_{\text{i + 1}}$, respectively, assuming as center the peak of the X-ray surface brightness map and out to approximately r₅₀₀ (see G19, for details). G19 fitted the i-th spectrum with the plasma emission code apec (Astrophysical Plasma Emission Code⁴; Smith et al. 2001), as implemented in the software for the spectral analysis XSPEC (Arnaud 1996). For our study, the main results of each spectral analysis are the normalization of the thermal emission Norm and the spectroscopic temperature T_sp. In the model-data comparison, we excluded the bins covering a region partly or totally within 50 kpc from the cluster center, where the presence of a central brightest cluster galaxy and the interplay between gas cooling and heating, not accounted for in our models, could be relevant.

We adopted a Cartesian reference system with origin in the cluster center, where r = (x, y, z) is the position vector in the 3D space, such that the x-axis is along the line of sight. In the apec emission code, Norm is a function of the 3D density distribution of the plasma. In the i-th cylindrical shell, with radii ${\widehat{R}}_{\text{i}}$ and ${\widehat{R}}_{\text{i + 1}}$, and having the x-axis as symmetry axis, we can write it as

$$\begin{array}{c}\hfill Nor{m}_{\mathrm{i}}:=C{\displaystyle {\int }_0^{2\pi }{\int }_{{\widehat{R}}_{\mathrm{i}}}^{{\widehat{R}}_{\mathrm{i}+1}}}\left(2{\displaystyle {\int }_0^{\infty }}{n}_{\mathrm{H}}{n}_{\mathrm{e}}\mathrm{d}x\right)\widehat{R}\mathrm{d}\widehat{R}\mathrm{d}\widehat{\varphi },\end{array}$$(1)

where $\widehat{\varphi }$ is the azimuthal angle in the plane of the sky, n_e(r) and n_H(r) are the electron number density of the ICM and the equivalent number density of the hydrogen nuclei, respectively, and C = 10⁻¹⁴/[4πD_ang²(1+z₀)²], with D_ang and z₀ the angular distance and systemic redshift of the cluster, respectively, is the normalization factor adopted in XSPEC. In our analysis we compared our models with the quantity

$$\begin{array}{c}\hfill {\mathcal{N}}_{\mathrm{i}}=\frac{Nor{m}_{\mathrm{i}}}{\pi [{({\widehat{R}}_{\mathrm{i}+1}/[\mathrm{arcmin}])}^2-{({\widehat{R}}_{\mathrm{i}}/[\mathrm{arcmin}])}^2]},\end{array}$$(2)

whose values are provided by G19.

2.2. ICM thermodynamic properties from the SZ effect

We present here the properties of the X-COP clusters as inferred from the SZ effect signal in the Planck maps by G19 and reported in the X-COP website.

G19 analyzed the Compton y-parameter maps, under the assumption that the relation between the “SZ pressure”, P_SZ, and the y-parameter in Eq. (3) of G19 holds. Radial profiles were extracted in N_bins,SZ circular bins, and then deconvolved by the Planck PSF and geometrically deprojected under the assumption of spherical symmetry. It follows that they recovered the SZ pressure as the average in the j-th spherical shell with radii r_j and r_{j + 1}, P_SZ,j, where j = {1, …, N_bins,SZ} (see Sect. 2.5 of G19, for details). In this analysis, a covariant matrix accounts for the errors and the cross-correlations between the values of P_SZ,j at all j. Given that the diagonal terms of the covariance matrix are dominant, we neglected the off-diagonal terms and considered only the errors from the diagonal terms in the model-data comparison. Given that the sizes of the three central radial bins in P_SZ,j, as widely discussed in G19, are smaller than the PSF of Planck, we excluded them from the model-data comparison. The used P_SZ data range from ∼0.5r₅₀₀ out to ∼3r₅₀₀ in all the clusters.

2.3. Cluster mass distribution

The spherically averaged mass density profile of the observed clusters and those formed in cosmological simulations are routinely described by the virial mass M₂₀₀ and by the concentration c₂₀₀ = r₂₀₀/r₋₂, where r₋₂ is the radius at which the logarithmic slope of the density profile is −2.

We took the distribution of M₂₀₀ measured from WL analysis, when available, as the reference estimate of the cluster mass. Specifically, we took from Col. 4 of Table A.2 of Herbonnet et al. (2020) the average and the 1σ scatter of WL mass estimates for the clusters A85, A1785, A2142, and Zw1215. These values defined our a priori normal distribution of M₂₀₀ for each of these clusters. For the remaining clusters, we assumed the a priori normal distribution of M₂₀₀ as follows. We took the median of the M₂₀₀ from Table 2 in Eckert et al. (2019), who estimated the M₂₀₀ distributions of each cluster of the X-COP sample using a correction for the hydrostatic masses based on the assumption of the universal baryon fraction. We conservatively adopted as the standard deviation (normalized to the median M₂₀₀) σ = (σ_ubf²+σ_WL²)^1/2, where σ_ubf is (for each cluster) the average between the upper and lower relative uncertainties on M₂₀₀ (taking the data from Table 2 of Eckert et al. 2019), and σ_WL is the average relative uncertainty on WL estimates of M₂₀₀ for the aforementioned five clusters with WL measurements.

We assumed a Gaussian prior distribution of lnc₂₀₀. For each cluster the median lnc₂₀₀ is obtained from the median M₂₀₀, using the mass-concentration relation in Eq. (5) of Ragagnin et al. (2021), with coefficients given by their Eq. (7) evaluated for the cosmological parameters adopted in this work. Following Ragagnin et al. (2021), for all clusters we assumed 0.38 as the standard deviation of the prior distribution of lnc₂₀₀.

3. Models for a pressure-supported rotating ICM

We present here two families of axisymmetric cluster models with rotating ICM that we applied to our sample of clusters: in one family the ICM is assumed to have no turbulence, while in the other the turbulence of the ICM is accounted for. Models without turbulence, though idealized, are interesting, because they allow us to estimate the maximum room for rotation in the ICM, and could be not too unrealistic, if the negligible turbulence support, measured in the central regions of the Perseus cluster and A2029 (Hitomi Collaboration 2016; Miller et al. 2025), is confirmed at larger distances from the center and in other clusters.

Both families belong to the class of models presented by B24 (see also Bianconi et al. 2013; Nipoti et al. 2015), in which a pressure-supported and rotating ICM has a composite polytropic distribution (e.g., Curry & McKee 2000) and is in equilibrium in the gravitational potential of a dark-matter halo (the self-gravity of the gas is neglected). While B24 considered also oblate and prolate dark halos, here, for the sake of simplicity, we limited ourselves to the case of spherical dark halos: in particular, we adopted a spherical Navarro-Frenk-White (NFW; Navarro et al. 1996) gravitational potential, which is fully determined by the parameters M₂₀₀ and c₂₀₀ (see, e.g., Sect. 2.2 of B24). We generalized the formalism of B24 by interpreting the gas pressure p as

$$\begin{array}{c}\hfill p=n{k}_{\mathrm{B}}{T}_{\mathrm{gas}}+p_{\mathrm{turb}},\end{array}$$(3)

where n is the gas number density and p_turb is the turbulent pressure, which we parameterized as p_turb = α_turbp, with α_turb = α_turb(p) a dimensionless quantity in the range 0 ≤ α_turb < 1, representing the fraction of turbulent pressure support. The ICM temperature is thus given by

$$\begin{array}{c}\hfill T_{\mathrm{gas}}=(1-{\alpha }_{\mathrm{turb}})\frac{p}{n{k}_{\mathrm{B}}}·\end{array}$$(4)

3.1. Intrinsic properties of the models

Here we describe the intrinsic properties of the models adopting a cylindrical coordinate system (R, ϕ, z). Given that the gas distribution is barotropic, the gas rotation velocity u_ϕ depends only on R: as in B24, we assumed

$$\begin{array}{c}\hfill u_{\varphi }(R)=4u_{\mathrm{peak}}\frac{S}{{(1+S)}^2},\end{array}$$(5)

where S ≡ R/R_peak, R_peak is the radius at which u_ϕ is maximum and⁵u_peak ≡ u_ϕ(R_peak). This functional form of u_ϕ(R) is relatively flexible and is able to reproduce the average rotation velocity profiles of the ICM found in cosmological simulations for clusters with masses comparable to those of X-COP clusters (see B24).

In the composite polytropic models here adopted, the gas number density distribution is

$$\begin{array}{c}\hfill n(R,z)=n_{\star }{[1-\frac{\gamma \prime -1}{\gamma \prime }\frac{n_{\star }\mu m_{\mathrm{p}}}{p_{\star }}\mathrm{\Delta }{\mathrm{\Phi }}_{\mathrm{eff}}(R,z)]}^{\frac{1}{\gamma \prime -1}},\end{array}$$(6)

and the pressure distribution is

$$\begin{array}{c}\hfill p(R,z)=p_{\star }{\left(\frac{n(R,z)}{n_{\star }}\right)}^{{\gamma }^{\prime }},\end{array}$$(7)

with γ′=γ′_IN and γ′=γ′_OUT in the inner and outer regions of the cluster, respectively. The effective potential Φ_eff, depending on M₂₀₀, c₂₀₀, u_peak and R_peak, is defined in Eqs. (24) and (25) of Bartalesi et al. (2024). In the above equations, n_⋆ = n(R_break, 0) and p_⋆ = p(R_break, 0), with R_break the radius in the equatorial plane where the value of the polytropic index changes; μ is the mean molecular weight and m_p is the proton mass. We assumed a position-independent metallicity of 0.3 times the solar value, for which μ = 0.6 and n/n_e = 1.94, or n_e = 1.17n_H in a fully ionized plasma. For models both with and without turbulence, we found it convenient to define the parameter T_⋆ ≡ p_⋆/(n_⋆k_B), which is the gas temperature at (R_break, 0) when α_turb = 0 and can be interpreted as an equivalent gas temperature at (R_break, 0) when α_turb ≠ 0.

In models with turbulent ICM (0 ≤ α_turb < 1) we adopted the following law to describe the dependence of α_turb on p:

$$\begin{array}{c}\hfill {\alpha }_{\mathrm{turb}}(p)=({\alpha }_{\inf }-{\alpha }_0)\frac{ln[1+p/(\xi p_{\star })]}{p/(\xi p_{\star })}+{\alpha }_0,\end{array}$$(8)

where ξ is a dimensionless parameter. In Eq. (8), α₀ and α_inf are the asymptotic values of α_turb for p ≫ ξp_⋆ and p ≪ ξp_⋆, respectively (the transition between the two regimes occurs where the pressure has values around ξp_⋆). Given that, for a pressure-supported cluster, the pressure decreases outward, α_turb increases outward if α_inf > α₀. The contribution to the ICM support against cluster gravity from turbulent motions is routinely measured as a spherically averaged α_turb profile in clusters formed in cosmological simulations. Using a multi-scale filtering technique (see Vazza et al. 2012), which is believed to isolate the uncorrelated component of the velocity field of the ICM, Angelinelli et al. (2020, hereafter A20) accurately estimated the profiles of α_turb for a sample of 68 clusters formed in a cosmological, non-radiative simulation without active galactic nucleus or stellar feedback. They, then, presented a functional form in their Eq. (11) that successfully reproduces the outward increasing radial dependence of the median α_turb profile of these clusters. The functional form of α_turb presented in Eq. (8) is thought to generalize that of A20 to nonspherical analyses, in particular to axisymmetric systems for the scope of this work.

In summary, the family of models without turbulence (α_turb = 0) has 9 free parameters: γ′_IN, γ′_OUT, n_⋆, T_⋆, R_break, M₂₀₀, c₂₀₀, R_peak, and u_peak. The family of models with turbulence has 12 free parameters: α_inf, α₀ and ξ, in addition to the nine parameters of the other family.

3.2. From models to observational data

We define here the quantities, obtained from the intrinsic properties of our models (see Sect. 3.1), that can be compared with the observational data detailed in Sect. 2. We use the Cartesian reference system and the radial grids ${\widehat{R}}_{\text{i}}$ and r_j introduced in Sect. 2.

In the apec emission code, T_sp is formally the temperature of an isothermal emitting plasma (see Sect. 2.1). Given that the ICM in our model is multi-temperature, we adopted, as an approximation of T_sp, the spectroscopic-like temperature T_sl (see Mazzotta et al. 2004, for details) so that we could compare the T_sl derived from our models to T_sp. T_sl in the i-th cylindrical shell, with radii ${\widehat{R}}_{\text{i}}$ and ${\widehat{R}}_{\text{i + 1}}$, is

$$\begin{array}{c}\hfill T_{\mathrm{sl},\mathrm{i}}=\frac{{\displaystyle {\int }_0^{2\pi }{\int }_{{\widehat{R}}_{\mathrm{i}}}^{{\widehat{R}}_{\mathrm{i}+1}}{\int }_0^{\infty }}{n}_{\mathrm{H}}{n}_{\mathrm{e}}{T}_{\mathrm{gas}}^{1/4}\widehat{R}\mathrm{d}x\mathrm{d}\widehat{R}\mathrm{d}\widehat{\varphi }}{{\displaystyle {\int }_0^{2\pi }{\int }_{{\widehat{R}}_{\mathrm{i}}}^{{\widehat{R}}_{\mathrm{i}+1}}{\int }_0^{\infty }}{n}_{\mathrm{H}}{n}_{\mathrm{e}}{T}_{\mathrm{gas}}^{-3/4}\widehat{R}\mathrm{d}x\mathrm{d}\widehat{R}\mathrm{d}\widehat{\varphi }},\end{array}$$(9)

where T_gas is the intrinsic temperature of the ICM. The SZ pressure, obtained from the intrinsic properties of the ICM in our models, is assumed to be the average of the thermal pressure distribution of the ICM in the j-th spherical shell, with radii r_j and r_{j + 1},

$$\begin{array}{c}\hfill P_{\mathrm{SZ},\mathrm{j}}=\frac{3{\displaystyle {\int }_{r_{\mathrm{j}}}^{r_{\mathrm{j}+1}}}{n}_{\mathrm{e}}{k}_{\mathrm{B}}{T}_{\mathrm{gas}}{r}^2\mathrm{d}r}{r_{\mathrm{j}+1}^3-r_{\mathrm{j}}^3},\end{array}$$(10)

where r is the spherical radius and k_B the Boltzmann constant.

In the comparison between the models and the observed dataset, we assumed that our axisymmetric models are observed edge-on, that is, with line of sight x orthogonal to the symmetry and rotation axis z. As we have seen in Sect. 2, the observational data we want to compare with are measurements of quantities in either circular annuli or in spherical shells centered in the cluster center. Though we could perform analogous measurements for our models starting from Eqs. (1), (9), and (10), in order to save computational time in the model-data comparison, we approximated such measurements with estimates obtained by considering only properties of the models in the equatorial plane (z = 0). The details of this approximation, which turns out to be sufficiently good for our purposes, are given in Appendix D. In particular, imposing that the gas density is zero when |x|> r_NSZ or |y|> r_NSZ, where r_NSZ is the outer radius of the outermost spherical shell, we applied the approximations (D.3) and (D.4) to our models. Thus, Eqs. (1) and (9) become, respectively,

$$\begin{array}{c}\hfill Nor{m}_{\mathrm{i}}=4\pi C{\displaystyle {\int }_{{\widehat{R}}_{\mathrm{i}}}^{{\widehat{R}}_{\mathrm{i}+1}}}y\mathrm{d}y\left({\displaystyle {\int }_0^{r_{\mathrm{NSZ}}}}\mathrm{d}x{n}_{\mathrm{H}}(R,0)n_{\mathrm{e}}(R,0)\right)\end{array}$$(11)

and

$$\begin{array}{c}\hfill T_{\mathrm{sl},\mathrm{i}}=\frac{{\displaystyle {\int }_{{\widehat{R}}_{\mathrm{i}}}^{{\widehat{R}}_{\mathrm{i}+1}}}y\mathrm{d}y\left({\displaystyle {\int }_0^{r_{\mathrm{NSZ}}}}\mathrm{d}x{n}_{\mathrm{H}}(R,0)n_{\mathrm{e}}(R,0)T_{\mathrm{gas}}^{1/4}(R,0)\right)}{{\displaystyle {\int }_{{\widehat{R}}_{\mathrm{i}}}^{{\widehat{R}}_{\mathrm{i}+1}}}y\mathrm{d}y\left({\displaystyle {\int }_0^{r_{\mathrm{NSZ}}}}\mathrm{d}x{n}_{\mathrm{H}}(R,0)n_{\mathrm{e}}(R,0)T_{\mathrm{gas}}^{-3/4}(R,0)\right)},\end{array}$$(12)

with $R=\sqrt{x^2+y^2}$; Eq. (10) becomes

$$\begin{array}{c}\hfill P_{\mathrm{SZ},\mathrm{j}}=\eta \frac{3{\displaystyle {\int }_{r_{\mathrm{j}}}^{r_{\mathrm{j}+1}}}{y}^2\mathrm{d}y{n}_{\mathrm{e}}(|y|,0)k_{\mathrm{B}}{T}_{\mathrm{gas}}(|y|,0)}{r_{\mathrm{j}+1}^3-r_{\mathrm{j}}^3},\end{array}$$(13)

where η is a parameter accounting for the possible systematic offset between the electron pressure of the ICM and that measured from the Compton y-parameter maps (see Appendix A).

4. The statistical method

For each cluster of our sample we computed both α_turb = 0 and α_turb ≥ 0 models. Taking mass distributions based on available estimates (see Sect. 2.3), these models, presented in Sect. 3, are tuned to reproduce the radial profiles of 𝒩, T_sp and P_SZ as resolved in the XMM-Newton and Planck maps (see Sects. 2.1 and 2.2) via a Markov chain Monte Carlo (MCMC) method. This section details the statistical method: we present the form of the likelihood and the prior distributions in Sects. 4.1 and 4.2, respectively, and the method used to estimate the posterior distribution of the model in Sect. 4.3. The details of the MCMC algorithm are described in Appendix B.

4.1. Likelihood

Using the models described in Sect. 3, we reconstructed the intrinsic properties of the ICM (density, turbulent and thermal pressures, and rotation speed) and the gravitational potential in which it is in equilibrium. From the model of the intrinsic quantities of the ICM, we evaluated the radial profile of 𝒩, T_sl and P_SZ as described in Sect. 3.2 (we recall that the normalization of the P_SZ profile is regulated by η, which is one of the free parameters of the MCMC). The statistical errors, σ_stat, on the observed 𝒩, T_sp and P_SZ are obtained as the average between the upper and lower 1σ errors reported in the data analyses. We considered three parameters, ϵ_N, ϵ_T, and ϵ_SZ, to account for a systematic contribution, ϵσ_stat with ϵ = {ϵ_N, ϵ_T, ϵ_SZ}, to the uncertainty in the values of 𝒩, T_sp and P_SZ, respectively. ϵ_N, ϵ_T and ϵ_SZ are assumed to have the same values in all the radial bins. Specifically, we assumed each datum of 𝒩, T_sp and P_SZ to be randomly generated from a normal distribution with a median equal to the corresponding 𝒩, T_sl and P_SZ in the model and a standard deviation equal to $\sqrt{1+{ϵ}^2}{\sigma }_{\text{stat}}$. It follows that the generic likelihood of each datum Q_m, im, ℒ(Q_m, im), where Q_m with m = {1, 2, 3} indicates 𝒩, T_sp, or P_SZ, respectively, and i_m the radial bin in which Q_m was extracted, is a normal distribution with the same median and standard deviation as described above. The radial bins used in the fitting are defined in Sects. 2.1 and 2.2.

4.2. Prior distributions

Let us consider the marginal prior of the generic parameter θ_h, π(θ_h), where h = {1, …, H} with H the number of the parameters in the fitting. We found it convenient to use ν_⋆ ≡ log[n_⋆/(10⁻³ cm⁻³)] as the MCMC parameter instead of n_⋆. In addition to the parameters of the models (see Sect. 3.1), for the joint analysis of the X-ray and SZ data we introduced the four parameters η, ϵ_N, ϵ_T and ϵ_SZ (see Sect. 4.1), so for the α_turb = 0 models we have H = 13 parameters (γ′_IN, γ′_OUT, ν_⋆, T_⋆, R_break, log ϵ_N, log ϵ_T, u_peak, R_peak, η, log ϵ_SZ, M₂₀₀, and lnc₂₀₀), while for the α_turb ≥ 0 models we have H = 16 parameters (log ξ, α_inf and α₀, in addition to those of the α_turb = 0 models).

We assumed π(M₂₀₀) and π(lnc₂₀₀) to be Gaussian, with the medians and standard deviations reported for each cluster in Table 1 (see also Sect. 2.3). For all the other parameters we assumed a uniform prior distribution, with upper and lower bounds reported in Table 2. The choice of the upper and lower bounds of the parameters regulating the turbulent support (log ξ, α₀ and α_inf) requires a brief comment. The parameter log ξ sets the distance from the center where there is a transition between an inner region (p ≫ ξp_⋆) in which α_turb → α₀ and an outer region (p ≪ ξp_⋆) in which α_turb → α_inf. We adopted a prior on log ξ uniform in the range [−1.5, 0]. This is motivated by the expectation (based on the properties of the SRM model from B24) that values of log ξ in the range [−1.5, 0] map the radial range (0.15−1)r₂₀₀, which ensures the transition does not occur within the cluster core (log ξ = 0 sets the transition in the α_turb profile at R ≈ R_break in the equatorial plane). For our models applied to the X-COP clusters, a turbulent velocity dispersion higher than σ_RGS = 500 km s⁻¹ at a distance smaller than ≈150 kpc would violate the 90% upper limit on the broadening of the X-ray emitting lines from RGS data. We thus set the upper bound of π(α₀) to α_0,RGS = σ_RGS²/[σ_RGS²+k_BT₀/(μm_p)], where T₀ is the best-fit observed spectroscopic temperature in the innermost annulus of each cluster (taken from G19). α_inf accounts for the turbulent support at large radii, which is observationally unconstrained. We thus assume 0.95 as upper bound of π(α_inf), which corresponds to dominant turbulent pressure in the outskirts (see Sect. 3.1). We set the lower bounds of π(α₀) and π(α_inf) to 0, to conservatively include the possibility of negligible turbulence.

4.3. Posterior distributions

According to Bayes’ theorem, the posterior distribution of the parameter vector θ = (θ₁, …, θ_H), 𝒫(θ), is obtained from $ln\mathcal{P}={\displaystyle \sum _{h=1}^H}ln[\pi ({\theta }_h)]+{\displaystyle \sum _{m=1}^3\sum _{i_m}}ln[\mathcal{L}(Q_{i_m})]$, where π(θ_h) and ℒ(Q_m) are described in Sects. 4.1 and 4.2, respectively. The MCMC algorithm individually maximizes ln𝒫 of each cluster. When it converges to a stationary posterior, we obtain the discrete posterior consisting of a sample of the parameter vector (see Appendix B for details). The relative frequency of the values of the parameters in this sample can be visualized in corner plots, which we present for the parameters M₂₀₀, lnc₂₀₀, u_peak, R_peak and α₀ (if α_turb ≠ 0) for all clusters in Appendix E.

We chose to characterize the marginal posterior distributions using the 68% HPDI as a credible interval (see Sect. 1), which is very similar to the interval between the 16th and 84th percentiles for symmetric distributions, but tends to be more meaningful than the 16th and 84th percentiles for highly skewed distributions. For the sake of conciseness, in the following we refer to the lower and upper bounds of the 68% HPDI credible interval on the marginal posterior distribution on a parameter, simply as the lower and upper limits on that parameter.

For each parameter vector belonging to our sampling of the posterior distribution (see above), we computed a few relevant quantities for each cluster, needed for the interpretation of the results presented in Sect. 5. In particular, we derived the profiles of 𝒩, T_sl and P_SZ as described in Sect. 3.2, and the rotation curve using Eq. (5). The 1D velocity dispersion of the gas particles at a given point in the equatorial plane is given by ${\sigma }_{\text{gas,1D}}(R,0)=\sqrt{p(R,0)/\left[\mu m_{\text{p}}n(R,0)\right]}$. As a measure of the dynamical importance of rotation, we computed, as a function of radius R, the ratio between the gas rotation velocity u_ϕ(R) and the gas velocity dispersion σ_gas,1D(R, 0) in the equatorial plane, u_ϕ/σ_gas,1D. Given that we are interested in the effect of the turbulence when fitting the thermodynamic data with a rotating ICM model, in the α_turb ≥ 0 models, we also derived the α_turb(p) profile through Eq. (8) once we had obtained the pressure profile (see Sect. 3.1). The 1D turbulent velocity dispersion in the equatorial plane is then given by ${\sigma }_{\text{turb,1D}}(R,0)=\sqrt{{\alpha }_{\text{turb}}(p)}{\sigma }_{\text{gas,1D}}(R,0)$, where p = p(R, 0).

In Appendix E we show for each cluster the profiles of 𝒩, T_sl, and P_SZ for the α_turb = 0 model. The 𝒩, T_sl and P_SZ profiles of the α_turb ≠ 0 models, not shown, are virtually indistinguishable from the corresponding profiles of the α_turb = 0 models. We also show for each cluster the u_ϕ profile of the α_turb = 0 model, and the u_ϕ/σ_gas,1D and α_turb profiles of the α_turb ≥ 0 model.

5. Results

We present here the results obtained by applying our models, presented in Sect. 3, to the observed thermodynamic profiles of the X-COP objects through the statistical method described in Sect. 4, focusing in particular on the radial distribution of the ICM rotation and turbulence.

5.1. Models of rotating ICM without turbulence

Our models are able to closely reproduce the observed properties of the ICM (Sect. 2) of all the clusters of our sample from a radius comparable to the brightest cluster galaxy size out to ∼3r₅₀₀, the outermost radial bin in the PlanckP_SZ profile (see Appendix E). The residuals between the median of the model and the data are, as a rule, not larger than 6% in the radial profiles of 𝒩, T_sp and P_SZ (the only exceptions are A644 and A2319, for which the residuals are between 6% and 10%). The median value of ϵ_N inferred in the sample is 4.5 (the median ϵ_N in each cluster ranges from 1, in Zw1215, to 10, in A85). The values of ϵ_T and ϵ_SZ are lower than 1 in all the clusters of the sample, apart from the ϵ_T in A644, A1644, and A2319, whose medians range from 2 to 4. The median errors inferred as $\sqrt{1+{ϵ}^2}{\sigma }_{\text{stat}}$ are plotted in the thermodynamic profiles in Appendix E.

For all the clusters, the posterior distribution 𝒫 of u_peak is unimodal, with a credible interval that encloses the mode (see the corner plots in Appendix E). In most cases, the range between the lower and upper limits of u_peak is wide, implying that the uncertainty on the inferred values of u_peak is large. The posterior lower limit of u_peak is very close to the lower bound of π(u_peak), except for three clusters (A1795, A2255, and A2319).

Considering our full sample of objects, we characterized the distribution of the profiles of any quantity, say q(R), by defining 1σ lower and upper limits of as follows. At any radius R in the equatorial plane, we define the 1σ lower limit as the median of the 16th percentiles of the distribution of q evaluated at R in each cluster model. The 1σ upper limit is defined in the same way, but by taking the 84% percentiles.

Figure 1 summarizes the rotation properties of our models of the clusters of our sample. In the left panel, we show the median u_ϕ/σ_gas,1D profiles of the clusters (see Sect. 4.3). A common trend emerges: this population of clusters exhibits a slightly outward-increasing median profile of u_ϕ/σ_gas,1D. However, the lower limit of the population (lower bound of the gray band in the plot) has everywhere u_ϕ/σ_gas,1D < 0.2, corresponding to ρu_ϕ²/p < 0.04, which indicates a negligible rotation support. In the right panel, where we show the median profiles of the rotation velocity in the cluster models, the lower limit of the population (lower bound of the gray band in the plot) corresponds to velocity everywhere below 100 km s⁻¹. The population upper limit of the u_ϕ/σ_gas,1D profile (upper bound of the gray band in the left panel) significantly increases outward up to nearly 0.9 (corresponding to a rotation velocity of ≈600 km s⁻¹; see the right panel).

Taking the difference between the population upper and lower limits as a measure of the 1σ scatter, we find that the scatter in u_ϕ/σ_gas,1D and in u_ϕ is up to 0.7 and 500 km s⁻¹, respectively, implying that the determination of the properties of the population in our models is very uncertain. Three clusters have median u_ϕ/σ_gas,1D and u_ϕ profiles significantly higher than the 1σ upper limit of the population: A1795 and A2319 have high u_ϕ/σ_gas,1D and u_ϕ in the outskirts (with u_ϕ higher than 1000 km s⁻¹), while A2255 at intermediate radii (with u_ϕ ≈ 500 km s⁻¹).

We conclude that, if turbulence in the ICM is negligible throughout the cluster (which is an hypothesis consistent with the currently available measurements in the central regions of the Perseus cluster and A2029; Hitomi Collaboration 2016; Miller et al. 2025), the current data of the X-COP clusters leave room for an important rotation in the ICM (especially in the three clusters A1795, A2319 and A2255), even if in most cases they are also consistent within the 68% credible interval with negligible rotation of the ICM.

In the next section we consider the case in which the contribution of ICM turbulence is accounted for in the models.

5.2. Models of rotating ICM with turbulence

When applying the α_turb ≥ 0 model to the clusters of our sample, we find that the posterior distributions of the parameters γ′_IN, γ′_OUT, ν_⋆, R_break, log ϵ_N, log ϵ_T, η, log ϵ_SZ, M₂₀₀ and lnc₂₀₀ are very similar to the case of the α_turb = 0 model. We find instead significant differences in the posterior distributions of the parameters T_⋆, u_peak and R_peak. This confirms the expectation that the X-ray and SZ data alone are insufficient to distinguish between the two forms of nonthermal support considered in our study. However, the functional forms adopted for the radial profiles of rotation and turbulence play a pivotal role in separating the rotational contribution to nonthermal support from that from turbulence. In all the clusters, the upper limit of u_peak for α_turb ≥ 0 tends to be lower than for α_turb = 0, which indicates that rotation is less important when the ICM is turbulent. This is expected, because turbulence contributes to the pressure support against gravity, thus reducing the room for rotational support.

Figure 2 summarizes the properties of the α_turb ≥ 0 models for the clusters of our sample. As shown in the upper-left panel, when α_turb ≠ 0 the population upper limit on u_ϕ/σ_gas,1D (upper bound of the gray band in the plot) is ≈0.6 (corresponding to a rotation speed of ≈500 km s⁻¹; see the upper-right panel), smaller than ≈0.9 found when α_turb = 0, but still indicating room for significant rotation support, even in the presence of turbulence. Conversely, as shown in the lower panels of Fig. 2, the 1σ scatter of the turbulent-to-total-pressure ratio (α_turb, in its left panel) and 1D turbulence velocity dispersion (σ_turb,1D, in its right panel) in the population ranges from ≈0.05 and 200 km s⁻¹ up to 0.4 and 500 km s⁻¹, characterized by a very significant rise of population upper limit of α_turb in the outskirts.

It is particularly interesting to compare α_turb = 0 and α_turb ≥ 0 models of the clusters whose α_turb = 0 model shows significant rotation, namely A1795, A2319, and A2255. The models of A1795 and A2319 have similar properties. The α_turb = 0 models of both objects have a posterior lower limit of u_peak around 850 km s⁻¹, much higher than the lower bound of the prior of u_peak (see the corner plot in Figs. 4 and 7 of Appendix E, respectively). This corresponds to a median u_ϕ/σ_gas,1D exceeding 1 in the outskirts and thus significantly deviating from the average properties of the population (see the right panel of Fig. 1). Instead, when α_turb ≥ 0, the models of A1795 and A2319 have posterior lower limit of u_peak coincident with the lower bound of π(u_peak) (50 km s⁻¹). This translates into a median u_ϕ/σ_gas,1D profile, with a value below 0.4 in the outskirts, that everywhere falls within the 1σ interval of the population (see the upper-left panel of Fig. 2). The lower-left panel of Fig. 2 shows the median turbulent-to-total pressure ratio profile α_turb(R) of our clusters. A1795 and A2319 have α_turb significantly increasing outward, with α_turb > 0.5 in the outskirts, standing out from the 1σ band of the considered cluster population. We can interpret this as a significant shift of the nonthermal support from rotation to turbulence in the outer regions.

In the α_turb ≥ 0 model, A2255 exhibits the highest median u_ϕ/σ_gas,1D profile in a region close to the center, making it the most promising candidate for detecting rotation in the ICM with XRISM/Resolve. The credible interval of 𝒫(u_peak) in the α_turb ≥ 0 model consists of two separate intervals corresponding to relatively low and high values of u_peak, centered at 100 km s⁻¹ and 600 km s⁻¹, respectively. 𝒫(u_peak) is relatively flat, from the lower bound of the low-u_peak credible interval (that is the lower bound of the prior) to the upper bound of the high-u_peak credible interval, which is approximately equal to the 1σ upper limit of u_peak in the α_turb = 0 model. Given the relatively high rotation velocity of the ICM in the inner or intermediate regions, this cluster warrants further investigation through a dedicated XRISM pointing in its center. In Sect. 5.4 we propose an observational strategy with XRISM/Resolve to validate the rotation peak of 500 km s⁻¹ at a distance of 200 kpc from the center.

In Appendix C we compare the flattening of the X-ray surface brightness distribution measured in our models with the measurements of the X-COP clusters. We found that the α_turb = 0 models (being more flattened) tend to perform better than the α_turb ≥ 0 models in reproducing the observed distribution of the axial ratios estimated from the X-ray maps. This might indicate that the ICM turbulence support is not as high as inferred with our analysis. However, it must be stressed that in our models, where the gravitational potential is spherically symmetric, only the rotation determines the flattening of the X-ray surface brightness distribution. The flattening of the ICM can also be influenced by the shape of the dark-matter halo, which can deviate significantly from spherical symmetry (e.g., Allgood et al. 2006). This is widely discussed in B24, where models of the rotating ICM in equilibrium in both prolate and oblate halos were studied.

5.3. Comparison with the published studies on ICM motions

To our knowledge, the present work is the first attempt to reproduce the temperature, density and pressure profiles of the observed clusters with models with rotating and turbulent ICM, consistent with the available mass estimates. Here we compare our results with published works presenting constraints on the contribution of nonthermal pressure support from observed datasets, where the nonthermal component is totally ascribed to turbulence, and on the turbulent support in clusters formed from cosmological hydrodynamic simulations.

Ettori & Eckert (2022, hereafter E22) present constraints on the contribution of a nonthermal component to the total pressure, which is assumed to be in hydrostatic equilibrium in the gravitational potential of the X-COP clusters. The nonthermal pressure is parametrized as a power-law dependence upon the gas density with the normalization and the slope as free parameters. These free parameters are constrained by requiring that the gas mass fraction at an overdensity of 500 and 200 ρ_crit matches a “universal” value, which was derived from cosmological simulations, accounting for the baryonic mass fraction and subtracting a statistical contribution for stellar mass (see also Eckert et al. 2019). By construction, the model is sensitive to the regions where a gas mass fraction measurement (based on the initial assumption that the total pressure support is only thermal) is available, and it can only account for nonthermal pressure support in systems where the observed gas fraction exceeds the universal value.

For the sake of simplicity, we compare our α_turb = 0 rotation velocity profile, which contains the entire information on the nonthermal support, with the 1D turbulent velocity dispersion profiles obtained in E22 for the same cluster (see their Eq. 21). In the regions where the gas mass fractions are available (r ∼ r₅₀₀ − r₂₀₀), the 1D turbulent velocity dispersion profiles of E22 for most clusters lies within the range between the 16th and 84th percentiles of our reconstructed rotation velocity profiles. However, in A1795 the 1D turbulent velocity dispersion profile from E22 is significantly lower than the 16th percentile of our rotation velocity profile. Although a detailed investigation of the cause of this discrepancy is beyond the scope of this comparison, we suggest a potential inconsistency in the mass estimates⁶ rather than issues related to the different modeling. Overall, in E22, A2255 and A2142 show a central 1D velocity dispersion of about 400 km s⁻¹, and are therefore promising candidates for a detectable ICM rotation. While A2142 was found with low nonthermal support in the central regions in our model (see the left panel of Fig. 1), we confirm the high level of nonthermal support near the center of A2255.

In B24, we compared the rotation curves of our models with those recovered from some cosmologically simulated clusters. Given that the inferred constraint on the turbulence in our models impacts the inferred rotation, we investigate here how our constraints on the turbulent component of the ICM velocity field relate to the predictions of A20, who accurately recovered the turbulent part of the ICM velocity field in their simulated sample (see Sect. 3.1, for details). In the lower-left panel of Fig. 2, we overlay the median profile of α_turb from A20. This average profile is well in agreement with our median results, lying within the interval between the 16th and 84th percentiles of the distribution of the α_turb profiles in the X-COP population from ∼50 kpc out to ∼3r₅₀₀, although local tensions in specific cases are expected and observed (see the plots in Appendix E).

Fig. 1. Equatorial plane velocity profiles of the ICM in the αturb = 0 models of the clusters of our sample. Left panel: Ratio between the gas rotation velocity and velocity dispersion (measured in the equatorial plane) as a function of the cylindrical radius normalized to the virial radius, r200 = 3M200/[800πρcrit]1/3, computed using the median value of M200 for each cluster. Right panel: Gas rotation velocity as a function of the cylindrical radius in physical units. In both panels the dashed lines are the median profiles of the individual clusters (as indicated in the legend), whereas the gray band represents the interval between the 16th and 84th percentiles of the distribution of the profiles of the entire sample of clusters.

Fig. 2. Equatorial plane rotation and turbulence profiles of the ICM in the αturb ≠ 0 models of the clusters of our sample. Upper panels: Same as Fig. 1 but for models with αturb ≥ 0. Lower left panel: Median turbulent-to-total pressure ratio, αturb, as a function of the cylindrical radius normalized to the virial radius, r200 (see the caption of Fig. 1), for the same models as in the upper panels. Lower right panel: 1D turbulent velocity dispersion as a function of the cylindrical radius in physical units for the same models as in the upper panels. The color coding and the meaning of the gray bands are the same as in Fig. 1. In the lower right panel, we overplot the αturb profile measured in simulations by A20 (solid black line; the profile was taken from their Eq. (11), with a0, a1, and a2 equal to their best-fit values, as reported in the first row of their Table 1, and assuming r200, m/r200 ≈ 1.67).

5.4. Detecting the ICM rotation with XRISM: Abell 2255

In this section we study the capability of high-resolution (eV-scale) X-ray spectroscopy to detect and characterize the rotation of the ICM. In particular, the X-ray calorimeter Resolve on board the XRISM mission⁷ is the only currently available instrument that can offer such a spectral resolution; however, it has a limited collecting area and FoV, meaning it can only efficiently capture the brightest regions. In B24, we studied the significance of the measurement of the rotation through the estimates of the centroid of the iron emission lines, the impact of the line broadening on these estimates and its cross-correlations with the other fitting parameters in the analyses of a set of XRISM/Resolve mock spectra generated from three realistic (but generic) rotating ICM models. Here, we present a feasibility study applied to A2255, the most interesting target because of its highest median rotation velocity (≈500 km s⁻¹ as inferred from the α_turb ≠ 0 model) near the center in the considered sample of objects.

A2255 is a disturbed, non-cool-core system at redshift z₀ = 0.0809, settled in the North Ecliptic Pole SuperCluster. It shows a flat X-ray surface-brightness profile in its core and an X-ray morphology extended along the E–W axis (Sakelliou & Ponman 2006), which could indicate the presence of rotation. Spectral analysis of XMM-Newton data shows that the temperature distribution presents azimuthal variations. These characteristics were interpreted as due to a shock moving along the E–W direction with Mach number ∼1.24 and with a core crossing occurred ∼0.15 Gyr ago (Sakelliou & Ponman 2006). Also, some properties inferred from multiwavelength observations, such as the fact that the X-ray centroid of the cluster does not coincide with the position of the brightest galaxy, are consistent with the merger scenario (Burns et al. 1995). While there is currently no solid evidence of a radio counterpart of E–W shock, many relics were detected in the outskirts (some seen in projection), and a radio halo was detected in the center (Pizzo et al. 2008; Pizzo & de Bruyn 2009). Recent deep radio observations of A2255 suggested that very composite mechanisms, including seed particles injected by radio galaxies and spread in the ICM by turbulent motions and weak shocks, are at work at the same time to generate the radio halo (Botteon et al. 2020, 2022). Interestingly, Roettiger & Flores (2000), using a hydrodynamical N-body simulation, found that off-axis mergers are an efficient mechanism for transferring angular momentum to the ICM, which suggests that large-scale rotation of the ICM could be present in post-merging objects (see also Roettiger et al. 1998).

We simulate a pointing with XRISM/Resolve centered on the core of A2255 and including a rotation speed of 500 km s⁻¹. When we extract a spectrum from the entire FoV of Resolve, consisting of 6 × 6 pixels covering an area ≈3 × 3 arcmin², and affected by a PSF with a half-power diameter of about 1.3 arcmin (≈130 kpc, at the redshift of A2255), both turbulent and rotational motions would contribute to the measured broadening of the X-ray emission lines. To associate the rotation with the shifting of the centroids of the X-ray emitting lines and decouple it from the turbulent motions in the ICM, we consider the following observational strategy. We assume that, if the ICM rotation is present in the Resolve exposure, it can be resolved by looking at the different redshift measured in two complementary spectra extracted from the single observation. To evaluate this, we divide the FoV into two complementary sectors of 3 × 6 pixels. We simulate the ICM emission using the thermal plasma emission model bapec⁸ with a 6.2 keV temperature, the value spectroscopically measured with an XMM exposure of the central regions (see G19) and abundance equal to 0.3 times the solar abundance in Anders & Grevesse (1989) (see Ghizzardi et al. 2021, for details). The normalization is set to the fraction (∼20%) of the 0.5−2 keV band luminosity estimated within the hydrostatic r₅₀₀ that falls within the Resolve FoV using a β model (Cavaliere & Fusco-Femiano 1978) with parameters the scale radius r_c = 0.15r₅₀₀ and the exponent β = 0.65. Given that we are considering half of the Resolve FoV, the normalization of the bapec model is defined as half of the estimated flux. In the simulated spectrum, the nominal redshift z₀ is modified as $z_{\mathrm{u}}=(1+z_0)\sqrt{(1+u_{\text{peak}}/c)/(1-u_{\text{peak}}/c)}-1$ (Roncarelli et al. 2018), where c is the light speed. Given that the peak of the rotation lies within the 3 × 6 sector, in this equation we considered the peak rotation velocity, u_peak (positive for a receding ICM). For z₀ = 0.0809 and u_peak ∼ ±500 km s⁻¹, z_u = 0.0791 and 0.0827, respectively. The most prominent emission line complex, key to accurately recovering the redshift of the ICM, is that of iron (Fe XXVI) at rest-frame energy equal to 6.7 keV. The centroid of it should be shifted, in the observer’s frame, to (6.188, 6.199, 6.209) keV for a rotation of +500, 0, −500 km/s, respectively. These values correspond to differences in the locations of the Fe XXVI lines of ΔE ≈ 10 − 21 eV that cannot be resolved with CCD-like spectral resolution (like for XMM), but are well above the Resolve spectral resolution and accuracy of about 4.8 and 0.2 eV, respectively. We include a broadening of the X-ray emission lines equal to the median turbulent velocity dispersion (≈300 km s⁻¹) inferred from our model. All the spectra are absorbed by the expected Galactic column density of 2.5 × 10²⁰ cm⁻². This parameter is kept fixed in the following spectral analysis. Using the responses for the closed Gate-Valve configuration, and adding an instrumental and a cosmic background and a Poissonian noise, we simulate in XSPEC the 1.8−10 keV spectrum of the rotating ICM, as shown in Fig. 3. Importantly, the instrumental background is significantly lower than the thermal model, and the cosmic background adds a relatively small contribution to the flux of the 6.1−6.2 keV emission line complex. We proceed by fitting the broadband (1.8−10 keV) spectrum, fixing the best-fit value of the temperature, and fitting the region around the Fe XXV–XXVI complex again. We obtain that a minimum exposure of 100 ks is required to reach accuracy and precision of < 1% on the measurement of the redshift, allowing us to identify the centroid of the redshifted emission lines at more than 5σ level of confidence for a rotation greater than ±400 km s⁻¹ (ΔE ≈ 17 eV), and at > 3σ for a rotation down to ±200 km/s (ΔE ≈ 8.2 eV). The temperature, line broadening and metallicity are consistently measured with an accuracy on the order of 10−20%. The PSF is not an issue when determining the centroid of the emitting lines in two different sectors. Moreover, based on the analysis of the spectral simulations presented in B24, no significant cross-correlation between the bapec parameters affects our reconstruction, which assures us that it is possible to find evidence of a rotating ICM in A2255.

Fig. 3. Simulated spectrum (black crosses), with an exposure of 100 ks, in half of the Resolve FoV with the best-fit bapec model overplotted. The thermal model only is shown with the green line, the instrumental background in blue, and the thermal model plus the instrumental and cosmic background in red. The most prominent lines are the Fe XXV and Fe XXVI complexes in the range 6−7 keV.

6. Conclusions

We have presented equilibrium models with a rotating ICM for 11 galaxy clusters from the X-COP sample. In the models, based on those of B24, the gas has an axisymmetric composite polytropic distribution, in equilibrium in the gravitational potential of a spherical NFW dark-matter halo. For each cluster we considered both the case of a turbulent ICM (α_turb ≥ 0) and the case of an ICM with negligible turbulence (α_turb = 0). The profile of rotation velocity and the distribution of turbulent velocity dispersion are described with flexible functional forms that are able to reproduce the radial profiles obtained from the spherically averaged analyses of clusters formed in cosmological simulations.

We considered the observed profiles of the normalization of the thermal continuum (which scales with the gas density squared), of the spectroscopic gas temperature, and of the SZ-derived pressure as obtained by G19 using XMM-Newton and Planck measurements. We constrained the cluster gravitational potentials using estimates of the virial mass inferred either from the WL analysis by Herbonnet et al. (2020) or under the assumption of the “universal” value of the baryon fraction from Eckert et al. (2019). For each cluster in our sample, we used an MCMC method to infer the posterior distributions of the parameters of the models using observational data on the thermodynamic profiles and the aforementioned constraints on the gravitational potential. In the comparison with the observational datasets, we assumed that we were observing our cluster models with a line of sight orthogonal to the rotation and symmetry axis.

Rotation of the ICM has not been detected so far, but rotation-induced shifts in the centroids of X-ray emission lines could be measured with the currently operating XRISM/Resolve X-ray spectrometer. To assess the feasibility of such rotation-velocity measurements, for our rotating-ICM model of the cluster Abell 2255, we computed and analyzed a mock spectrum of a central exposure with XRISM/Resolve. Our main findings are the following:

• Our models closely reproduce the observed radial profiles of the thermodynamic quantities of the ICM with median residuals of ≈5%.

• According to our α_turb = 0 models, there is room for non-negligible rotation in the ICM of massive clusters. In our sample of clusters, the median rotation-velocity-to-velocity-dispersion ratio, u_ϕ/σ_gas,1D, increases from ≈0.1 at R = 0.1r₂₀₀ to ≈0.4 at R = r₂₀₀. Significant deviation from this trend is found in Abell 2255, with a median u_ϕ/σ_gas,1D of ≈0.5 (u_ϕ ≈ 600 km s⁻¹) at 0.1r₂₀₀.

• The room for ICM rotation is also significant when the contribution of turbulence to the pressure support is accounted for. On average, the α_turb ≥ 0 models of the clusters of our sample have outward-increasing median u_ϕ/σ_gas,1D profiles going from ≈0.1 at 0.1r₂₀₀ to ≈0.3 at r₂₀₀. A2255 has a median u_ϕ/σ_gas,1D of ≈0.4 (u_ϕ ≈ 500 km s⁻¹) at 0.1r₂₀₀, making it a very promising candidate for detectable rotation in the ICM.

• We show that in A2255 a rotation velocity of 400 km s⁻¹ at a distance of 100 kpc from the center (as predicted by our model) can be detected at more than a 5σ level of confidence with a 100 ks central exposure with XRISM/Resolve.

The analysis presented in this paper has some limitations. For instance, we neglected the effect of bulk motions, any ICM inflows and outflows, and the relative motion between the ICM and the dark-matter distributions, commonly known as sloshing. Moreover, in the presence of sharp temperature variations along the line of sight, which might be the case in mergers, the estimate of the temperature based on the spectroscopic-like temperature may be inaccurate. This can bias the reconstruction of the cluster mass profile in our model and by extension its rotation curve. Moreover, in our model we are working under the assumption that the peak of the ICM emission coincides with the minimum of the cluster gravitational potential. Any offset between the two in real clusters can bias the mass profile reconstruction.

Nonetheless, as reported in Sect. 2, by construction the X-COP clusters lack evidence of large-scale morphological disturbances (for instance, as induced by major mergers and coherent bulk motions) both in their X-ray emission maps and in the spectral analysis, and are the closest representation of a symmetric ICM in equilibrium with the underlying potential in the nearby Universe (Ghirardini et al. 2019). This guarantees an accurate reconstruction of the gas density and temperature distribution in these systems through an azimuthally averaged analysis, which also included a treatment of some possible residual systematic effects that were considered as an additive term in the error budget and radially constant. However, disturbances limited to local regions, on scales of a few hundred kiloparsecs or less, are observed: for instance, the centroid of the X-ray emission in A2255 differs from the X-ray peak significantly (the value of the centroid shift is the second highest one, just below what is measured in the well-known merging object A2319; see Table 2 in Eckert et al. 2022b). This effect is partially mitigated by using large radial bins in our analysis, and by excluding the central regions where this offset should have the largest impact on the reconstructed mass profile.

Additional sources of systematics can come from the inaccuracy or non-universality of the “universal value” of the baryon fraction (see Sect. 5.1 of Eckert et al. 2019, for a discussion). This could bias the marginal priors on M₂₀₀ and c₂₀₀ and, consequently, the inferred rotation curve of the ICM. However, we conservatively added in quadrature the typical WL error to that of these mass estimates.

In parallel to the rotation of the ICM, the rotation of the cluster galactic component has also been an object of study. Hwang & Lee (2007) searched for the presence of rotation in the galactic components of a sample of 56 clusters. They identified 12 clusters exhibiting observational characteristics (such as a rotation-velocity-to-velocity-dispersion-ratio between 0.5 and 0.7 and a significant velocity gradient across the entire plane of the sky) that a dynamically significant rotation would induce. However, none of these clusters are included in the X-COP sample. Stronger indications of the rotation of member galaxies have been found in the inner regions of the clusters Abell S1013, MACS J1206.2−0847 and A2107, with rotation-velocity-to-velocity-dispersion ratios of ≈0.15 (Ferrami et al. 2023) and ≈0.6 (Song et al. 2018), respectively. It is worth mentioning that the ICM can have a non-negligible rotational support even in the absence of dynamically unimportant rotational motions in the galactic component, as shown in the simulated cluster discussed in Roettiger & Flores (2000), where the collisional ICM significantly rotates whereas the collisionless galactic and dark matter components do not. Robust determination of the angular momentum in the galactic and hot-gas components of galaxy clusters, from optical and X-ray spectroscopy, respectively, would provide further constraints on models that describe the formation and evolution of cosmic structures.

Improved spectral (and spatial) resolution in the X-ray band, as expected with the proposed European (NewAthena⁹) and Chinese (HUBS¹⁰) missions, will open a new era in the study of the ICM velocity field (see Roncarelli et al. 2018; Clerc et al. 2019; Zhang et al. 2024). At the same time, more sophisticated models will be needed to reproduce, simultaneously, the radial profiles of the thermodynamic quantities and the shifting and broadening of the emitting lines of the ICM as resolved spatially and spectroscopically in the X-ray bands. These models will provide robust assessments of the contribution of the bulk motions to the energy budget of the ICM and of the underlying mass distribution of galaxy clusters, with important implications for our knowledge of the thermal balance of the ICM and cosmological structure formation.

Data availability

Appendix E is available on the Zenodo repository at the link https://doi.org/10.5281/zenodo.15037889

Acknowledgments

We thank the anonymous referee for very useful comments that helped improve the article. We thank Dominique Eckert for helpful discussions on data analysis. T.B. acknowledges funding from the European Union NextGenerationEU. S.E. acknowledges the financial contribution from the contracts Prin-MUR 2022 supported by Next Generation EU (n.20227RNLY3 The concordance cosmological model: stress-tests with galaxy clusters), ASI-INAF Athena 2019-27-HH.0, “Attività di Studio per la comunità scientifica di Astrofisica delle Alte Energie e Fisica Astroparticellare” (Accordo Attuativo ASI-INAF n. 2017-14-H.0), and from the European Union’s Horizon 2020 Programme under the AHEAD2020 project (grant agreement n. 871158).

References

Appendix A: Offset between X-ray and SZ data

We present here the method used to account for the possible systematic offset between X-ray and SZ data. The joint fitting to X-ray and SZ data was performed under the assumption of a nuisance parameter (e.g., G19) η = P_SZ/(n_ek_BT), where n_e and T are inferred from the fitting to 𝒩 and T_sp data (see Sect. 2.1). For each cluster, we inferred the marginal posterior of η: the corresponding 16th-84th percentile intervals are shown in the lower panel of Fig. A.1. In the upper panel of Fig. A.1, we plot the distribution of η in the entire sample¹¹ obtained for the α_turb ≥ 0 models: we inferred a median value of 1.03 and a 1σ scatter¹² of 0.05. The result for the α_turb = 0 model is similar, with a median value of 1.02 and a 1σ scatter of 0.06. η = 1 is the expected value when the clusters have random orientations in the plane of the sky and random shapes of the 3D ICM distribution (e.g., Kozmanyan et al. 2019). The fact that in both models η = 1 lies in the interval between the 16th and 84th percentiles of the distribution of η inferred in the entire sample is an indication of the goodness of our joint analysis.

Fig. A.1. Upper panel: Probability density function (pdf) of the systematic offset parameter, η, in our cluster sample. The gray band ranges between the 16th and 84th percentile of η in the entire sample. Bottom panel: Median η of each cluster (as indicated in the vertical axes) according to the αturb ≥ 0 model, with error bars corresponding to the 16th and 84th percentiles. For reference, η = 1 (no offset) is indicated with a red vertical line.

Appendix B: Set-up of the Markov chain Monte Carlo method

Using the emcee¹³ sampler (Foreman-Mackey et al. 2013), we ran the Markov chain and drew the samples of the parameter vector for each cluster model from the posterior 𝒫. Empirically, we noticed that a combination of Nelson et al. (2013) and Goodman & Weare (2010) moves, as implemented in the emcee sampler, reduces the autocorrelation between successive iterations and optimizes the exploration of the parameter space. We ran the MCMC for 60000 iterations with 500 chains. It is standard to consider the end of the burn-in phase when 𝒫 is found to be stationary as the iteration number increases. The burn-in phase is longer when α_turb ≥ 0 than when α_turb = 0, but as a rule it lasts less than 5000 steps. In few cases the number of steps of the burn-in phase is higher: ≈15000 for both models of A2319. To be conservative, in all cases we discarded the first 40000 steps. We then applied a thinning of 700 steps to build the posterior sample.

Appendix C: Flattening of the X-ray surface-brightness distributions

Most X-ray surface brightness distributions of the X-COP clusters as resolved in the XMM-Newton maps exhibit clear evidence of a deviation from circular symmetry. Because of the presence of rotation, our models of the X-COP clusters have oblate axisymmetric density distributions, even if the gravitational potential is assumed spherically symmetric. As we assumed that the models are observed edge-on, their surface brightness distributions are flattened. Here we compare the flattening of the models with that observed in the XMM maps.

For the models we evaluated the surface-brightness distribution in each point of the plane of the sky (y, z) as

$$\begin{array}{c}\hfill \mathrm{\Sigma }(y,z)=2\mathrm{\Lambda }{\displaystyle {\int }_{|y|}^{r_{\mathrm{NSZ}}}}\frac{n_{\mathrm{i}}(R,z)n_{\mathrm{e}}(R,z)R\mathrm{d}R}{\sqrt{R^2-y^2}},\end{array}$$(C.1)

where n_i = n − n_e and Λ is the cooling function, which we assume to be independent of temperature (and then of position), as appropriate when considering the soft X-ray band typically used for imaging. We limited ourselves to the region in the plane of the sky, defined by 1 kpc $\le |y|\le {\widehat{R}}_{\text{NX}}$ and 1 kpc $\le |z|\le {\widehat{R}}_{\text{NX}}$, where ${\widehat{R}}_{\text{NX}}$ is the outer radius of the N_bins,X-th radial bin in the plane of the sky (see Sect. 2.1).

Following B24 (see their Sect. 4.2), we use as measure the flattening of the X-ray 2D surface brightness distribution the axial ratio of isophotes, ζ, computed from the eigenvalues of the 2 × 2 matrix containing the centroids in y, z and yz weighted over the surface brightness. For each cluster model, we randomly extracted 100 parameter vectors from the posterior sample and computed ζ for each parameter vector.

In the upper panel of Fig. C.1, we compare the distribution of ζ in our models to those directly measured in the real clusters. The ζ distribution of the α_turb = 0 models (with median 0.91) is more similar to the observed distribution of ζ (with median 0.90) than that of the α_turb ≥ 0 models (with median 0.97). However, we note in the α_turb = 0 distribution a non-negligible tail at ζ < 0.75 with no counterpart in the observed distribution. The lower panel of Fig. C.1, where we compare the ζ measurements in the models individually with the observed ζ of each cluster, shows that the models of A1795 and A2319 perform better when α_turb ≥ 0 than when α_turb = 0, while for the other clusters the intervals between the 16th and 84th percentiles in the two models have substantial overlap. In A644 and A2142, neither the α_turb = 0 nor α_turb ≥ 0 models reproduce the observed ζ at the 1σ confidence level. This might indicate that, in addition to the rotation of the ICM, other properties not accounted for by our models, such as deviations from spherical symmetry of the gravitational potential, non-negligibly contribute to the observed flattenings.

Fig. C.1. Observed and model axial ratio, ζ, of the X-ray surface-brightness distributions of the clusters of our sample. Upper panel: Observed distribution of ζ (black) compared to the distributions of ζ of the αturb = 0 (red) and αturb ≥ 0 (blue) models. Lower panel: Observed (squares) ζ compared with the median values of ζ obtained for each cluster (as indicated in the vertical axis) for the αturb = 0 (circles) and αturb ≥ 0 (triangles) models. The error bars indicate the 16th and 84th percentiles measured in our models. The observed data points have very small 1σ uncertainties (δζ/ζ ∼ 10−3), which are not shown in the plot.

Appendix D: Geometrical approximations in the model-data comparison

Let us consider a generic quantity Q(R, z) that depends on the intrinsic properties of the ICM in an axisymmetric cluster model represented in cylindrical coordinates (R, ϕ, z). In addition to the cylindrical coordinate system, we adopt also a Cartesian system with the same z axis, and with x and y related to R and ϕ by $R=\sqrt{x^2+y^2}$ and ϕ = arctan(x/y). Taking x parallel to the line of sight, we also define a system of polar coordinates in the plane of the sky $(\widehat{R},\widehat{\varphi })$, related to y and z by $\widehat{R}=\sqrt{y^2+z^2}$ and $\widehat{\varphi }=arctan(z/y)$. For comparison with the observations we need to evaluate integrals of Q over a cylindrical shell with symmetry axis along the line of sight x, and with inner and outer radii ${\widehat{R}}_{\text{i}}$ and ${\widehat{R}}_{\text{i + 1}}$, respectively:

$$\begin{array}{c}\hfill P_{Q,i}={\displaystyle {\int }_0^{2\pi }{\int }_{{\widehat{R}}_{\mathrm{i}}}^{{\widehat{R}}_{\mathrm{i}+1}}}f(\widehat{R},\widehat{\varphi })\widehat{R}\mathrm{d}\widehat{R}\mathrm{d}\widehat{\varphi },\end{array}$$(D.1)

where $f(\widehat{R},\widehat{\varphi })\equiv 2{\displaystyle {\int }_0^{\infty }}Q(R,z)\text{d}x$. We also need to evaluate integrals over spherical shells with inner and outer radii r_j and r_{j + 1}, respectively:

$$\begin{array}{c}\hfill I_{Q,j}=\frac{{\displaystyle {\int }_{r_{\mathrm{j}}}^{r_{\mathrm{j}+1}}}Q(R,z)r^2\mathrm{d}r}{{\displaystyle {\int }_{r_{\mathrm{j}}}^{r_{\mathrm{j}+1}}}{r}^2\mathrm{d}r}.\end{array}$$(D.2)

Given that our models are neither spherically symmetric nor circularly symmetric in the plane of the sky, the numerical calculation of Eqs. (D.1) and (D.2) is computationally expensive and would represent the bottleneck of our MCMC algorithm. We thus chose to replace Eqs. (D.1) and (D.2) with less time-consuming approximations.

We limited the evaluation of Q to z = 0 plane. Assuming $f(\widehat{R},\widehat{\varphi })\approx f(\widehat{R},0)$ and Q(R, z)≈Q(R, 0), Eqs. (D.1) and (D.2) become

$$\begin{array}{c}\hfill P_{Q,i}\approx 2\pi {\displaystyle {\int }_{{\widehat{R}}_{\mathrm{i}}}^{{\widehat{R}}_{\mathrm{i}+1}}}f(|y|,0)y\mathrm{d}y\end{array}$$(D.3)

and

$$\begin{array}{c}\hfill I_{Q,j}\approx \frac{3{\displaystyle {\int }_{r_{\mathrm{j}}}^{r_{\mathrm{j}+1}}}Q(|y|,0)y^2\mathrm{d}y}{r_{\mathrm{j}+1}^3-r_{\mathrm{j}}^3},\end{array}$$(D.4)

respectively.

All Tables

All Figures

Fig. 1. Equatorial plane velocity profiles of the ICM in the αturb = 0 models of the clusters of our sample. Left panel: Ratio between the gas rotation velocity and velocity dispersion (measured in the equatorial plane) as a function of the cylindrical radius normalized to the virial radius, r200 = 3M200/[800πρcrit]1/3, computed using the median value of M200 for each cluster. Right panel: Gas rotation velocity as a function of the cylindrical radius in physical units. In both panels the dashed lines are the median profiles of the individual clusters (as indicated in the legend), whereas the gray band represents the interval between the 16th and 84th percentiles of the distribution of the profiles of the entire sample of clusters.

In the text

Fig. 2. Equatorial plane rotation and turbulence profiles of the ICM in the αturb ≠ 0 models of the clusters of our sample. Upper panels: Same as Fig. 1 but for models with αturb ≥ 0. Lower left panel: Median turbulent-to-total pressure ratio, αturb, as a function of the cylindrical radius normalized to the virial radius, r200 (see the caption of Fig. 1), for the same models as in the upper panels. Lower right panel: 1D turbulent velocity dispersion as a function of the cylindrical radius in physical units for the same models as in the upper panels. The color coding and the meaning of the gray bands are the same as in Fig. 1. In the lower right panel, we overplot the αturb profile measured in simulations by A20 (solid black line; the profile was taken from their Eq. (11), with a0, a1, and a2 equal to their best-fit values, as reported in the first row of their Table 1, and assuming r200, m/r200 ≈ 1.67).

In the text

	Fig. 3. Simulated spectrum (black crosses), with an exposure of 100 ks, in half of the Resolve FoV with the best-fit bapec model overplotted. The thermal model only is shown with the green line, the instrumental background in blue, and the thermal model plus the instrumental and cosmic background in red. The most prominent lines are the Fe XXV and Fe XXVI complexes in the range 6−7 keV.
In the text

Fig. A.1. Upper panel: Probability density function (pdf) of the systematic offset parameter, η, in our cluster sample. The gray band ranges between the 16th and 84th percentile of η in the entire sample. Bottom panel: Median η of each cluster (as indicated in the vertical axes) according to the αturb ≥ 0 model, with error bars corresponding to the 16th and 84th percentiles. For reference, η = 1 (no offset) is indicated with a red vertical line.

In the text

Fig. C.1. Observed and model axial ratio, ζ, of the X-ray surface-brightness distributions of the clusters of our sample. Upper panel: Observed distribution of ζ (black) compared to the distributions of ζ of the αturb = 0 (red) and αturb ≥ 0 (blue) models. Lower panel: Observed (squares) ζ compared with the median values of ζ obtained for each cluster (as indicated in the vertical axis) for the αturb = 0 (circles) and αturb ≥ 0 (triangles) models. The error bars indicate the 16th and 84th percentiles measured in our models. The observed data points have very small 1σ uncertainties (δζ/ζ ∼ 10−3), which are not shown in the plot.

In the text