© The Authors 2026

1 Introduction

It has been clearly established that globular clusters (GCs) are not as simple as previously thought, in their shape, their dynamics, and in their internal stellar populations. The first evidence of GCs differing from spherically symmetric objects was reported by Pease & Shapley (1917). Using photographs taken by the 60-inch reflector at Mount Wilson, the authors manually counted the number of stars in different sectors. They concluded that several clusters, among them NGC 6205, had an elongated shape and were elliptical rather than spherical. Since then, many studies have focused on quantifying the ellipticity of GCs (Frenk & Fall 1982; Geyer et al. 1983; White & Shawl 1987; Chen & Chen 2010), unravelling a flattening in a significant number of clusters. These findings have indicated the importance of understanding the origin of the deviation from spherical symmetry of GCs. Determining whether this deviation arises from internal dynamics, such as rotation or velocity anisotropy, or from external tidal effects from the Milky Way (MW), will allow how GCs form, evolve, and interact with their environment to be discovered.

Several studies have suggested rotation as the main driver of the shape of GCs. For example, the flattening of NGC 5139 (Meylan & Mayor 1985; Pancino et al. 2003), NGC 104 (Mayor et al. 1984; Meylan & Mayor 1985; Bellini et al. 2017), NGC 7089 (Pryor et al. 1986), and NGC 5904 (Lanzoni et al. 2018a) has been attributed mainly to rotation. Using image processing of Palomar and SRC Sky Survey material, White & Shawl (1987, hereafter WS87) studied the flattening of 100 Galactic GCs. They found that 32% of their sample showed an axial ratio inferior to 0.9. By comparing the orientation of the semi-major axis of the GCs with the galactic centre, they concluded that there was no evidence of tidal effects from the galaxy shaping GCs and that rotation or velocity anisotropies were more likely to contribute to the shape of the cluster. Fabricius et al. (2014) found a tight correlation between the ellipticity of a cluster and its central rotation in 11 MW GCs. Using 3D kinematics data, Dalessandro et al. (2024) found a correlation between their newly introduced α parameter, probing the relative strength of the rotation signal over the disordered motion, and their best-fit ellipticity values, supporting the conclusion that flattening is linked to rotation. Kamann et al. (2018) also reported a trend with more elliptical clusters corresponding to higher rotation. They highlighted that other mechanisms, such as velocity anisotropy or tidal effects, could blur this correlation. This last point has also been emphasised by Sollima et al. (2019) as a potential reason for the absence of correlation between ellipticity and rotation. Using ellipticities from Chen & Chen (2010, hereafter CC10) and rotations derived from Gaia DR2 on 62 MW GCs, they did not find any correlation between these two parameters nor between the isophotal minor axis and rotation axis. This conclusion is also supported by the results from Bianchini et al. (2018), who computed the V/σ parameter as the ratio between the peak of the rotation that they measured and the central velocity dispersion from Baumgardt & Hilker (2018). Bianchini et al. (2018) did not find any significant correlation with ellipticity when using the values from WS87. In contrast, Cruz Reyes & Anderson (2024, hereafter CR24) noticed an alignment between the semi-minor axis of NGC 5139, NGC 104, and NGC 6341 and their rotation axis, suggesting rotation as the driving mechanism of their shape. Rotation has also been detected in the GC NGC 1904, and it is likely behind the slight flattening observed in the direction perpendicular to the rotation axis (Leanza et al. 2022).

The authors of CC10 studied the shape and orientation of 116 Galactic GCs, and they found that GCs close to the Galactic bulge exhibiting a high flattening tend to have their semi-major axis pointing towards the Galactic bulge. They concluded that tidal effects from the external Galactic potential might impact the shape of GCs located in the Galactic bulge. The authors of CC10 compared their results with those of WS87 and found significant discrepancies. They investigated the outliers in further detail to understand the possible cause of these differences and proposed two main interpretations. The first one is linked to the data. CC10 used infrared observations from 2MASS covering a broader spatial extent than the optical data used by WS87. They highlighted the possibility that dust extinction could significantly alter the recovered shape of the ellipticity when using optical observations more sensitive than infrared. The spatial extent of the available data can also play a significant role, as the dense centre of the cluster is expected to be rounder due to more stellar interactions than the outer regions of the cluster (Binney & Tremaine 2008). The second source of discrepancy can be linked to the different methods. In WS87, the authors used surface photometry to compute the shape of GCs, which can be severely impacted by bright stars, while the authors of CC10 used a probabilistic star-counting technique, which could be affected by the spacing of the grid. The dependency of ellipticity measurements on the method used had already been noted by Fall & Frenk (1985), corroborating the point from CC10 that bright stars might dominate the surface brightness maps and therefore lead to less reliable results than the probabilistic star-counting method used in CC10 and introduced in Chen et al. (2004).

Fabricius et al. (2014) and Cruz Reyes & Anderson (2024) both used an eigenvector analysis of the spatial stellar distribution, with the ellipticity and position angle directly related to the eigenvalues and eigenvectors of the covariance matrix of the stellar positions. Such a method is straightforward and mathematically simple. However, it faces issues when dealing with outliers or for datasets with a low number of available stars that can severely impact the results (Jolliffe 2002).

This brief overview of the investigation of the flattening of GCs highlights the importance of providing a robust method capable of delivering consistent results across diverse datasets. Current approaches often depend strongly on the adopted technique, the data quality, and the number of stars available, which can lead to inconsistent ellipticity estimates. In this context, a method that remains robust even for a limited sample of stars is crucial. Such an approach will be particularly relevant for future studies focusing on the shape of multiple stellar populations (MSPs) in GCs, for instance. The origin of these MSPs remains a puzzle, and understanding their flattening patterns could provide valuable information to confirm or rule out scenarios on their formation. However, these analyses are typically restricted to a small number of stars, mostly along the red giant branch (RGB). They require a method to measure the ellipticity that is both reliable and minimally sensitive to sample size. Providing a consistent and reproducible framework for measuring GC shapes will therefore not only clarify the role of rotation and other factors in cluster morphology but also open the way to investigating possible differences in the spatial structure of MSPs.

In this article, we present a robust method that we used to compute the shape of 29 Galactic GCs. In Sect. 2 we describe the photometric data available. In Sect. 3, after testing two commonly used methods, we present our robust method to recover the ellipticity of GCs, which performs well on small datasets and is robust to the presence of outliers. We present our results and explore the potential causes for the flattening of GCs in Sect. 4. Finally, we present our conclusions in Sect. 5.

2 Photometric data

In this study, we focus our analysis on RGB stars when computing ellipticities. This choice was motivated by the need for a stellar population with high and uniform photometric completeness across the entire field of view, including in the cluster cores, where crowding can affect the completeness.

We acknowledge that this selection significantly reduces the size of our stellar sample and excludes fainter stars that could enhance the statistical robustness of our measurements. However, the robust method presented in Sect. 3.4 is corrected to perform well on small datasets (see Appendix A). Dedicated tests conducted on a representative subset of clusters (see Appendix B) indicate that ellipticity estimates using all stars might suffer from photometric incompleteness. In contrast, focusing our analysis on RGB stars yields stable, unbiased measurements.

