© The Authors 2026

1. Introduction

Galaxy clusters are powerful cosmological probes, providing competitive constraints on the growth of structure in the late-time Universe (for a recent review, see Clerc & Finoguenov 2023). The utility of clusters as cosmological probes hinges on accurately relating their observable properties–such as X-ray luminosity, gas mass, and optical richness–to their total mass. This is achieved through scaling relations. However, these relations exhibit significant intrinsic scatter, a substantial fraction of which is attributed to the complex merger history and dynamical state of clusters (Kravtsov et al. 2006). Clusters grow hierarchically through mergers and accretion (Frenk et al. 1999), and these events can disrupt the thermal equilibrium of the intracluster medium (ICM), alter the gravitational potential, and induce transient features in the measured ICM properties (Poole et al. 2007; ZuHone et al. 2011).

A critical manifestation of this dynamical activity is the presence of substructure. In the optical regime, substructure reveals itself as infalling groups of galaxies within a cluster’s halo (Dressler & Shectman 1984). In the X-ray regime, it is evidenced by asymmetries, shocks, and cold fronts in the hot gas (Markevitch & Vikhlinin 2002), as well as direct detections of the intragroup medium of infalling groups (Haines et al. 2018). Identifying and characterizing this substructure is therefore not merely a taxonomic exercise but a crucial step in understanding the assembly history of clusters and in controlling the systematic biases that currently limit the precision of cosmological constraints (Pratt et al. 2019). This focus on substructure is motivated by a critical challenge in cluster cosmology: current studies are systematics-limited. As shown by Clerc & Finoguenov (2023), cosmological constraints become inconsistent when using X-ray luminosity versus optical richness, indicating a significant bias. Damsted et al. (2024) identified that a primary driver of this discrepancy is the correlation between optical richness and the large-scale environment, potentially due to variations in cluster formation time affecting the richness–mass relation.

The advent of large, uniform surveys has revolutionized this field. In the optical, algorithms like redMaPPer (Rykoff et al. 2014) applied to the Legacy Surveys (Dey et al. 2019) have cataloged millions of clusters and their member galaxies. In the X-ray, the eROSITA (extended ROentgen Survey with an Imaging Telescope Array) instrument aboard the Spectrum-Roentgen-Gamma (SRG) mission (Predehl et al. 2021) is providing the first all-sky X-ray survey since ROSAT, yielding a vast catalog of clusters with detailed morphological information (Bulbul et al. 2024; Sanders et al. 2025). This multiwavelength dataset enables, for the first time, a statistically robust investigation of substructure across a large population of clusters.

We leveraged these new datasets to perform a comprehensive study of substructure in galaxy clusters. We applied the hierarchical density-based clustering algorithm HDBSCAN (Campello et al. 2013; McInnes et al. 2017) to the redMaPPer galaxy cluster catalog to identify substructure based on the spatial distribution of probabilistically assigned member galaxies. We then cross-matched these optical clusters with their X-ray counterparts in the first eROSITA All-Sky Survey (eRASS1; Merloni et al. 2024) and analyzed their X-ray morphological properties using the largest set of morphological parameters compiled to date (Sanders et al. 2025). Finally, we investigated the impact of substructure on the scaling relation between X-ray luminosity and optical richness.

The paper is structured as follows: In Sect. 2 we describe the substructure analysis of the redMaPPer catalog and present the analysis of the matched eROSITA clusters. In Sect. 3 we present the results of our analysis and discuss their implications. We conclude in Sect. 4. Throughout the paper we assume a flat Λ cold dark matter cosmology with H₀ = 70 km s⁻¹ Mpc⁻¹ and Ω_m = 0.3. Uncertainties are quoted for the 68% confidence level.

2. Data and methods

2.1. redMaPPer catalog

The main sample for this study comprises 4,029,875 optical clusters identified in the DESI Legacy Imaging Surveys (Kluge et al. 2024) at redshifts z < 1 using the redMaPPer algorithm (Rykoff et al. 2014). In this catalog, each galaxy is assigned a probabilistic membership, and a cluster’s total richness (λ) is calculated as

$$\begin{array}{c}\hfill \lambda =\frac{S}{1-f}·{\displaystyle \sum p_{\text{mem}}},\end{array}$$(1)

where S is a scale factor that corrects for the survey-limiting magnitude, f is the masking fraction, and ∑p_mem is the sum of membership probabilities. Since partial cluster coverage can bias substructure detection–a primary focus of this study–we imposed quality cuts of S < 1.3 and f < 0.2 to select only well-surveyed systems. While this reduces the final sample size by 13%, it ensures the reliability of our subsequent substructure analysis, for which completeness in absolute cluster number is less critical than data quality for the retained objects.

2.2. HDBSCAN

We used the hierarchical density-based clustering algorithm HDBSCAN (McInnes et al. 2017), originally developed by Campello et al. (2013), to identify substructure within the redMaPPer cluster catalog. HDBSCAN is run separately for each redMaPPer cluster, and the input data to the algorithm contains only the sky positions of the member galaxies in the projected plane–meaning no color, magnitude, or redshift information is used. The primary parameter affecting the clustering results is the minimum cluster size. We adopted a value of ten galaxies for the minimum subcluster size, a choice informed by the efficiency results of SPIDERS’ (SPectroscopic IDentification of eROSITA Sources) follow-up of redMaPPer members (Clerc et al. 2016) and considerations of the redMaPPer completeness corrections (Rykoff et al. 2014). The results of Clerc et al. (2016) indicate that optical characterization of low-richness systems can be highly unreliable, mainly due to the increased chance of projection effects.

