© The Authors 2025

1. Introduction

Galaxy clusters are the largest and most massive gravitationally bound objects observed in the Universe. Their large fraction of dark matter (DM) and highly ionised plasma (the intracluster medium, or ICM; often referred to as ‘the gas’) makes them an ideal workbench for both observational cosmology (e.g., see Reiprich & Böhringer 2002; Allen et al. 2011) and X-ray astronomy (Sarazin 1988). Even though they contribute only a small percentage to the total cluster mass, galaxies remain a valuable tool for investigating the dynamical state of these systems.

Both the gas and the galactic components of most relaxed clusters are traditionally considered to be fully pressure-supported, in the sense that their kinetic energy is taken to be in the form of random motions of their constituents. The ICM is frequently assumed to be in hydrostatic equilibrium, which provides a convenient model to derive the total mass of the cluster (for a review, see Ettori et al. 2013. Similarly, galaxies in clusters are generally considered (as a collisionless component) to be in a quasi-stationary state of dynamical equilibrium with no or negligible systematic motions (see, e.g., Carlberg et al. 1997). In practice, such a picture largely derives from the lack of detailed spectroscopic studies, much like how in the case of globular clusters non-rotating King models remained the standard paradigm until deeper kinematical studies became available (see Djorgovski & Meylan 1994; Zocchi et al. 2012; Bianchini et al. 2013). It follows that a convincing detection of systematic rotation of galaxy clusters would significantly change our current picture of these systems, pushing forward the current understanding of their formation, evolution, and the role of the interaction with their surroundings. Cluster rotation could have different origins, such as tidal torques (Peebles 1969) or off-axis mergers (Ricker 1998; Roettiger & Flores 2000). In other words, the existence of non-thermal dynamics in these systems might call for a significant revision of the standard cosmological framework, as a definitive understanding of the presence of angular momentum in the large-scale structure of the Universe still remains elusive.

There is a growing effort to quantify and characterise the various non-thermal pressure contributions of the ICM, such as turbulent and bulk motions, including rotation (Bianconi et al. 2013; Gianfagna et al. 2021; Ettori & Eckert 2022; Bartalesi et al. 2024, 2025). Correspondingly, a limited number of studies have hinted at the possible existence of non-negligible systematic rotation of the galactic component (Hwang & Lee 2007; Tovmassian 2015), and mostly within the central regions of some specific clusters (Materne & Hopp 1983; Ferrami et al. 2023). An interesting example is Abell 2107, for which there is some evidence of rotation for both the ICM (Liu & Tozzi 2019) and the galaxies (Oegerle & Hill 1992; Kalinkov et al. 2005; Song et al. 2018).

Furthermore, most of the past investigations that looked for rotation of the galactic component have relied (in the best cases) on only a few tens of data points with measured kinematics and have primarily focussed on the asymmetries of the velocity distribution in the relevant position-velocity diagrams (Smith 1936; Rood et al. 1972; Dressler 1981). However, in recent years, the significant increase in the availability of high-quality kinematic data for galaxies in clusters has encouraged the use of more accurate and robust rotation-detecting methods (e.g., see Hwang & Lee 2007; Ferrami et al. 2023). This approach has already been proven to be successful in the parallel field of globular clusters, where a similar improvement in the quantity and quality of the data made it clear that many of these structures are partially rotation-supported (see, e.g., Bianchini et al. 2013; Kacharov et al. 2014; Leanza et al. 2022).

As a collisionless component, it is believed that galaxies in clusters have reached a state of partial relaxation via the mechanism of violent relaxation (see Lynden-Bell 1967). While this process is not necessarily associated with the production of rotational motion, it leaves the core regions of the cluster generally more relaxed than the outer regions. Combined with the new developments on galaxy and globular clusters mentioned above, this suggests that when present, systematic rotation may be expected to be differential, with the central region typically rotating rigidly and the outskirts possibly not rotating at all.

This paper has two main objectives. First, we look for evidence of systematic rotation in the galactic components of a sample of Abell clusters; this involves the use of a well-known algorithm that evaluates the rotational velocity and the projected position angle of the rotation axis. Second, we aim to clarify and quantify the ‘low-number bias’ that we mentioned above and assess its impact on the estimates of rotational support as a function of the number of kinematic data points considered.

This paper is organised as follows. The sample of clusters with the relevant photometric and spectroscopic data is presented in Section 2. In Section 3, we outline the methods of investigation used for both the search of systematic rotation and the evaluation of the related low-number bias. Then, in Section 4, we show the results of the analysis followed by their discussion in Section 5. The final conclusions and remarks are reviewed in Section 6.

2. Data samples

In this paper, the analysis and comparison of our results with those of previous studies is carried out on two different datasets, collected with different instruments.

2.1. Sample for rotation detection

The search of rotation signatures is mainly conducted with the sample of clusters presented in D’Eugenio et al. (2015). They discarded the clusters at redshift z >  0.1 and those with richness ℛ <  2 from the Abell et al. (1989) catalogue, resulting in a sample of 20 clusters. From D’Eugenio et al. (2015), we have extracted the following properties, shown in Table 1.

• The co-ordinates of the cluster centre, which hereafter we assume to be those of their brightest central galaxy. When the second brightest galaxy falls within 0.5 mag from the brightest one, the cluster centre is taken to be the mid-point between the two.

• The estimated value of the cluster redshift, z.

• The angular distance from the cluster centre, R (in arcmin), that corresponds to 1.5 Mpc at the cluster redshift, calculated with the Planck cosmology (H₀ = 67 km s⁻¹ Mpc⁻¹, Ω_m = 0.32, Ω_Λ = 0.68, see Planck Collaboration VI 2020). This radius is employed as a standard scale within which we look for rotation. We opted for this particular value because, for typical clusters, it manages to balance the two effects mentioned at the end of Sect. 1. A larger radius would likely dilute the rotation signal, whereas a smaller one would reduce the number of member galaxies, preventing a more statistically robust analysis. Note that these radii are calculated under the assumption that the clusters do not have any significant peculiar motion with respect to the Hubble flow.

The optical and spectroscopic data of this sample were taken from SDSS DR18 (Almeida et al. 2023). For each cluster, we selected the right ascension and declination of the galaxies within the angular radius, R. We then extracted only the most reliable spectroscopic redshifts (zWarning = 0), listed with the ‘GALAXY’ spectroscopic type. Membership selection was made by performing a Gaussian fit to the redshift distribution in the whereabouts of the z listed in Table 1. Note that fast-rotating clusters might show a strongly bi-peaked redshift distribution, but no major deviations from the usual Gaussian distribution of redshifts are expected for slow rotators. This allows us to adopt the described fitting procedure even in the presence of a relatively low rotational velocity.

We then obtained a new estimate of the cluster redshift (z_cl) as the mean of the Gaussian curve. Therefore, the fitted standard deviation (σ_z) can be converted into the global velocity dispersion of the galactic component with

$$\begin{array}{c}\hfill {\sigma }_{\mathrm{fit}}=\frac{c}{1+z_{\mathrm{cl}}}{\sigma }_z.\end{array}$$(1)

Finally, we selected member galaxies by means of a 3σ_fit-clipping of the line-of-sight velocity relative to the cluster mean, i.e., a galaxy with measured redshift z_g was considered a member if its line-of-sight velocity, v, was

$$\begin{array}{c}\hfill |v|=|z_{\mathrm{g}}-z_{\mathrm{cl}}|c<3{\sigma }_{\mathrm{fit}}.\end{array}$$(2)

For the sample for rotation detection, the fitted cluster redshift, the final total number of member galaxies, and the cluster global velocity dispersion are shown in the second, third, and fourth columns of Table 2. Note that the sample for rotation detection is reduced to 17 clusters because SDSS does not have spectroscopic data for A16, A1186, and A2256.

2.2. Bias-test sample

Additional tests of the low-number bias have been performed on a different dataset, consisting of three rich and more distant (z ∼ 0.4) clusters (Abell S1063, MACS J1206.2−0847, and Abell 370) with high-quality data, which previous studies have suggested the possibility of rotation for and identified as good candidates for this type of study. The dynamics of Abell S1063 (z ≈ 0.35) and MACS J1206.2−0847 (hereafter MACS J1206; z ≈ 0.44) have been studied extensively by Ferrami et al. (2023). The data primarily come from the CLASH-VLT program (Rosati et al. 2014). It builds on the Cluster Lensing And Supernova survey with Hubble (Postman et al. 2012), a multi-wavelength photometric survey that observed 25 massive galaxy clusters and on a spectroscopic campaign using the VIMOS instrument (Le Fèvre et al. 2003) on the Very Large Telescope. This campaign was further enhanced by observations from the MUSE (Multi-Unit Spectroscopic Explorer, Bacon et al. 2010) spectrograph, providing comprehensive data on the kinematics of galaxies within the innermost regions of the clusters (∼2 arcmin² for Abell S1063 and ∼2.6 arcmin² for MACS J1206; see, respectively, Caminha et al. 2016, 2017). After membership selection, the final dataset consists of 1215 galaxies for Abell S1063 and 658 galaxies for MACS J1206, each covering an area of 25 × 25 arcmin² (or, equivalently, of 7.7 × 7.7 Mpc² and 8.9 × 8.9 Mpc², respectively).

For Abell 370 (z ≈ 0.38), we used the data sample assembled by Lagattuta et al. (2022), who studied the structure and kinematics of this strong-lensing cluster. This catalogue consists of 416 galaxy redshifts, spanning both the central region and the outskirts, targeted with the MUSE spectrograph on multiple visits. The initial spectroscopic dataset is based on a narrower field, covering a total area of 4 arcmin² around the identified centre (Lagattuta et al. 2019). Lagattuta et al. (2022) included shallower MUSE observations covering more external regions, resulting in a final central coverage of an area of 14 arcmin² (or 1.46 Mpc²) around the core. The photometric data of Abell 370 is composed of high-resolution images mainly from the Hubble Frontier Fields programme (HFF, Lotz et al. 2017) for the cluster core, integrated with less deep measurements on the outskirts from the Beyond Ultra-deep Frontier Fields And Legacy Observations (BUFFALO, Steinhardt et al. 2020) project, with coverage of nearly 6 × 6 arcmin² (1.9 × 1.9 Mpc²). The basic properties of this second cluster sample are recorded in the bottom part of Table 1.

3. Methods

3.1. Rotation measurements

The method we used to estimate the position angle of the rotation axis and the rotational velocity has been employed in several previous studies, primarily on the available stellar kinematical data of globular clusters (see, e.g., Cote et al. 1995; Bianchini et al. 2013; Kacharov et al. 2014; Leanza et al. 2022, 2023). It operates as follows. The initial step is to subtract the mean velocity of the cluster from the line-of-sight velocity of the member galaxies. Then, the projected image of the cluster is divided into two halves by a line passing through the centre at a given position angle, θ_i. For each of the two resulting sub-samples, the average line-of-sight velocity, ⟨v⟩_i, is calculated and the difference between them is then plotted against θ_i. Subsequently, the position angle is increased by 10° (so that θ_i + 1 = θ_i + 10°), and the same procedure is repeated until the entire 360° are covered. Finally, the resulting pattern is fitted with a sine function so that the angle at which the difference between the mean velocities reaches its maximum corresponds to the position angle of the rotation axis, while the amplitude of the sine provides an estimate of a representative value of the systematic rotation velocity. An example of this procedure is shown in Fig. 1. (In passing we note that, in the left frame of the figure, the contours are not circular, because of the co-ordinates adopted in the map).

Fig. 1. Example of the fitting procedure used to search for systematic rotation of the galactic component in clusters of galaxies. Left: Projected spatial distribution of galaxies of Abell 2199, colour-coded according to their line-of-sight velocity relative to the cluster mean (v). The dashed contours represent three different cuts within which we performed the test, in order to monitor whether the velocity amplitude would increase or decrease and how the axis position angle changes with radius. Right: Corresponding sine fits, from which we extracted the parameters related to the largest scale (R) reported in Table 2.

This method has been applied to the galactic component of Abell S1063 and MACS J1206 by Ferrami et al. (2023). Similar methods are described and used in the literature; in particular, they have been used to study the possible presence of systematic rotation in clusters of galaxies (see, e.g., Dressler 1981; Hwang & Lee 2007). We decided to test the methods used in Hwang & Lee (2007) and Ferrami et al. (2023) on a set of simplified simulated galaxy clusters, for which the rotational velocity, v_rot, and velocity dispersion, σ, were known in advance, to determine which method is statistically more accurate in retrieving the correct rotation parameters as a function of v_rot/σ. The specifics of the simulation and the results are presented in Appendix A. It turns out that the method presented above is consistently better at recovering the projected orientation of the rotation axis, and it is probably preferable when considering clusters with a relatively low rotational velocity (v_rot/σ ≲ 0.30).

3.2. Low-number bias

In carrying out our investigations we found that, somewhat unexpectedly, systems with smaller numbers of kinematical data points often appeared to exhibit a stronger signal of rotation, and thus we became aware that we might have to face a bias problem. Note that, as is outlined in Section 1, there are two main reasons why some systems may have limited kinematic measurements: either due to insufficient spectroscopic observations, or because the analysis is restricted to the cluster core. Here, new observations might not be feasible simply because of a lack of suitable targets in the region considered. To investigate the extent to which the number of member galaxies influences the fitted rotation parameters, we repeatedly applied the rotation detection algorithm described above to randomly selected subsets of N member galaxies of the sample of 17 clusters listed in the upper part of Table 1. Starting from the total sample of N_tot galaxies, we progressively reduced the sample size by randomly removing ten galaxies at a time, and extracted v_rot 5000 times for each value of N. Under the assumption that the value fitted when N = N_tot is consistent with the one calculated with infinite statistics (that is to say, the true rotation velocity), we studied the trend of the average rotation velocity to velocity dispersion ratio as a function of N.

A simple model has been found to provide an approximate estimate of the bias introduced by the low-number statistics; for a derivation of the following Eq. (3) and for the definition of the correlation coefficient, the reader is referred to Appendix B. If a cluster is rotating, the difference between the mean velocities calculated on each side of the rotation axis projected on the sky will be non-zero. Let us call Δv half this mean velocity difference, as measured in a sample of N galaxies. For N → ∞ this value approaches v_rot. Similarly, we can define σ_v as the velocity dispersion calculated for a sample of N galaxies, and σ its limit for large N. Then, the ratio ⟨|Δv|⟩/σ_v as a function of N is

$$\begin{array}{c}\hfill \frac{⟨|\mathrm{\Delta }v|⟩}{{\sigma }_v}=[\frac{v_{\mathrm{rot}}}{\sigma }\left(\frac{v_{\mathrm{rot}}}{\sigma }\sqrt{\frac{N}{\widehat{\rho }}}\right)+\sqrt{\frac{\widehat{\rho }}{\pi N}}exp(-\frac{v_{\mathrm{rot}}^2}{{\sigma }^2}\frac{N}{\widehat{\rho }})]g(N),\end{array}$$(3)

where $\widehat{\rho }=1/(1+\rho )$ depends on the correlation coefficient, ρ, between the samples on each side of the rotation axis, and

$$\begin{array}{c}\hfill g(N)=\sqrt{\frac{N-1}{2}}\frac{\mathrm{\Gamma }((N-1)/2)}{\mathrm{\Gamma }(N/2)}.\end{array}$$(4)

Equation (3) tells us that ⟨|Δv|⟩/σ_v as a function of N depends on two parameters: the limit ratio, v_rot/σ, and the correlation coefficient, ρ. As was previously mentioned, the first parameter is specific to each cluster and is approximated using the value obtained when all the galaxies are considered. In the case of the second parameter, under the assumption that the correlation coefficient is the same for each cluster, we calibrate its value in the following section using the eight clusters for which significant evidence of systematic rotation has been found.

4. Results

4.1. Rotation detection and model calibration

As was noted at the beginning of Sect. 3.1, Fig. 1 illustrates the method that we applied to derive the rotation parameters listed in Table 2. The same figure also clarifies what happens in the application of the method when, for a given cluster with a given set of data, different cuts are performed in order to consider smaller and smaller regions. In the case shown in Fig. 1, the cuts have been made at 1R, 0.6R, and 0.2R. In the absence of the low-number bias, this procedure might provide indications as to the radial dependence of rotation. The low-number bias tends to confuse the picture, because the different cuts define datasets with different values of N (see comments at the beginning of Sect. 3.2). In this respect, one should note that the present method leads to an estimate of the typical rotation speed, as alinear velocity; for rigid rotation, which is constant angular velocity, the typical rotation speed decreases linearly with decreasing radius. This combination of factors prevents an unambiguous interpretation of the trend shown by the three curves on the right panel of Fig. 1.

In 10 of the 17 clusters in the sample for rotation detection, the fitted parameter values indicate a significant rotation support over random motions (v_rot/σ ≳ 0.15, see Table 2). However, a high value of this parameter does not necessarily guarantee that the cluster is actually rotating. In fact, a similar signal may be caused by various factors, including the following:

1. The cluster could consist of two (or more) separate clusters at slightly different angular positions in the sky. This is the case for systems exhibiting a distinct sub-clustering of galaxies, often associated with a strongly two-horned distribution of redshifts.

2. The cluster may be the result of a recent (or ongoing) merging event and may not have reached a quasi-equilibrium state yet. Whether or not off-axis mergers could mimic the presence of systematic rotation, the search for angular momentum in clusters of galaxies is best carried out in systems that appear to have already settled in a quasi-equilibrium state. In this respect, X-ray data of the ICM offers a particularly useful probe of the dynamical state of clusters, with unrelaxed systems displaying evidence of disturbances such as shock fronts, filaments, or irregular isophotes.

3. The dataset could contain an insufficient number of data points with measured kinematics, and therefore suffer from the low-number bias mentioned in Sect. 3.2.

Among our sub-sample of ten clusters for which rotation is indicated, for the following analysis we discarded the ones that previous investigations had shown to consist of close but distinct substructures: this is the case for A1035 (see Kopylov & Kopylova 2007) and A2151 (see Monteiro-Oliveira et al. 2022). This left us with eight clusters with a fairly relaxed optical morphology and whose redshift distribution exhibits indications of systematic rotation (A1367, A1650, A2029, A2065, A2142, A2199, A2255, and A2670). Interestingly, Bartalesi et al. (2025) found that the current X-ray data leave room for significant systematic rotation of the ICM in three clusters of their dataset, the most promising candidate being A2255. Furthermore, three of the eight clusters selected above are considered to be typical, cool-core, relaxed clusters with regular X-ray morphology and no evidence of recent merging events: A1650 (Govoni et al. 2009), A2029 (see, e.g, Buote & Tsai 1996), and A2199 (see, e.g., Markevitch et al. 1999). For the second objective of this paper – that is, to quantify the low-number bias – we consider the eight clusters with a sufficiently regular optical morphology and significant evidence of systematic rotation.

Figure 2 shows the resulting trends of the average ⟨|Δv|⟩/σ_v when gradually removing galaxies from the eight selected clusters (5000 resamplings for each N). The rotation parameter, normalised to its limit value when using all the galaxies, shows a steep increase as the number of kinematic data points is reduced. We also plot the related analytical curves (Eq. (3)) and the corresponding residuals, Δ. The optimal value of ρ was determined by minimising the sum of the clusters’ mean square errors. For this optimisation process, the case with N = 10 was excluded, as the hypotheses at the basis of the simple model described in Appendix B may not be applicable in this context (for example, the Gaussianity of the distributions and the equal number of the galaxies on the two sides of the axis, see Appendix B). The final best-fit value of the correlation coefficient turned out to be ρ_opt = −0.788. A more general approach is taken in Appendix C, in which both the limit ratio, v_rot/σ, and the parameter, ρ, are fitted simultaneously. This different procedure confirms that the latter coefficient remains approximately the same for each cluster. This justifies our original assumption that all clusters exhibit a similar correlation coefficient when considering the two sub-samples of galaxies.

Fig. 2. Upper panel: Low-number bias for the eight clusters selected in this paper compared to the related analytical curve with the best-fit ρ = −0.788. To illustrate the statistical uncertainty without overcrowding the graph, we only show the A2199 one-sigma error band forreference. Lower panel: Residuals, Δ, between the data points and the analytical curves. In the legend, we use Nt = Ntot, vr = vrot, and σ calculated as the standard deviation of the full set of galaxies (this is what ultimately causes the difference between the values of vr/σ reported here and the vrot/σfit in Table 2, where σfit is extracted from the Gaussian fit of the redshift distribution). The case of N = 10 was excluded from the analysis, as undesired additional systematic effects are expected to appear with such a small number of galaxies.

As is suggested by Fig. 2, there appears to be an inverse relation between the magnitude of the rotation signal and the level of the bias, in the sense that the bias is low when the rotation signal is high. This is clarified further by Fig. 3, where the analytical curves as a function of N are plotted for different fixed true ratios, v_rot/σ, in the limit N → ∞. As the value of v_rot/σ increases, the influence of N on the measured ratio diminishes. This is due to the fact that the folded Gaussian distribution of velocities undergoes only minimal alterations at its tails (see Appendix B). Furthermore, the values of the eight clusters employed in this study, along with those identified by Hwang & Lee (2007), are superimposed on the curves. These last clusters exhibit minimal bias, as is evidenced by their high rotational velocities.

Fig. 3. Different impact of the low-number bias between fast- and slow-rotating clusters. Because of their stricter selection criteria, Hwang & Lee (2007) found clusters whose rotational support is so high that they can be almost assumed to be bias-free, even though most of them have only a few tens of data points with measured kinematics.

4.2. Bias testing with radial cuts

Finally, we tried to determine whether the bias affects the measurements in a similar way when the galaxies were removed from the projected distribution with progressively smaller radial cuts, rather than being randomly selected across the entire extent of the cluster. We conducted this test by using the three massive clusters with a large number of identified cluster members (Abell S1063, MACS J1206, and Abell 370) of the bias-test sample, introduced in Sect.2.2. As was noted earlier in Subsection 4.1, in relation to the interpretation of Fig. 1, an increase in |Δv|/σ_v when focussing on the central regions could be explained by a real rotation signal rather than the low-number bias, even though for rigid rotation the representative value of rotation decreases when the radius of the disc is decreased.

Another important issue that needs to be addressed is the fact that removing galaxies with consecutive radial cuts produces only one realisation of the cluster for each N, whereas in Fig. 2 each point is the average of 5000 resamplings. Therefore, the specific distribution of galaxies in each cluster introduces a large fluctuation in the |Δv|/σ_v trend. In order to mitigate this effect, we adopted the following procedure.

1. To be consistent with the previous tests, we focussed our attention on 1.5 Mpc as the outermost scale (R, see Table 1), containing N_tot galaxies.

2. For every fixed value of N increasing in steps of 10 up to N_tot, we calculated the cut-off radius containing the N galaxies.

3. For each N, corresponding to a radial cut, we made 500 different resamplings by randomly removing 40% of the galaxies. After trying several removal percentages, we settled, somewhat arbitrarily, on the value of 40%, as this represents a good compromise between removing too few galaxies, which would not mitigate the |Δv|/σ_v fluctuation, and removing too many galaxies, which would significantly alter the overall appearance of the cluster. We then calculated for each resampling its |Δv|/σ_v and took their average. Note that the resulting ⟨|Δv|⟩/σ_v may be an overestimate, as we expect the low-number bias to come into play when using a sub-sample of galaxies.

4. Finally, given this undesired effect introduced in the previous step, we corrected for the low-number bias using the analytical model with ρ = −0.788 calibrated in Sect. 4.1.

Figure 4 shows the final ⟨|Δv|⟩/σ_v trend for Abell S1063, MACS J1206, and Abell 370, after applying the procedure described above. In the first two cases, the measured ratios of rotation velocity over dispersion appear to roughly align with the trend predicted by our model of the bias. This supports the idea that, in these two clusters, if there is any systematic rotation of the member galaxies, it might be a negligible component of their dynamics. At variance with these two clusters, Abell 370 shows a significant deviation from the predicted trend of the bias in the central regions; this apparently provides compelling evidence that the cluster core has a substantial amount of rotationalsupport.

Fig. 4. Rotation velocity over velocity dispersion trends when galaxies are removed by successive radial cuts in the three clusters with a large number of identified members with measured line-of-sight velocities: Abell S1063, MACS J1206, and Abell 370 (in the legend, we use Nt = Ntot and vr = vrot). As in Fig. 2, we display only the AS1063 one-sigma error band to illustrate resampling uncertainty without overcrowding the figure. While for the first two the low-number bias model reproduces the observed trend reasonably well, A370 deviates significantly in the central region – convincing evidence of systematic rotation in its core. The vertical lines indicate the number of galaxies within a 0.5 Mpc radius.

5. Discussion

The extent to which the low-number bias can affect the correct recovery of the rotation parameters varies from case to case, as it depends mainly on the rotational velocity and the number of galaxies used to extract it. Using Eq. (3) and the definition of relative error

$$\begin{array}{c}\hfill \epsilon =\frac{⟨|\mathrm{\Delta }v|⟩/{\sigma }_v}{v_{\mathrm{rot}}/\sigma }-1,\end{array}$$(5)

we can study this dependence by finding the minimum number of galaxies, N, required to get an error on the order of ε in the case of a true ratio, v_rot/σ. The results are shown in Fig. 5. We find that galaxy clusters with v_rot/σ ≳ 0.25 are essentially bias-free when using around 100 galaxies or more; for the 0.15 ≲ v_rot/σ ≲ 0.20 range, it would be preferable to use more than 120 data points to limit the error to ε ≲ 0.1. Finally, it appears that slow-to-non-rotators (v_rot/σ ≲ 0.15) require a significantly larger number of galaxies (up to more than 300) to validate the rotation signal.

Fig. 5. Minimum number of galaxies needed to achieve a relative error, ε, in the case of a true ratio of vrot/σ.

The low-number bias might be studied and should also play a role in the context of globular clusters, in those cases in which the search for rotation is carried out with similar procedures. A real case example of how this bias can impact rotation measurements is studied in Appendix D. However, for the study of star motions in globular clusters other sources of information are often available (in particular, proper motions for Galactic globular clusters).

A final point worth discussing is how the position angle (PA) varies when we remove galaxies, either by bootstrapping or by radial cuts. As can be expected, if there is indeed a systematic redshift distribution of the cluster member galaxies, the PA should remain fairly centred on a single value, even if we do not use the entire dataset. A narrower distribution of PAs is therefore a sign of a more stable and thus more reliable rotationmeasure. We have investigated this fact by keeping track of the difference between the 84^th and 16^th percentiles (P₈₄ − P₁₆) of the distribution of PAs extracted from the random subsets of cluster galaxies used for the characterisation of the low-number bias (see Section 3.2, 5000 resamplings for each N). Representative results for rotating and non-rotating clusters are shown in Fig. 6. This shows that in both cases the confidence interval containing 68% of the PAs becomes larger as fewer galaxies are used. However, the non-rotating cluster A168 shows a much steeper increase, suggesting the absence of a preferred direction of the PA, and hence the lack of a clear systematic distribution ofredshifts.

Fig. 6. Difference between the 84th and 16th percentile of PAs for each N of a non-rotating (Abell 168) cluster and a rotating (Abell 2199) one (in the legend, we use Nt = Ntot, vr = vrot, and σ calculated as the standard deviation of the full set of galaxies).

We can try to use this result to strengthen the evidence of rotation in the core of Abell 370 that we argued for when looking at the deviation from the model of the bias in Fig. 4. Using the same bootstrapping technique we introduced in Sect. 4.2, we plot in Fig. 7 the median values of the fitted PAs in the four innermost radii (corresponding to N = 20, 30, 40, 50, where the deviation in Fig. 4 is most noticeable). To assess the stability of these PA measurements, we also show the interval between the 16^th and 84^th percentiles for each N, providing a measure of their statistical uncertainty. Interestingly, A370 PAs are spread over a much smaller region than the other two clusters, with a maximum angular separation between the medians of 16°, while for Abell S1063 the separation is 52° and for MACS J1206 it is 127°. We can further characterise this fact by considering the minimum and maximum P₈₄ − P₁₆ interval for the three clusters, which are, respectively, 57° (N = 20) and 98° (N = 40) for Abell S1063, 81° (N = 50) and 136° (N = 30) for MACS J1206, and 24° (N = 20) and 36° (N = 50) for A370. Its narrow interval of PAs apparently confirms that the innermost region of the galactic component of Abell 370 indeed has a non-negligible amount of systematic rotation.

Fig. 7. Medians of the fitted PAs for the four smallest values of N (the radial co-ordinate) for the three clusters in the bias-test sample. The coloured arcs represent the intervals between the 16th and 84th percentiles of PAs.

6. Conclusions

The first objective of this study has been to investigate the presence of systematic rotation of the galactic component in a sample of 17 rich and nearby (z <  0.1) Abell clusters. To achieve this, we applied a standard algorithm to extract the rotation parameters (v_rot, PA) by fitting the difference in the average line-of-sight velocity of galaxies, Δv(θ), between the two regions defined by an axis at position angle θ, as a function of θ itself. After rejecting the clusters that exhibit clear evidence of subclustering in the optical data, we are left with eight clusters that show a sufficiently high rotational support over the velocity dispersion (v_rot/σ ≳  0.15). Three of them, A1650 (v_rot ≈ 176 km/s, PA ≈109°, v_rot/σ = 0.23), A2029 (v_rot ≈ 229 km/s, PA ≈58°, v_rot/σ = 0.18), and A2199 (v_rot ≈ 154 km/s, PA ≈78°, v_rot/σ = 0.21), do not show any evidence of a recent merging event in X-ray data and previous literature considers them to be typical, cool-core, relaxed galaxy clusters.

We also investigated whether the number of kinematic data points has an effect on the determination of the rotation parameters. Naively, one might assume that when using a reduced sample of member galaxies, the fitted rotational amplitude and axis position angle would still provide accurate, though less precise, estimates. Our findings indicate that this is not always the case. Using the full set of eight tentatively rotating clusters with a sufficiently regular optical morphology, we find a clear tendency to overestimate the amplitude of the rotation signal when the dataset is reduced by multiple random resampling at increasingly small values of N (see Fig. 2). The effect also depends on the rotational velocity, since the faster the rotation the smaller the number of galaxies required to minimise the bias. To quantify this relation, we introduced a simple analytical model (see Eq. (3) and Appendix B) and used the eight clusters to find the best-fit value for the free correlation parameter, ρ, associated with the model. The optimal value for this parameter was extracted by least square minimisation and was found to be ρ_opt ≈ −0.79.

The final model was used to determine the minimum value of N required to limit the fractional increment, ε, of the measured ratio over its limit value, v_rot/σ (see Eq. (5)). Our results indicate that a reliable (ε ≲ 0.15) detection of non-negligible rotation (v_rot/σ ≳  0.15) should utilise at least 120 data points with measured kinematics..

Finally, we tested whether our model could account for the increase in the v_rot/σ parameter when focussing on the central regions of the clusters. To do this, we used high-quality MUSE data from three rich clusters and ran the rotation detection algorithm described above on subsets of galaxies obtained from concentric radial cuts. While for two of them (Abell S1063 and MACS J1206) the data points remain fairly centred around the analytical model curve of the bias, the innermost values of Abell 370 deviate significantly from what the model predicts (Fig. 4). Thus, we argue that this deviation confirms the presence of a significant rotation of the core of Abell 370, supported by the fact that the rotation axis PAs for these regions in this cluster are contained within an interval much narrower than those for the other two clusters (see Fig. 7).

Acknowledgments

The authors thank the anonymous referee for useful comments that improved the presentation of the work. We acknowledge financial support through grants MIUR2017 WSCC32 and MIUR2020 SKSTHZ. We thank P. Rosati and A. Mercurio for the spectroscopic datasets of Abell S1063 and MACS J1206 and for useful discussions.

References

Appendix A: Performance of the two rotation detecting algorithms

In order to evaluate the efficiency of the rotation-detecting algorithms, it would be useful to have a number of rotating cluster-projected images of which the kinematics is known. Since this is not available for real clusters, we must refer to simulated images of clusters of galaxies. We adopt a statistical approach to describe the galactic component by means of a continuous, non-negative, phase-space distribution function f(⃗x, ⃗v). Here, we do not aim at constructing a physically justified model for the galactic component, but only at providing a reasonable representation of the distribution of galaxies in phase space when some systematic rotation is present.

Clusters are bound and quasi-stationary structures and, despite its low collisionality (with respect to galaxy-galaxy scattering; see, e.g., Bahcall 1977), the galactic component may be considered to have reached an equilibrium state, at least in the central regions. This picture provides arguments in favour of using a distribution function that closely resembles the Maxwell-Boltzmann f(E)=Aexp[−aE], where $E=\frac{v^2}{2}+\mathrm{\Phi }$ is the single-galaxy energy per unit mass and $a=\frac{1}{{\sigma }^2}$, with σ² the velocity dispersion, assumed here to be constant. The simplest way to introduce rotation is to replace the energy per unit mass with the Jacobi integral

$$\begin{array}{c}\hfill H=E-J·\omega ,\end{array}$$(A.1)

where ⃗ω and ⃗J have, respectively, the dimensions of an angular velocity and an angular momentum per unit mass. This quantity represents the single-galaxy energy per unit mass in the frame of reference rotating with a constant angular velocity vector ⃗ω (see Landau & Lifshitz 1980), and, following the Jeans theorem (Jeans 1915), it is a suitable integral of motion in the case of axisymmetric potentials. This leaves us with the following distribution function:

$$\begin{array}{c}\hfill f=Aexp[-\frac{1}{{\sigma }^2}(\frac{v^2}{2}-J·\omega +\mathrm{\Phi })].\end{array}$$(A.2)

Although rigid rotation is better motivated (at least for the core region of the cluster), to test the accuracy of the rotation-detecting algorithms we decided to use a rotation law of the form ω(R)=v_rot/R, where R is the distance from the rotation axis, so that the rotational velocity would remain constant. The apparent violation of the Jeans theorem is not an issue in this context, because we are not interested in the evolution of the galactic component and we only wish to check the quality of the algorithms used, in the presence of systematic motions in the galactic component. Finally, to complete the description of the spatial distribution of galaxies, we model the total gravitational potential with a NFW profile (Navarro et al. 1997):

$$\begin{array}{c}\hfill \mathrm{\Phi }=-{\mathrm{\Phi }}_0\frac{r_s}{r}log(1+\frac{r}{r_s}).\end{array}$$(A.3)

To produce simulated images of a rotating cluster, we have developed a Monte Carlo simulation that utilises the Metropolis algorithm (Metropolis et al. 1953) to sample 250 phase-space points (the three dimensional position and velocity vectors of the galaxies) of 500 clusters, from the simple distribution function in Eq.(A.2). In this way, even when we narrow our field of view to match the scale radius R, the number of data points drops to a typical value of the order of 120, which is still large enough to be considered unbiased, regardless of the value of v_rot/σ (see Fig.5).

We tested the two most used rotation-detecting methods. The first one is that described in Sect.3: here, for simplicity, we will refer to it as AVD (Average Velocity Difference). The second one was introduced by Hwang & Lee (2007), here referred to as HL: it employs a fitting procedure for the i-th observed galaxy with radial velocity v_i and projected angular position θ_i, using the function

$$\begin{array}{c}\hfill v_{p,i}(v_{\mathrm{rot}},{\theta }_{\mathrm{ax}})=v_{\mathrm{rot}}sin({\theta }_i-{\theta }_{\mathrm{ax}}).\end{array}$$(A.4)

The optimal fitting parameters v_rot and θ_ax are extracted through the least-square minimisation of

$$\begin{array}{c}\hfill {\chi }^2={\displaystyle \sum _i}{(v_{p,i}-v_i)}^2.\end{array}$$(A.5)

The results are shown in Fig. A.1. It is important to note that the performance of the two methods depends on the adopted distribution function: in particular, the results may vary for a different angular velocity law ω(R), a different gravitational potential and even for a different number of galaxies. For this reason, we argue that the main conclusions of this test are that the AVD method is likely to be more appropriate for slow rotating clusters; for higher rotational velocities, the HL method is better. In addition, the AVD method appears to be more accurate in recovering the correct axis orientation.

Fig. A.1. Root mean square errors of vrot/σ and projected axis position angle extracted with the two methods described in the main text.

Appendix B: Derivation of Equation (3)

Let us call the sample of galaxies on one side of a projected axis passing through the cluster centre A and the sample on the other side B. The starting assumptions are that the distributions of velocities of A and B are Gaussian with different mean but with the same variance, which is the same as the variance of the total sample. This can be written as

$$\begin{array}{c}\hfill v_{A/B}\sim \mathcal{N}({\mu }_{A/B},{\sigma }_{A/B}^2).\end{array}$$(B.1)

If the axis that divides the two samples is the rotation axis, then we expect μ_A = −μ_B and σ_A = σ_B. Furthermore we know that since we have already subtracted the cluster member velocities from the Hubble flow, the overall distribution of velocities is

$$\begin{array}{c}\hfill v\sim \mathcal{N}({\mu }_A,{\sigma }_A^2)+\mathcal{N}({\mu }_B,{\sigma }_B^2)=\mathcal{N}(0,2{\sigma }_A^2(1+\rho ))=\mathcal{N}(0,{\sigma }^2),\end{array}$$(B.2)

where ρ is the Pearson correlation coefficient. The Pearson correlation coefficient is defined as

$$\begin{array}{cc}\hfill \rho & =\frac{\mathrm{Cov}(v_A,v_B)}{{\sigma }_A{\sigma }_B}=\frac{\mathbb{E}\left[(v_A-{\mu }_A)(v_B-{\mu }_B)\right]}{{\sigma }_A{\sigma }_B}\hfill \\ \hfill & =\frac{{\displaystyle {\int }_0^{R_{\max }}{\int }_0^{R_{\max }}}(v_A(R_A)-\mu )(v_B(R_B)-\mu )\overline{\mathrm{\Sigma }}(R_A,R_B)\mathrm{d}{R}_A\mathrm{d}{R}_B}{{\sigma }_A{\sigma }_B},\hfill \end{array}$$(B.3)

where $\overline{\mathrm{\Sigma }}(R_A,R_B)$ is the joint probability density of drawing two galaxies at R_A and R_B within a maximum radius R_max. The correlation coefficient is a dimensionless and bounded (−1 ≤ ρ ≤ 1) measure of the covariance between two samples of galaxy line-of-sight velocities drawn from opposite sides of the rotation axis. The value of ρ incorporates the effect of slight asymmetries in the surface number density and the effect of the random motion component of the member galaxy velocity distribution via the structure of the joint probability density. We expect the velocity measurements at opposite sides of the rotation axis to be anticorrelated, therefore −1 <  ρ <  0 and we can then define $\widehat{\rho }=1/(1+\rho )$. We consider ρ to be a summary statistics of the structure and dynamics of the cluster members within a maximum projected radius, and we estimate its value from a fit of a set of clusters in Sect.4.1.

Considering now a finite sample of cluster members, the Central Limit Theorem (CLT) leads to

$$\begin{array}{c}\hfill ⟨v_{A/B}⟩\sim \mathcal{N}({\mu }_{A/B},{\sigma }_A^2\frac{2}{N})=\mathcal{N}({\mu }_{A/B},{\sigma }^2\frac{\widehat{\rho }}{N}),\end{array}$$(B.4)

assuming that the cluster contains N galaxies, half in A and half in B.

The halved difference between the mean velocities in A and B is distributed as

$$\begin{array}{c}\hfill \frac{⟨v_A⟩-⟨v_B⟩}{2}=\mathrm{\Delta }v\sim \mathcal{N}(v_{\mathrm{rot}},{\sigma }^2\frac{\widehat{\rho }}{2N}),\end{array}$$(B.5)

where v_rot = (μ_A − μ_B)/2.

Using the properties of the folded normal distribution, the mean of the absolute value of Eq.(B.5) is

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill ⟨|\mathrm{\Delta }v|⟩& =\langle \left|\mathcal{N}(v_{\mathrm{rot}},{\sigma }^2\frac{\widehat{\rho }}{2N})\right|\rangle \hfill \\ \hfill & =\sigma \sqrt{\frac{\widehat{\rho }}{\pi N}}exp(-\frac{v_{\mathrm{rot}}^2}{{\sigma }^2}\frac{N}{\widehat{\rho }})+v_{\mathrm{rot}}\left(\frac{v_{\mathrm{rot}}}{\sigma }\sqrt{\frac{N}{\widehat{\rho }}}\right).\hfill \end{array}\end{array}$$(B.6)

Noting that the standard deviation $\stackrel{\sim }{\sigma }$ of a finite sample extracted from a Gaussian distribution with variance σ² is

$$\begin{array}{c}\hfill \stackrel{\sim }{\sigma }=\sigma \sqrt{\frac{2}{N-1}}\frac{\mathrm{\Gamma }(N/2)}{\mathrm{\Gamma }((N-1)/2)},\end{array}$$(B.7)

we obtain Eq.(3).

Appendix C: Simultaneous fit of v_rot/σ and ρ

Fig. C.1. Low-number bias trends of the eight selected clusters in the paper as a function of the number N of galaxies included in the analysis. The continuous lines are our model best fit curves to the data, where both vrot/σ and ρ are free parameters. Their fitted values are reported in the legend, alongside the total number of galaxies Nt = Ntot.

A simultaneous fit of both v_rot/σ and ρ allows us to retrieve a new estimate of the limit rotation velocity over velocity dispersion ratio and a distinct correlation coefficient for each cluster, whose value we previously imposed to be a constant. The resulting fits are shown in Fig. C.1. Note that we previously constrained the curves to 1 for N = N_tot because we assumed this was the ’true’ value of v_rot/σ, obtained using the maximum available statistics. In practice, removing this constraint is equivalent to treating the N = N_tot case as an approximation of the ideal limit with infinite sample size. Consequently, the values of |Δv|/σ_v at N = N_tot (referred to as v_rot/σ in Fig. 2) are in this case also affected by the low-number bias, and thus we expect them to be slightly larger than the fitted values of v_rot/σ in Fig. C.1. However, as with any other N, the size of the bias in the measured value of |Δv|/σ_v at N = N_tot depends both on the fitted v_rot/σ ratio and N_tot. This explains why the limit ratio of a relatively fast-rotating, rich cluster, such as A2199, has barely changed, whereas that of the slower A1367 has decreased more, despite its N_tot being almost as large as that of A2199.

Although the optimal values of ρ are slightly larger than those obtained using the first method, they are remarkably consistent with each other, with an average value of −0.759 ± 0.015. This confirms that this coefficient remains approximately the same for every cluster. As already mentioned in Appendix B, its numerical value is directly linked to the structural and kinematical properties of the clusters through the joint probability density function $\overline{\mathrm{\Sigma }}(R_A,R_B)$. However, it is not easy to justify and construct a distribution function that can reproduce our estimate of ρ, so, for the moment, the origin of its numerical value remains unclear.

Appendix D: Low-number bias on a sample of globular clusters

As noted above, the low-number bias is not exclusive to galaxy clusters, but it is of interest in all cases where the search for rotation is carried out using line-of-sight velocities alone.

For instance, this is the case when studying the kinematics of globular clusters for which star proper motions are not known or not included in the analysis. Here, we consider a recent paper by Petralia et al. (2024), in which they analyse the kinematics of a sample of 21 globular clusters from the APOGEE-2 survey and select the clusters exhibiting signatures of systematic rotation. The principal method they use is essentially identical to the one used in this paper, with a slight change in the notation: N_* = N_tot and A_fit = 2|Δv|. Their findings reveal that 21 of the 23 globular clusters included in their sample show a signature of systematic rotation (we should note, however, that they adopt a looser definition of rotating cluster).

The results on the bias discussed in this paper suggest that a reassessment of the fitted rotational velocity of some of the globular clusters in Petralia et al. (2024) might be appropriate. In particular, Table D.1 highlights the globular clusters that appear most likely to be affected by the presence of the bias. Specifically, our model suggests that the signal observed in the three globular clusters with |Δv|/σ_v ≤ 0.10 (NGC 6171, NGC 6273, and NGC 6752) could be largely attributed to the low-number bias.

All Tables

All Figures

Fig. 1. Example of the fitting procedure used to search for systematic rotation of the galactic component in clusters of galaxies. Left: Projected spatial distribution of galaxies of Abell 2199, colour-coded according to their line-of-sight velocity relative to the cluster mean (v). The dashed contours represent three different cuts within which we performed the test, in order to monitor whether the velocity amplitude would increase or decrease and how the axis position angle changes with radius. Right: Corresponding sine fits, from which we extracted the parameters related to the largest scale (R) reported in Table 2.

In the text

Fig. 2. Upper panel: Low-number bias for the eight clusters selected in this paper compared to the related analytical curve with the best-fit ρ = −0.788. To illustrate the statistical uncertainty without overcrowding the graph, we only show the A2199 one-sigma error band forreference. Lower panel: Residuals, Δ, between the data points and the analytical curves. In the legend, we use Nt = Ntot, vr = vrot, and σ calculated as the standard deviation of the full set of galaxies (this is what ultimately causes the difference between the values of vr/σ reported here and the vrot/σfit in Table 2, where σfit is extracted from the Gaussian fit of the redshift distribution). The case of N = 10 was excluded from the analysis, as undesired additional systematic effects are expected to appear with such a small number of galaxies.

In the text

	Fig. 3. Different impact of the low-number bias between fast- and slow-rotating clusters. Because of their stricter selection criteria, Hwang & Lee (2007) found clusters whose rotational support is so high that they can be almost assumed to be bias-free, even though most of them have only a few tens of data points with measured kinematics.
In the text

Fig. 4. Rotation velocity over velocity dispersion trends when galaxies are removed by successive radial cuts in the three clusters with a large number of identified members with measured line-of-sight velocities: Abell S1063, MACS J1206, and Abell 370 (in the legend, we use Nt = Ntot and vr = vrot). As in Fig. 2, we display only the AS1063 one-sigma error band to illustrate resampling uncertainty without overcrowding the figure. While for the first two the low-number bias model reproduces the observed trend reasonably well, A370 deviates significantly in the central region – convincing evidence of systematic rotation in its core. The vertical lines indicate the number of galaxies within a 0.5 Mpc radius.

In the text

	Fig. 5. Minimum number of galaxies needed to achieve a relative error, ε, in the case of a true ratio of vrot/σ.
In the text

	Fig. 6. Difference between the 84th and 16th percentile of PAs for each N of a non-rotating (Abell 168) cluster and a rotating (Abell 2199) one (in the legend, we use Nt = Ntot, vr = vrot, and σ calculated as the standard deviation of the full set of galaxies).
In the text

	Fig. 7. Medians of the fitted PAs for the four smallest values of N (the radial co-ordinate) for the three clusters in the bias-test sample. The coloured arcs represent the intervals between the 16th and 84th percentiles of PAs.
In the text

	Fig. A.1. Root mean square errors of vrot/σ and projected axis position angle extracted with the two methods described in the main text.
In the text

	Fig. C.1. Low-number bias trends of the eight selected clusters in the paper as a function of the number N of galaxies included in the analysis. The continuous lines are our model best fit curves to the data, where both vrot/σ and ρ are free parameters. Their fitted values are reported in the legend, alongside the total number of galaxies Nt = Ntot.
In the text