The Stetson ground-based homogeneous catalogue (Stetson et al. 2019) is one of the most extensive photometric studies publicly available. It provides homogeneous photometry in the U, B, V, R, and I filters for 48 GCs in the MW. The Hubble Space Telescope UV Globular Cluster Survey (HUGS) catalogue (Piotto et al. 2015; Nardiello et al. 2018) targets the inner regions of each cluster, providing photometry in the F275W, F336W, F438W, F606W and F814W filters for 56 GCs in the MW. We combined both catalogues in order to produce a wide-field view of each cluster in our sample, using the positions of all RGB stars in both catalogues to compute the ellipticity profiles of the GCs. The method for combining and cleaning the photometric catalogues is described in this section.

In order to merge the HST photometry with the ground-based catalogue in such a way that allows for detailed analysis of the overall shape of each cluster, the spatial completeness of the individual photometric catalogues was calculated. In this method, outlined in depth in Sect. 2.6 of Leitinger et al. (2023), the square HST field of view was reshaped in order to smoothly transition into the ground-based photometry. To do this, we first took the original, full catalogue of HST photometry and defined annuli, spaced by 1^′′, from the centre of the cluster to the outer regions. Within each annulus, 360 artificial points were equally distributed. If an artificial point was within proximity to a real star in the photometry, that artificial point was considered to be within a spatially complete region. The first annulus to reach 50% spatial completeness defined the outer edge of the HST field, effectively changing the square field of view to more of a circular shape. When the ground-based photometry was then added to the HST photometry, we removed any stars in the ground-based photometry which overlapped with the HST field, ensuring a smooth transition between the HST and ground-based fields of view.

Cleaning the two catalogues included removing non-member stars that did not belong to the cluster. More specifically, non-member stars were removed with the use of proper motions provided by Gaia DR3 (Gaia Collaboration 2016, 2021), following a method closely aligning with the one outlined in Sect. 2.4 of Leitinger et al. (2023). The HST catalogues include a ‘membership probability’ column which we used to clean the data of non-members, using a threshold of >75%. The majority of stars in the final sample of each cluster fall in the 95−100% range with <30 stars per cluster showing <95% membership probability. We cross-matched the combined HST and Stetson catalogue with Gaia DR3 using TOPCAT¹. Differential reddening correction was then completed for both the HST and ground-based photometry, following the method described in Sect. 2.2 of Leitinger et al. (2023), using reddening maps from Pancino et al. (2024). To ensure the correctness of the match, we derived approximate V and I magnitudes from Gaia G magnitudes to investigate whether a roughly one-to-one correlation existed for the matched stars. Then, we used a χ² test with the right ascension (RA) and declination (DEC) proper motions components taken from Vasiliev & Baumgardt (2021) and selected stars with χ²<5, which provided a reasonable balance between rejecting non-members and retaining likely cluster members. Stars with unreliable photometry in both catalogues were removed based on the sharp quality indicator, while the ground-based photometry from Stetson et al. (2019) was additionally cleaned using the χ quality indicator, as outlined in Sect. 2.3 of Leitinger et al. (2023). Red giant branch stars were isolated through a combination of polynomial fitting and N−σ clipping of the CMD in various photometric combinations outlined in Sect. 2.5 of Leitinger et al. (2023). During this process, we removed horizontal branch (HB) and asymptotic giant branch (AGB) stars, as well as stars which were located too far from the RGB (i.e. blue stragglers or excessively red outliers). In terms of the photometric completeness of the two individual photometric samples, we assume 100% completeness for HST stars at magnitudes brighter than the sub-giant branch as specified by Anderson et al. (2008), while the Stetson catalogue is complete between magnitudes of 12<V<19 mag across all radii, according to Stetson et al. (2019). The final sample was then comprised of member RGB stars with >50% spatial completeness and reliable photometry.

The final catalogue spans a broad spatial range: extending from 2 half-light radii in the least-covered clusters, to 23.6 half-light radii in the most well-covered ones, with a median coverage of 7.1 half-light radii. This coverage was estimated using the distance of the most distant star from the cluster centre in each case, and the half-light radius from Harris (1996). The original catalogue from Leitinger et al. (2023,2025) provides data for 30 GCs. Here, we have excluded NGC 6752 due to significant gaps in the ground-based observations from Stetson et al. (2019), which would greatly affect the ellipticity results.

3 Methods

As mentioned in Sect. 1, different methods to measure the shape of GCs can yield different results. In this section, we investigate two methods to compute the ellipticity of GCs. In the first method, described in Sect. 3.1, we smooth the spatial distribution of stars in the GCs using a two-dimensional kernel density, and fit ellipses to different contour levels; this approach makes it possible to extract an ellipticity profile as a function of the distance to the cluster centre, and is similar to the one employed by Stetson et al. (2019). In the second method, presented in Sect. 3.2, we investigate the data-driven principal component analysis (PCA) method that identifies the major and minor axes of the clusters based on the variance of the data. Each method has its advantages, assumptions, and limitations, which we further discuss in Sect. 3.3, before presenting our robust method in Sect. 3.4.

For all the methods, our inputs are stellar positions with respect to the centre of the cluster. To avoid introducing any projection effect, we converted the on-sky coordinates in Cartesian coordinates, using the formula specified in Eq. (1) from van de Ven et al. (2006).

3.1 Method 1: Fitting of density contours

We first used kernel density estimation (KDE; Scott 2015) to estimate the density of stars in the Cartesian space. This non-parametric method creates a smooth and continuous function representing the density of stars. Then, we defined a set of contour levels ranging from 10% to 90% of the maximum density value, and we fitted the coefficients A, B, C, D, G, and H representing an ellipse described by the general conic form $$F(x,y)=A{x}^2+Bxy+C{y}^2+Dx+Gy+H=0$$(1) to the stellar positions using a least-squares method (Halir & Flusser 1998). Then, we derived the ellipse’s semi-major axis, semi-minor axis, and position angle from the coefficients.

3.2 Method 2: Principal component analysis

Principal component analysis is a method that aims to reduce the number of dimensions in large datasets in order to extract the most important information (Pearson 1901). In our case, the PCA is useful for finding the axes with the most information encoded.

The dataset of stellar positions is represented by the matrix D: $$\mathbf{D}=\left(\begin{array}{cc}{x}_1& y_1\\ x_2& y_2\\ ⋮& ⋮\\ x_n& y_n\end{array}\right).$$(2) Each row (x_i, y_i) represents the position of a star in a 2D Cartesian coordinate system with the origin in the centre of the cluster. We implemented PCA using the singular value decomposition (SVD), which decomposes D into three matrices: $$\mathbf{D}=\mathbf{U}\Sigma {\mathbf{V}}^T.$$(3) In this decomposition U} is an N × N orthogonal matrix, where N is the length of the data matrix D, in our case the number of stars in our dataset. Σ is a diagonal matrix with the singular values σ₁ and σ₂, which are related to the square roots of the eigenvalues of the covariance matrix of the data. V}^T is a 2 × 2 orthogonal matrix, where each column represents the principal direction (eigenvector) of the data. The singular values σ₁ and σ₂ represent the spread of the data along each principal component. The square of each singular value gives the variance along each principal component, which is proportional to the eigenvalues of the covariance matrix. This decomposition thus provides the principal directions of the data and their associated variances.