A key feature of HDBSCAN is its treatment of low-density data points (i.e., galaxies in our case) as noise. This effectively cleans our catalog of galaxies in sparsely populated regions that are unlikely to be genuine cluster members. The parameter min_samples determines the conservatism of the clustering, where a higher value results in more points being classified as noise. For our analysis, we set this to the minimum value, min_samples = 1 galaxy. The parameter α, which additionally affects the conservatism of the clustering, was set to 0.05. This value works as a distance scaling parameter; thus, increasing α reduces the influence of local density variations, while decreasing it increases the sensitivity to smaller density fluctuations.

HDBSCAN does not natively facilitate the return of single-component clusters, and configuring it to do so can introduce a bias inherent to the cluster extraction method, known as the Excess of Mass algorithm (Müller & Sawitzki 1991; Campello et al. 2013). The algorithm computes a hierarchy of clusters and selects the most stable ones. Here, the term “mass” refers to a mathematical quantity related to the stability of the clusters during the clustering procedure, and it is not linked to any physical mass of the astrophysical systems. To mitigate this bias, we implemented a two-step analysis procedure:

1. HDBSCAN was first executed with a configuration that prohibits single-component clusters. This results in the identification of either multi-subcluster structures or the classification of all points as noise.

2. For clusters where no substructure is identified in the first step, we reran the algorithm with parameters adjusted to enable the identification of single-component clusters.

Using the same min_cluster_size = 10 for this second run would be overly conservative. Therefore, for these singleton clusters, we used the original redMaPPer richness, λ, as an adaptive minimum cluster size, justified by their prior validation as genuine clusters. All other clustering parameters remain unchanged.

This two-step procedure ensures a consistent cleaning process across both the substructure-rich and no-substructure samples. Following the clustering, the cluster richness is recalculated according to Eq. (1), with any galaxies designated as noise by HDBSCAN assigned a membership probability p_mem = 0. For comparison, we also present results from a control analysis where single-component clusters are identified in a single step from the outset.

2.3. Performance verification

We validated the performance of the HDBSCAN clustering algorithm by applying it to a mock catalog derived from the Uchuu simulations (Comparat et al. 2025). Our mock sample consisted of a cylindrical volume projected along the line of sight, containing 210 clusters and 2524 individual dark matter halos.

To quantify the algorithm’s performance, we calculated two metrics for each simulated cluster with at least ten member galaxies: the fraction of galaxies that are correctly grouped and the fraction that are misplaced. With our chosen parameters, HDBSCAN correctly groups 77.1% of galaxies, on average. Of the remaining 22.9%, 8.5% are incorrectly labeled as noise, while 14.4% are erroneously assigned as members of other clusters.

Our analysis was performed solely in the 2D projected sky plane. The clustering performance could potentially be improved by executing the algorithm in 3D comoving space, an approach employed in spectroscopic studies that combine methods like mclust with friends-of-friends (Tempel et al. 2017). However, this requires spectroscopic redshift information and is therefore reserved for future work.

2.4. X-ray properties

The same catalog of redMaPPer clusters was used to identify extended X-ray sources as galaxy clusters in eRASS1 (Kluge et al. 2024; Merloni et al. 2024), with detailed follow-up studies providing extensive morphological and physical characterizations (Sanders et al. 2025). This established framework enables us to demonstrate the practical utility of our substructure analysis.

The eROSITA telescope, the primary instrument on the SRG mission, is designed to perform an all-sky survey in the X-ray band. A catalog of 12 247 clusters detected by eROSITA is presented in Bulbul et al. (2024). Sanders et al. (2025) conducted a comprehensive analysis of the morphological properties of these clusters, producing the largest uniformly derived set of X-ray morphological parameters to date. This set includes established parameters from previous studies, such as central density (e.g., Ghirardini et al. 2022), cuspiness (inner density slope; Vikhlinin et al. 2007), concentration (e.g., Santos et al. 2008), fit-peak offset, power ratios (Buote & Tsai 1995), the Gini coefficient (Abraham et al. 2003; Lotz et al. 2004), photon asymmetry (Nurgaliev et al. 2013), centroid shift, and ellipticity (Ghirardini et al. 2022). The study also introduces novel parameters: the slosh parameter (H) and the multipole magnitudes M₁–M₄.

In addition to these individual parameters, Sanders et al. (2025) define two composite disturbance scores:

• D_shape: A combined disturbance score based on shape parameters, including ellipticity (ϵ), slosh (H), multipole magnitudes M₁–M₄, and the fit-peak offset (F).

• D_Comb: This score incorporates the same parameters as D_shape but also includes the concentration (c_80–800).

In both cases, a value of 1 indicates a maximally disturbed morphology. All parameters are summarized in Table A.1. As some parameters are sensitive to the definition of the cluster center, values are computed using both the best-fitting centroid and the peak of X-ray emission; the latter are denoted with an asterisk (*).

2.5. Optical and X-ray counterparts

We identified optical and X-ray counterparts by cross-matching the redMaPPer and eROSITA catalogs, considering all pairs within R₅₀₀ and with a maximum redshift difference of 0.02(1 + z). A total of 314 X-ray sources are associated with multiple optical clusters, and 81 optical clusters have multiple associated X-ray sources; these multiple associations will be discussed separately in Sect.3.3. The remaining sample consists of 6421 unique X-ray sources and optical cluster pairs.

As anticipated, not every X-ray cluster has an identified optical counterpart due to differences in sky coverage between the Legacy Surveys and eROSITA. We performed substructure analysis on the successfully matched counterparts using the method outlined in Sect. 2.2.

3. Results and their implications

3.1. Substructure detection in eROSITA clusters

By configuring HDBSCAN to permit single-component clusters in an initial single-step run, we identified 6120 eROSITA clusters without substructure and only 301 clusters exhibiting between 2 and 7 subclusters. This significant imbalance indicates a strong bias toward classifying systems as single-component, likely due to the algorithm’s inherent conservatism when forced to return singleton clusters. The resulting small sample of clusters with substructure is statistically insufficient for meaningful further analysis. Consequently, we excluded the results from this initial, biased clustering from the remainder of this work.

We instead based our analysis on the two-step clustering procedure described previously. This method identified a total of 2159 eROSITA clusters with substructure (containing between 2 and 11 subclusters) and 4262 clusters without substructure.

A number of these clusters were previously identified as merging systems by Wen et al. (2024), who analyzed a sample of 338,841 clusters with spectroscopic redshifts from a parent catalog of 1.58 million clusters (Wen & Han 2024). Their method identified overdensities in the stellar mass distribution for clusters up to z ∼ 1.5, yielding 39 282 partner systems associated with 33 126 main clusters. Within this sample, 184 of our substructure clusters and 200 of our non-substructure clusters are present among their main clusters. Furthermore, 16 of our substructure clusters and 39 without substructure are cataloged as partner systems in their work.

Focusing on massive clusters (λ > 30), Wen et al. (2024) identified 7845 systems with merging subclusters. This subset includes 352 clusters where we detect substructure and 183 where we do not. From the same massive cluster sample, they also identified 3446 post-collision mergers, which encompass 243 of our substructure clusters and 118 without substructure.

Overall, the level of agreement between our substructure detections and the independent results of Wen et al. (2024) is considerable. The observed differences can be partially attributed to fundamental methodological distinctions, particularly in the initial assignment of cluster membership. The redMaPPer algorithm, which underpins our catalog, limits membership assignment to within the estimated virial radius (R_200c). This approach effectively filters out very early-stage mergers whose signatures may lie predominantly outside this boundary. Therefore, our results specifically characterize the role of substructure within the virial radius, as defined by redMaPPer.

As a result of our analysis, we refined the richness estimates for all clusters and introduced a new measurement for the substructure richness, which we designate as the infalling richness, λ_infalling. We defined the host cluster to be the subcluster returned by HDBSCAN with the highest richness, denoted λ_host. We applied the same corrections used for the total cluster richness, λ, as defined in Sect. 2.2.

Numerical simulations indicate that the substructure mass function exhibits a consistent shape when normalized by the host halo mass. This property allowed us to assess the completeness of our substructure detection method. By substituting mass with richness and applying the mass-richness relation by McClintock et al. (2019), we could compare our empirical results against theoretical expectations.

Figure 1 presents the distribution of mass ratios between the infalling galaxy groups and their respective host clusters (μ = M_infalling/M_host) for different bins of total cluster richness. We compared our measurements to the unevolved subhalo mass function of Giocoli et al. (2008):

$$\begin{array}{c}\hfill \frac{\mathit{dN}}{dln(M_{\text{infalling}}/M_{\text{total}})}=N_0{x}^{-\alpha }exp(-6.283x^3),x=\left|\frac{M_{\text{infalling}}}{\alpha M_{\text{total}}}\right|,\end{array}$$(2)

Fig. 1. Distribution of mass ratios, μ = Minfalling/Mhost, between the infalling groups and their host clusters in different bins of total cluster richness. The dashed purple and yellow lines show the unevolved subhalo mass functions of Giocoli et al. (2008) and Jiang & van den Bosch (2014), respectively.

where M_infalling is the mass of the infalling subcluster, M_host is the host mass, M_total is the total mass of the descendant cluster, α = 0.8 and N₀ = 0.21. Here, we assume that M_total = M_host + M_infalling. Additionally, we plotted the subhalo mass function of Jiang & van den Bosch (2014), which adds an extra power-law term to the fitting function of Giocoli et al. (2008).

The form and steepness of the mass-ratio distribution of detected infalling subclusters reproduce the unevolved subhalo mass function of Jiang & van den Bosch (2014) reasonably well, especially at high richness and high μ values. At low values of μ, the distribution reflects incompleteness caused by a finite detection threshold. For clusters with λ ≥ 20, our substructure detections are complete down to μ > 0.4, while for λ ≥ 40, completeness extends to μ > 0.25. At the lowest richnesses, only nearly equal-mass mergers (μ > 0.8) are reliably detected.

The redshift-richness distribution of our sample, shown in Fig. 2, indicates that incompleteness is primarily a function of low richness, with no strong redshift dependence. Since our substructure measurements are highly incomplete for λ < 20, we limited the subsequent analysis to λ ≥ 20. We introduced two bins of total cluster richness and additionally created subsamples where we impose a lower limit of μ > 0.35, at which the fraction of clusters exhibiting substructure is approximately 25% in both richness bins. For completeness, we first present the results without the μ cut and then with the cut applied. Sample sizes for each subset are provided in Fig. 3. Single-component clusters are more prevalent at lower richness, while the occurrence of substructure increases with richness due to the finite substructure detection threshold. This trend underscores the importance of the detection threshold: although previous studies report that 40 − 70% of clusters show substructure (Schuecker et al. 2001; Smith et al. 2005; Wen & Han 2013), our results demonstrate that this fraction is highly sensitive to the limits imposed on substructure richness and its fractional contribution to the total cluster richness.

Fig. 2. Redshift vs. richness for the full sample. Clusters with substructure are shown in black, and clusters without in red.

Fig. 3. Number of substructure and no-substructure clusters for the full λ ≥ 20 sample and the two richness bins.

3.2. Influence of substructure on X-ray morphology