In PCA, the eigenvector associated with the largest eigenvalue (or the largest singular value σ₁) is the direction along which the variance is maximised. This direction is identified as the major axis of the data distribution. Conversely, the eigenvector associated with the smallest eigenvalue (or the smallest singular value σ₂) corresponds to the minor axis. The ellipticity and the position angle are then computed as $$e=1-\frac{{\sigma }_2}{{\sigma }_1}$$(4) $$\Phi =arctan\left(\frac{v_{1y}}{v_{1x}}\right),$$(5) where v_1x and v_1y are the components of the eigenvector associated with the largest singular value (or eigenvalue), indicating the direction of maximum variance.

Fig. 1 Recovered ellipticity as a function of the true ellipticity using the KDE method (upper left), the PCA method (upper right), the robust PCA method (lower left), and the robust PCA method after bias correction (lower right). The points are coloured according to the number of stars in the mock GCs.

3.3 Robustness and accuracy of the methods

In this section we test the two different methods described above on mock GCs while varying three parameters: the number of stars (N), the true ellipticity (e_true), and the true angle of the ellipse (Φ_true), defined as the angle between the x-axis and the semi-major axis of the ellipse. We also tested the robustness of the methods with respect to the presence of outlier stars. These tests are particularly important to ensuring that the selected method provides accurate and meaningful information on the global ellipticity of GCs.

We generated mock positions of stars following a King profile (King 1966). Then, we applied a transformation to the dataset to impose an ellipticity and an angle. Each test was realised 1000 times, on slightly different mock GCs with the same imposed ellipticity and position angle. We fixed the seed of the random generator of stars to make our tests reproducible and comparable. The mean value and the error bars on the parameters were obtained by averaging and taking the standard deviation of the realisations.

In the first test that we ran, our aim was to quantify the impact of the number of stars on the error on the ellipticity for different values of the true ellipticity. To do so, we uniformly sampled 300 times the number of stars N between 50 and 1500, the true ellipticity e_true between 0.01 and 0.3, the true angle Φ_true between 0^° and 180^°, and the core radius r_c between 0.5 arcmin and 15 arcmin. We used these parameters as inputs to mimic the distribution of stars in the GCs from our sample. We present the results in the upper panel of Fig. 1, where we plot the recovered ellipticity e_rec as a function of the true ellipticity e_true. The points are coloured according to the number of stars. In the left panel, the ellipticity was recovered using the KDE method, and in the right panel using the PCA method. In both cases, two biases are present. First, the number of stars strongly impact the error on the ellipticity. When N<300 (dark blue points), both methods overestimate the true ellipticity, visible as the large number of points located above the black identity line in Fig. 1. The other source of bias comes from small ellipticities. If the cluster is close to spherical, both methods tend to overestimate the ellipticity. This is visible as larger errors for smaller e_true in Fig. 1. The position angle of the ellipse, not shown here, is always accurately recovered, except for very small values of the ellipticity, where the position angle is not defined.

The second test that we ran was a robustness check. We manually added ten outlier stars to mock GCs containing 1000 stars (1% of outliers), and inspect the mean value of the ellipticity and position angle recovered using the KDE and PCA methods. We present the results in Fig. 2 for mock GCs with e_true=0.15 and Φ_true=45^° (top panel), and with e_true=0.25 and Φ_true=60^° (bottom panel). The values for N and e_true have been chosen to avoid including any of the biases presented above. We ran this test with two types of outliers. The first type of outliers were added in the direction of the semi-major axis (yellow points), and the outliers of the second type were added in the direction of the semi-minor axis (pink points). In the first case, adding 1% of outlier stars has a strong impact on the recovered ellipticity regardless of the method used and of the value of the true parameters. The ellipticity was overestimated due to the presence of these outliers along the semi-major axis. In the second case, when outliers are added along the semi-minor axis, the recovered position angle was shifted by 10^° to 90^°, and the ellipticity was underestimated by 0.05 to 0.2, depending on the method. This means that outliers distributed perpendicular to the intrinsic flattening tend to rotate the recovered major axis and make the cluster appear rounder than it actually is.

These tests on mock data suggest the necessity of developing a robust and accurate method to measure the ellipticity in GCs. This will be very important when measuring the ellipticity of clusters where the number of stars available is low, for example, when studying the shape of multiple stellar populations in GCs. One of the main advantages of the PCA method is its speed. Therefore, we decided to adapt the PCA to make it robust to outliers and to correct for the biases present in small datasets and for round clusters.

Fig. 2 Impact of outliers on the recovered ellipticity and position angle for the KDE (circles), PCA (squares), and robust PCA (diamonds) methods for two sets of true parameters. The colour of the points indicates the location of the outliers. The points in pink with a positive x and negative y correspond approximately to the direction of the semi-minor axis, and the points in orange with a positive x and positive y correspond approximately to the direction of the semi-major axis of the ellipse. Outliers were added within fixed quadrants relative to the cluster centre, not strictly along the ellipse axes, which leads to small asymmetries in the recovered angles when the ellipse position angle changes. The dashed green line represents the true value of the ellipticity and angle, taken to 0.15 and 45° (top panel) and to 0.25 and 60° (bottom panel). We note that for the Robust PCA method, only one point is visible, thus highlighting the robustness of the method. Regardless of the direction where outliers are added, the ellipticity and position angle are accurately recovered.

3.4 Robust PCA method and bias correction

The robust PCA uses a distinct approach to compute the covariance matrix, relying on the minimum covariance determinant method (MCD; Rousseeuw 1984). This method identifies a subset of the data with the smallest possible covariance determinant, which measures the spread of the distribution. By minimizing this determinant, the method selects the subset of data that is most tightly clustered, thus reducing the influence of outliers on the final covariance estimate. We used the scikit-learn class MinCovDet² to compute the covariance matrix using the MCD method. Then, as done in Sect. 3.2, we computed the eigenvalues and eigenvectors of the covariance matrix and estimated the ellipticity and angle of the data.

To test the robustness of this method, we performed the same test as in Sect. 3.3. We added the ellipticities recovered using the robust PCA in Fig. 2 as diamonds. The recovered ellipticities and angles are almost superimposed on the lines indicating the true values of these parameters. This agreement holds regardless of the directions in which outliers were added. This method is thus very robust to the presence of outliers. However, we find that it is still sensitive to the number of stars, as attested by the lower-left panel of Fig. 1, which shows the recovered ellipticity e_rec as a function of the true ellipticity e_true. In the figure, points are coloured according to the number of stars. Therefore, another correction should be applied. We constructed a bias correction function using a set of mock data to correct for the systematic biases linked to small datasets (N ≃ 100) and to the overestimation of the ellipticity. More information on the bias correction we performed is available in Appendix A. The correction of the bias function is visible in the lower-right panel of Fig. 1. After applying the correction, most of the points fall close to the 1: 1 black line. The bias linked to a small number of stars is efficiently tackled.

4 Ellipticity of 29 MW GCs

To compute the global ellipticity and position angle, we randomly sampled 90% of all the stars in each GC 1000 times. We then used the mean value and standard deviation of the ellipticities and position angles obtained for these realisations as measurements and corresponding error estimates for these quantities.

4.1 Comparison with literature

In this section we present the ellipticity and position angle values obtained for the 29 Galactic GCs in our sample by using the robust PCA method described in Sect. 3.4. For clusters with N_stars<1500, we further applied a correction model to correct for the biases linked to small datasets (see Appendix A). Then, we compare these values with those from WS87, CC10, and CR24.

In our analysis, we refer to GCs as round or spherical when they have ellipticities of e<0.05 and as flattened or elliptical otherwise. According to this criterion, among our 29 GCs, we found 16 flattened GCs.

The results of our analysis are presented in Table C.1, where we provide the values of ellipticity e, together with ellipticity measured by WS87, CC10, and CR24 (e_WS87, e_CC10, and e_CR24, respectively). We compare our results with WS87, CC10, and CR24 in the left, middle, and right panels of Fig. 3, respectively. Overall, our values show better agreement with WS87, with reasonable consistency within 3 σ, except for five GCs: NGC 6981, NGC 6838, NGC 6366, NGC 6656, and NGC 7089. Particularly striking is the situation of NGC 6838, NGC 6656, and NGC 6366. The first appears to be round in WS87, but we find it to be very flattened, and for the second and third, the situation is reversed. The analysis by CC10 resulted in higher ellipticity values in general, which differ from our findings. This discrepancy may arise because their measurements are more sensitive to bright stars and the outer regions of globular clusters, where flattening tends to be more pronounced. The case of NGC 6838 is once again particularly interesting, as CC10 found it to be very flattened, while WS87 found it to be round. Our measurement appears to be somewhat in between the estimates of the two works. We discuss these results further in Sect. 4.2.2.

Finally, CR24 computed the ellipticity of 163 globular clusters using the on-sky distribution of cluster members and the PCA method. We compare our results with theirs in the right panel of Fig. 3. Our results differ for the majority of clusters we have in common, as we find higher values of ellipticity for 14 out of 26 clusters. These discrepancies can be linked to two points. First, CR24 used the right ascension and declination of stars for the PCA analysis, while we used the Cartesian coordinates. Computing the covariance matrix from stellar positions in RA and Dec will indubitably generate projection effects and errors in the distance estimate, which will propagate to the measurement of ellipticity. Second, as demonstrated in Sect. 3.3, only a robust implementation of the PCA can prevent the impacts caused by sensitivity to outliers. So if outliers are present in the data used by CR24, this will impact the resulting ellipticities. This issue also applies to methods based on fitting density contours, although the impact of outliers is generally less visible than in the case of PCA (see Fig. 2).

Fig. 3 Ellipticities of the 29 Galactic GCs in our sample measured using the robust PCA method compared with the estimates from the literature provided by White & Shawl (1987, left), Chen & Chen (2010, middle), and Cruz Reyes & Anderson (2024, right). The solid lines indicate equal ellipticities. Each GC is represented by a coloured point (see the legend at the top), and error bars are indicated.

4.2 The role of relaxation time, rotation, velocity anisotropy, and tides

In this section we explore the factors that can influence the shapes of GCs. The first quantity that we investigated is the relaxation time. It can be defined as the timescale on which encounters between stars, leading to energy exchange, produce major changes in the structure of the cluster, without disturbing its dynamical equilibrium (Meylan & Heggie 1997). Our results are shown in Fig. 4, where we put into relation the ellipticity that we measure for the 29 MW GCs in our sample and the relaxation time at the half-mass radius taken from Harris (1996)³. The GCs are coloured according to the time of the last disc crossing from Pancino et al. (2024). A trend is visible in the figure: clusters with longer relaxation times (dynamically young) are, in general, more elliptical, and those with short relaxation times (dynamically old) are rounder. An evident exception is NGC 6838, which is very flattened (e=0.14), even if it appears to be dynamically old. As mentioned in Sect. 4.1, its high flattening might be related to its recent crossing through the Galactic disc.

These findings indicate the importance of understanding the origin of the deviation from spherical symmetry of GCs. To identify the primary driver behind the observed flattening in GCs (i.e. whether it is internal rotation, velocity anisotropy, or tidal effects), we first considered the role of rotation.

The V/σ diagram is a tool frequently used to study elliptical galaxies (Cappellari et al. 2007; Emsellem et al. 2011; Brough et al. 2017). The relation between projected rotation signature and flattening is derived from the tensor virial theorem (Binney 1978) and for an edge-on system is given by $${\left(\frac{V}{\sigma }\right)}_{ed}=k_1\sqrt{\frac{e}{1-k_{2}e}},$$(6) where the subscript e d indicates the edge-on case, V is the maximum rotation velocity along the line-of sight, σ is the central velocity dispersion, e is the ellipticity, and k₁ and k₂ are constants that account for the intrinsic properties of the system, including its geometry and kinematics. In our analysis, we assume k₁=0.831 and k₂=0.896 (Cappellari et al. 2007; Bianchini et al. 2013).

Such a tool has already been applied to GCs to understand the link between flattening and rotation or anisotropy. Using this diagram with ellipticities from WS87, Fabricius et al. (2014) concluded that rotation is the driving factor of the ellipticities, as higher flattening exhibits higher V/σ. However, their analysis did not consider projection effects, which have a significant impact when i<30^°. For example, Bianchini et al. (2013) investigated three clusters: ω− Cen, NGC 104 and NGC 7078. By placing them in the V/σ diagram and correcting for their inclination, they showed that the positions of all three clusters were shifted, demonstrating the impact of projection effects. They highlighted that deviation from the line of isotropic oblate rotators could be a hint that several mechanisms are acting together, such as inclination, differential rotation, and pressure anisotropy. Kamann et al. (2018) found a similar trend in this diagram; increasing flattening with higher rotational signature, though with a lower significance, using ellipticities from WS87.

The quantities presented in Eq. (6) need to be corrected to account for the inclination of the GCs. This can be done as follows (Binney & Tremaine 1987): $${\left(\frac{V}{\sigma }\right)}_{ed}={\left(\frac{V}{\sigma }\right)}_{obs}\frac{\sqrt{1-\delta {cos}^{2}i}}{sini}$$(7) $$e_{\mathrm{intr}}=1-\sqrt{1+e(e-2)/{sin}^{2}i}$$(8) $$\delta =1-\frac{1+{(V/\sigma )}^2}{[1-\alpha {(V/\sigma )}^2]\Omega \left(e_{\mathrm{intr}}\right)}$$(9) where i represents the inclination angle, δ the anisotropy and Ω(e_intr) is a function of the cluster’s intrinsic ellipticity e_intr or inherent shape, excluding projection effects, provided in Eq. (16) of Cappellari et al. (2007). α is a dimensionless parameter depending on the streaming velocity and the stellar density.

Similarly to Cappellari et al. (2007), we adopt a fixed value of α=0.15. For the small ellipticities and rotation strengths relevant to our clusters, varying α has a negligible effect on the results (see Fig. 2 of Binney 2005).

In our analysis, we use the amplitude A of the best-fit rotation curve as a proxy for the projected rotational velocity V, as commonly adopted in past studies (Bianchini et al. 2013; Kamann et al. 2018). We present our results in Fig. 5, where we plot our ellipticity values and the order-to-random motion parameter A_LOS/σ from Leitinger et al. (2025). Three clusters, NGC 5053, NGC 6934, and NGC 6981, lack a A_LOS/σ measurement and are therefore not included in the plot. The yellow line indicates the theoretical relation for an isotropic oblate rotator viewed edge-on (Eq. (6)). The other colourful lines show the relation between V/σ and e when i ≠ 90^°, as described by Eqs. (7), (8), and (9). Each point is coloured according to the inclination angle of the GC. Inclinations in Leitinger et al. (2025) are reported over 0−180^°, but for our analysis we adopt the standard 0−90^° convention⁴. Leitinger et al. (2025) inclination angles are reported in the last column of Table C.1. Therefore, each GC should lie in the corresponding coloured line if they behave as isotropic oblate rigid rotators. While most of the clusters follow the line corresponding to their inclination in the framework of the isotropic oblate rotator, a few clusters deviate from their expected location. To further analyse these results, we divided the clusters into three groups. In the first group, we gather the GCs with a flattening coherent with their rotational signature (Sect. 4.2.1). In the second group discussed in Sect. 4.2.2, we include the clusters more flattened than expected. Finally, the third group gathers clusters less flattened than expected (Sect. 4.2.3). All clusters have uncertainties associated with their inclination angles, also reported in Table C.1. These uncertainties should be taken into account in the interpretation of the results, as large uncertainties can lead to a wide range of values for sin i in Eqs. (7) and (8), resulting in variations in the intrinsic rotation strength and flattening, therefore challenging the interpretation of the results. The most affected clusters are those with low inclination and large uncertainties, such as NGC 6838.