We compared the distributions of morphological parameters between clusters with and without identified substructure using the Anderson-Darling (AD) test. Additionally, we used the Kolmogorov-Smirnov (KS) test and a permutation test (shuffle test) with the mean and median as test statistics. The resulting probability values (p-values) indicate the likelihood that the two samples are drawn from identical distributions. Since we tested multiple hypotheses on the same dataset, we had an increased likelihood of false positives. This effect can be addressed using a Bonferroni correction; however, the results would be overly conservative since many of the morphological parameters are correlated with each other. We were careful in drawing any conclusions from the results.

In the full λ ≥ 20 sample, without any richness binning or cuts in μ, our analysis reveals statistically significant differences (p < 0.05) in all morphological parameters between the two groups except slosh (H) and multipole magnitudes M₁, M₃, and M₄. The combined disturbance scores, D_comb and D_shape, show that clusters with substructure are generally more disturbed. Their cumulative distribution functions indicate that a larger fraction of these clusters exhibit higher disturbance values at any given threshold. The AD test yields p-values of 1.00  ×  10⁻⁴ and 0.049 for D_comb and D_shape, respectively. The KS test yields a p-value of 5.97 × 10⁻¹² for D_comb but fails to find a significant difference between the D_shape distributions, indicating a much stronger distinction in the D_comb parameter.

This pronounced difference in D_comb is primarily driven by large disparities in the concentration parameter c_{80 − 800}, where clusters with substructure display significantly flatter surface brightness profiles. Concentration is defined as

$$\begin{array}{c}\hfill c_{80-800}=log\frac{I_X(80\text{kpc})}{I_X(800\text{kpc})},\end{array}$$(3)

where I_X(r) is the integrated surface brightness within radius r. A flat profile yields a value near −2, while a steep, peaked profile approaches 0. In order to account for the effects of the point spread function and background, Sanders et al. (2025) compute the concentration using their cluster model, rather than the observed surface brightness. Cluster models are computed with the MultiBand Projector in 2D (MBProj2D) tool (Bulbul et al. 2024). Sanders et al. (2025) use monochromatic 0.2–2.3 keV band images as the input to MBProj2D. The concentration is then evaluated from the model before point spread function convolution and background addition. An AD test on the c_{80 − 800} distributions returns a p-value of 1.00 × 10⁻⁴, a KS yields a p-value of 3.21 × 10⁻¹⁷, and the permutation test gives p-values of 2.0 × 10⁻⁴ for both the mean and median, robustly confirming the significance of this result. The distribution of c_{80 − 800} is shown in Fig. 4.

Fig. 4. Cumulative distributions of the concentration c80 − 800 for clusters with λ ≥ 20.

The permutation test does not find significant differences in the means of the photon asymmetry, A_phot and A_phot^*, distributions, while the median shows a clear difference with p-values of 1.00 × 10⁻⁴ for both distributions. This is caused by the long left tail of the no-substructure sample due to a couple of extreme outliers. Clusters 1eRASS J105803.8-041038 and 1eRASS J030104.7-554410 have A_phot values of −2935.484 and −665.066, while the median is ∼0.71. These outliers increase the variance of the distribution and reduce the robustness of the mean as a test statistic. This showcases the importance of using multiple statistical tests when comparing the distributions. Next, we extended our comparison of morphological parameter distributions to different richness ranges.

In the high-richness bin (λ > 40), AD tests yield p-values of 9.00 × 10⁻⁴ for D_comb and 0.016 for D_shape, indicating significant differences between clusters with and without substructure. Beyond the combined disturbance scores, we find significant differences in several other parameters as described in the following paragraph. These include central densities, inner density slope, concentration, fit-peak offset, centroid shift, and the multipole magnitude M₄.

As expected for disturbed systems, clusters with substructure also exhibit flatter generalized density profiles, as indicated by the inner density slope parameters α and α₅₀, and flatter surface brightness profiles, as measured by the concentration parameter c_{80 − 800}. In contrast, we find no significant difference in the c₅₀₀ distributions, which are computed within a fraction of R₅₀₀ rather than a fixed physical aperture. Further supporting their disturbed nature, clusters with substructure show higher centroid shift (ω) values, indicating greater spatial variability in the centroid of emission, and higher M₄ multipole magnitudes, suggesting more pronounced large-scale asymmetries. No significant differences are detected in the distributions of several other morphological parameters, including the power ratios P₁₀–P₃₀, Gini coefficient, photon asymmetry, ellipticity (ϵ), slosh parameter (H), and multipole magnitudes (M₁–M₃). This suggests that the presence of substructure is most evident in a specific subset of morphological indicators.

In the low-richness bin (20 < λ < 40), differences in shape-based morphological parameters are less pronounced. We find no significant differences in Gini coefficient, ellipticity, slosh, multipole magnitudes M₁–M₄, or D_shape. The only shape parameter showing a significant difference is the fit-peak offset (AD p-value = 8.00 × 10⁻³), with substructure clusters exhibiting larger offsets between the best-fitting center and the X-ray emission peak. Significant differences are also found in non-shape parameters: central densities (lower at fixed R₅₀₀ fraction but higher at fixed 50 kpc radius), inner density slope (flatter for substructure clusters), concentration (flatter profiles in c_{80 − 800}, c_{80 − 800}^*, c₅₀₀, and c₅₀₀^*), power ratios P₁₀–P₄₀, photon asymmetry, centroid shift, and D_comb.

A surprising result in this richness bin is that power ratios P₁₀–P₄₀ suggest clusters without substructure are less spherical than those with substructure. We note that the difference is minimal, and if we were to use a Bonferroni correction it would not be statistically significant. Therefore, we refrained from drawing any conclusions from this seemingly puzzling result.