Fig. 4 Ellipticity (e) of 29 Galactic GCs obtained using the robust PCA method against the relaxation time at the half-mass radius, log th, from Harris (1996). The points are coloured according to the last crossing time through the Galactic disc, log tc, taken from Pancino et al. (2024). The name of each cluster is indicated close to the corresponding point. Both times are given in years. Dynamically young clusters tend to be more elliptical.

Fig. 5 Rotation strength, taken from Leitinger et al. (2025), as a function of the flattening we measured. The names of GCs are indicated close to the corresponding points and coloured according to the inclination angle of the cluster, also taken from Leitinger et al. (2025). Low inclinations are plotted in blue while high inclinations are coloured in yellow. The uncertainties associated with the inclination angles are not shown in this figure. Instead, we report them in the last column of Table C.1. The dashed lines indicate the relation for isotropic oblate rotators viewed at different inclination angles.

Fig. 6 Comparison of the position angle of the rotation axis from Leitinger et al. (2025) and the position angle measured in this work. The black line indicates the identity. In this case, this represents an alignment between the semi-minor axis of the ellipse and the rotation axis of the cluster. Here, we have only included the clusters classified as flattened and with an available position angle for the rotation axis. Some clusters are classified as non-rotating in Leitinger et al. (2025), but they still have a measured rotation strength and position angle for the rotation axis. They are plotted with a star symbol.

4.2.1 Group 1: GCs with a flattening coherent with their rotational signature

The majority of the GCs appear to have the ellipticity and rotation strength consistent with what is expected in the isotropic oblate rotator framework, given their inclination. This is the case of the spherical and low or non-rotating clusters: NGC 288, NGC 6218, NGC 6254, NGC 6366, NGC 6809, and NGC 7099.

Among the GCs with a flattening coherent with their rotational signature, eleven are elliptical: NGC 1261, NGC 3201, NGC 4590, NGC 5024, NGC 5272, NGC 5286, NGC 5986, NGC 6121, NGC 6205, NGC 6341, and NGC 7078. Their ellipticity is in good agreement with the one expected in the isotropic oblate rotator framework, given their rotational signature, indicating the major role of rotation in shaping these clusters. As further confirmation that rotation is the main driver of the flattening in these clusters, we analyse the orientation of their flattening axis and their rotation axis. In Fig. 6, we show the position angles of the ellipse we measure, Φ, and compare them to the rotation axis position angles, θ₀, from Leitinger et al. (2025). To compare our measured ellipse position angles with the rotation axis position angle, we account for the difference in reference frames. Our PCA-derived angles are measured counter-clockwise from the west direction, whereas the angles in Leitinger et al. (2025) are measured counter-clockwise from north. To bring both measurements to the same reference, we apply a 90^° shift to their angles. After this adjustment, we wrapped the rotation angles to the [0^°., 180^°) interval. We also subtracted 90^° from the ellipse’s position angle to plot the semi-minor axis instead of the semi-major axis, allowing for a direct comparison with the rotation angle. Here, we only consider the clusters flagged as elliptical (e>0.05). Non-rotating clusters with a measured A/σ in Leitinger et al. (2025) are plotted with a star symbol. Their rotation angle might be unreliable. The black line indicates the identity: a cluster falling on this line would suggest that its flattening could be due to rotation, as the semi-minor axis of the cluster aligns with the axis of rotation. Five clusters show position angles consistent with the identity line within 1 σ, suggesting a strong alignment between the rotation and flattening axes: NGC 3201, NGC 5986, NGC 6205, NGC 6341, and NGC 7078. That is also the case of NGC 2808 and NGC 104, but their position in the V/σ diagram prevents us from constructing a definitive conclusion. Their case is examined in more detail in Sects. 4.2.2 and 4.2.3. We found that the geometric position angles and rotational position angles are consistent within 2 σ for the clusters NGC 5024, NGC 5286, and NGC 5904. We conclude that the flattening in these eight clusters is mainly due to rotation. Finally, the position angle of the rotation axis and the semi-minor axis of the ellipse differ from the identity for NGC 5272. However, this cluster still lies along the isotropic oblate rotator relation in the A/σ−e diagram, suggesting that its flattening is primarily driven by rotation, and that the observed misalignments may result from projection effects rather than anisotropy.

4.2.2 Group 2: Clusters that are more flattened than expected

Three clusters are located well below their expected position in Fig. 5, with a much lower rotation strength (A/σ) than expected for their ellipticity: NGC 2808, NGC 4833, and NGC 6838. Many studies have confirmed rotation in NGC 2808 (Bellazzini et al. 2012; Lardo et al. 2015; Jeffreson et al. 2017; Kamann et al. 2018; Sollima et al. 2019; Martens et al. 2023). Gaia DR2 observations have revealed a population of extra-tidal candidate stars. These stars show a clear over-density in the trailing direction of the cluster’s orbit, a feature commonly associated with tidal stripping and disc-shocking events (Kundu et al. 2021). Deep DECam imaging further revealed a substantial population of stars beyond the King tidal radius, as well as the presence of extended tidal tails with complex morphology (Carballo-Bello et al. 2017). These findings suggest that NGC 2808 resides in a dynamically rich region of the Galaxy and has likely undergone significant tidal interactions, which may have contributed to its enhanced flattening. As we find an excellent agreement between the ellipse minor axis and the rotation axis of the cluster in Fig. 6, we conclude that the high flattening of NGC 2808 is likely a combination of internal rotation and tidal interactions.

The cluster NGC 6838 is a non-rotating cluster (Bianchini et al. 2018; Vasiliev 2019; Leitinger et al. 2025). It is dynamically old, and therefore isotropy throughout the cluster is expected as a natural consequence of dynamical relaxation. We find it to have a significant flattening. This suggests that additional mechanisms such as tidal interactions play a role in determining its shape. Chen & Chen (2010) suggested that the high flattening of NGC 6838 might be related to its recent crossing through the Galactic disc some 16 Myr ago (Vande Putte & Cropper 2008), much shorter than the relaxation time of the cluster, which is 270 Myr (Harris 1996). More recent estimates of the time since the last galactic disc crossing have been computed by Pancino et al. (2024), using proper motion based on Gaia DR3. They found a last disc crossing time of 31.2 Myr for NGC 6838, almost twice the value from Vande Putte & Cropper (2008). If the cluster was affected by disc shocking during this passage, we might also expect some stripped stars to remain in its vicinity. Although Ibata et al. (2021) did not list NGC 6838 among clusters with detected extra-tidal structures, this is likely due to its position at low Galactic latitude, where stream detectability is more challenging. Disc shocking could also disperse stars more diffusely, preventing coherent structures such as stellar streams from forming. Chen & Chen (2010) identify some extra-tidal member stars with a clumpier spatial distribution, and Ferrone et al. (2023) predict that NGC 6838 should experience substantial mass loss due to repeated disc crossings. Given its present low mass (∼ 5 × 10⁴ M_⊙; Baumgardt & Hilker 2018), it is plausible that this process has already removed a significant fraction of its original stellar content. We therefore consider the remnant tidal effects from disc crossings to be a compelling explanation for the apparent flattening of NGC 6838, supported by its orbital history, low mass, and indirect evidence of past mass loss.

The case of NGC 4833 is particularly interesting. We find a flattening near 0.14, in good agreement with CC10, but the cluster is non-rotating (Lardo et al. 2015; Bianchini et al. 2018; Vasiliev 2019; Leitinger et al. 2025). So rotation cannot account for its ellipticity. This cluster has a last crossing time through the Galactic disc among the shortest in our sample, as shown by its dark blue colour in Fig. 4. A reason for its high flattening could be tidal effects from this last disc crossing. However, the results on the ellipticity using all stars are greatly different for this specific cluster, as reported in Appendix B. Dedicated simulations would be valuable to further investigate this complex case.

While most clusters have uncertainties on their inclination angle which do not significantly impact the results and their interpretation, NGC 4833 and NGC 6838 both exhibit large uncertainties: the inclination angle of NGC 4833 could vary between 35^° and 87^° (almost edge-on), and between 4^° (face-on) and 64^° for NGC 6838. This wide range results in a large variation of sin i, leading to significant potential variations in the intrinsic flattening and rotational strength. Therefore, the interpretation for these clusters should be considered with caution due to these large uncertainties.

4.2.3 Group 3: Clusters that are less flattened than expected

pt The clusters NGC 104, NGC 1851, NGC 5904, NGC 6101, NGC 6656, and NGC 7089 are located above the isotropic rotator line corresponding to their inclination, as shown in Fig. 5, suggesting that they are less flattened than expected given their rotation. This might suggest that another (or more) mechanisms contribute to the cluster shape and counteract the effect of rotation.

The two clusters NGC 104 and NGC 5904 both lie close to the edge-on line, above their expected location. Inspecting their position in Fig. 6, the semi-minor axis of NGC 104 is perfectly aligned with its rotation angle. It is also the case of NGC 5904 within 1.5 σ. This suggests rotation as the primary driver of the flattening in these two clusters, with a secondary mechanism partially counteracting this effect. NGC 104 exhibits radial anisotropy in the outer regions of the cluster, defined as the intrinsic random motion between stars, as opposed to the anisotropy induced by the fast rotation of the cluster (Milone et al. 2018; Cordoni et al. 2020; Libralato et al. 2022; Leitinger et al. 2025). This leads us to conclude that the flattening of NGC 104 could be linked to both rotation and radial anisotropy.

We note that our results for this cluster are in very good agreement with Bianchini et al. (2013). NGC 5904 exhibits either isotropy (Libralato et al. 2022) or mild tangential anisotropy (Leitinger et al. 2025). Its most recent passage through the Galactic disc is among the oldest of the globular clusters in our sample. Neither anisotropy nor tidal effects seem to account for the shape of the cluster. Given that its position in the V/σ diagram is consistent within two sigma with the isotropic oblate rotator line at the corresponding inclination angle, and that the semi-minor axis aligns closely with the rotation axis, we conclude that rotation is the primary driver of the cluster’s flattening.

The cluster NGC 1851 is a dynamically old cluster. Several studies have detected rotation within the cluster (Bellazzini et al. 2012; Lardo et al. 2015; Kamann et al. 2018; Martens et al. 2023; Petralia et al. 2024; Leitinger et al. 2025), while others have not found evidence of it (Bianchini et al. 2018; Vasiliev 2019; Sollima et al. 2019). The cluster is characterised as isotropic (Leitinger et al. 2025). Given that the cluster’s position in Fig. 5 agrees within 2−σ with the expected location in the context of the isotropic oblate rotator, and that discrepancies have been found regarding the rotation pattern, we conclude that there is no conflict with our results.

Leitinger et al. (2025) detected for the first time a low rotation signal in NGC 6101. They also highlighted the presence of radial anisotropy in the outer regions of the cluster. The cluster exhibits peculiar dynamics, characterised by the absence of mass segregation (Dalessandro et al. 2015) and a moderate likelihood of hosting a large population of black holes (Peuten et al. 2016; Askar et al. 2018; Weatherford et al. 2020). Tidal tails have also been detected (Ibata et al. 2021). The large error bars associated with the rotation strength and the ellipticity, combined with a complex internal dynamics and tidal interaction events, lead us to conclude that several parameters are likely to contribute to the shape of NGC 6101.

Notably, NGC 6656 is a well-studied GC (Lane et al. 2010; Kamann et al. 2018; Bianchini et al. 2018; Vasiliev 2019; Sollima et al. 2019; Leitinger et al. 2025). It exhibits a high rotational signature and relatively little ellipticity. We note that both WS87 and CC10 found an ellipticity greater than 0.1 for this cluster, which would shift its position to the right in Fig. 5, making it compatible with the results expected in the framework of the isotropic oblate rotator. The last galactic disc crossing time has been estimated to be 2.7 Myr (Pancino et al. 2024), placing NGC 6656 to be the latest GC in our sample crossing the Galactic disc. This cluster is a good candidate in the search for tidal tails (Kunder et al. 2014). It also hosts groups of stars with differences in heavy elements abundances (Marino et al. 2009; Mucciarelli et al. 2015), suggesting a potential extragalactic origin or that it is a remnant of a disrupted dwarf galaxy (Da Costa et al. 2009). However, according to its kinematics, NGC 6656 is expected to be an in-situ cluster born in the MW main disc (Massari et al. 2019). In this context, the complex dynamics and history of NGC 6656 could explain its deviation from the isotropic oblate rotator framework exhibited in Fig. 5.

Significant rotation has also been detected in NGC 7089 (Kimmig et al. 2015; Kamann et al. 2018; Bianchini et al. 2018; Vasiliev 2019; Sollima et al. 2019; Martens et al. 2023), as well as tidal tails (Grillmair 2022). The cluster exhibits mild radial anisotropy in the outer regions (Leitinger et al. 2025). A possible explanation for its low flattening and high rotation could be the preferential stripping of stars on radial orbits during tidal interactions, as attested by the presence of tidal tails. This could lead to a rounder appearance of the cluster, despite the significant internal rotation.

Fig. 7 Ellipticity computed using the robust PCA method with all stars (x-axis) and with all stars within 3rh (y-axis). Only four clusters show ellipticities incompatible within their respective error bars.

4.3 Impact of radial coverage