The p-values for D_shape from AD, KS, and shuffle tests across richness bins are plotted in Fig. 5. Additionally, the D_shape and D_comb distributions for the two richness bins are shown in Figs. 6 and 7. These figures clearly show how the differences in the shape parameters are less pronounced in the low-richness bin. Applying a substructure richness ratio cut (μ > 0.35) does not alter our conclusions regarding disturbance parameters, and we find no significant differences between the full substructure sample and the μ > 0.35 subset, confirming that incompleteness at low μ does not bias our results. Figures 6 and 7 also include the combined disturbance score distributions for the μ > 0.35 subsamples. These figures demonstrate that the μ > 0.35 cut has virtually no effect as the difference is barely noticeable in the plots.

Fig. 5. p-values for the comparison of DShape distributions for the full λ ≥ 20 sample and the two richness bins. Results from the AD tests are shown in red, the KS tests in cyan, and the permutation tests in purple (mean) and black (median). The horizontal dashed and dotted black lines indicate the 0.05 and 0.01 significance levels, respectively.

Fig. 6. Cumulative distributions of the combined disturbance score Dshape for the substructure (black), substructure μ > 0.35 (blue), and no-substructure (red) groups for two richness bins.

Fig. 7. Same as Fig. 6 but for Dcomb.

All p-values for the statistical tests on the full λ ≥ 20 sample and the two richness bins are summarized in Tables A.2, A.3, and A.4. Overall, clusters with substructure are consistently more disturbed across all richness levels.

3.3. Multiple associated sources

Of the 92 optical clusters with multiple associated X-ray sources, we identify substructure in 48 clusters. Based on visual inspection, the substructure identified by HDBSCAN aligns well with the positions of the eROSITA clusters in 37 of these 48 cases. In the remaining 11 cases, the optically identified substructure does not correspond to the locations of the X-ray sources. Among the 37 clusters where the splitting is successful, we find that 19 contain more subcomponents than the number of associated X-ray sources. We detect no substructure in 44 optical clusters that have multiple X-ray counterparts.

For the 314 X-ray clusters linked to multiple optical clusters, there are 632 optical cluster associations in total. Among these, we find substructure in 315 optical clusters and no substructure in the remaining 317.

3.4. Scaling relations

To investigate the impact of substructure on scaling relations, we analyzed the correlation between optical richness and X-ray luminosity. We used the X-ray luminosity in the 0.2–2.3 keV rest-frame band, expressed as ln(L_xE(z)⁻¹ erg⁻¹ s), and the richness as ln(λ/E_z), where E_z² = Ω_m(1 + z)³ + Ω_Λ accounts for the evolution of the scaling relations. The richness was normalized by its median value of 27.9.

We performed parameter estimation using the LINMIX¹ routine from Kelly (2007), which implements a Bayesian linear regression model that accounts for measurement uncertainties in both variables. The linear relation is given by

$$\begin{array}{c}\hfill \eta =\alpha +\beta \xi +ϵ,\end{array}$$(4)

where α is the intercept, β is the slope, and ϵ represents the intrinsic scatter, assumed to be normally distributed with mean zero and variance σ². The actual values (η, ξ) are related to the observed values by

$$\begin{array}{cc}\hfill x_i& =\xi +x_{\text{err}}\hfill \\ \hfill y_i& =\eta +y_{\text{err}},\hfill \end{array}$$(5)

where x_err and y_err are the errors in the data x_i and y_i. We ran 100 000 Markov chain Monte Carlo iterations, discarding the first 1000 steps as burn-in. Reported values correspond to the medians of the posterior distributions.

Results from the LINMIX fits are summarized in Tables 1 and 2. We fit the relation for the full sample, as well as separately for clusters with and without substructure.

Our sample covers a wide range in richness and redshift. To ensure high completeness and minimize contamination, we limited our analysis to the region where the eRASS1 sensitivity is high, characterized by a redshift-dependent richness completeness limit λ > z × 100, analogous to previous ROSAT All-Sky Survey studies (Klein et al. 2019; Finoguenov et al. 2020). Consequently, low-richness clusters are only well sampled at low redshifts. To disentangle redshift and richness trends, we analyzed the scaling relations in bins of both richness and redshift.

Given that substructure measurements are highly incomplete for λ < 20, we present results only for λ ≥ 20. The relations for these clusters are shown in Figs. 8 and 9. The scatter in the relation is comparable between clusters with and without substructure. We note some differences in the slope and normalization; however, the values compensate for each other, making the relation comparable between the samples. We check for further redshift and richness trends below. Since Damsted et al. (2023) report large scatter in scaling relations at z < 0.2, we first tested how substructure contributes to these results.

Fig. 8. X-ray luminosity (LX) vs. richness (λ) for clusters with λ ≥ 20. Clusters with substructure are marked with black circles, and clusters without as red circles. The solid black (solid cyan) line represents the best-fit scaling relation for clusters with (without) substructure. The dashed blue line shows the relation fitted to the full sample (λ ≥ 20). All fits were obtained using the LINMIX Bayesian regression method, which accounts for measurement errors in both variables.

Fig. 9. Effect of substructure on the X-ray luminosity–richness scaling relation. Contours indicate the 68% and 95% confidence regions of the posterior distributions for the scaling relation parameters. Clusters with substructure are represented in black, and clusters without in red. The blue contours and line correspond to the full sample. The parameters α, β, and σ denote the intercept, slope, and intrinsic scatter of the relation ln(LX) = α + βln(λ)+ϵ, where ϵ ∼ 𝒩(0, σ2).

Figure 10 shows redshift-dependent trends for two richness bins, separately for clusters with and without substructure. Clusters without substructure exhibit very mild evolution in scatter with redshift, while those with substructure show strong evolution in both richness bins. We also observe a higher normalization in the scaling relations for clusters with substructure, with the difference being comparable to the scatter in that subsample. Notably, this effect is only pronounced at low redshifts. If interpreted in the context of the growing importance of cool cores at low redshifts (McDonald et al. 2013), the excess luminosity in merging clusters could be naturally explained, along with the dominance of cool cores at low z. One may further speculate whether survival times for low-entropy components increase at low redshifts, possibly due to stronger magnetic fields, which can be further verified by corresponding predictions for redshift evolution of the nonthermal radio emission from galaxy clusters (see van Weeren et al. 2019 for a review).

Fig. 10. Effect of substructure on the redshift and the richness-dependent trends on the X-ray luminosity vs. richness scaling relation. Left panel: Clusters with substructure. Right panel: Clusters without substructure. Blue shows z < 0.2 and red z ≥ 0.2. Solid lines show λ ≥ 40 samples, and dashed lines 20 ≤ λ < 40 samples.

The richness distributions–including original redMaPPer values and our updated estimates (as described in Sect. 2.2) for clusters with and without substructure–are shown in Fig. 11. Clusters with substructure are slightly cleaner after HDBSCAN processing than those without, which partly contributes to the higher normalization but does not fully explain it, especially the redshift dependence observed in Fig. 10. Fitting the relation using the original redMaPPer richness, we find that using our updated HDBSCAN-based richness reduces the scatter by 7% for clusters with λ ≥ 20.

Fig. 11. Comparison of the redMaPPer and updated λ distributions for clusters with and without substructure.

Additionally, we find that scatter is a strong function of richness. This is evident in Fig. 12, which shows no-substructure clusters at z < 0.2 in three richness bins, including the range 10 ≤ λ < 20 (which is well sampled at low z). A clear decrease in scatter with increasing richness is observed. The other parameters show strong correlation and the values are consistent within 2σ. We have marginal evidence that the relation is steeper for low-mass clusters.

Fig. 12. Effect of richness on the X-ray luminosity vs. richness scaling relation. Blue shows the λ ≥ 40 clusters, black 20 ≤ λ < 40 clusters, and red 10 ≤ λ < 20 clusters. All clusters are without substructure.

Our results are not affected by the limitation λ > z × 100. We can observe the same redshift and richness-dependent trends in Fig. 13 where the limitation was not applied.

Fig. 13. Effect of substructure on the redshift and richness-dependent trends on the Lx − λ scaling relation without the λ > z × 100 limitation. Left panel: Clusters with substructure. Right panel: Clusters without substructure. Red shows z ≥ 0.3, black 0.15 ≤ z < 0.3, and blue z < 0.15. Solid lines show λ ≥ 40 samples, and dashed lines 20 ≤ λ < 40 samples.

We also fit the scaling relation to the subsamples of multiple associated sources described in Sect. 3.3. For the subsample of multiple X-ray sources (i.e., optical clusters matched to more than one eROSITA detection), we included only those systems where the substructure identified by HDBSCAN aligns well with the multiple X-ray sources and the λ > z × 100 condition is fulfilled. This left us with 41 subclusters. Of these, 34 are at z < 0.2. The mean subcluster richness of this sample is $\overline{\lambda }=26.3$. Here, the subcluster richness is not limited to λ ≥ 20. We find that these clusters exhibit a relation comparable to that of the general substructure sample with similar richness and redshift.

For the subsample of multiple optical clusters (i.e., X-ray sources matched to more than one optical cluster), we summed the richness values (λ) of all associated optical clusters within a single X-ray source to form a total richness estimate. This combined richness is used in the scaling relation fit. The mean richness and redshift of the subsample are $\overline{\lambda }=71.4$ and $\overline{z}=0.31$. We find that this subsample exhibits a relation comparable to the general substructure sample of massive clusters (λ ≥ 40), with slightly increased scatter for a given redshift.

4. Conclusions

Based on our comprehensive analysis of substructure in galaxy clusters using HDBSCAN on the redMaPPer catalog and its correlation with X-ray morphology and scaling relations, we draw the following conclusions:

• Our two-step HDBSCAN procedure successfully identifies substructure in a significant fraction of clusters, revealing 2159 clusters with substructure (containing between 2 and 11 subclusters) and 4262 without. The method effectively cleans galaxy membership assignments, reducing scatter in richness estimates. The mass ratio μ = M_infalling/M_host is a crucial parameter, with completeness limits of μ > 0.4 for λ ≥ 20 and μ > 0.25 for λ ≥ 40.

• Our substructure detections show considerable agreement with independent catalogs of merging clusters (Wen et al. 2024), validating our methodology. Differences arise primarily from the initial cluster membership definition in redMaPPer, which is limited to the virial radius (R_200c) and may filter out very early-stage mergers. Thus, our results characterize substructure within the virialized region of clusters.

• Clusters with substructure are unequivocally more disturbed. They exhibit significantly different distributions in key morphological parameters, particularly flatter surface brightness profiles (lower concentration c_{80 − 800}), higher centroid shifts (ω), and larger multipole moments (M₄). The combined disturbance score (D_comb), which incorporates concentration, shows a far stronger distinction than the shape-based D_shape score. This confirms that substructure is a primary driver of morphological disturbance in clusters.

The presence of substructure significantly affects the L_X − λ scaling relation. We find:

• Clusters with substructure exhibit a higher normalization in the L_X − λ relation at low redshifts (z < 0.2). This is consistent with the survival of cool cores from infalling subclusters, which boost the X-ray luminosity, and aligns with the known increased prevalence of cool cores at low redshifts.