In the results presented in Sect. 4.2, the ellipticity has been computed using all the RGB stars available in our sample. This results in an inhomogeneous radial coverage between clusters, with some clusters being covered up to 2 r_h (e.g. NGC 6121 and NGC 6656) and others up to 19 r_h (e.g. NGC 1261 and NGC 2808). To evaluate the impact of this inhomogeneous radial coverage, we restricted the data to 3r_h and removed all the stars with a distance to the cluster centre superior to 3r_h. We note that such a cut can also introduce its own bias towards lower ellipticity, as the centres of clusters lose their initial conditions faster than their outer regions. We plot in Fig. 7 the total ellipticity e_tot, as reported in the results from Sect. 4.2, versus the ellipticity within 3r_h. We remove NGC 6121 and NGC 6656 from the sample, as their radial coverage extends only to 2r_h. We highlight the names of the clusters for which the ellipticity results deviate by more than 1.5 σ, indicating a significant discrepancy beyond expected uncertainties. NGC 104, NGC 1261, NGC 6341 all have an ellipticity within 3r_h significantly lower than when using all the stars. The two latest are among the clusters with the most extensive radial coverage. This difference could be attributed to the fact that the inner regions of clusters tend to relax faster than their outer parts. Consequently, restricting the analysis to within 3r_h excludes the outer, more elliptical regions, leading to a lower overall measured ellipticity. Regarding NGC 104, the analysis from Sect. 4.2.3 revealed that a mechanism is likely counteracting the effect of rotation. The ellipticity within 3 r_h leads to a similar conclusion.

The cluster NGC 5986 is the only cluster with a statistically significant higher ellipticity within 3 r_h. This is peculiar, as we expect the centre of a cluster to relax faster and the outer parts to be more subject to tidal effects, potentially making them more elliptical. Lanzoni et al. (2018b) detected peculiar kinematics for this cluster, with a solid-body rotation pattern, and a velocity dispersion profile flattening after one half-mass radius. In this scenario, the inner regions of the cluster could preserve the effects of rotation more strongly. When stars at large radii are included in the analysis, their isotropic kinematics could dilute the rotation-induced flattening in the central parts, resulting in a lower ellipticity. Such a result remains compatible with the analysis from Sect. 4.2.1, which shows that NGC 5986’s high flattening is likely due to rotation. We therefore conclude that the interpretation of the ellipticity results discussed in Sect. 4.2 remains robust when adopting a homogeneous radial coverage between clusters.

5 Conclusion

We have used photometric ground-based observations (Stetson et al. 2019) combined with HST data (Nardiello et al. 2018; Piotto et al. 2015) to model the flattening of 29 Galactic GCs. To do this, we implemented and investigated two methods that provide a measurement of the flattening and of the position angle of the elliptical shape representing the cluster. The first method relies on a kernel density estimate of the distribution of stars and an ellipse fitting to different isodensity contours. This approach has the advantage of providing an ellipticity profile and is particularly interesting for clusters for which we expect a variation of the flattening with the radial distance. However, the method is sensitive to outliers, which can considerably alter the results. The second method we investigated is PCA, which captures the direction along which the variance of the data is minimum and maximum. Though PCA provides results very quickly, it is sensitive to outliers. The two methods that we tested are not robust and not accurate, as they systematically overestimate small ellipticities, leading to spurious flattening measurements. Therefore, we further developed and improved the fastest of these methods, PCA, to make it robust to the presence of outliers and correct for the bias linked to small datasets. After carefully testing it on mock data, we applied the robust PCA to the 29 MW GCs in our sample.

We reported deviation from spherical symmetry for 16 GCs. Using the V/σ diagram, a tool frequently used to link the rotation of elliptical galaxies to their ellipticities, we further analysed the dominant mechanism shaping GCs. For NGC 3201, NGC 5024, NGC 5286, NGC 5904, NGC 5986, NGC 6205, NGC 6341, and NGC 7078, their location in the V/σ diagram combined with a good alignment between the semi-minor axis and the rotation axis lead us to conclude that rotation is the main driver of the flattening in these eight GCs.

The clusters NGC 104 and NGC 5904 also exhibit a good alignment between their semi-minor axis and rotation axis, reinforcing the idea that rotation plays a dominant role in shaping these clusters, even though their flattening is lower than expected from their rotational signature in the isotropic oblate rotator framework. For other clusters, such as NGC 6838 and NGC 6656, tidal interactions from the galaxy and/or complex formation history are a plausible explanation for their apparent shape.

We note that the PCA method we developed is both robust and accurate, but it does not account for possible radial variations in ellipticity. Capturing such trends would require extending the method to include the radial binning of stars and constructing ellipticity profiles. Although feasible, this extension requires careful consideration to ensure the reliability of the recovered profiles, which goes beyond the scope of this study. To ensure that our results and interpretations remain valid using the different radial coverage available, varying between 1.9 r_h and more than 20 r_h for some clusters, we computed the ellipticity using all the stars within a 3 half-light radius for 27 clusters with a radial coverage superior to this limit. Apart from a few clusters discussed in Sect. 4.3, which require careful considerations, we find that the interpretation and discussion of our results remain largely valid.

While the isotropic oblate rotator model is a useful way to compare ellipticity and rotation, it was first developed for elliptical galaxies, which usually have enough rotation for the model to apply. Most non-rotating GCs also match the framework quite well, showing low flattening (as expected). An exception is NGC 6838, which we discussed extensively in Sect. 4.2.2. Clusters exhibiting anisotropy also tend to agree within the error-bars, except for NGC 104, whose radial anisotropy likely explains its offset from the relation. This means the framework can still be a good first-order guide for GCs, but deviations could be expected for clusters exhibiting anisotropy.

Obtaining robust flattening measurements and rotation strengths of a significant fraction of MW GCs would be interesting for testing of the isotropic oblate rotator framework and to better constrain the causes of the flattening. From a modelling point of view, it would be crucial to take into account the kinematic complexity of these systems, including, for example, rotation and pressure anisotropy in dynamical models of GCs. Moreover, a careful exploration of the time evolution of these systems is needed to identify which mechanisms play a significant role in shaping the properties of present-day GCs.

In this first paper of a two-part series, our main goal was to present the method and general results. In a follow-up paper, we will apply this method to study ellipticity differences among multiple stellar populations in GCs, where the available data are limited, hence the need for a method that performs reliably when using small datasets (Fréour et al., in prep.).

Data availability

The Python codes used to obtain the results presented in this article are available upon request.

Acknowledgements

The authors are very grateful for the referee’s careful reading of the manuscript and their valuable comments on the issue of photometric incompleteness and inhomogeneous radial coverage, which significantly improved the quality of this paper. The authors would like to thank Francisco Aros for the helpful comments on the manuscript and Paolo Bianchini for the discussions. EL acknowledges support from the ERC Consolidator Grant funding scheme (project ASTEROCHRONOMETRY, https://www.asterochronometry.eu, G.A. no. 772293). EP and EL acknowledge funding by the European Union (ERC-2022-AdG, “StarDance: the non-canonical evolution of stars in clusters”, Grant Agreement 101093572, PI: E. Pancino). Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Research Council. Neither the European Union nor the granting authority can be held responsible for them. This research has made use of NASA’s Astrophysics Data System (ADS) and of the SIMBAD database operated at the Strasbourg astronomical Data Center (CDS), Strasbourg (France).

References

Appendix A Bias correction

Our aim is to construct a bias correction function using a set of mock data, to correct for the systematic biases linked to small datasets (N ≃ 100) and to the overestimation of the ellipticity for clusters close to spherical. To do this, we proceeded as follows. We fixed the angle to Φ_true=45^° because, contrary to the ellipticity, changing its value does not impact the results. We generated mock data of GCs with true ellipticity varying from e_true=0.01 to e_true=0.2 with steps of 0.01, and a total number of stars N going from N=40 to N=1500 in steps of 20. For each pair of parameters (e_true, N), we created 500 different GCs using a random distribution of stars following a King (1962) profile. We fix the core radius of the clusters to the median value of our sample r_c=0.6 arcmin, taken from Harris (1996). The truncation radius is also fixed (r_t=17.5 arcmin) and computed from the median value of the concentration of the clusters, also taken from Harris (1996). Using the robust PCA algorithm, we computed the ellipticity e_rec as the mean value of the ellipticity of the 500 realisations, and the associated error by taking the standard deviation. The problem is now a classical regression problem where we seek to find the true ellipticity (target) of an elliptical distribution of stars given a total number of stars and the (biased) recovered ellipticity (features). Several methods can be used to solve such a problem. Among the ones we tried are the polynomial regression (where the relationship between the input variables and the output is a n-degree polynomial) and the random forest regression (which creates multiple decision trees, each learning slightly different patterns from the data by using random subsets of data and features, and then averages their predictions to produce a final output). We find that the random forest regression brings a better correction than the polynomial regression. We also fit a regression model to the residuals to make the correction more accurate, using a Gradient Boosting Regressor, better at capturing complex relationships in data. Figure A.1, representing the true ellipticity versus the ellipticity recovered, illustrates the results of this correction. The blue dots indicate the value of the ellipticity recovered by the robust PCA before correction, while the orange circles represent the value of the ellipticity after bias correction.

Fig. A.1 Ellipticity recovered using the robust PCA method before (blue dots) and after (orange dots) correction for the bias due to a small number of stars. The grey area shows the uncertainties in the recovered ellipticity. The dashed line indicates the true ellipticity. The mock GCs were simulated using 300 stars.

Appendix B Impact of photometric incompleteness on the results

In this article we use RGB stars only to compute the ellipticity. For most clusters, the full stellar samples are affected by some degree of photometric incompleteness, particularly at the faint end (due to detection limits).

We evaluate the impact of this effect on four representative clusters spanning different sample sizes and ellipticities: NGC 6838 (>100,000 stars in the full sample, e=0.107 ± 0.014), NGC 6254 (14,000 stars, very round e=0.025 ± 0.004), NGC 4833 (7,000 stars, round) and NGC 7089 (4,000 stars, e=0.078 ± 0.006). In Fig. B.1, we plot the distribution of stars as a function of V magnitude and define approximate completeness thresholds for each cluster (in green). We clean the sample by removing stars with V magnitudes above this limit. Then, we run the robust PCA code and compute the ellipticity as described in Sect. 4. The results are presented in Table B.

The table summarises the ellipticity measurements for four representative clusters using three different stellar samples: the full sample (N_all), the magnitude-limited sample after applying a faint-end cut (N_vcut), and the RGB-only sample (N_RGBs). For NGC 6254, NGC 6838, and NGC 7089, the ellipticity values remain broadly consistent across the datasets, with minor variations.

Interestingly, NGC 4833 shows a notably higher ellipticity when considering only RGB stars (e_RGBs=0.14 ± 0.02) compared to both the full sample (e_all=0.045 ± 0.006) and the magnitude-limited sample (e_Vcut=0.037 ± 0.011), as shown in Table B. Given that the RGB sample is well defined and the ellipticity measurement is robust against small-number statistics, this difference likely reflects other factors. This discrepancy could stem from crowding or radial-dependent incompleteness in the central regions, which are not fully accounted for by the faint-end magnitude cut.

Fig. B.1 Number of sources per magnitude bin (log scale) for four clusters: NGC 4833 (upper-left), NGC 6254 (upper-right), NGC 6838 (lower-left), and NGC 7089 (lower-right). In each panel, the dashed green lines represent the faint-end threshold cut.

Appendix C Table of ellipticity results

All Tables

All Figures

	Fig. 1 Recovered ellipticity as a function of the true ellipticity using the KDE method (upper left), the PCA method (upper right), the robust PCA method (lower left), and the robust PCA method after bias correction (lower right). The points are coloured according to the number of stars in the mock GCs.
In the text

Fig. 2 Impact of outliers on the recovered ellipticity and position angle for the KDE (circles), PCA (squares), and robust PCA (diamonds) methods for two sets of true parameters. The colour of the points indicates the location of the outliers. The points in pink with a positive x and negative y correspond approximately to the direction of the semi-minor axis, and the points in orange with a positive x and positive y correspond approximately to the direction of the semi-major axis of the ellipse. Outliers were added within fixed quadrants relative to the cluster centre, not strictly along the ellipse axes, which leads to small asymmetries in the recovered angles when the ellipse position angle changes. The dashed green line represents the true value of the ellipticity and angle, taken to 0.15 and 45° (top panel) and to 0.25 and 60° (bottom panel). We note that for the Robust PCA method, only one point is visible, thus highlighting the robustness of the method. Regardless of the direction where outliers are added, the ellipticity and position angle are accurately recovered.

In the text

	Fig. 3 Ellipticities of the 29 Galactic GCs in our sample measured using the robust PCA method compared with the estimates from the literature provided by White & Shawl (1987, left), Chen & Chen (2010, middle), and Cruz Reyes & Anderson (2024, right). The solid lines indicate equal ellipticities. Each GC is represented by a coloured point (see the legend at the top), and error bars are indicated.
In the text

Fig. 4 Ellipticity (e) of 29 Galactic GCs obtained using the robust PCA method against the relaxation time at the half-mass radius, log th, from Harris (1996). The points are coloured according to the last crossing time through the Galactic disc, log tc, taken from Pancino et al. (2024). The name of each cluster is indicated close to the corresponding point. Both times are given in years. Dynamically young clusters tend to be more elliptical.

In the text

Fig. 5 Rotation strength, taken from Leitinger et al. (2025), as a function of the flattening we measured. The names of GCs are indicated close to the corresponding points and coloured according to the inclination angle of the cluster, also taken from Leitinger et al. (2025). Low inclinations are plotted in blue while high inclinations are coloured in yellow. The uncertainties associated with the inclination angles are not shown in this figure. Instead, we report them in the last column of Table C.1. The dashed lines indicate the relation for isotropic oblate rotators viewed at different inclination angles.

In the text

Fig. 6 Comparison of the position angle of the rotation axis from Leitinger et al. (2025) and the position angle measured in this work. The black line indicates the identity. In this case, this represents an alignment between the semi-minor axis of the ellipse and the rotation axis of the cluster. Here, we have only included the clusters classified as flattened and with an available position angle for the rotation axis. Some clusters are classified as non-rotating in Leitinger et al. (2025), but they still have a measured rotation strength and position angle for the rotation axis. They are plotted with a star symbol.

In the text

	Fig. 7 Ellipticity computed using the robust PCA method with all stars (x-axis) and with all stars within 3rh (y-axis). Only four clusters show ellipticities incompatible within their respective error bars.
In the text

	Fig. A.1 Ellipticity recovered using the robust PCA method before (blue dots) and after (orange dots) correction for the bias due to a small number of stars. The grey area shows the uncertainties in the recovered ellipticity. The dashed line indicates the true ellipticity. The mock GCs were simulated using 300 stars.
In the text

	Fig. B.1 Number of sources per magnitude bin (log scale) for four clusters: NGC 4833 (upper-left), NGC 6254 (upper-right), NGC 6838 (lower-left), and NGC 7089 (lower-right). In each panel, the dashed green lines represent the faint-end threshold cut.
In the text