• Substructure introduces greater intrinsic scatter into the scaling relation at low redshifts. The survival of cool cores from infalling subclusters produces highly luminous systems, where at times multiple cool cores can even be present. Occasionally, the cool cores are destroyed in the mergers, which leads to under-luminous systems. The presence of both over- and under-luminous systems increases the scatter.

• The effect of substructure on scaling relations is strongly redshift-dependent. The increased scatter and normalization are dominant at low redshifts, while at higher redshifts the differences between clusters with and without substructure become less distinct, likely due to the declining prevalence of stable cool cores.

• The ability to detect substructure and its impact is a strong function of both cluster richness and redshift. Substructure is more common in high-richness systems. Furthermore, morphological parameters themselves correlate with richness; lower-richness clusters are generally less spherical and symmetrical regardless of substructure.

In summary, internal substructure is a fundamental property of galaxy clusters that must be accounted for in cosmological studies. It is a major source of scatter in scaling relations and a significant contributor to morphological disturbance, particularly at low redshifts. Correcting for its effects, through methods like the one presented here, is essential for achieving the precision required in next-generation cluster cosmology. Future efforts incorporating spectroscopic data for 3D clustering and larger samples from upcoming surveys will further refine our understanding of substructure and its role in cluster evolution.

Data availability

The catalogs of the identified redMaPPer substructure are available at the CDS via https://cdsarc.cds.unistra.fr/viz-bin/cat/J/A+A/710/A57.

Acknowledgments

We thank the referee for the insightful comments, which helped us to improve this paper. The work is based on the Legacy Survey data. The Legacy Surveys consist of three individual and complementary projects: the Dark Energy Camera Legacy Survey (DECaLS; Proposal ID #2014B-0404; PIs: David Schlegel and Arjun Dey), the Beijing-Arizona Sky Survey (BASS; NOAO Prop. ID #2015A-0801; PIs: Zhou Xu and Xiaohui Fan), and the Mayall z-band Legacy Survey (MzLS; Prop. ID #2016A-0453; PI: Arjun Dey). DECaLS, BASS and MzLS together include data obtained, respectively, at the Blanco telescope, Cerro Tololo Inter-American Observatory, NSF’s NOIRLab; the Bok telescope, Steward Observatory, University of Arizona; and the Mayall telescope, Kitt Peak National Observatory, NOIRLab. Pipeline processing and analyses of the data were supported by NOIRLab and the Lawrence Berkeley National Laboratory (LBNL). The Legacy Surveys project is honored to be permitted to conduct astronomical research on Iolkam Du’ag (Kitt Peak), a mountain with particular significance to the Tohono O’odham Nation. NOIRLab is operated by the Association of Universities for Research in Astronomy (AURA) under a cooperative agreement with the National Science Foundation. LBNL is managed by the Regents of the University of California under contract to the U.S. Department of Energy. This project used data obtained with the Dark Energy Camera (DECam), which was constructed by the Dark Energy Survey (DES) collaboration. Funding for the DES Projects has been provided by the U.S. Department of Energy, the U.S. National Science Foundation, the Ministry of Science and Education of Spain, the Science and Technology Facilities Council of the United Kingdom, the Higher Education Funding Council for England, the National Center for Supercomputing Applications at the University of Illinois at Urbana-Champaign, the Kavli Institute of Cosmological Physics at the University of Chicago, Center for Cosmology and Astro-Particle Physics at the Ohio State University, the Mitchell Institute for Fundamental Physics and Astronomy at Texas A&M University, Financiadora de Estudos e Projetos, Fundacao Carlos Chagas Filho de Amparo, Financiadora de Estudos e Projetos, Fundacao Carlos Chagas Filho de Amparo a Pesquisa do Estado do Rio de Janeiro, Conselho Nacional de Desenvolvimento Cientifico e Tecnologico and the Ministerio da Ciencia, Tecnologia e Inovacao, the Deutsche Forschungsgemeinschaft and the Collaborating Institutions in the Dark Energy Survey. The Collaborating Institutions are Argonne National Laboratory, the University of California at Santa Cruz, the University of Cambridge, Centro de Investigaciones Energeticas, Medioambientales y Tecnologicas-Madrid, the University of Chicago, University College London, the DES-Brazil Consortium, the University of Edinburgh, the Eidgenossische Technische Hochschule (ETH) Zurich, Fermi National Accelerator Laboratory, the University of Illinois at Urbana-Champaign, the Institut de Ciencies de l’Espai (IEEC/CSIC), the Institut de Fisica d’Altes Energies, Lawrence Berkeley National Laboratory, the Ludwig Maximilians Universitat Munchen and the associated Excellence Cluster Universe, the University of Michigan, NSF’s NOIRLab, the University of Nottingham, the Ohio State University, the University of Pennsylvania, the University of Portsmouth, SLAC National Accelerator Laboratory, Stanford University, the University of Sussex, and Texas A&M University. BASS is a key project of the Telescope Access Program (TAP), which has been funded by the National Astronomical Observatories of China, the Chinese Academy of Sciences (the Strategic Priority Research Program “The Emergence of Cosmological Structures” Grant # XDB09000000), and the Special Fund for Astronomy from the Ministry of Finance. The BASS is also supported by the External Cooperation Program of Chinese Academy of Sciences (Grant # 114A11KYSB20160057), and Chinese National Natural Science Foundation (Grant # 12120101003, # 11433005). The Legacy Survey team makes use of data products from the Near-Earth Object Wide-field Infrared Survey Explorer (NEOWISE), which is a project of the Jet Propulsion Laboratory/California Institute of Technology. NEOWISE is funded by the National Aeronautics and Space Administration. The Legacy Surveys imaging of the DESI footprint is supported by the Director, Office of Science, Office of High Energy Physics of the U.S. Department of Energy under Contract No. DE-AC02-05CH1123, by the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility under the same contract; and by the U.S. National Science Foundation, Division of Astronomical Sciences under Contract No. AST-0950945 to NOAO. This work is based on data from eROSITA, the soft X-ray instrument aboard SRG, a joint Russian-German science mission supported by the Russian Space Agency (Roskosmos), in the interests of the Russian Academy of Sciences represented by its Space Research Institute (IKI), and the Deutsches Zentrum für Luft- und Raumfahrt (DLR). The SRG spacecraft was built by Lavochkin Association (NPOL) and its subcontractors, and is operated by NPOL with support from the Max Planck Institute for Extraterrestrial Physics (MPE). The development and construction of the eROSITA X-ray instrument was led by MPE, with contributions from the Dr. Karl Remeis Observatory Bamberg & ECAP (FAU Erlangen-Nuernberg), the University of Hamburg Observatory, the Leibniz Institute for Astrophysics Potsdam (AIP), and the Institute for Astronomy and Astrophysics of the University of Tübingen, with the support of DLR and the Max Planck Society. The Argelander Institute for Astronomy of the University of Bonn and the Ludwig Maximilians Universität München also participated in the science preparation for eROSITA.

References

Appendix A: Tables

Table A.1 summarizes all of the morphological parameters from Sanders et al. (2025). Tables A.2-A.4 contain all of the p-values from the performed statistical tests. The statistical tests were performed using the SciPy (Virtanen et al. 2020) modules for Python. Due to the finite number of permutations, the p-values provided by the AD and permutation tests are capped at 1.0 × 10⁻⁴, while the KS-test has no such limitation. Figure A.1 shows the p-values from the AD and KS tests plotted against each other. Here entries with KS p-values lower than the AD-test cap of 1.0 × 10⁻⁴ are left out in order to make a reasonable comparison. We find that the tests are generally in good agreement.

Fig. A.1. p-values from the AD and KS tests plotted against each other.

All Tables

All Figures

	Fig. 1. Distribution of mass ratios, μ = Minfalling/Mhost, between the infalling groups and their host clusters in different bins of total cluster richness. The dashed purple and yellow lines show the unevolved subhalo mass functions of Giocoli et al. (2008) and Jiang & van den Bosch (2014), respectively.
In the text

	Fig. 2. Redshift vs. richness for the full sample. Clusters with substructure are shown in black, and clusters without in red.
In the text

	Fig. 3. Number of substructure and no-substructure clusters for the full λ ≥ 20 sample and the two richness bins.
In the text

	Fig. 4. Cumulative distributions of the concentration c80 − 800 for clusters with λ ≥ 20.
In the text

	Fig. 5. p-values for the comparison of DShape distributions for the full λ ≥ 20 sample and the two richness bins. Results from the AD tests are shown in red, the KS tests in cyan, and the permutation tests in purple (mean) and black (median). The horizontal dashed and dotted black lines indicate the 0.05 and 0.01 significance levels, respectively.
In the text

	Fig. 6. Cumulative distributions of the combined disturbance score Dshape for the substructure (black), substructure μ > 0.35 (blue), and no-substructure (red) groups for two richness bins.
In the text

	Fig. 7. Same as Fig. 6 but for Dcomb.
In the text

Fig. 8. X-ray luminosity (LX) vs. richness (λ) for clusters with λ ≥ 20. Clusters with substructure are marked with black circles, and clusters without as red circles. The solid black (solid cyan) line represents the best-fit scaling relation for clusters with (without) substructure. The dashed blue line shows the relation fitted to the full sample (λ ≥ 20). All fits were obtained using the LINMIX Bayesian regression method, which accounts for measurement errors in both variables.

In the text

Fig. 9. Effect of substructure on the X-ray luminosity–richness scaling relation. Contours indicate the 68% and 95% confidence regions of the posterior distributions for the scaling relation parameters. Clusters with substructure are represented in black, and clusters without in red. The blue contours and line correspond to the full sample. The parameters α, β, and σ denote the intercept, slope, and intrinsic scatter of the relation ln(LX) = α + βln(λ)+ϵ, where ϵ ∼ 𝒩(0, σ2).

In the text

	Fig. 10. Effect of substructure on the redshift and the richness-dependent trends on the X-ray luminosity vs. richness scaling relation. Left panel: Clusters with substructure. Right panel: Clusters without substructure. Blue shows z < 0.2 and red z ≥ 0.2. Solid lines show λ ≥ 40 samples, and dashed lines 20 ≤ λ < 40 samples.
In the text

	Fig. 11. Comparison of the redMaPPer and updated λ distributions for clusters with and without substructure.
In the text

	Fig. 12. Effect of richness on the X-ray luminosity vs. richness scaling relation. Blue shows the λ ≥ 40 clusters, black 20 ≤ λ < 40 clusters, and red 10 ≤ λ < 20 clusters. All clusters are without substructure.
In the text

	Fig. 13. Effect of substructure on the redshift and richness-dependent trends on the Lx − λ scaling relation without the λ > z × 100 limitation. Left panel: Clusters with substructure. Right panel: Clusters without substructure. Red shows z ≥ 0.3, black 0.15 ≤ z < 0.3, and blue z < 0.15. Solid lines show λ ≥ 40 samples, and dashed lines 20 ≤ λ < 40 samples.
In the text

	Fig. A.1. p-values from the AD and KS tests plotted against each other.
In the text