© ESO 2021

1. Introduction

The isotropy of the late Universe has been a question of great debate in recent decades and a conclusive answer has yet to be given. As precision cosmology enters a new era with numerous experiments covering the full electromagnetic spectrum, the underlying assumption of isotropy for widely adopted cosmological models has to be scrutinized as well. The importance of new, independent tests of high precision cannot be understated. A possible departure of isotropy in the local Universe could have major implications for nearly all aspects of extragalactic astronomy.

Galaxy clusters, the largest gravitationally bound objects in the Universe, can be of great service for this purpose. As a result of the multiple physical processes taking place within them, different components of clusters can be observed almost throughout the full electromagnetic spectrum (e.g., Allen et al. 2011). This provides us with several possible cosmological applications for these objects.

Such a possible application, for instance, comes from the so-called scaling relations of galaxy clusters (e.g., Kaiser 1986). These are simply the correlations between the many cluster properties and can usually be described by simple power-law forms. Some of the measured cluster properties depend on the assumed values of the cosmological parameters (e.g., X-ray luminosity), while others do not (e.g., temperature). Utilizing scaling relations between properties of these two categories can provide us with valuable insights into different aspects of cosmology.

More specifically, the cosmic isotropy can be investigated using such methods. In Migkas & Reiprich (2018) and Migkas et al. (2020) (hereafter M18 and M20, respectively), we performed such a test with very intriguing results. We studied the isotropy of the X-ray luminosity–temperature (L_X − T) relation, which we used as a potential tracer for the isotropy of the expansion of the local Universe. In M20 we combined the extremely expanded HIghest X-ray FLUx Galaxy Cluster Sample (eeHIFLUGCS, Reiprich 2017; Pacaud et al., in prep.) with other independent samples and detected a ∼4.5σ anisotropy toward the Galactic coordinates (l, b)∼(305°, −20°). The fact that several other studies using different cosmological probes and independent methods find anisotropies toward similar sky patches makes these findings even more interesting.

Multiple possible systematics were tested separately as potential explanations for the apparent anisotropies, but the tension could not be sufficiently explained by any such test. Therefore, the anisotropy of the L_X − T relation seems to be attributed to an underlying, physical reason. There are three predominant phenomena that could create this: unaccounted X-ray absorption, bulk flows (BFs), and Hubble expansion anisotropies. Firstly, the existence of yet undiscovered excess X-ray absorption effects could bias our estimates. The performed tests in M20 however showed that this is very unlikely to explain the apparent anisotropies. Further investigation is nonetheless needed to acquire a better understanding of these possibilities.

Secondly, coherent motions of galaxies and galaxy clusters over large scales, BFs, could also be the cause of the observed cluster anisotropies. The objects within a BF have a peculiar velocity component toward a similar direction due to the gravitational attraction of a larger mass concentration such as a supercluster. These nonrandom peculiar velocities are imprinted in the observed redshifts. If not taken into account, they can result in a biased estimation of the redshift-based distances of the clusters and eventually their other properties (e.g., L_X). The necessary BF amplitude and scale to wash away the observed cluster anisotropies by far surpasses ΛCDM expectations, which predicts that such motions should not be present at comoving scales of ≳200 Mpc (see references below). If such a motion is confirmed, a major revision of the large-scale structure formation models might be needed. These BFs have been extensively studied in the past using various methods (e.g., Lauer & Postman 1994; Hudson et al. 2004; Kashlinsky et al. 2008, 2010; Colin et al. 2011; Osborne et al. 2011; Feindt et al. 2013; Appleby et al. 2015; Carrick et al. 2015; Hoffman et al. 2015; Scrimgeour et al. 2016; Watkins & Feldman 2015; Peery et al. 2018; Qin et al. 2019). However, most of these studies used galaxy samples, which suffer from the limited scale out to which they can be constructed. No past study has used cluster scaling relations to investigate possible BF signals to our knowledge. A more in-depth analysis to determine if such effects are the origin of the observed cluster anisotropies is thus necessary.

The vast majority of cluster studies ignore the effects of BFs in the observed redshifts and assume the peculiar velocities are randomly distributed. The frequent use of heliocentric redshifts in local cluster studies, instead of cosmic microwave background (CMB) frame redshifts, might also amplify the introduced bias from BFs. Hence, such discovered motions could strongly distort the results for most cluster studies along with their cosmological applications.

The third possible explanation for our results is an anisotropy in the Hubble expansion and the redshift–distance conversion. To explain the observed anisotropies obtained in M20 solely by a spatial variation of H₀, we would need a ∼10% variation between (l, b)∼(305°, −20°) and the rest of the sky. This scenario would contradict the assumption of the cosmological principle and the isotropy of the Universe, which has a prominent role in the standard cosmological model. This result was derived in M20 based on relatively low-z data, where the median redshift of all the 842 used clusters is z ∼ 0.17. Hence, for now, it is not possible to distinguish between a possibly primordial anisotropy or a relatively local anistropy purely from galaxy clusters. New physics that interfere with the directionality of the expansion rate would be needed in that case, assuming of course that no other underlying issues cause the cluster anisotropies. Many recent studies have tackled this question with contradictory results (e.g., Bolejko et al. 2016; Bengaly et al. 2018; Colin et al. 2019; Soltis et al. 2019; Andrade et al. 2019; Hu et al. 2020; Fosalba & Gaztañaga 2021; Salehi et al. 2020; Secrest et al. 2021).

The interplay of the different possible phenomena and the effects they may have on cluster measurements can make it hard to identify the exact origin of the anisotropies. Additionally, a combination of more than one phenomena might affect the cluster measurements because the existence of a dark gas cloud, for instance, does not exclude the simultaneous existence of a large BF. In order to provide a conclusive answer to the problem, other tests with the same cluster samples must be utilized. These new tests should have somewhat different sensitivities than the L_X − T relation. This way, we can cross check if the anisotropies also appear in tests sensitive only to X-ray effects, BFs, or cosmological anisotropies, for instance.

In this work, alongside L_X and T, we also measure and use the total integrated Compton parameter Y_SZ, the half-light radii of the clusters R, and the infrared luminosity L_BCG of the brightest galaxy of each cluster (BCGs). This allows us to study ten different scaling relations between these properties and test their directional dependence. Based on these results, we try to fully investigate the nature of the observed anisotropies and provide conclusive results. Throughout this paper, unless otherwise stated, we use a ΛCDM cosmology with H₀ = 70 km s⁻¹ Mpc⁻¹, Ω_m = 0.3 and Ω_Λ = 0.7 to constrain the scaling relation parameters.

The paper is organized as follows: In Sect. 2 we describe the measurements of the used samples. In Sect. 3 we describe the modeling of the scaling relations and the statistical methods and procedures used to constrain their directional behavior. In Sect. 4, the full sky best-fit parameters for the ten scaling relations are presented. In Sect. 5, we study the anisotropic behavior of scaling relations, which are sensitive only to X-ray absorption effects and not cosmological factors. In Sect. 6, the focal point of this work is presented, namely, the anisotropy of scaling relations that are sensitive to cosmological anisotropies and BFs. In Sects. 8 and 7, the possible systematic biases are discussed, and the comparison of our results to those from isotropic Monte Carlo simulations is presented. Finally, in Sects. 9 and 10 the discussion and conclusions of this work are given.

2. Sample and measurements

As a general basis, we use the eeHIFLUGCS sample. The only exception is the use of the full Meta-Catalog of X-ray detected Clusters of galaxies (MCXC, Piffaretti et al. 2011) for the L_X − Y_SZ relation. This is done to maximize the available clusters since both L_X and Y_SZ are available beyond eeHIFLUGCS. From eeHIFLUGCS, a slightly different cluster subsample is used for different scaling relations, where all available measurements for both quantities that enter the scaling relation are considered; the vast majority of the subsamples naturally overlap. For the four scaling relations that include T, we use the same sample as in M20. In a nutshell, the sample includes 313 galaxy clusters, the vast majority of which are part of the eeHIFLUGCS sample. It is a relatively low redshift sample with median z = 0.075 (z ∈ [0.004, 0.45]), which includes objects with an X-ray flux of f_{X, 0.1−2.4 keV} ≥ 5 × 10⁻¹² erg s⁻¹ cm^{−2 1}. The Galactic plane region (b ≤ |20° |) and the regions around the Virgo cluster and the Magellanic clouds were masked and no clusters are considered from there. The spatial distribution of the sample in the rest of the sky is very homogeneous. The X-ray luminosity L_X, temperature T, redshift z, R₅₀₀, and metallicities of the core Z_core and 0.2−0.5 R₅₀₀ annulus Z_out are obtained as described in M20. The only change we make in the sample is the spectral fit of NGC 5846. When Z_out is left free to vary it results to unrealistically large values, also affecting the measured T and making NGC 5846 a strong outlier in the L_X − T plane. Thus we repeat the spectral fitting with a fixed value of Z_out = 0.400 Z_⊙ (the median of the sample), which returns T = 0.927 ± 0.013 keV.

For the rest of the scaling relations, the number of clusters used depends on the availability of each measurement, which is described below and summarized in Table 2. Finally, we also use the ASCA Cluster Catalog (ACC, Horner 2001), after excluding all the common clusters with eeHIFLUGCS and after further cleaning its X-ray luminosities, as described in M18 and M20.

2.1. Total integrated Compton parameter Y_5R500

Galaxy clusters can be observed in the submillimeter regime through the thermal Sunyaev–Zeldovich effect (tSZ). The tSZ is a spectral distortion of the CMB toward the sky positions of galaxy clusters caused by inverse Compton scattering of CMB photons by the hot electrons of the intracluster plasma. This causes a decrement (increment) in the observed temperature of the CMB at frequencies below (above) ∼217 GHz. The amplitude of this spectral distortion is proportional to the Comptonization parameter y defined as

(1)

where σ_T is the Thompson cross section, m_e is the electron rest mass, c is the speed of light, k_B is the Boltzmann constant, n_e the electron number density, T_e the electron temperature, and l.o.s. denotes the line of sight toward the cluster. The product of k_Bn_e(r)T_e(r) represents the electron pressure profile P_e(r). We adopted the form of the latter from the widely used generalized Navarro–Frenk–White (GNFW) pressure profile (Nagai et al. 2007) with the “universal” parameter values obtained by Arnaud et al. (2010). The required values of R₅₀₀ and z were taken from M20 (when available) or MCXC.

A common technique for the extraction of the tSZ signal from multifrequency data sets are matched multifilters (MMFs). These MMFs allowed us to construct a series of optimal spatial filters that are built from the spectral energy distribution (SED) of the tSZ together with a spatial template computed from the expected pressure profile of clusters. Matched multifilters have been widely used in blind cluster searches by the ACT, SPT, and Planck collaborations (Hasselfield et al. 2013; Hilton et al. 2021; Bleem et al. 2015, 2020; Planck Collaboration XXIX 2014; Planck Collaboration XXVII 2016) and allow for the estimation of the integrated Comptonization parameter Y_5R500, where Y_5R500 = ∫y dΩ, where Ω is the chosen solid angle.

We estimated Y_5R500 from the latest version of the 100, 143, 217, 353, 545, and 857 GHz all-sky maps delivered by the Planck High Frequency Instrument (HFI; Planck Collaboration III 2020), which were published during the final data release in 2018 (R3.00; Planck Collaboration I 2020). The Y_5R500 values were extracted from Planck data with MMFs that are built using the relativistic tSZ spectrum (e.g., Wright 1979; Itoh et al. 1998; Chluba et al. 2012) and have zero response to the kinematic SZ (kSZ) spectrum. We also tested values of Y_5R500 obtained using standard MMFs derived from the non-relativistic tSZ spectrum without additional kSZ removal and there were completely negligible changes in the results. Further details of the processes used to extract the Y_5R500 values are presented in detail in Erler et al. (2019, hereafter E19) and in Appendix A.1.

We measured Y_5R500 for all the 1743 entries in MCXC. We derived a signal-to-noise ratio (S/N) of > 1 for 1472 clusters, S/N > 2 for 1094 clusters, S/N > 3 for 746 clusters, and S/N > 4.5 for 460 clusters. The latter matches the threshold set by the Second Planck Sunyaev–Zeldovich Source Catalog (PSZ2, Planck Collaboration XXVII 2016, hereafter A16). We compared our derived Y_5R500 values with those from PSZ2 in Appendix A.1. For our work, we only kept clusters with S/N > 2 to avoid using most measurements dominated by random noise, but without losing too many clusters. This S/N threshold is a safe choice since we already know galaxy clusters exist at these sky positions through their X-ray detection. We always check if increasing the S/N threshold alters our results.

Finally, the cluster quantity that enters the cluster scaling relations is the total integrated Comptonization parameter Y_SZ, which is given by

(2)

where D_A is the angular diameter distance of the cluster in kpc. The dependence of Y_SZ on the cosmological parameters therefore enters through D_A.

2.2. Half-light radius R

The apparent angular size R_app in the sky within which half of the total X-ray emission of a cluster is encompassed, is a direct observable. By using a cosmological model, this observable (measured in arcmin) can be converted to a physical size of a cluster, the half-light radius R, where . The size of a cluster correlates with many other cluster properties and can be used to construct scaling relations. We measured R_app for all eeHIFLUGCS clusters, and additionally for all MCXC clusters with f_{X, 0.1−2.4 keV} ≥ 4.5 × 10⁻¹² erg s⁻¹ cm⁻². This led to 438 measurements.

To measure R_app the ROSAT All-Sky Survey (RASS) maps were used. The count-rate growth curves were extracted for all eeHIFLUGCS clusters, and additionally for all MCXC clusters with f_{X, 0.1−2.4 keV} ≥ 4.5 × 10⁻¹² erg s⁻¹ cm⁻². By a combination of applied iterative algorithms and visual validation by six different astronomers, the plateau of the count-rate growth curves (i.e., the boundaries of cluster X-ray emission) were determined. The radii R_app were then found, and converted to R using the default cosmology. The exact details of this process will be presented in Pacaud et al. (in prep.).

The point spread function (PSF) of the XRT/PSPC imager of ROSAT varies with the off-axis angle and the photon energies. In general, it is ≤1.5′. To avoid strong biases due to PSF smearing, we excluded all the clusters with R_app ≤ 2′. This left us with 418 cluster measurements with a median R_app = 3.76′. Residual PSF smearing effects are still expected to affect the measurements. Since the scaling relations including R are currently inconclusive (see Appendix C.1), we neglect these effects for now. In future work, any PSF effects will be fully taken into account. eROSITA will also be able to provide R_app values that are rather insensitive to PSF effects because of its better spatial resolution.

2.3. Near-infrared BCG luminosity L_BCG

The BCGs were found for the 387 clusters of the eeHIFLUGCS sample. To determine the BCG for all clusters, we used the optical and near-infrared (NIR) data from the SDSS (York et al. 2000), Pan-STARRS (Kaiser et al. 2002, 2010), VST ATLAS (Shanks et al. 2015), DES (Abbott et al. 2018), 2MASS (Skrutskie et al. 2006) and WISE (Wright et al. 2010) catalogs. The redshifts of the galaxies were either taken from the SDSS catalog or the NASA/IPAC Extragalactic Database (NED)². All the galaxy magnitudes were corrected for Galactic extinction and the proper k-correction was applied. The exact details on the BCG selection are described in Appendix A.2.

For this work, we use the magnitudes coming from the 2MASS catalog for two reasons. First and foremost, 2MASS returns the largest number of available BCGs for our sample. Out of the 387 clusters of eeHIFLUGCS, we detected the BCG in 2MASS for 331 of them. Secondly, the infrared BCG luminosities are not strongly sensitive to extinction effects, minimizing the potential risk of unaccounted absorption biases. We exclude all BCGs with z < 0.03, as they appear to be systematically overluminous based on all L_BCG scaling relations, indicating a flattening of the L_BCG relations at low cluster masses. This flattening was also observed by Bharadwaj et al. (2014).

Additionally, the redshift evolution of the BCG luminosity versus the other quantities is unknown. As a consequence of the large scatter of the L_BCG scaling relation as shown later in the paper, the applied evolution cannot be left free to vary simultaneously with the other scaling relation parameters because no reliable constraints are obtained. Subsequently, we also opt to exclude all clusters with z ≥ 0.15. This allows us to ignore any evolution during the model fitting. Even if this added a small bias in our estimates, it would not be expected to affect the anisotropy analysis because the cluster redshift distributions of different sky regions are similar (see M20). As such, any potential bias would cancel out when comparing cluster subsamples from different sky regions. These criteria eventually leave us with 244 clusters with 2MASS L_BCG measurements with a median redshift of 0.069.

2.4. X-ray determined N_H, Xray

We determined the total hydrogen column density N_H,Xray using the X-ray spectra of the 313 clusters from M20. The exact details of our methodology can be found in Appendix A.3. Overall, we were able to obtain a safe estimation of N_H,Xray for 156 clusters. Several systematics might creep in during the whole process, hence we approach the analysis done with these measurements conservatively.

2.5. X-ray luminosity L_X and redshift z for clusters not included in M20

For the 1430 MCXC clusters not included in the M20 sample, we started from their L_X values given in MCXC, which are corrected for the absorption traced by the neutral hydrogen column density only. We further corrected these values to account for the total hydrogen absorption based on the Willingale et al. (2013, hereafter W13) N_Htot values. The procedure is exactly the same as that followed in M20. The only difference here is the use of fixed T = 5 keV and Z = 0.4 Z_⊙ values for all clusters in the XSPEC (Arnaud 1996) apec⋅phabs model because we do not have spectral measurements for these clusters. The redshift values were adopted from MCXC.

2.6. ACC sample

We measured Y_SZ for the ACC sample (Horner 2001) as well, using the same procedure as for our sample. We excluded all the clusters already included in any of our different subsamples. We also excluded the 55 clusters with Y_SZS/N < 2. This results in 168 clusters with X-ray luminosity and temperature values; 113 of these clusters have a Y_SZ measurement with S/N > 2. Cross-checking the results of a completely independent cluster sample with the results from our sample is crucial to understand the origin of the observed anisotropies (e.g., to exclude that sample selection effects may bias our results). The properties of the ACC sample as we use it are already given in detail by M18 and M20. In a nutshell, the sample mainly consists of massive galaxy clusters, spanning across z ∈ [0.009, 0.839], with a median redshift of 0.226 (for the 113 clusters with Y_SZ measurements). The temperatures are obtained by a single-thermal model for the whole cluster, while the X-ray luminosities L_bol are given in the bolometric band within the R₂₀₀ of the cluster (measured within 0.5–2 keV and within a “significance” radius, and then extrapolated). All measurements are performed by Horner (2001) and the ASCA telescope, while we corrected the L_bol for the total absorption similarly to the L_X corrections of our sample. The necessary R₅₀₀ values to measure Y_SZ were obtained using the mass-temperature scaling relation of Reichert et al. (2011). The latter mostly uses XMM-Newton-derived temperatures, which might differ from the ASCA temperatures in general. Thus, we compared the ACC temperatures with our temperatures for the common clusters between the two samples and applied the necessary calibration factors before calculating R₅₀₀.

3. Scaling relations

We study ten cluster scaling relations in total, namely the L_X − Y_SZ, L_X − T, Y_SZ − T, L_X − L_BCG, L_BCG − T, R − Y_SZ, R − L_X, R − L_BCG, R − T, and Y_SZ − L_BCG relations. The first nine relations are sensitive to either additionally needed X-ray absorption corrections (1), BFs (2), or possible cosmological anisotropies (3). The Y_SZ − L_BCG cannot trace any of the above effects (4), as explained in Appendix C.3. An observed anisotropy could then point to systematics in the measurements or methodology, which may affect the other nine scaling relations of interest.

In principle, it is very challenging to distinguish cases (2) and (3) since their effects on the observed scaling relations are similar. One way to distinguish between the two is to perform a tomography analysis, analyzing redshift shells individually and see if the anisotropies persist at all scales. Evidently, the power with which each scaling relation traces a possible effect varies, since it depends on the scatter of the relation and the exact way each effect might intervene with each measurement. The information of what effect each scaling relation can detect is given in Table 1. The exact explanation on how a scaling relation detects (or not) an anisotropy origin is given in Sects. 5 and 6.

The form of the studied scaling relations between some measured cluster quantities Y and X is given by

(3)

where C_Y is the calibration term for the Y quantity, the term E(z)=[Ω_m(1 + z)³ + Ω_Λ]^1/2 accounts for the redshift evolution of the Y − X relation, γ_YX is the power index of this term, A_YX is the normalization of the relation, C_X is the calibration term for the X quantity, and B_YX is the slope of the relation. The calibration terms C_Y and C_X are taken to be close to the median values of Y and X, respectively. These values are shown in Table 2, together with the assumed (self-similar) values of γ_YX.

3.1. Linear regression

Similar to M20, to constrain the scaling relation parameters we perform a linear regression in the logarithmic space using a χ² minimization procedure. We consider two separate cases.

The first case is when we assume to know the universal, isotropic cosmological parameters (or when they do not matter owing to canceling out between the two parts of the scaling relation) and wish to constrain the normalization, slope, and scatter of a scaling relation. We then constrain the desired parameters by minimizing the expression

(4)

where N is the number of clusters, , ³, σ_log Y and σ_log X are the Gaussian logarithmic uncertainties for Y and X, respectively (derived as in M20), and σ_int, YX is the intrinsic scatter of the Y − X relation with respect to Y in orders of magnitude (dex). Following Maughan (2007), σ_int, YX is iteratively increased and added in quadrature as an extra uncertainty term to every data point until the reduced ⁴. σ_tot, YX is the total scatter, equal to the average value of the denominator of Eq. (4) for all considered clusters.

The second case only applies when there is a strong cosmological dependency on the best-fit scaling relations (i.e., all the T scaling relations). In this case we assume the normalization of the scaling relation to be known and direction-independent. This is a reasonable assumption since this quantity is associated with intrinsic properties of the clusters and as such there is no obvious reason why it should spatially vary. In this case, the free parameters we wish to constrain are the Hubble constant H₀ or the BF amplitude and direction u_BF, while the slope is treated as a free, nuisance parameter. To do so, we minimize the following equation:

(5)

where D is either the luminosity distance D_L (e.g., for L_X − T) or the angular diameter distance D_A (e.g., for Y_SZ − T); D_obs is the observed distance given the measurements Y and X and their exact scaling relation Y − X; D_th is the theoretically expected distance based on the cosmological parameters (e.g., H₀); the redshift z and the existing BF u_BF; σ_Di, obs is the statistical uncertainty of the observed distance, which is a function of the measurement uncertainties of Y and X; and σ_int,D is the intrinsic scatter of the relation in megaparsec units. The observed distance D_obs enters every scaling relation differently. Generally, it is given by

(6)

In this equation, D_H0 = 70 is the distance found for H₀ = 70 km s⁻¹ Mpc⁻¹, which enters in the calculation of Y. The use of D_H0 = 70 in Eq. (6) cancels out this cosmological contribution to Y. For this, we need for L_X − T, Y_SZ − T, R − T, and L_BCG − T, respectively. This leaves us with X-ray flux (in the rest frame of the cluster), apparent size in arcmin, Y_5R500, and BCG flux respectively. Moreover, D_th is given by

(7)

where the (1 + z)^±1 factor depends on if we are considering D_L or D_A. Also, z is the cosmological redshift, which is given by (e.g., Harrison 1974; Dai et al. 2011; Feindt et al. 2013; Springob et al. 2014), where z_obs is the observed redshift (converted to the CMB rest frame), u_BF is the amplitude of the BF, and ϕ is the angular distance between a cluster and the BF direction. We should note that the sign of in front of u_BF is changed from “−” to “+” to avoid confusion with negative velocities.

We should stress that when searching for cosmological anisotropies, the calculated H₀ variations express relative differences between regions. The absolute H₀ values cannot be constrained by cluster scaling relations only since a fiducial H₀ value was initially assumed to calibrate the relation. Hence, the H₀ anisotropy range always extends around the initial H₀ choice. Apparent H₀ fluctuations could also mirror other underlying cosmological effects we are not yet aware of.

3.2. Bulk flow detection

We follow two different methods to estimate the best-fit BF. Firstly, we fit the scaling relations to all the clusters, adding a BF component in our fit, as in Eqs. (5) and (7). The free parameters are the direction and amplitude of the BF along with the intrinsic scatter. The parameters A_YX and B_YX are left free to vary within their 1σ uncertainties of the considered (sub)sample, when no BF is applied. The redshift of every cluster (and thus D_th) is affected by the BF, depending on the angle between the direction of the cluster and the BF direction. Since the spatial distribution of clusters is nearly uniform, any applied BF does not significantly affect the best-fit A_YX and B_YX, but it does affect the scatter. This procedure essentially minimizes the average Y residuals around each scaling relation; that is, we get for a smaller σ_int, YX; the procedure is thus labeled the “minimum residuals” (MR) method.

The second method we follow to detect the best-fit BF is that of the “minimum anisotropies” (MA) method. When we add a BF component directly toward the most anisotropic region, other regions start altering their behavior and might appear anisotropic. Absolute apparent isotropy cannot then be achieved. Therefore, the amplitude and direction of the BFs are found so that the final, overall anisotropy signal of the studied relation is minimized. For this, we repeat the anisotropy sky scanning (Sects. 3.4 and 3.5) for a different BF every time. This procedure is computationally expensive, and thus we reduce the number of bootstrap resamplings to 500 when estimating the uncertainties of the BF characteristics.

We consider independent redshift shells (redshift tomography) and cumulative redshift bins to constrain the BF motions with both the MR and MA methods. It is important to note that flux-limited samples suffer from certain biases when determining BFs. Probably the most important bias is the nonuniform radial distribution of objects. If this is not taken into account, the BFs found for spherical volumes with a large radius are strongly affected by low-z objects and do not clearly reflect the larger-scale motions. This can be partially accounted for if the iterative redshift steps with which we increase the spherical volumes encompass similar numbers of clusters (e.g., Peery et al. 2018, and references therein). That way, the contribution of every scale to the overall BF signal is averaged out. In this work, to minimize the biased contribution of low-z systems to the larger volumes, we consider redshift bins with similar cluster numbers. Finally, the size and redshift range of our sample did not allow us to use the kSZ signal of the clusters to detect any BFs (see Sect. 9).

3.3. Parameter uncertainties

To estimate the uncertainties of the best-fit parameters, we use a bootstrap resampling method with replacement. We randomly draw 5000 resamplings from the studied (sub)sample of clusters with the same size, simultaneously constrain their best-fit parameters, and obtain the final distribution of the latter. The quoted parameter uncertainties refer to the 68.3% credible interval defined by the positive and negative 34.1_th percentile of the distributions with respect to the best-fit value of the full (sub)sample. Each parameter value distribution is considered separately from the rest⁵. This method provides more conservative and robust parameter uncertainties than the approach followed in M20, since it depends only on the true, random variation of the studied statistic and not on analytical expressions (e.g., Δχ² limits).

3.4. Detection of anisotropies, parameter sky maps, and direction uncertainties

The methodology followed in this work is described in detail in Sect. 4.3 of M20. In a nutshell, we consider cones of various radii (θ = 60° −90°) and point these cones toward every possible direction in the sky (with a resolution of ≤1°). Each time, we consider only the clusters within each cone and we obtain the best-fit parameters with their uncertainties, as described in the two previous sections. We then create a color-coded full sky map based on the best-fit parameters of every direction.

There are only two differences with the method followed in M20. Firstly, in this work we leave the slope B_YX free to vary, instead of fixing it to its best-fit value for the full sample. That way, the parameter of interest is marginalized over B_YX. In M20 we demonstrated that this choice does not significantly alter our results compared to the case in which B_YX is kept fixed, but it constitutes a more conservative approach and thus we make it the default. Secondly, we slightly change the applied statistical weight of the fitted clusters, which depends on their angular distance from the center of the cone. In M20 we divide the statistical uncertainties σ_log Y and σ_log X with a normalized cosine factor that shifted from 1 (center of cone) to 0 (edge of cone), despite the actual radius of the cone. Under certain conditions however⁶, this method can slightly overestimate the final statistical significance of the observed anisotropies. Even if this does not have an effect on the conclusions of M20, in this work we choose to follow a more conservative approach, dividing the uncertainties with the cosθ₁ term, where θ₁ is the angular separation of a cluster from the center of the cone we consider. This is motivated by the fact that the effects of a BF or an H₀ anisotropy to the measured cluster distance scale with cosθ₁. Also, this results in the same weighting as in M20 for θ = 90° and in a weaker weighting for θ < 90°.

To provide anisotropy direction uncertainties, we perform bootstrap resamplings similar to those described in Sect. 3.3. For every sample used to create a sky map, we create 1000 bootstrap resamplings and perform the sky scanning again for each. The 68.3% limits of the posterior distribution of the maximum anisotropy direction are reported as the 1σ limits.

3.5. Statistical significance of anisotropies

We wish to assess the statistical significance of the observed differences of the behavior of the scaling relations toward different directions. Briefly, we obtain the best-fit values of the fitted parameters together with their uncertainties for every direction (cone) in the sky. We then identify the region that shares the most statistically significant anisotropy from the rest of the sky (similar to a dipole anisotropy)⁷. We assess the statistical significance of the deviation between them by

(8)

where p_1, 2 are the best-fit values for the two independent subsamples and σ_p1, 2 are their uncertainties derived by bootstrapping. Finally, the statistical significance (sigma) maps are color coded based on the observed anisotropy level between every region and the rest of the sky.

3.6. Monte Carlo simulations

To further validate the effectiveness of our methodology and the statistical significance of the observed anisotropies, we perform Monte Carlo (MC) simulations. We create isotropic cluster samples to which we apply the same procedure as in the real data. That way we estimate the frequency with which artificial anisotropies would be detected in an isotropic Universe, and compare this to the real data. More details are described in Sect. 7.

3.7. Summary of statistical improvements compared to M20

In this section, we summarize the improvements in the statistical analysis of this work compared to M20: (1) The parameter uncertainties at every stage of this work are found by bootstrap resampling instead of Δχ² limits. (2) During the sky scanning for identifying anisotropies, all the nuisance parameters (e.g., the slope B) are left free to vary, instead of fixing them to their best-fit values. (3) The statistical weighting of clusters during the anisotropy searching is relaxed to avoid creating any artificial anisotropies. (4) Uncertainties of the anisotropy directions are provided. (5) Monte Carlo simulations are carried out to further assess the statistical significance of the anisotropies.

4. General behavior of the ten scaling relations

As a first step we constrain the overall behavior of the observed scaling relations when considering all the available data from across the sky. The effects of selection biases are discussed in Sect. 8.7. The overview of the best-fit results of all the scaling relations is given in Table 2, while the scaling relations themselves are plotted in Fig. 1.

4.1. The L_X − T relation

The full analysis of the L_X − T relation for our sample is presented in detail in M20. The only changes compared to the M20 results are the slightly changed T for NGC 5846 (see Sect. 2), and the use of bootstrap resampling for estimating the parameter uncertainties. Accounting for these changes, the best-fit values for the L_X − T relation remain fully consistent with M20 and with previous studies. A detailed discussion is available in Sect. 5.1 of M20, while a more recent work confirming our results can be found in Lovisari et al. (2020).

The L_bol − T results for ACC are shown in M18 and M20. The only difference in this work is that we perform the χ²-minimization in the L_bol axis, using bootstrap for estimating the parameter uncertainties. Since the bolometric X-ray luminosity is used for ACC (within R₂₀₀), both A_LT and B_LT are larger than the results of our sample. Also, σ_int is ∼38% larger for ACC.

4.2. The L_X − Y_SZ relation

The L_X − Y_SZ scaling relation has been studied in the past (e.g., Morandi et al. 2007; Planck Collaboration X 2011; Planck Collaboration XI 2011; De Martino & Atrio-Barandela 2016; Ettori et al. 2020; Pratt & Bregman 2020) mostly using Y_SZ from Planck and L_X from ROSAT data and the MCXC catalog. Both of these quantities are efficient proxies of the total cluster mass; they also scale with each other. Their scatter with respect to mass is mildly correlated (e.g., Nagarajan et al. 2019). This results in L_X − Y_SZ having the lowest scatter among all the scaling relations used in this study.

As mentioned in Sect. 2.1, we measured Y_SZ for 1095 MCXC clusters with S/N > 2. In studying the L_X − Y_SZ relation for these objects, we see that there are significant systematic differences between cluster subsamples based on their physical properties. For example, clusters with low N_Htot or high z tend to be significantly fainter on average than clusters with high N_Htot or low z, respectively. Surprisingly, the same behavior persists even when the original MCXC L_X values and the Y_SZ values from PSZ2 are used. As we increase the S/N threshold, these inconsistencies slowly fade out. More details on these effects can be found in Appendix D. Owing to this, we chose to apply a low S/N ≥ 4.5 threshold to the Y_SZ values (same as in the PSZ2 sample), leaving us with 460 clusters with a median z ∼ 0.14. The clear benefits of this choice are that clusters with different properties now show fully consistent L_X − Y_SZ solutions and no systematic behaviors are observed. Additionally, the intrinsic scatter of the L_X − Y_SZ relation decreases drastically compared to cases with lower S/N thresholds, allowing us to put precise constraints on the best-fit parameters and the possibly observed anisotropies of the relation. For these 460 clusters, the best-fit values are in full agreement with past studies within the uncertainties. The scatter we obtain however is lower than most past studies, most probably resulting from the use of the same R₅₀₀ between L_X and Y_SZ in our analysis. The observed slope B_LY ∼ 0.93 is slightly larger than the self-similar prediction of B_LY = 0.8.

For ACC, no trends for the Y_SZ residuals are observed for S/N > 2, with any of the cluster parameters. Therefore all 113 clusters can be used. The slope lies again close to unity, while the scatter is ∼32% larger than when our sample is used. Nevertheless, the L_X − Y_SZ scatter for ACC is sufficiently small to allow for precise constraints. We should note that for S/N > 4.5, the scatter of ACC is similar to our sample, but this cut would leave us with only 67 clusters; this number is not sufficient for our purposes.

4.3. The L_X − L_BCG relation

The L_X − L_BCG scaling relation has not been extensively studied in the past. Mittal et al. (2009) and Bharadwaj et al. (2014) used 64 and 85 low-z clusters, respectively, to compare the L_X − L_BCG relation between cool- and non-cool-core (CC; NCC) clusters, finding mild differences mostly for the slope. Furnell et al. (2018) also studied the correlation of the stellar mass of BCGs (which is proportional to its luminosity) with the L_X of the cluster. We use significantly more clusters to study the L_X − L_BCG relation, namely 244. The best-fit slope of our analysis is less steep than the best-fit slope found by Bharadwaj et al. (2014); in contrast to this work the latter authors used the bolometric X-ray luminosity. The scatter we obtain is considerably smaller than theirs, although still significantly large. This is the scaling relation with the largest scatter out of the ten we examined. The L_X residuals do not show any systematic behavior as a function of cluster properties.

4.4. The R − L_X relation

The relation between R and L_X for galaxy clusters has not been investigated before to our knowledge. In this work we use the 418 clusters with both measurements available to constrain this relation. The redshift evolution of R − L_X is unknown, thus we attempt to constrain it since the high number of clusters and the low scatter allowed us to do so. We leave γ_RLX free to vary, simultaneously with the rest of the parameters. This results in γ_RLX = −2.15 ± 1.51. The implied evolution is not statistically significant since the limited redshift range of the sample does not allow us to constrain it more efficiently. The same evolution as in the R₅₀₀ − L_X relation (self-similar prediction of γ ∼ 1.3) is not expected since a redshift evolution between the half-light radius R and R₅₀₀ is also expected (e.g., owing to the time-varying, CC cluster fraction).

Since γ_RLX = −2.15 describes our data best, we fix γ_RLX to this value for the rest of our analysis. After performing the fit for all 418 clusters, we note that clusters at z < 0.01 are systematically down-scattered compared to the rest. To avoid any biases during our anisotropy analysis, we excluded the latter clusters. The final subsample used consists of the remaining 413 clusters. The R − L_X scatter residuals appear to be randomly distributed with respect to z, RASS exposure time, and N_Htot. However, a non-negligible correlation is observed between the residuals and the apparent half-light radius R_app. The clusters with the lowest R_app appear to be down-scattered in the R − L_X plane, and vice versa. A mild correlation of the R residuals is also observed with the offset between the X-ray peak and the BCG position. The latter can be used as a tracer of the dynamical state of the cluster, which correlates with the existence (or not) of a CC in the center of the cluster (e.g., Hudson et al. 2010; Zitrin et al. 2012; Rossetti et al. 2016; Lopes et al. 2018). Cool-core clusters are expected to strongly bias the scaling relations involving R, since they emit most of their X-ray photons from near their centers. As a result, the half-light radius is lower compared with non-CC clusters at a fixed mass. More details on that can be found in Appendix C.1.

Finally, the slope is lower than the self-similar prediction for the R₅₀₀ − L_X relation (B ∼ 0.33), while there is a moderate scatter. Surprisingly, the latter is the largest observed between all R scaling relations.

4.5. The Y_SZ − T relation

The Y_SZ − T relation has been previously studied by several authors, for example, Morandi et al. (2007), Planck Collaboration XI (2011), Bender et al. (2016), and Ettori et al. (2020) among others. This relation has never been studied before with such a large number of clusters as used in this work. Since the T measurement is needed, we use the sample from M20. We retrieved 260 clusters with Y_SZ measurements with S/N > 2 from the M20 sample. No systematic differences in the Y_SZ − T relation are observed for different cluster subsamples with different properties. The Y_SZ residuals remain consistent with zero with increasing T, z, N_Htot, and other cluster parameters, while A_YT, B_YT and σ_int, YT also stay constant with an increasing S/N cut. These results clearly indicate that the applied S/N threshold does not introduce any strong biases to the Y_SZ − T relation (see Appendix D.2 for more details). The best-fit parameter values that we obtain for these clusters are in line with previous findings. The value of the slope agrees with the self-similar prediction (B_YT = 2.5), while the scatter is lower than that for the L_X − T relation. Finally, a single power law perfectly describes the relation since a change in A_YT, all and B_YT, all is not observed for a changing low T cut.

When using the ACC sample, we obtain a similar scatter with our sample, but with a slightly steeper slope and larger normalization. This is because temperatures are measured for ACC considering the entire cluster, leading to generally lower T than our sample (where T is measured within 0.2−0.5 R₅₀₀).

4.6. The Y_SZ − L_BCG relation

The Y_SZ − L_BCG scaling relation has not been studied in the past, and it is constrained for the first time in this work. The two quantities are expected to scale with each other since they both scale with the cluster mass. We have measurements for 214 clusters with S/N > 2 for the Y_SZ measurement and the applied redshift limits for the BCGs. When we performed the fit, we did not detect any strong systematic behavior of the residuals as a function of cluster parameters, with the exception of Y_SZ S/N, with high Y_SZ S/N clusters tending to be up-scattered. This behavior persists even with an increasing S/N threshold. Although the best-fit B_YLBCG stays unchanged for different S/N cuts, the A_YLBCG, varies by ∼3.4σ for S/N > 4.5 (Fig. D.1). Since the residuals are a function of S/N in any case, and since the Y_SZ − L_BCG relation offers only limited insight into the anisotropy analysis, we adopt S/N > 2, but suggest caution because of the aforementioned dependence of A_YLBCG. The slope lies close to linearity, while the scatter is ∼10% smaller than the L_X − L_BCG scaling relation.

4.7. The R − Y_SZ relation

A relation between Y_SZ and the X-ray isophotal radius R is expected to exist since both quantities scale with cluster mass. Such a relation has not been observationally constrained however until now. For this work, both quantities were measured for 347 clusters. The redshift evolution of the relation is not known, however owing to the small scatter this value can be obtained observationally, similarly to the R − L_X relation. We find that the scatter is minimized for γ_RYSZ = −2.72 ± 1.41. Expectingly, the result is similar to γ_RLX. We note that we look for the evolution describing our data best and not necessarily for the true evolution. Fixing γ_RYSZ to its best-fit value, we repeat the fitting for the rest of the parameters. The residuals of R behave in exactly the same way as for R − L_X. While no trend is seen with respect to z, the systematic behaviors as functions of R_app and of the X-ray-BCG offset persist. The observed scatter is the smallest scatter between all R scaling relations, ∼20% smaller than the R − L_X scatter. Finally, the slope of R − Y_SZ is also the flattest compared to the other scaling relations.

4.8. The L_BCG − T relation

Similar to the L_X − L_BCG relation, the L_BCG − T relation has not been the focus of many studies in the past. Bharadwaj et al. (2014) studied the scaling of L_BCG with the mass of clusters. The latter is however obtained through a cluster mass–temperature scaling relation. For the L_BCG − T relation the exact redshift evolution is also not known, and constraining this value from the existing data is not trivial. This is because of the large scatter of the L_BCG − T relation and the other simultaneously fit parameters. As a result, from the 259 clusters for which we measured T and L_BCG, we again use only the 196 of the clusters with z < 0.15 to safely ignore any existing redshift evolution of the relation. The scatter of the L_BCG − T relation is considerably large, namely ∼70% larger than the L_X − T scatter when the latter is minimized with respect to T. The scaling between the two quantities however is clear, although with a less steep slope than the reversed L_X − T scaling (∼2.8). The slope has also a larger relative uncertainty (∼10%) than the normalization (∼5%).

4.9. The R − L_BCG relation

The R − L_BCG is another scaling relation that is presented for the first time in this work. Both measurements are available for 243 clusters, after the previously described L_BCG redshift cuts. The intrinsic scatter of the relation is similar to the other R relations, while the slope is slightly larger than the R − L_X and R − Y_SZ relations. The same behavior for the R residuals is observed as in the R − L_X and R − Y_SZ scaling relations.

4.10. The R − T relation

The R − T relation has not been studied extensively in the past, although some studies with both observations and simulations were performed two decades ago (Mohr & Evrard 1997; Mohr et al. 2000; Verde et al. 2001). The relation was used as a cosmological probe as well by Mohr et al. (2000) to constrain Ω_m and Ω_Λ. These authors however used only a few tens of clusters to constrain the relation. In this work, we used 308 clusters for which both R and T were measured. Mohr et al. (2000) argued that there is no redshift evolution for the R − T relation. However, this conclusion specifically depends on their methodology for measuring R and cannot be adopted for our method. Therefore, we choose to constrain the redshift evolution from our data, leaving it free to vary. We obtain γ_RT = −1.98 ± 1.42, similar to the other R relations. Fixing γ_RT = −1.98, we obtain the final best-fit results. The slope is the largest (∼0.57) among all R scaling relations. The same R residual behavior as for the other R scaling relations is observed.

5. Anisotropies due to unaccounted X-ray absorption effects

In this section we study the scaling relations, whose observed anisotropies could be caused purely by previously unknown soft X-ray effects, such as extra absorption from “dark”, metal-rich, gas and dust clouds. If the anisotropies observed in the L_X − T relation by M20 were due to such effects, we should detect the same anisotropies in the scaling relations studied in this section. We report the most anisotropic directions and the statistical significance of the observed tension and quantify the amount X-ray absorbing material that should exist to fully explain the discrepancy.

5.1. L_X − Y_SZ anisotropies

5.1.1. Our sample

The anisotropies of the L_X − Y_SZ relation are the focal point of the search for hidden X-ray effects that we were not aware of in the past. This relation exhibits the lowest scatter and the largest number of clusters among all scaling relations studied in this work. This allows for precise pinpointing of its anisotropies. Even more importantly, L_X − Y_SZ is almost completely insensitive to any spatial H₀ variations, since both quantities depend on cosmological parameters in the same way. The two quantities also scale almost linearly with each other. Thus, if we changed H₀, no significant change in the best-fit A_LYB_LY would be observed. Analytically, as and , their ratio (considering their best-fit B_LY = 0.928) would be

(9)

Thus, a ∼15% spatial variation of H₀ would cause a nondetectable ∼2% variation in A_LY. Moreover, the L_X − Y_SZ relation is very insensitive to BFs as well. Based on the above calculation, if a BF of ∼1000 km s⁻¹ existed at z ∼ 0.05 toward a sky region, this would only lead to a < 2% increase in A_LY of this region.

Therefore, any statistically significant anisotropies in the L_X − Y_SZ relation should mainly be caused by unaccounted X-ray absorption effects acting on L_X. To scan the sky, we adopt a θ = 60° cone. This returns at least 72 clusters for each cone, which, considering the very low scatter of the relation, are sufficient to robustly constrain A_LY. The variation and significance maps are shown in Fig. 2. We detect the most anisotropic sky region toward (l, b)=(118°, +7°), which shares a 3.5σ anisotropy with the rest of the sky. The 100 clusters within this region appear to be 12 ± 3% fainter on average than the rest of the sky. The extra N_Htot needed to explain this discrepancy is ∼3.9 ± 1.1 × 10²⁰ cm⁻² (the uncertainties were symmetrized). If the assumed hydrogen quantity based on W13 was the true quantity, then the metallicity of the absorbing material toward that direction would need to be Z ∼ 1.54 ± 0.16 Z_⊙ (currently assumed to be Z = Z_⊙). Considering that the specific region lies close to the Galactic plane, this metallicity value does not seem unlikely. The A_LY and the anisotropy significance maps of L_X − Y_SZ are represented in Fig. 2. We can see that there are no anomalously bright regions. This indicates that there are no regions with significantly lower-than-solar metallicities of the Galactic material. Assuming availability of a much larger number of clusters with L_X and Y_SZ measurements, ideally extending to low Galactic latitudes, this scaling relation could potentially be used as a new probe of the ISM metallicity.

It should be stressed that in M20, the most anisotropic direction for the L_X − T relation was found to be (l, b)∼(300°, −20°). Based on the above test, this region does not show any signs of extra, previously unaccounted absorption. Adding up to the numerous tests done in M20, this further supports the hypothesis that the observed L_X − T anisotropies are not caused by unmodeled Galactic effects.

5.1.2. Joint analysis of the L_X − Y_SZ relation for our sample and ACC

If the L_X − Y_SZ anisotropies seen in our sample originate by the unaccounted absorption effects of a yet undiscovered mass (or a higher interstellar gas and dust metallicity), then we should obtain similar results for the ACC sample. Before extrapolated to L_bol, the flux of the ACC clusters was initially measured within the 0.5–2 keV energy range, and thus is sensitive to X-ray absorption effects. Therefore, jointly analyzing the two samples should provide us with better insight into any possible X-ray absorption issues.

To combine the results of the two independent cluster samples, we perform a joint likelihood analysis. The applied method is that followed in M20 (Sect. 8.2), where we determine the overall apparent H₀ variation. In this work, the joint parameter is the L_X − Y_SZ normalization for every region over the normalization of the full sample (A_LY/A_LY, all), marginalized over the slope. We extract the posterior likelihood of A_LY/A_LY, all for every sky region, for both samples. Then by multiplying the two posterior likelihoods, we obtain the combined, final one.

Performing the joint analysis, we find that the most anisotropic region lies toward (l, b)=(122°, +8°), where the clusters appear to be fainter in average by 11 ± 5% than the rest of the sky. Our sample dominates the joint fitting owing to the much higher number of clusters and lower scatter. However, since ACC does not show a strongly deviating behavior toward that region, the statistical significance of the anisotropy drops to just 2.1σ. To fully explain the mild tension, an undetected N_Htot ∼ 3.5 ± 1.7 × 10²⁰ cm⁻² would be needed, or alternatively a Z ∼ 1.49 ± 0.23 Z_⊙ for the already-detected Galactic gas and dust.

As such, the tension could be attributed to chance and not necessarily to an unaccounted X-ray absorption on top of the already applied absorption correction. This is also indicated by the MC simulations later on. The normalization and sigma maps of the L_X − Y_SZ anisotropies can be found in Fig. 2.

5.2. L_X − L_BCG anisotropies

Following the same reasoning as for L_X − Y_SZ, the L_X − L_BCG scaling relation cannot detect cosmological anisotropies or BFs, since both quantities depend on D_L and the slope is close to unity. Also, L_BCG is not expected to suffer from any excess extinction, since the near-infrared K_s filter of 2MASS shows a nearly negligible sensitivity on extinction (Schlafly & Finkbeiner 2011), and the L_BCG values were already corrected for the known extinction effects. Consequently, the only origin of any observed (statistically significant) anisotropies should be a stronger true X-ray absorption than the adopted value, affecting L_X. We adopt a θ = 75° cone to scan the sky so each cone contains ≥70 objects. Fewer clusters per cone might lead to strong cosmic and sample variance effects, which can result in overestimated anisotropic signals. Empirically, we conclude that 70 is a sufficient number of clusters to (mostly) avoid such effects. The anisotropy maps of L_X − L_BCG are shown in the bottom panel of Fig. 2.

The most anisotropic region turns out again to have the lowest A_LXLBCG, toward (l, b)=(171°, −22°). It appears to be 23 ± 14% dimmer than the rest of the sky. Its statistical significance however does not overcome 1.7σ, and therefore the relation is consistent with being statistically isotropic. This might be due to the large scatter and parameter uncertainties of the relation and not necessarily due to the lack of anisotropy-inducing effects. The reported direction is also 57° away from the faintest direction found by the joint L_X − Y_SZ analysis, although within ≤1.5σ. The necessary excess N_Htot to explain this mild discrepancy in the L_X − L_BCG relation is ∼8.7 ± 5.1 × 10²⁰ cm⁻². Alternatively, a metal abundance of Z ∼ 2.4 ± 0.8 Z_⊙ of the already-detected hydrogen cloud would also alleviate this small tension. Thus, the L_X − L_BCG relation does not show any indications of previously unknown X-ray absorption, although the large scatter of the relation limits the confidence of our conclusions. Finally, it is noteworthy that for once more, no apparent anisotropy exists toward (l, b)∼(300°, −20°).

5.3. Comparison between N_H, Xray and N_Htot

As a final test for detecting potential excess absorption effects, we compare the X-ray determined N_H,Xray with the N_Htot value given in W13, used in our default analysis. If a region showed a systematically larger N_H,Xray, it could indicate an extra, previously unaccounted X-ray absorption taking place toward there. We perform this comparison for the 156 clusters left after the cuts we applied, which are described in Sect. 2.4. The best-fiT relation is cm⁻². The X-ray based values are systematically higher than N_Htot, except for the high N_Htot range where the two values converge. This behavior marginally agrees with the findings of Schellenberger et al. (2015), who used ChandraT measurements and the same abundance table as in this work. However, this N_Htot discrepancy should not be taken at face value, since the overall comparison is biased by the exclusion of clusters with large N_H,Xray uncertainties. As discussed in Appendix A.3, these are mostly low N_H,Xray clusters lying below the equality line since the N_H,Xray measurement for these clusters is very challenging owing to the lower spectral cut at 0.7 keV. This selection effect would, therefore, tend to flatten the slope. Thus, this systematic up-scatter of the remaining N_H,Xray values, does not necessarily mean that the true X-ray absorption is higher than previously thought, but it is probably the result of unaccounted systematics. To avoid such issues as much as possible, we simply compare the results of a region against the rest of the sky to evaluate if this region deviates more from the W13 values compared to the rest of the sky. This assumes that any systematics should not be direction-dependent. It should be borne in mind that this is an approximate, complementary test as a result of its limitations and not a stand-alone check for excess X-ray absorption.

Firstly, we consider the region within 45° around (l, b)=(281°, −16°). This is the most anisotropic and faintest region of the L_X − T relation as found in M20, when only our sample is considered. For the 16 clusters lying within this region We find that the best-fiT relation is cm⁻². As also shown in the top panel of Fig. 3, the behavior of this region is completely consistent with the rest of the sky. Assuming that the low N_Htot clusters have a larger, true N_H,Xray, then the clusters of the tested region would be less biased, since they have a larger median N_Htot, compared to the rest of the sky. Consequently, there is no indication that an untraced, excess X-ray absorption is the cause behind the L_X − T anisotropies.

The second region we consider is within 45° around (l, b)=(122° , + 8°). This region shows mild indications of uncalibrated true X-ray absorption that differs from the W13 values. Using its 11 clusters, we find cm⁻². The results (Fig. 3) are again consistent with the rest of the sky and do not reveal any signs of a biased applied X-ray absorption correction.

5.4. Overall conclusion on possible X-ray biasing effects on the anisotropy studies

The main purpose of the analysis in this section was to address if the strong observed L_X − T anisotropy found in M20 is caused by unaccounted for soft X-ray absorption that is not traced by the N_Htot values of W13. For instance, a yet undiscovered gas and dust cloud, or some galactic dust with oversolar metallicities, could cause such an effect. Additionally, we also wished to discover if there is such a region elsewhere in the sky.

All the different tests we performed, mostly independent from each other, did not show any signs of possible absorption biases toward (l, b)∼(300°, −20°). Future tests, particularly with the eROSITA All-Sky Survey (eRASS), will reveal more information on that topic. For now though, we can safely conclude that the anisotropies found in M20 are not the result of any biases in the applied soft X-ray absorption correction.

If a mysterious hard X-ray absorbing material existed that did not affect soft X-rays, this would have an effect only on the measured T. However, this would lead to underestimated T values, which would cause the opposite anisotropic behavior compared to the M20 results. So we also reject this hypothetical scenario.

Finally, using the L_X − Y_SZ relation for our sample and ACC, we identified a region, (l, b)∼(122° , + 8°), which might imply that an extra soft X-ray absorption correction is needed for the clusters there. However, the statistical significance of this result is only ∼2σ, while the L_X − L_BCG and the N_H,Xray tests did not reveal any deviations toward that region. Future work will again give a clearer answer, but for now we conclude that no extra correction is needed for the L_X values of the clusters lying within this region.

6. Anisotropies due to bulk flows or H₀ variations

We have established that strong biases in the X-ray measurements related to previously uncalibrated absorption do not seem to exist, particularly toward (l, b)∼(300°, −20°). The anisotropies observed in M20 thus cannot arise from such issues. Two more possible origins of such anisotropies are BF motions or spatial variations of cosmological parameters. The first adds a systematic, non-negligible, local velocity component on the measured redshift of the objects toward a particular direction. If not accounted for, it leads to a systematic over- or underestimation of the distances of the clusters, and hence of their cosmology-dependent quantities. The same is true if the real underlying values of the cosmological parameters, for example, H₀, vary from region to region. This spatial H₀ variation in that case only has to extend up to the redshift range of our samples, as the result of a local, unknown effect, while a convergence to isotropy at larger scales would be perfectly consistent with our results. We proceed to search for anisotropies in scaling relations that are sensitive to these two phenomena, and we attempt to quantify the BF or the needed H₀ variation to explain the observations.

6.1. The L_X − T relation

Our inference of L_X would be strongly affected by either BFs or a cosmological anisotropy through the assumed luminosity distance (). At the same time, our T measurement would remain relatively unchanged. This allows us to predict the L_X values of the clusters across the sky based on their T and the globally calibrated L_X − T and attribute any directionally systematic deviations to BFs or H₀ anisotropies.

6.1.1. Our sample

The anisotropies of the L_X − T relation are extensively discussed in M20. In this section we update our results based on the new statistical methods we follow, whose differences with M20 are described in Sect. 3.7. Based on the observed scatter and the number of available clusters, we consider θ = 75° cones to scan the sky. The maximum anisotropy is detected toward (l, b)=(274°±43°, −9°±32°) (120 clusters), deviating from the rest of the sky by 19 ± 7% at a 2.8σ level. The decreased statistical significance of the anisotropy compared to the M20 results (3.64σ there) is due to the more conservative parameter uncertainty evaluation, which is now based on bootstrap resampling with marginalization over the slope, and the weaker statistical weighting of the clusters close to the center of the considered cones. Also, we do not compare the two most extreme regions with opposite behavior, but only the most extreme region with the rest of the sky. The A_LT/A_LT, all and the sigma maps are given in Fig. 4. As expected, they are not significantly different than those obtained in M20. The most noticeable difference is the lack of the particularly bright region close to the Galactic center. Now the brightest region is located ∼150° away from the faintest, indicating an almost dipolar anisotropy.

Cosmological anisotropies and BFs. We assume that the cause of the observed tension is an anisotropic H₀ value in a Universe without any BFs. We would need H₀ = 66.0 ± 1.7 km s⁻¹ Mpc⁻¹ toward (l, b)∼(274°, −9°), and H₀ = 72.1 ± 1.4 km s⁻¹ Mpc⁻¹ for the rest of the sky. The respective Hubble diagrams are compared in the top panel of Fig. 5.

We now assume that a BF motion is the sole origin of the observed anisotropies, where H₀ is isotropic. Based on the MR method, we find a BF of u_BF = 980 ± 300 km s⁻¹ toward (l, b)=(315°±34°, −10°±20°) for the full sample. This result is dominated however by lower z clusters, since there are only 26 clusters with z > 0.2. For the redshift range z ∈ [0, 0.06], we obtain u_BF = 1100 ± 410 km s⁻¹ toward (l, b)=(318°±37°, −5°±23°). We retrieve a similar BF for clusters within 270 Mpc as for the full sample. Both the direction and amplitude of the BF, as well as its statistical significance, stay within the uncertainties as we consider iteratively larger volumes. The detailed results are given in Table 4. For the concentric redshift bins z ∈ [0.06, 0.12] and z ∈ [0.12, 0.3], we obtain u_BF = 1170 ± 400 km s⁻¹ toward (l, b)=(262°±52° , + 2°±26°), and u_BF = 1040 ± 570 km s⁻¹ toward (l, b)=(253°±60°, −18°±31°),  respectively. The BF direction is consistent within 1σ between all redshift shells, not showing any convergence to zero even at ≳500 Mpc. The limited number of clusters beyond ∼500 Mpc poses a challenge however for the precise pinpointing of the BF at larger scales. The results are in tension with ΛCDM, which predicts much smaller BFs at scales of ≳200 Mpc (e.g., Li et al. 2012; Carrick et al. 2015; Qin et al. 2019).

Using the MA method, we find a BF of u_BF = 600 ± 260 km s⁻¹ toward (l, b)=(298°±25°, −21°±18°) for the full sample. When this BF is applied to our data, the L_X − T relation is consistent with isotropy within 1.4σ based on the usual sky scanning. Unfortunately, for the MA method we need to consider broader redshift bins than with the MR method to have enough data available per sky patch. For the redshift bins z ∈ [0, 0.09] we find u_BF = 690 ± 300 km s⁻¹ toward (l, b)=(268°±31°, −5°±23°), while iteratively increasing the cosmic volume does not significantly affect this result. For z > 0.09, while the anisotropy toward (l, b)∼(270°, −25°) persists at a ∼2σ level, the search for a BF is inconclusive owing to the limited number of clusters, which leads to large uncertainties. To get an idea of the possible BF signal in the MA method purely at larger scales, we can exclude local clusters (≤300 Mpc, z < 0.067)⁸ and only consider the 170 clusters at larger distances (z > 0.067). For these, we obtain u_BF = 620 ± 310 km s⁻¹ toward (l, b)=(293°±39°, −12°±27°). Although there is some overlap with the z ∈ [0, 0.09] results, these findings serve as a hint for BFs extending to larger scales.

Both methods reveal large BFs and there are no signs of fading when larger volumes or different redshift bins are considered. The direction appears to be relatively consistent between both methods, well within the 1σ uncertainties. The main difference is that the MA method returns a smaller amplitude for the BF, which is consistent within 1σ though with the MR method. Finally, the L_X − T anisotropies are not subject to a specific redshift bin but consistently extend throughout the z range.

6.1.2. Joint analysis of the L_X − T relation for our sample and ACC

As shown in M20 (and in Appendix B.2), ACC shows a similar anisotropic behavior as our sample, even though it is completely independent of the latter. We perform a joint likelihood analysis of the two independent samples using the 481 individual clusters included in the samples. We express the apparent anisotropies in terms of H₀, similar to Fig. 23 of M20. The results are plotted in Fig. 4. The free parameter H₀ seems to vary within ∼65−76 km s⁻¹ Mpc⁻¹. The most anisotropic region is found toward (l, b)=(284°, −4°). The best fit of this region is H₀ = 66.2 ± 1.6 km s⁻¹ Mpc⁻¹, while for the rest of the sky we gets H₀ = 72.7 ± 1.5 km s⁻¹ Mpc⁻¹. While the posterior H₀ range is identical to that found in M20, the significance of the anisotropies do not exceed 3σ, owing to the more conservative methodology followed, along with the non-use of XCS-DR1, in this work. The peak anisotropy is close to the result from our sample since it dominates the joint fit, which is also revealed because the σ level compared to our sample alone increases only slightly. There is a second, weaker peak on the maximum anisotropy position of ACC.

In terms of a BF motion, there is no meaningful way to combine the two independent data sets in an analytical way similar to the H₀ analysis, since any BF has meaning only within a certain z range. The redshift distribution of the two samples differ, however, and in every given redshift shell, one data set dominates over the other.

6.2. The Y_SZ − T relation

The anisotropies of the Y_SZ − T relation are presented in this work for the first time. The quantity Y_SZ strongly depends on the angular diameter distance (), which is affected by BFs and possible spatial changes of the cosmological parameters. At the same time, T is independent of these effects and thus the same reasoning as in the L_X − T relation applies. There are three advantages to the Y_SZ − T relation. Firstly, the scatter is clearly smaller than the L_X − T relation allowing for more precise constraints. Secondly, Y_SZ is unaffected by absorption issues, while T (determined using spectra with photon energies of > 0.7 keV) only has a weak dependence on uncalibrated absorption effects with much less severe effects than L_X. Thus, in practice no Y_SZ − T anisotropies can occur from unaccounted Galactic absorption issues. Finally, owing to the applied redshift evolution of Y_SZ − T, the latter is more sensitive to BFs than L_X − T.

6.2.1. Our sample

As a consequence of the low observed scatter and the 260 clusters with quality Y_SZ measurements, we consider θ = 60° cones to scan the sky. The maximum anisotropy is detected toward (l, b)=(268°±34°, −16°±29°) (57 clusters), deviating from the rest of the sky by 27 ± 7% at a 4.1σ level. The direction agrees remarkably with the results from the L_X − T relation. The amplitude of the anisotropy is larger, owing to the narrower cones, which limit the contamination of unaffected clusters in each cone. The statistical significance of the anisotropies is also considerably larger as a result of the smaller Y_SZ − T scatter and the narrower scanning cones.

This result strongly demonstrates that the observed L_X − T anisotropies are not due to X-ray absorption issues. The A_YT/A_YT, all and the sigma maps are shown in Fig. 6. We should note that there is a mild correlation between the L_X and Y_SZ scatter compared to T due to the physical state of galaxy clusters. If the origin of the anisotropies was sample-related (e.g., a surprisingly strong archival bias), we would expect to see similar anisotropies in L_X − T and Y_SZ − T. However, this would still not explain the large statistical significance of the Y_SZ − T anisotropies or the fact that we see a similar effect in ACC. Nonetheless, we further explore this possibility later in the paper (Sects. 7 and 8.6), and confirm that this is not the reason for the agreement of the two relations. Finally, if we use the Y_SZ values from PSZ2 instead, we obtain similar anisotropic results (Appendix D.3).

Cosmological anisotropies and BFs. We now investigate the necessary H₀ variations to fully explain the observed anisotropies. We would need H₀ = 62.3 ± 2.2 km s⁻¹ Mpc⁻¹ toward (l, b)∼(268°, −16°) and H₀ = 72.1 ± 0.9 km s⁻¹ Mpc⁻¹ for the rest of the sky. The result is consistent within 1.4σ with the H₀ value obtained from the joint L_X − T analysis. The corresponding Hubble diagram is shown in the bottom panel of Fig. 5.

The observed anisotropies cannot be caused by a local spatial variation of Ω_m. The Y_SZ − T relation is rather insensitive to such changes because of the low z of the clusters and the opposite effect that an Ω_m change would have on the Y_SZ value, and in the redshift evolution of the relation. However, the anisotropies appear to be even stronger now compared to the L_X − T relation.

Assuming the apparent anisotropies are due to a BF affecting the entire sample, using the MR method we find u_BF = 950 ± 340 km s⁻¹ toward (l, b)=(263°±39°, −22°±20°). For the redshift bin z ∈ [0, 0.07], we obtain u_BF = 1060 ± 390 km s⁻¹ toward (l, b)=(254°±42°, −17°±19°), similar to the results from the full sample. By gradually expanding the redshift range, the best-fit results remain within the uncertainties without an amplitude decay. For the z ∈ [0.07, 0.12] bin, we obtain u_BF = 840 ± 490 toward (l, b)=(312°±61°, −34°±31°). The BF drifts by ∼50°, but remains within the (large) uncertainties and within the marginal anisotropic region. The amplitude is also slightly decreased, but still well above the ΛCDM prediction for such scales. For the z > 0.12 clusters we find u_BF = 1110 ± 670 km s⁻¹ toward (l, b)=(321°±103°, −42°±45°). There is some indication that the BF persists toward a similar direction up to scales larger than ∼500 Mpc, however the poor constraining power of the z > 0.12 subsample does not allow for robust conclusions.

From the MA method, we obtain u_BF = 960 ± 290 km s⁻¹ toward (l, b)=(267°±22°, −28°±16°), which is completely consistent with the MR method. Within z ≤ 0.09, the constrained BF is u_BF = 1200 ± 350 km s⁻¹ toward (l, b)=(254°±22°, −18°±16°). For z > 0.09, the amplitude is reduced (u_BF = 720 ± 380 km s⁻¹), while its direction (l, b)=(242°±87°, −13°±22°) is rather uncertain (pointing however toward a similar sky patch).

We note that the Y_SZ − T anisotropies reveal similar BFs than the L_X − T case. In both cases, the MR and MA methods agree on the full sample and the low z regime, while the results are uncertain for higher redshifts. The scale out to which the apparent BF extends and its amplitude by far exceed the ΛCDM expectations.

6.2.2. Joint analysis of the Y_SZ − T relation for our sample and ACC

The ACC sample shows a very similar anisotropic Y_SZ − T behavior to our sample (Appendix B.3). We combine the two independent samples and their 376 different clusters to jointly constrain the apparent H₀ anisotropies with the Y_SZ − T relation. The H₀ spatial variation with the sigma maps are given in the bottom panel of Fig. 6. The combined maximum anisotropy direction is found toward (l, b)=(276°±26°, −14°±20°), with H₀ = 63.4 ± 2.5 km s⁻¹ Mpc⁻¹, at a 4.3σ tension with the rest of the sky (13 ± 3%). The statistical significance is slightly increased compared to our sample alone, while the direction is mostly determined by our sample. The obtained anisotropies remarkably agree with the joint L_X − T results. Finally, a dipole form of the anisotropy is apparent in the sigma maps.

6.3. The L_BCG − T relation

The final scaling relation that can potentially trace cosmological anisotropies and BFs is the L_BCG − T relation. Absorption effects are rather irrelevant for this relation, as explained in Sect. 5.2. Furthermore, the latter strongly depends on the luminosity distance (). The large scatter of the relation and the fewer number of clusters compared to L_X − T and Y_SZ − T constitute the main disadvantages of L_BCG − T . As a result of these limitations, θ = 90° scanning cones are considered to include the maximum possible number of clusters per cone.

Despite the aforementioned disadvantages, the L_BCG − T relation can offer additional insights on the observed anisotropies. A 1.9σ anisotropy is detected toward (l, b)=(257°±55°, −12°±38°), where the BCGs appear to be 19 ± 10% fainter than the rest of the sky. The normalization anisotropy map is shown in the bottom panel of Fig. D.6. This mild tension does not provide sufficient statistical evidence for a deviation of isotropy. However, the agreement of the direction with the L_X − T and Y_SZ − T results offers additional confirmation of the existence of the physical phenomenon causing this.

Cosmological anisotropies and BFs. In terms of H₀, we obtain H₀ = 66.2 ± 2.5 km s⁻¹ Mpc⁻¹ for that direction and H₀ = 72.6 ± 2.3 km s⁻¹ Mpc⁻¹ for the opposite hemisphere. The H₀ values are consistent with the other scaling relations for similar directions.

To search for any BFs, we can only consider the full sample, owing to the restricted 0.03 < z < 0.15 range, and the few available clusters. For the MR method, we obtain u_BF = 580 ± 370 km s⁻¹ toward (l, b)=(293°±50° , + 2°±29°). With the MA method, we find similar results as well, shown in detail in Table 4. Even though the direction is rather uncertain, the BF results are consistent with those found in L_X − T and Y_SZ − T.

6.4. Combined anisotropies of L_X − T, Y_SZ − T, and L_BCG − T

We now combine the information from the L_X − T, Y_SZ − T, and L_BCG − T anisotropies, for our sample and ACC, into one single H₀ anisotropy map, ignoring the peculiar velocities of clusters within the CMB reference frame. The joint analysis procedure is the same as before, where the different H₀ posterior likelihoods for every region from every scaling relation are combined.

6.4.1. Limitations of joint analysis

This method shows certain limitations that need to be kept in mind before we present the joint results. Firstly, it assumes that the underlying effect causing the anisotropies does not affect T, but only L_X, Y_SZ, and L_BCG. This is of course true for cosmological anisotropies and BFs. This also requires that there are no unaccounted X-ray absorption effects that bias the T measurement, which should be the case according to the analysis presented in Sect. 5. If a strong temperature outlier exists due to chance, it would not significantly affect the results by propagating to all scaling relations. The effect of the outlier would be buried under the average behavior of the rest of the clusters in its region, if no reason for an anisotropy exists. The posterior H₀ likelihood from the entire region is then combined with the same results from the other scaling relations.

Secondly, a moderately correlated scatter exists between L_X and Y_SZ, which relates to the physical state of galaxy clusters. We discuss this in detail in Sect. 8.6. When the average properties of clusters are similar from sky region to sky region (and any anisotropies arise purely from cosmological effects), this correlated scatter is not expected to bias our combined results. If, however, a strongly inhomogeneous spatial distribution of cluster populations (e.g., CC and non-CC clusters) exist, then the correlation between L_X and Y_SZ artificially boosts the statistical significance of the observed anisotropy. Nevertheless, we account for this effect in the MC simulations later on (Sect. 7) and show that the obtained statistical significance of the anisotropies remains rather unchanged.

Thirdly, if our methodology suffers from a systematic bias in the estimation of the statistical significance, this would affect all three scaling relations. Since their results are combined, this bias would be amplified in the final, joint anisotropy estimation. However, in Sects. 7 and 8.3 we quantify this bias and take it into account, showing that it has no significant effect on the final conclusions. Bearing in mind the above, we proceed to combine all the available information together to construct the final map for the apparent H₀ spatial variation.

6.4.2. Apparent H₀ anisotropy from joint analysis

When L_X − T, Y_SZ − T, and L_BCG − T results for both our sample and ACC are combined, we obtain H₀ = 66.5 ± 1.0 km s⁻¹ Mpc⁻¹ toward (l, b)=(273°±40°, −11°±27°), and H₀ = 72.8 ± 0.6 km s⁻¹ Mpc⁻¹ for the rest of the sky. There is a 5.4σ tension between the two values, demonstrating the high statistical significance of the detected anisotropies.

The small H₀ uncertainties are not surprising if we consider that five nearly independent results (that generally agree) have been combined into one map and as such, there is plethora of information for every sky patch. Moreover, the normalization of the scaling relations (which relates to the calibration of true cluster distances) is assumed to be perfectly known, as discussed in Sect. 3.1. This removes a dominant source of H₀ uncertainty from which probes that attempt to put absolute constraints on H₀ suffer. However, it is not needed for the relative spatial differences to which we are interested in. The joint H₀ map together with the statistical significance map are shown in Fig. 7. The observed anisotropy forms a dipole.

7. Comparison with isotropic Monte Carlo simulations

The well-established statistical methods used up to now provide reliable parameter uncertainties and a robust estimation of the rareness of the observed anisotropies. Nonetheless, biases we are not aware of, can still be present in our analysis and lead to overestimated anisotropy signals. The same is true for cosmic variance. To further investigate this, we need to apply our analysis to simulated isotropic MC simulations. For the L_X − T, Y_SZ − T, and L_BCG − T scaling relations, we create 10 000 simulated isotropic samples similar to the real samples. Analyzing these samples with the same procedures as in the main analysis, we test the frequency with which anisotropies equal or larger than those observed in the real data are retrieved. We also calculate the frequency with which the directions of the anisotropies for different scaling relations are randomly found to be as close as in the real data.

7.1. Constructing the isotropic simulated samples

To build every isotropic MC realization, we create the same number of clusters as in the real sample (for both our sample and ACC). We start from the L_X − T relation. We keep the coordinates, redshifts, and T (together with σ_log T) of the simulated clusters fixed to the real values. This is done to incorporate all the possible effects that could create anisotropies, including the spatial distribution of the real clusters. Based on the best-fit scaling relation for the full, real sample, we calculate the predicted L_X value. We then add a random offset to the latter, drawn from a log-normal distribution with a standard deviation of . Since low-mass systems might exhibit a larger intrinsic scatter, we consider three different values of . We divide the sample in three equally sized subsamples according to their T value, and for every subsample we constrain σ_int, LT ⁹ given the best-fit A_LT and B_LT for the entire sample. Based on the T value of each simulated cluster, we then use the respective σ_int, LT to draw the random offset of the L_X value. We ensure that the posterior distribution of the best-fit values of the simulated samples follow the input values for every parameter.

After that, we need to simulate Y_SZ for every cluster with S/N > 2 by taking into account the correlated scatter of Y_SZ and L_X with T. We first predict Y_SZ for every cluster based on the best-fit Y_SZ − T relation and T. Then, we further add the expected scatter of Y_SZ based on the observed best-fit correlation in Sect. 8.6 and the random, simulated L_X scatter from before. On top of that, further noise is added on the Y_SZ value based on the scatter of this correlation, which is drawn from a log-normal distribution combining the statistical uncertainties of the observed L_X and Y_SZ residuals and the intrinsic scatter. We constrain the best-fit A_YT, B_YT, and σ_intr, YT for all 10 000 samples and confirm their distribution follow the values from the real samples.

A simulated, isotropic L_BCG value is also drawn for the 196 clusters with 0.03 < z < 0.15 in the same way as for L_X − T. The infrared BCG luminosity does not show any correlation in its scatter with L_X and Y_SZ, and therefore no extra procedures are needed. We finally repeat the L_bol − T and Y_SZ − T simulations for ACC as well.

7.2. Results

We firstly explore our sample and ACC separately for the available scaling relations. At the end, we combine the results to estimate the overall probability that our entire findings are due to chance. The histograms of the maximum anisotropy detected in the MC samples for every scaling relation, together with the results from the real samples, are shown in Fig. 8.

7.2.1. Our sample

For the L_X − T relation, we find that, using our methodology, there is a 1.6% probability (p = 0.016) to observe a ≥2.8σ anisotropy within an isotropic Universe. The statistical significance of the two methods does not vary significantly.

For the Y_SZ − T relation, the chance to randomly observe a ≥4.1σ anisotropy in an isotropic Universe is p = 0.011. This is significantly more probable than implied by standard analysis and demonstrates how apparent anisotropies might be present as a consequence of noise and cosmic variance. The effect of the latter are enhanced by the narrower adopted scanning cone for Y_SZ − T than for L_X − T, which results in fewer clusters per cone. Nevertheless, the probability of such results to occur randomly is still very low and does not meaningfully change with increasing cone size. For example, when we adopt a θ = 75° scanning cone, the default analysis returns a 3.1σ anisotropy, while the MC simulations return p = 0.013, which is in closer agreement than before. The effects of cosmic and sample variance are sufficiently taken into account by the MC isotropic simulations. Finally, the statistical significance is in any case expected to drop slightly by increasing the cone width, since more clusters that are less affected by the anisotropy are included per cone.

If we consider the probability that the observed L_X − T and Y_SZ − T anisotropies emerge simultaneously for a simulated sample, then the probability dramatically drops to p = 0.001. On top of that, we need to consider the direction agreement of the anisotropies of the real samples, which is within 10°. Out of the 10 000 isotropic MC samples, only one exhibits this level of anisotropy for both scaling relations and within such a narrow cone (p = 10⁻⁴). It is important to stress again that the correlation between the L_X and Y_SZ values has been taken into account as explained before.

For the L_BCG − T relation, a ≥1.9σ anisotropy is detected for 42% of the samples, mostly owing to the large scatter. Although this result alone is consistent with isotropy, when combined with the L_X − T and Y_SZ − T results, it returns p = 4.2 × 10⁻⁴, without accounting for the three similar anisotropy directions. To account for that also, we first find the maximum angular separation θ among the most anisotropic directions of the three scaling relations. We then see how often all three anisotropy directions lie within θ. This is done for all simulated samples. We then multiply this probability with the probability that comes purely from the amplitude of the anisotropies (i.e., p = 4.2 × 10⁻⁴). Thus, we obtain p = 3.5 × 10⁻⁶. To sum up, from our sample alone, there is a 1 in ∼286 000 probability that an observer would obtain our results owing to statistical noise, cosmic and sample variance, or scatter correlation between L_X, Y_SZ, and L_BCG.

7.2.2. Joint probability

Before we assess the overall probability of our results considering ACC as well, we wish to test ACC alone. For the L_bol − T relation, the observed anisotropy amplitude is detected in the isotropic samples with p = 0.26; this anisotropy amplitude is much less rare than implied by the default methods and is consistent with isotropy. The observed behavior of the Y_SZ − T relation yields p = 0.055 alone, which is more statistically significant than L_bol − T. The joint probability to observe such anisotropies simultaneously in a sample is p = 0.028. If they were also to be separated by only < 8° as in the real data, then the probability significantly drops to p = 0.003. There are strong indications of an existing anisotropy from the ACC alone, although not enough to completely exclude the scenario of a random event.

Since the two samples are completely independent, the overall probability of our results is given by multiplying the distinct probabilities from the two samples. This yields p = 1.2 × 10⁻⁸. Furthermore, we need to account for the agreement in the anisotropy direction between our sample and ACC. The maximum angular separation between the results of the two samples, for any two scaling relations, is 50°. Since the two samples are independent, we would expect a uniform distribution in the angular separation of their simulated maximum anisotropy directions. Thus, a ≤50° separation would have a ∼28% chance to occur.

Finally, combining all the available information, we find that our results are practically impossible to occur randomly within an isotropic Universe without an underlying physical cause, since the probability for this to happen is p = 3.4 × 10^{−9 10}.

8. Possible systematics

Several possible biases could jeopardize the interpretation of the observed cluster anisotropies. A large number of them, including numerous X-ray and sample-related issues were tested and discussed for the L_X − T relation in M20. Here we explore some additional effects that might undermine the significance of our results.

8.1. Cool-core and morphologically relaxed clusters

Cool-core clusters have a strong central peak in their surface brightness profile, and they are known to be intrinsically brighter in X-rays than NCC clusters (e.g., Mittal et al. 2011). This bias propagates to scaling relations when the two physical quantities are not similarly affected by the (N)CC nature of clusters. Previous studies have found differences for such scaling relations between CC and NCC, or morphologically relaxed (i.e., regular) and disturbed clusters (e.g., Maughan 2007; Pratt et al. 2009; Zhang et al. 2017; Maughan et al. 2012; Bharadwaj et al. 2015; Lovisari et al. 2020). Environmental effects have been also found to mildly correlate with some cluster properties (M18, Manolopoulou et al. 2021).

For classifying clusters as CC or NCC, we need a robust proxy of the dynamical state of the clusters, such as the central cooling time, shape of the surface brightness profile, and concentration parameter. Future work will soon provide such information for eeHIFLUGCS, which will allow for a more precise calibration of our scaling relations. For now, we use the offset between the X-ray peak and the BCG position (XBO) to categorize clusters as morphologically relaxed or disturbed. This categorization is not the same as CC and NCC, however there is a rather strong correlation between them. Specifically, the XBO has been shown to approximately correlate with the existence of a CC in the center of the cluster (e.g., Hudson et al. 2010; Zitrin et al. 2012; Rossetti et al. 2016; Lopes et al. 2018). We consider the clusters with XBO < 0.01 R₅₀₀ as morphologically relaxed and possibly CC. We also consider those with XBO > 0.08 R₅₀₀ as disturbed, possibly NCC clusters. Each of these subsamples constitutes ∼30% of our sample.

8.1.1. Spatial distribution of relaxed and disturbed clusters

For any anisotropy study, this bias is only relevant if the spatial distribution of morphologically relaxed and disturbed clusters is not relatively uniform. As a consequence of the homogeneous selection of our sample however, we see only mild, random variations in the fraction of such clusters across the sky (Fig. 9). There is a small excess of relaxed clusters toward (l, b)∼(60° , + 30°). This mild distribution imbalance strongly correlates with the directional behavior of the R scaling relations (see Appendix C.1 for more details). Finally, it is important to stress that the (l, b)∼(280°, −15°) region shows an average behavior in its cluster population, and hence no relevant bias is expected there (as is further shown later).

8.1.2. Lack of bias in L_X − T, Y_SZ − T, and L_BCG − T

Relaxed and disturbed clusters do not show any meaningful difference in their Y_SZ − T and L_BCG − T normalization, as shown in Fig. 10. Therefore, the possible CC bias is irrelevant for these two scaling relations. For L_X − T, relaxed clusters appear 26 ± 10% brighter than the disturbed clusters (right panel of Fig. 10), which is generally expected. Of course, all these clusters are not found in only one region, but they are distributed sparsely across the sky. Thus, their effects are hardly detectable over the possible cosmological effects. Even more importantly, there is only an average number of relaxed and disturbed clusters toward the L_X − T anisotropy region, hence no considerable bias is expected.

This was shown in M18, where using core-excised L_X did not have any effect on the L_X − T anisotropies of the HIFLUGCS sample (Reiprich & Böhringer 2002). Excluding the possibly relaxed clusters and considering only clusters within superclusters, which are more likely to be disturbed, did not affect the anisotropies in M18 either.

To provide further evidence that the dynamical state of clusters does not have an effect on the detected L_X − T anisotropies, we create 10⁵ randomly drawn bootstrap subsamples (same process as in M20), independent of direction. We investigate the correlation between the A_LT and the median XBO for every subsample. Figure 11 shows that there is only a very weak anticorrelation between the median XBO of the sample and A_LT (Pearson’s correlation coefficient r = −0.11). This simply shows that it is extremely improbable for enough relaxed clusters to be included in a subsample to make a noticeable difference in A_LT.

From all the above, it is evident that an inhomogeneous distribution of morphologically relaxed and disturbed clusters is not the reason behind the observed anisotropies.

8.2. Malmquist and Eddington biases

8.2.1. Malmquist bias

Flux-limited samples such as eeHIFLUGCS suffer from the so-called Malmquist bias (MB). Clusters that are intrinsically brighter than the others, are overrepresented in the sample, especially close to the flux limit. If not taken into account, this results in a biased estimation of the scaling parameters of the true, underlying cluster population (e.g., Hudson et al. 2010; Mittal et al. 2011; Eckert et al. 2011). This bias is expected to be stronger for scaling relations including L_X since the selection of the sample was conducted based on the X-ray flux of the clusters. At the same time, the effects of the MB should be relatively weaker for the other scaling relations. We wish to assess if the MB could affect the statistical significance of the detected anisotropies. We focus on the L_X − T relation since this is probably the most affected scaling relation and the one we used in M20.

If the MB influences all regions equally, then there is no effect on the detected anisotropies. This is a naive expectation since the MB is directly linked to the scatter of the L_X − T relation; for decreasing scatter, the effects of the MB are also weakened. The scatter of a relation is an intrinsic cluster property because it mainly depends on the physical state of these objects. In an isotropic Universe, there is no obvious reason why such an intrinsic cluster property would spatially vary. A risk factor that would make the MB relevant is if the cluster population differs from region to region (e.g., due to archival bias). For large enough, homogeneously selected samples such as our own, significant differences are not expected, but nevertheless we try to quantify the possible bias in our anisotropy estimates.

For L_X − T, the faintest, maximum anisotropy region lies toward (l, b)∼(274°, −9°), containing 120 clusters with σ_int = 58.9 ± 4.6%. The rest of the 193 clusters are 19% brighter with a 2.8σ tension and with σ_int = 50.9 ± 3.5%¹¹. The scatter of the faint region is actually larger than the rest of the sky at a 1.4σ level. According to Vikhlinin et al. (2009), the observed scatter estimate is also the true, underlying scatter for X-ray flux-limited samples.

The MB correction to be applied in the normalization of the L_X − T for low-z, flux-limited samples such as our own is (Vikhlinin et al. 2009). Consequently, for the (l, b)∼(274°, −9°) region we have MB_corr ≈ 0.594 ± 0.048, and for the rest of the sky MB_corr ≈ 0.678 ± 0.036. Therefore, there is a 84% chance that after the proper MB correction, the faint region will become even more anisotropic compared to the rest of the sky. Based on the MB_corr uncertainties, there is only a 0.007% (4σ) probability that the full amplitude of the L_X − T anisotropy is a result of the MB.

Similar to the effects of the scatter, we might expect that larger L_X uncertainties could lead to larger normalizations and vice versa. However, the faintest region again has a slightly larger average L_X uncertainty than the rest of the sky, while the median value is similar between the two subsamples.

This result strongly suggests that the MB is not the reason behind the L_X − T anisotropies, even though there is a small chance ( < 16%) that it might lead to a slight overestimation of the statistical significance of the findings. However, considering the even larger relative anisotropy of the Y_SZ − T relation and the weaker effects that MB has there, we conclude that the latter is probably irrelevant for our analysis.

8.2.2. Eddington bias

The Eddington bias refers to the fact that at large distances, only luminous, massive clusters exceed the flux limit of a survey. Thus, they are overrepresented at large z. This leads to a different L_X, Y_SZ, and T (i.e., mass) distributions between low and high z. Low- and high-mass systems might share slightly different scaling laws. If the z distributions of the compared cluster subsamples were significantly different, this could lead to a biased comparison and artificial anisotropies. However, in M20 we showed that the z and T distributions of the (l, b)∼(280°, −15°) are very similar to the rest of the sky, as expected for a homogeneously selected sample. This is further confirmed by our physical properties test in Sect. 8.4. Subsequently, the Eddington bias is rather irrelevant to our analysis.

8.3. Zone of avoidance bias

The fact that the direction of the observed anisotropies mostly lies close to the Galactic plane gap, or the zone of avoidance (ZoA; |b ≤ 20° |) is another issue that might raise concerns about the underlying origin of these anisotropies. In Sect. 5, we demonstrated that uncalibrated X-ray absorption issues are not the reason behind L_X − T anisotropies and definitely cannot explain the Y_SZ − T anisotropies. Subsequently, the only possible bias that might lead to the artificial detection of anisotropies close to the ZoA is the applied methodology or some unaccounted archival bias toward that region.

We use two independent tests to ensure that our applied methodology is not biased toward this region. Firstly, we repeat the sky scanning for the L_X − T, Y_SZ − T, and L_BCG − T relations of our sample, this time without applying any weighting on the cluster uncertainties; these weights were based on each the distance of each cluster from the center of each cone. As a result, the scanning algorithm does not take into account the spatial distribution of clusters and does not “see” the ZoA gap. The only information the scanning algorithm reads is the number of clusters in each cone. The cones close to the most anisotropic region happen to have average numbers of objects (Fig. E.2). Even if fewer clusters were included in the cones, in M20 we showed that the best-fit normalization is completely independent of the number of clusters.

When repeating the analysis for L_X − T, we obtain an anisotropy of 2.7σ toward (l, b)=(284°±44°, −4°±33°). The result is similar to the default analysis in terms of the anisotropy direction and amplitude. In a similar manner, for the Y_SZ − T relation we obtain a 3.9σ anisotropy toward (l, b)=(262°±33°, −22°±30°), again similar to the default case. The anisotropies of the L_BCG − T anisotropies are practically not affected. These results establish that the observed anisotropies in our data are not an artifact of the ZoA effect on our analysis. The results also show that the applied weighting does not significantly alter the anisotropy results.

The second test is twofold. Firstly, we consider the 10 000 isotropic MC simulated samples used in Sect. 7, which are drawn from the same distribution as the real sample. There the positions of the clusters are kept fixed and thus the ZoA remains empty. For the same scaling relations as before, we measure how many times the maximum anisotropy is found within the ZoA. We then compare this to the random expectation based on the fraction of the full sky area the ZoA covers, namely ∼33%. We find that ∼43% of the maximum anisotropy directions lie within the ZoA. This indicates that the applied statistical (distance) weighting of the clusters during the sky scanning can introduce a small bias in the direction of the detected anisotropies. However, this bias is very small and completely disappears if the procedure is repeated with uniform weighting. The real data anisotropies however do not drift in this case. Finally, in Sect. 7 we found that 1.6% and 1.1% of the isotropic simulated samples show larger anisotropies than the real data for the L_X − T and Y_SZ − T relations, respectively. The probability that such anisotropies are found within the ZoA for an isotropic simulated sample is 1% and 0.7%, respectively.

For the next part of this test, we randomly fill in the empty ZoA area of these 10 000 MC samples with ∼130−150 simulated clusters. These clusters have the same number density as the rest of the sky and the same T and scatter distributions. We again measure the level of the anisotropies and compare this level with the case in which the ZoA is excluded. That way we wish to see if excluding the clusters within the ZoA in the real data introduces a bias in our anisotropy estimates. We find that the amplitude of the anisotropies decreases by 14 ± 8%, which is expected owing to the larger number of isotropic data. Additionally, ∼36% of the most deviate directions are located within ZoA; this value is slightly decreased compared to the case when the ZoA clusters are excluded and is consistent with the random expectation. Hence, the gap of the ZoA can weakly affect the direction of the detected anisotropies, but not enough to compromise our results.

These tests combined strongly suggest that there is no significant bias in our methodology that favors and amplifies anisotropic signals close to the ZoA. The minor changes that can be caused are much smaller than the measured direction uncertainties. Relative plots for these tests can be found in Appendix E.1. Future surveys will offer a more robust cluster detection toward that region, which will help us pinpoint the anisotropies of the clusters slightly more accurately.

8.4. Anomalous combination of cluster properties

In M20 we investigated if a single average physical property of cluster subsamples can be associated with an anomalous behavior of A_LT. In this work, we wish to find out if such a behavior of A_LT (or A_YT) could also result from a certain combination of average cluster properties. To do so, we construct 10⁶ random bootstrap cluster subsamples of random size (10−60% of the total size of the sample), as done in M20. Except for the best-fit A_LT, B_LT and σ_int, LT, 12 more average properties of each subsample are also derived¹². We express A_LT as a function of all the parameters p_i (normalized by their sample mean p_mean), and their power-law index v_i, as follows:

(10)

The 1.132 term corresponds to the best-fit value for the full sample.

To understand which combination of properties could lead to a significantly altered A_LT, we need to constrain v_i. Knowing these, we can predict the best-fit A_LT based just on the physical properties of the clusters. To do so, we perform a Markov chain Monte Carlo (MCMC) fitting to the 10⁶ bootstrap subsamples. The exact details of the fitting process together with the relevant plots can be found in Appendix E.2. A_LT shows a negligible dependence on 11 out of the 14 parameters. There is a weak anticorrelation with the best-fit B_LT (v_B = −0.40 ± 0.22), a moderate anticorrelation with mean T (v_T = −0.87 ± 0.21), and a moderate correlation with mean z (v_z = +0.82 ± 0.19). In flux-limited samples and as a result of the Eddington bias, high-z clusters are also high-T clusters. As a result, the contribution of these two terms in the predicted A_LT balances out. In general, subsamples with local, hot clusters tend to appear fainter than average and vice versa.

To determine if the observed L_X − T anisotropies are caused by the physical cluster properties of the (l, b)∼(274°, −9°), we predict its A_LT using Eq. (10) and its mean cluster properties. We obtain A_LT ∼ 1.12 and A_LT ∼ 1.14 for the rest of the sky. The two results are similar and so we confirm that the strong apparent anisotropies are not caused by a possible bias as a result of an inhomogeneous distribution of different cluster populations. This result was expected, since in M20 we showed that the average physical cluster properties are similar across the sky. Repeating the analysis for Y_SZ − T returns similar results.

8.5. Temperature calibration issues

If the measured cluster temperatures suffered from calibration issues and were biased toward higher values within the (l, b)∼(280°, −10°) region, this would create apparent anisotropies in all scaling relations that include T. This would eventually lead to an overestimation of the statistical significance of the anisotropies, when we combine all the T scaling relations.

However, this scenario is highly unlikely for several reasons. Firstly, if T was overestimated toward (l, b)∼(280°, −10°), the same anisotropies should appear in the R − T scaling relation as well, but they do not (Appendix C.1). Secondly, there is no obvious reason why the ACC sample would show similar anisotropies to our sample, when its temperatures were measured two decades ago with a different X-ray telescope and data analysis process. Thirdly, clusters in that region were not observed within a certain, narrow time interval, but throughout many years. As such, any calibration issues could not have affected only these clusters, but it should have affected clusters from the entire sky. In that case, strong anisotropies would not have been observed because several clusters from one region are needed to show a systematic behavior in order for a statistically significant signal to appear.

One possible effect that could bias the directionality of T scaling relations is the fraction between clusters with measured T from Chandra (237 in total) or XMM-Newton (76). As shown in Schellenberger et al. (2015) and further confirmed in M20, T values measured with the two telescopes slightly (and systematically) differ. These values follow a clear, low-scatter relation however, which we used to convert XMM-Newton temperatures to the equivalent Chandra temperatures. Therefore, we should not expect any dependence of the results on the fraction of Chandra and XMM-Newton clusters. This is confirmed by the analysis in Sect. 8.4, where this fraction does not strongly correlate with the best-fit A_LT. Nevertheless, we test if the spatial variation of this fraction correlates with the observed anisotropies in L_X − T, Y_SZ − T, and L_BCG − T. We find that the (l, b)∼(280°, −10°) does not show any strong imbalance (Fig. E.4). There is an excess of Chandra clusters toward (l, b)∼(150° , + 40°), while an excess of XMM-Newton clusters is found toward (l, b)∼(334°, −50°). Moreover, when only Chandra clusters were used in M20, the same L_X − T anisotropies were observed.

Based on all the above, we can safely conclude that the observed anisotropies are not due to T measurement calibration issues. Finally, eRASS will provide cluster observations from a single instrument, analyzed within a strict time interval, and thus avoid such possible systematics.

8.6. Correlation of L_X and Y_SZ scatter with T

Past studies (i.e., Nagarajan et al. 2019, and references therein) have reported a positive correlation between the scatter of L_X and Y_SZ at fixed mass (or T equivalently). In this work we test this correlation with much larger samples as well as the correlation of the scatter between these two quantities and L_BCG. Considering the 260 clusters with both Y_SZ and L_X measurements, we confirm that their residuals with respect to T are correlated and have a Pearson’s correlation coefficient of r_corr = 0.668 ± 0.089. This is shown in Fig. 12. The best-fit line has a slope of ∼0.63, while the total scatter of the relation is ∼0.16 dex.

We should note that this correlation has no effect on the statistical significance values extracted for every scaling relation individually, although it can affect the combined statistical significance of the anisotropies as obtained in Sect. 6.4. In practice, this means that both Y_SZ − T and L_X − T would exhibit similarly biased anisotropies in the case of an inhomogeneous distribution of CC clusters or mergers for instance. Their results then could not be considered independent and combined as in Sect. 6.4. Our previous analysis showed that such a CC (or merger) bias is very unlikely to cause the observed anisotropies, since no correlation was found between the A_LT and A_YT and the average dynamical state of the clusters of a subsample. We fully explore the effect that this correlation has in our results, however, in the next section in which MC simulations are carried out.

Similar results are obtained for the ACC sample for both the correlation coefficient and the slope of the relation. Finally, the scatter of L_BCG with respect to T does not show any meaningful correlation with the scatter of L_X or Y_SZ.

8.7. Selection biases

As extensively argued in M20, selection biases are not expected to significantly affect our anisotropy results. We attempt to constrain the relative deviations between sky regions and thus we use the overall best-fit scaling relations as a frame of reference. The bias-corrected scaling relations are not a necessity for that as long as the cluster properties across the sky do not vary significantly. This was shown to be true in both M20 and in this work. Our results however are still directly comparable with similar literature results since many studies focus on the observed scaling relations. Nonetheless, we performed several tests to evaluate how some selection biases can affect the statistical significance of the anisotropies. The results from Sects. 8.1 and 8.2 support the notion that the observed anisotropies are not overestimated by selection biases. More tests are performed in Appendix D, yielding the same conclusion. The bias-corrected scaling relations will be presented in future work. Overall, the results of Sect. 8 further suggest that the observed anisotropies do not emerge as a result of unaccounted biases.

9. Discussion and conclusions

The study of the (an)isotropic behavior of galaxy cluster scaling relations opens a very promising window for testing the cosmological principle. Galaxy cluster samples are much more spatially uniform than SNIa samples, and they extent to much larger scales than the galaxy samples that are used to study BFs. The samples can also be observed in several wavelengths and allow us to measure many different properties of them. This provides us with plenty of nearly independent tests that trace cosmological phenomena and that can be carried out with the same sample of objects. In this work, we were able to perform numerous such tests that provided us with valuable insights into the apparent validity of the cosmological principle in the nearby Universe.

9.1. Known systematics do not explain the observed anisotropies

The L_X − T anisotropy toward (l, b)∼(280°, −15°) that was originally observed in M20 could in principle be due to uncalibrated X-ray effects, such as excess X-ray absorption toward this sky direction. With our new data and analysis strategy, especially from the L_X − Y_SZ relation, which shows no signs of extra X-ray absorption toward the region in question, we can now decisively exclude this possibility. The lack of excess X-ray absorption effects toward that region is further confirmed by the analysis of ACC and by the L_X − L_BCG relation of our sample. The comparison between the X-ray-determined N_Htot and the W13-based N_Htot also does not indicate such an effect. Furthermore, the same anisotropies persist in scaling relations that are insensitive to unaccounted X-ray absorption, such as Y_SZ − T and L_BCG − T. Thus, this is not the reason for the observed anisotropy.

Another possible solution to the observed anisotropies was the existence of biases related to the sample or the followed methodology, such as selection and archival biases, the effect of the ZoA gap, biases introduced by cluster morphology, temperature calibration issues, and correlations between cluster properties. Our current results clearly exclude these possibilities as well. We showed that the MB is more likely to underestimate the observed anisotropies rather than explain them. Also, it would not have such a strong effect for Y_SZ − T, L_BCG − T, and the ACC sample, which is not flux-limited. It was also shown that the fraction of morphologically relaxed and disturbed clusters does not spatially vary strongly enough to have any effect in the apparent anisotropy of the L_X − T relation. Even if it did, this would have no strong effect on the Y_SZ − T anisotropies or the results from the other samples. By using isotropic simulated data we further showed that the ZoA gap does not affect the amplitude of the anisotropies, while it has a very small effect in the direction of the latter. Even when we adjust our methodology to eliminate the effects of this gap, the results are the same. The isotropic MC simulations we performed verify the high (> 5σ) statistical significance of the observed anisotropies, when all the cluster scaling relations that are sensitive to cosmological phenomena are combined. In future work, anisotropic simulations, due to both H₀ anisotropies and BFs, will be employed to better understand the precision and accuracy with which our methodology can detect such phenomena.

Despite testing and rejecting the generally known biases in cluster analysis as the source of the observed anisotropies, the existence of currently unknown systematics cannot be excluded. Although it would be surprising for such a significant, direction-dependent, and yet undiscovered bias to exist, all possibilities should be thoroughly scrutinized. In case such a systematic is discovered in the future, its effects on past studies with similar cluster samples as our own should be investigated.

9.2. Cosmological phenomena behind the tension

For now, there are just two obvious remaining explanations for the persisting cluster scaling relation anisotropies. The first is a spatial variation of H₀ at a 9% level (assuming no BF motions) due to a primordial anisotropy or new physics at z ≲ 0.2 scales. The second is a BF motion of ∼900 km s⁻¹ extending to ≳500 Mpc scales (assuming H₀ is isotropic). Another way to interpret this is that the matter rest frame is not identical to the CMB rest frame within the investigated redshift range. All of these scenarios are in tension with ΛCDM and the generally adopted isotropic assumption after transitioning to the CMB rest frame. The results for every scaling relation that traces such phenomena are summarized in Tables 3 and 4. In Fig. 13, the amplitude and directions of the detected BFs are plotted for different redshift bins.

Unfortunately, the above explanations are not distinguishable with the current data sets and more high-z clusters are needed. The upcoming eRASS catalogs will provide us with thousands of such clusters. These distant objects will be of crucial importance when applying the techniques presented in this work, since we will be able to tell if, and at which scale, the cluster behavior converges to isotropy. The kSZ effect, which is caused by the peculiar motion of clusters in the CMB frame, can also be utilized to differentiate cosmological anisotropies from BFs. The limited number of clusters in our sample, together with the low redshifts and the large angular cluster sizes, did not allow us to extract useful information for this effect by stacking filtered kSZ maps of objects that presumably move toward a common direction. The reason for this lies in the large power of the primary CMB anisotropies at the scales that are relevant to our sample as well as the limited angular resolution and sensitivity provided by Planck. Next-generation SZ instruments will significantly improve over Planck’s sensitivity and angular resolution and allow for direct measurements of BFs using the kSZ effect.

9.2.1. Anisotropic Hubble expansion

In the case of an anisotropic Hubble expansion, a primordial anisotropy that extends to very large scales might be present. Such an anisotropy might correlate with the CMB dipole if the latter is not purely of kinematic origin. Many claims for the detection of such cosmological anisotropies have been made recently¹³. Fosalba & Gaztañaga (2021) find highly statistically significant anisotropies in the cosmological parameter constraints from the CMB using Planck toward a region consistent with our results. Secrest et al. (2021) also find a large-scale anisotropy on the matter distribution, rejecting cosmic isotropy and the solely kinematic interpretation of the CMB dipole at a 4.9σ level. If their result was solely caused by our local motion within the matter rest frame, this would correspond to a ∼800 km s⁻¹ velocity, similar to our BF results.

Similar large-scale anisotropies in the distribution of high-z radio and infrared sources have also been found (e.g., Tiwari et al. 2015; Colin et al. 2017; Bengaly et al. 2018; Rameez et al. 2018; Siewert et al. 2020), usually implying a ∼600−1500 km s⁻¹ motion compared to the matter rest frame. In other words, there are very strong indications that the matter rest frame differs from the CMB rest frame, which is ultimately assumed to be isotropic. If true, this would have crucial implications on the standard model of cosmology. Recently, several theoretical frameworks that accumulate such anisotropies were developed (e.g., Dąbrowski & Wagner 2020; Paliathanasis & Leon 2020; Das et al. 2021; Spallicci et al. 2021). On the other hand, results that do not show any strong evidence for departure for statistical isotropy have also been presented (e.g., Andrade et al. 2019; Bengaly et al. 2019). Future eRASS catalogs will also shed more light on this question since they will provide previously unmatched catalogs of millions of X-ray point sources, tracing the matter distribution out to very large scales.

In addition, much work has been done on testing the isotropy using SNIa. Colin et al. (2019) detect a 3.9σ anisotropy in the deceleration parameter using a maximum-likelihood analysis. In contrast, Soltis et al. (2019) find no evidence of departure from H₀ isotropy using a novel nonparametric methodology. Except for different SNIa samples, the two papers also follow different approaches in the treatment of the data. Many other studies (e.g., Chang & Lin 2015; Deng & Wei 2018; Sun & Wang 2019; Zhao et al. 2019; Salehi et al. 2020; Hu et al. 2020) also used SNIa combined with other probes, identifying an anisotropic direction in impressive agreement with our results. The statistical significance of their findings however is marginally consistent with an isotropic Universe. These authors also argue that the main problem of such SNIa studies is still the highly inhomogeneous distribution of the data across the sky, while the effect of peculiar velocity corrections has also been studied (e.g., Huterer 2020; Mohayaee et al. 2020).

Finally, other kinds of nonlocal tensions with the cosmological principle have recently been observed (Horvath et al. 2020; Shamir 2020). More novel and independent methods to constrain dark energy with galaxy clusters, and possibly test the cosmic isotropy in the future, were also developed recently (Korkidis et al. 2020; Pavlidou et al. 2020).

9.2.2. Large-scale bulk flows

The consistent observation of cluster anisotropies across several scaling relations could hint at the existence of large BFs that by far exceed the BF scales predicted by ΛCDM. In Fig. 13, a ≳600 km s⁻¹ motion is consistently detected with respect to the CMB rest frame; there are no signs of fading at ≳500 Mpc. The observed BF amplitude of spherical volumes of much larger radius might be overestimated since clusters at lower distances are expected to dominate the BF signal. However, this plays no role for the iterative redshift shells results at z > 0.12, which seem to hint that the large BF might persist farther (with a small statistical significance of ∼2σ). More high-z clusters are necessary to derive safe conclusions at these scales.

The ΛCDM model predicts negligible BF amplitudes at scales of ≳250 Mpc. Flow motions of similar, large amplitudes have been reported in the past for galaxy clusters (e.g., Lauer & Postman 1994; Hudson et al. 1999; Kashlinsky et al. 2008, 2010; Atrio-Barandela et al. 2015), while Osborne et al. (2011), Mody & Hajian (2012), and Planck Collaboration Int. XIII (2014) fail to see any large-scale BFs in the CMB kSZ data. Studies using thousands of galaxies consistently detect a BF toward the same direction as us, but with a ∼3−5 times smaller amplitude. However, it also seems to extend to scales larger than ∼300 Mpc, up to where the galaxy surveys are usable. For instance, Watkins et al. (2009) and Magoulas et al. (2016) find a ∼400 km s⁻¹ BF at ∼100−150 h⁻¹ Mpc scales toward (l, b)∼(295° , + 10°), which is moderately discrepant with the standard cosmological model. Carrick et al. (2015) and Boruah et al. (2020) find a ∼170 km s⁻¹ BF toward (l, b)∼(303° , + 3°) that extends to z > 0.067, but is not particularly inconsistent with the standard expectations. Lavaux et al. (2013) use the kSZ effect on galaxy halos and find a BF of ∼285 km s⁻¹ for 200 h⁻¹ Mpc scales. Watkins & Feldman (2015) also find a good agreement with ΛCDM at small scales, but a moderate disagreement at larger scales. Furthermore, Peery et al. (2018) find only a 2% chance for their detected BF at scales 150 h⁻¹ Mpc to occur within ΛCDM. At z < 0.07 scales. The reason behind the mild disagreement in the BF amplitude between our analysis and studies of peculiar motions of galaxies remains unknown. At larger scales, a direct comparison cannot be made yet owing to the lack of BF measurements based on galaxies.

Studies have tried to constrain the BF with SNIa data as well, which extend to much larger scales, but are highly inhomogeneous (Colin et al. 2011; Feindt et al. 2013; Appleby et al. 2015; Mathews et al. 2016; Salehi et al. 2021, and references therein). Such studies usually find similar results to the galaxy surveys, where mild BFs extend to larger-than-expected scales. However, they again stress the limitations of current SNIa samples for such studies and that conclusive answers cannot be given. Other studies have stressed that the usual Newtonian treatment of peculiar velocity fields can significantly underestimate the inferred BFs (e.g., Filippou & Tsagas 2021; Tsaprazi & Tsagas 2020). From all that, it is evident that it is still not very clear if the detected BF motions are entirely in agreement with the standard model and to which scale the BF persists.

Substantial, local matter inhomogeneities such as local voids or overdensities could contribute to giving rise to unexpected large BFs. If we are located away from the void center, this could cause the apparent expansion rate to be lower toward the more matter-dominated side. There have been claims for the existence of such large voids (or overdensities) that could create outflows out to ∼400 Mpc scales (Keenan et al. 2012; Rubart et al. 2014; Whitbourn & Shanks 2014; Shanks et al. 2019a,b; Tully et al. 2019; Böhringer et al. 2020; Kazantzidis & Perivolaropoulos 2020; Haslbauer et al. 2020, and references therein), and their effects on the measured H₀ have been studied. These scales are close to the median distance distribution of our cluster samples, and thus could affect our results. The predicted H₀ variation due to these voids is slightly lower than the H₀ variation we observe, but these values generally agree. Of course, the fact that these voids offer a possible (at least partial) explanation for the local cluster anisotropies does not alleviate the problem, since the existence of such large voids is in direct disagreement with the ΛCDM.

10. Summary

In this work, we applied a scrutinized test for the isotropy of the local Universe. We studied the anisotropy of ten galaxy cluster scaling relations, utilizing observations in X-rays, infrared, and submillimeter. Using the L_X − T, Y_SZ − T, and L_BCG − T scaling relations of eeHIFLUGCS and combining these with the completely independent ACC sample, we detected a ∼5.5σ anisotropy toward (l, b)∼(280°, −15°). This high statistical significance is further confirmed by applying our methods to isotropic Monte Carlo simulated samples. Considering the low median redshift of the cluster data (z ∼ 0.1), this direction agrees with a plethora of past studies of several probes. We robustly show that our data do not suffer from unknown X-ray absorption issues and our results do not originate from any kind of known biases. Moreover, it would be difficult for any typical bias to simultaneously explain the similar anisotropic results across nearly independent scaling relations in different wavebands. Nevertheless, more work is needed to further confirm our work is not prone to any currently unknown systematics.

If the observed anisotropies were the result of an anisotropic expansion rate, we would need a ∼9% spatial variation of H₀ to reconcile with the observations. Alternatively, we would need a BF motion of ∼900 km s⁻¹, possibly extending beyond 500 Mpc, to explain the obtained cluster anisotropies. As a consequence of the low redshift range of our samples, these two phenomena are currently inseparable for a typical observer. However, both of these scenarios are in tension with the standard assumptions in ΛCDM. Since these are currently the only available explanations for our observations, we are faced with a severe problem that needs to be solved. The future eRASS catalogs will help us understand these anisotropies better and determine if there is a scale of convergence with isotropy or if this anomaly extends to much larger scales.

---

Acknowledgments

We thank the anonymous referee for their constructive comments that helped us improve our manuscript. KM is a member of the Max-Planck International School for Astronomy and Astrophysics (IMPRS) and of the Bonn-Cologne Graduate School for Physics and Astronomy (BCGS), and thanks for their support. GS acknowledges support through NASA Chandra grant GO5-16137X. LL acknowledges financial contribution from the contracts ASI-INAF Athena 2019-27-HH.0, “Attività di Studio per la comunità scientifica di Astrofisica delle Alte Energie e Fisica Astroparticellare” (Accordo Attuativo ASI-INAF n. 2017-14-H.0), and from INAF “Call per interventi aggiuntivi a sostegno della ricerca di mainstream di INAF”. This research has made use of the NASA/IPAC Extragalactic Database (NED) which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration.

References

Appendix A: Measurements of Y_SZ, L_BCG, R, and N_H,Xray

A.1. Details on the Y_SZ measurement and comparison with the PSZ2 values

The MMFs algorithms are applied to the Planck data in steps that are analogous to those presented by A16 (Sect. 2.1). Fields with a size of 10° ×10° are extracted around the X-ray coordinates of each cluster at each of the Planck HFI bands. The Low Frequency Instrument (LFI) channels are excluded because of their low spatial resolution and lower sensitivity. The six Planck HFI maps for each cluster are then processed with a MMF, which filters and combines the maps into a single map. The HFI beams are approximated as two-dimensional Gaussians with a solid angle that is equivalent to that of the effective beams. The corresponding full width at half maximums (FWHMs) are taken from Table 12 of Planck Collaboration III (2020). The required SED of the tSZ effect was computed for both the non-relativistic case as well as with relativistic corrections (computed with SZpack; Chluba et al. 2012) as explained later. In all cases, the instrumental impact on the shape of the SED was taken into account by computing the bandpass-corrected spectra (see Eq. (13), Eq. (A1), and Table A1 in Erler et al. 2018).

Applying the filters and coadding the filtered maps yields a map of the deconvolved central Comptonization parameter in each pixel of the map. We determine the tSZ-measured center of each cluster by the position of the brightest pixel within a 15′ × 15′ box around the X-ray coordinates. Using the normalized cluster template, the value of is converted to Y_5R500 by

(A.1)

where y(θ) is the y-parameter profile of the cluster that has been normalized to unity, θ is the radial angular coordinate in the plane of the sky, and θ₅₀₀ is the apparent R₅₀₀ in arcmin. The uncertainties of θ₅₀₀ (which are not provided by MCXC) were ignored since their effect is negligible. The θ₅₀₀ values were calculated through the X-ray luminosity–mass relation of Arnaud et al. (2010) and the angular diameter distance to each cluster. The only meaningful contribution to the θ₅₀₀ uncertainties comes from the scatter of the scaling relation. However, owing to the shallow dependence of R₅₀₀ to L_X (), a fluctuation of L_X within the scatter would only cause a minor shift in θ₅₀₀ (see M20). Additionally, as shown in Fig. E.5, even shifts of 6% in R₅₀₀ would not cause a major change in Y_5R500, certainly less than the statistical uncertainties. Thus, it is reasonable to ignore any uncertainties of θ₅₀₀ for measuring Y_5R500.

A total of four different Y_5R500 values are extracted for every cluster. The first value is obtained via a MMF approach using the non-relativistic spectrum of the tSZ effect. This method is analogous to the MMF1 and MMF3 pipelines used by A16. The second value is obtained via a MMF approach using the relativistic spectrum of the tSZ effect computed using the M20 estimates for T. When the latter is not available, T is calculated through the L_X − T relation found in M20, where L_X is given by MCXC, after it has been corrected for the X-ray absorption based on the Willingale et al. (2013) values (following the same procedure as in M20). The third value is obtained via a constrained MMF approach (CMMF, E19) using the non-relativistic spectrum of the tSZ effect and the spectrum of the kSZ effect. The latter is used to remove any bias introduced by the kSZ effect. The fourth and last value is obtained identically to the third value, but using a relativistic spectrum of the tSZ effect.

The four Y_5R500 values do not differ significantly to each other ( < 7.2%), and any selection for the default Y_5R500 used in this work results in the same conclusions with only minimal numerical fluctuations. We choose as the default Y_5R500 the result of the fourth approach, where a CMMF approach is followed as described in E19, correcting for relativistic effects and attempting to remove any kSZ bias. The change in the Y_5R500 values due to the removal of the kSZ effect is usually much smaller than the 1σ uncertainties of the central values. This has nearly no effect in any BF detection from our analysis because the major contribution of BFs to Y_SZ comes from D_A and the (biased) redshift–distance conversion.

We compare the values of our default Y_5R500 with those from PSZ2 for 566 common clusters. The matching criteria were that the angular separation of the two given cluster centers should be ≤1.5° and the redshift difference Δz ≤ 0.01 between the two catalogs. Figure A.1 shows there is a linear relation between the two, which demonstrates the robustness of our method. Nevertheless, a systematic departure from the values reported by the Planck Collaboration is observed, especially in higher Y_5R500 values. This partially arises from the fact that A16 estimate the size of clusters directly from the SZ data by mapping the significance of the obtained values for Y_5R500 as a function of MMF filter templates of different cluster size, parameterized through the scale radius θ_s = θ₅₀₀/c₅₀₀, where c₅₀₀ is the concentration parameter of the cluster. However, for the vast majority of clusters, the θ_s − Y_5R500 plane is left largely unconstrained, for which reason the Planck Collaboration employs an XMM-Newton derived θ − Y_5R500 prior to break the degeneracy between the two parameters (e.g., Fig. 16 of A16). Additionally, small differences in the adopted centers might contribute to the observed scatter as well¹⁴. Although the data employed in this work are identical and our filtering is based on the same MMF technique as used by A16, differences in the reported values for Y_5R500 can therefore still appear. The comparison does not significantly change with the use of Y_5R500 values based on the other three E19 approaches we described in Appendix A.1.

A.2. Details on obtaining the near-infrared BCG luminosity L_BCG

The main criteria for a galaxy to be selected as a BCG were the following. Firstly, the galaxy must be the brightest in the r band of the optical regime or Ks and W1 in 2MASS and WISE of the NIR domain. Secondly, the galaxy should be located within a distance of R₅₀₀ from the X-ray center of the cluster. Both of these values are taken by the MCXC catalog. Thirdly, the redshift of the BCG must be within a certain range Δz of the cluster redshift. If the BCG redshift was estimated via photometric measurements, then Δz ≤ 0.02, whilst the allowed redshift difference for spectroscopic redshifts was Δz ≤ 0.01. If multiple BCG candidates were suggested for one cluster, the following criteria were considered. The BCG should possess a rather extended envelope and it should be selected by optical surveys over NIRs.

Two different corrections took place before selecting the BCG of every cluster. The galaxy magnitudes were corrected for galactic extinction, using the Schlegel maps (Schlegel et al. 1998) and assuming an extinction law (Fitzpatrick 1999) with R_V = 3.1. The photometric magnitudes were also corrected to the rest frame of individual galaxies by applying an appropriate k-correction. The k-correction code provided by Chilingarian & Zolotukhin (2011) was used for all data except for WISE because the coefficients for WISE filters are not included in this program. For WISE galaxies, we instead used the k-correction code from Blanton (Blanton & Roweis 2007).

A.3. Details on the X-ray determined N_H, Xray

To determine N_H,Xray, we used the exact same spectral fitting procedure as described in M20 when the X-ray redshift was fitted. This time the N_H,Xray is left free to vary instead. Its value was linked between the spectra from the 0 to 0.2 R₅₀₀ region and the 0.2 to 0.5 R₅₀₀ region. The redshift is kept fixed at its M20 value, while the temperature, metallicity and normalization values are left free to vary for each region separately. We used the 0.7–7.0 keV range to fit the spectra (to be consistent with the results from M20). The Asplund et al. (2009) abundance table was used for the fits.

We were able to constrain N_H,Xray for 213 out of 237 Chandra clusters and 27 out of 76 XMM-Newton clusters. When we performed the N_H,Xray − N_Htot comparison, we saw that the data residuals compared to the best-fit model behaved strongly as a function of the N_H,Xray measurement uncertainty. Clusters with large statistical uncertainties were significantly downscattered, increasing the scatter of the relation. These were mostly low N_H,Xray clusters. The 0.7 keV cut removes most of the power for detecting N_Htot (which shows up more strongly at lower energies). So when the N_Htot effects are already small, a constraint is nearly impossible. To sufficiently restrict the noise, we excluded the 35% of the sample with the largest (average) uncertainties, leaving us with 156 clusters with an N_H,Xray measurement uncertainty of < 25%. Although these cuts were chosen arbitrarily, they took place strictly before any anisotropy test to avoid any confirmation bias.

Appendix B: ACC results

In this section, we only use the ACC sample to study the anisotropic behavior of the L_bol − Y_SZ, L_bol − T, and Y_SZ − T relations. Owing to the much smaller number of clusters, we consider a θ = 90° scanning cones for all scaling relations, so at least 35 clusters lie within each cone.

B.1. The L_X − Y_SZ relation

We study the L_bol − Y_SZ relation for ACC alone. The results are shown in Fig. B.1.

The maximum anisotropy is found toward (l, b)=(24°±42° , + 22°±21°) at a 2.3σ level. These 59 clusters are 26 ± 11% dimmer than the rest. This level of anisotropy can be considered marginally consistent with an isotropic L_bol − Y_SZ relation. If however there was some extra X-ray absorption taking place toward that direction, it would correspond to an extra N_Htot ∼ 7.9 ± 3.4 × 10²⁰ cm⁻². This value is consistent with the L_X − Y_SZ result of our sample within 1.1σ. The maximum anisotropy region however is located 90° away from the L_X − Y_SZ direction. This could indicate that the normalization variation is due to statistical noise and does not reveal any unknown absorption. Alternatively, it could be attributed to the incapability of the θ = 90° cones to accurately pinpoint the direction where the hidden X-ray absorption takes place, especially since the direction uncertainties are considerably large. Based on all the above, it cannot be excluded that the mild anisotropies seen in the ACC sample are due to chance.

Finally, ACC further confirms our previous conclusion about the L_X − T anisotropies found in M20. There is no indication that they appear due to previously unaccounted X-ray effects, since the region toward (l, b)∼(300°, −20°) shows a completely consistent behavior with the rest of the sky.

B.2. The L_bol − T relation

For the L_bol − T relation of ACC, the maximum anisotropy is found toward (l, b)=(318°±45°, −9°±37°) in a 2.2σ tension with the rest of the sky¹⁵. This region appears fainter than the rest by ∼32 ± 13% on average. Its behavior is almost identical to the results of M20, as seen in Fig. B.1. However, the decreased statistical significance suggests that the L_bol − T relation for ACC alone is marginally consistent with isotropy. The angular separation between the most anisotropic regions of ACC and our sample is 40° and well within the uncertainties.

For ACC, we would need H₀ = 61.3 ± 4.2 km s⁻¹ Mpc⁻¹ toward (l, b)∼(318°, −9°) and H₀ = 72.2 ± 2.4 km s⁻¹ Mpc⁻¹ for the rest of the sky. The obtained H₀ value agrees within 1σ with the independent result from our sample.

For the BF explanation, applying the MR method to ACC, we find a BF of u_BF = 850 ± 410 km s⁻¹ toward (l, b)=(254°±52° , + 18°±27°). This direction is separated by 31° from the CMB dipole and by 67° from the maximum anisotropy region, although this difference is within 1σ. Owing to the limited number of data and the large scatter, we can only divide the sample into two independent redshift bins, z < 0.2 and z > 0.2. For the former we practically find the same BF as for the full sample, while for the latter we find no statistical significant evidence of a BF. However, the results are inconclusive owing to the limited number of available clusters in this sample and the subdominant effect BFs have at these redshifts compared to the intrinsic scatter.

For the MA method, we obtain u_BF = 810 ± 400 km s⁻¹ toward (l, b)=(324°±45°, −3°±32°). This direction is more in agreement with the maximum anisotropy direction of ACC, while the BF amplitude is similar. The sample size is not sufficient to consider individual redshift bins that would provide different insights than the full sample.

B.3. The Y_SZ − T relation

We repeat the Y_SZ − T anisotropy analysis with only the 113 ACC clusters with Y_SZS/N > 2 with a median z ∼ 0.22. The scatter is similar to our sample, but the limited number of clusters forces us to consider θ = 90° cones. This way we ensure we have at least 35 clusters in each cone (which inevitably leads to large uncertainties). The normalization and sigma maps are shown in Fig. B.1. We find an anisotropy of 3.2σ toward (l, b)=(311°±57°, −12°±39°), 42° away from the direction of our sample and well within the 1σ uncertainties. The relative difference of A_YT of this region compared to the rest of the sky is 35 ± 11%. The direction is also identical to that obtained with the L_bol − T relation with a larger statistical significance.

In terms of H₀, we obtain H₀ = 60.6 ± 3.6 km s⁻¹ Mpc⁻¹ toward the most anisotropic region and H₀ = 73.4 ± 1.9 km s⁻¹ Mpc⁻¹ for the opposite hemisphere. These values are consistent within < 1.4σ with the results of our sample and with the joint L_X − T analysis results.

For the BF scenario and the MR method, we obtain u_BF = 800 ± 400 km s⁻¹ toward (l, b)=(268°±42° , + 6°±29°) for the full sample. Once again for ACC, the direction is close to the CMB dipole, but slightly shifted compared to the maximum anisotropy direction. The BF has its usual amplitude, even though the used sample has a large median z ∼ 0.22. For z < 0.2 we obtain a slightly larger BF than for the full sample. For z > 0.2, the BF points toward (l, b)∼(320° , + 15°); however it is poorly constrained and not statistical significant. For the MA method, we find u_BF = 810 ± 370 km s⁻¹ toward (l, b)=(273°±38°, −11°±21°) for the full sample, which agrees with the MR method.

Appendix C: Anisotropies in the other scaling relations

The R − L_X, R − T, R − Y_SZ, R − L_BCG, and Y_SZ − L_BCG scaling relations cannot currently provide meaningful insights into the possible origin of the detected anisotropies. However, they may prove to be very useful for future tests with larger samples and with a better characterization of the dynamical state of the clusters. In this appendix, we discuss their results, limitations, and future potential.

C.1. Anisotropies of the R −  L_X, R − T, R − Y_SZ, and R −  L_BCG scaling relations

The scaling relations of the cluster effective radius are potentially very interesting once the systematic biases are properly handled and the sample sizes are further increased. For the R − L_X, R − Y_SZ, and R − L_BCG relations, both quantities depend on the cosmological parameters, but because of the rather flat slopes, the normalization weakly depends on the angular diameter distance (). Thus, a 15% spatial variation of H₀ would result in a 7−10% variation of A_RX. Much smaller effects are expected as a result of BFs. This level of anisotropy would not be detectable over random noise combined with the mild, existing systematic issues. This is the most crucial limitation of these relations. Future samples, such as the eRASS catalogs, will dramatically increase the sample sizes, which will reduce the random noise and allow the detection of ∼10% anisotropies.

For the R − T relation, there is an A_RT ∝ D_A dependency. For the same underlying effect, R − T shows half the normalization variation than L_X − T, Y_SZ − T, and L_BCG − T. Despite the weaker signal, this test can still be a valuable to detect cosmological anisotropies and BFs when applied to larger cluster samples, and when a more precise modeling of the systematics is feasible.

The quantity R is also rather insensitive to unaccounted X-ray absorption effects. The latter can only affect the determination of R₅₀₀, which is based on the absorption-affected flux. However, , so a 20% bias in L_X would only cause a < 4% bias in R₅₀₀, which is negligible. To confirm this, we remeasure R for all clusters, after changing the input N_Htot by ±50%. Almost all clusters (99.5% of the sample) show an R change of < 2.5%.

The main systematic effect compromising the R scaling relations is the effect of the CC clusters on the measurements. This was already discussed in Sect. 4.4 and is more extensively discussed here. The directional behavior of these relations clearly correlates with the dynamical state of the clusters. The latter affects their surface brightness profiles and consequently their half-light radii. For the R − T relation, the relaxed systems appear to be 45 ± 7% smaller (lower normalization) than the disturbed systems. This constitutes a 6.5σ deviation (Fig. C.1). It is clear then that the anisotropies possibly detected in these scaling relations are entirely driven by the slightly inhomogeneous spatial distribution of such clusters and not by cosmological phenomena. As discussed before, the R scaling relations (contrary to L_X − T and Y_SZ − T) are not sensitive enough to cosmological effects to overcome the bias coming from CC clusters.

When we attempt to scan the sky with a θ = 75° cone, we obtain the same sky pattern for all four scaling relations. The normalization anisotropy maps for R − L_X and R − Y_SZ are shown in Fig. C.2. The map for R − L_BCG is not shown since it is very similar to the other R maps. The smallest clusters are consistently found within 12° from (l, b)∼(64° , + 32°). This direction strongly correlates with the region in which the highest fraction of relaxed clusters in our sample is found, where there is also a lack of disturbed clusters, as shown in Fig. 9. The maximum anisotropy does not exceed 2.3σ for any relation, with a ≲14% variation of A_RX. As such, the effect of this bias is not extremely strong, but is sufficient to dominate over the effect of the possible cosmological anisotropies and BFs.

C.2. Calibration of R − T anisotropies for dynamical state of clusters

To study the underlying cosmological effects, we would need larger samples to further smooth out the sky distribution of different dynamical types of clusters. Also, we could attempt to calibrate the scaling relations for this systematic bias. This requires an observable proxy that effectively traces the existence of CC clusters. Unfortunately, such an observable is not currently available for our sample. Future work will soon provide a remarkable characterization of the core state of the eeHIFLUGCS clusters using numerous independent measurements. Proper calibration of the R relations will be then possible, followed by their application for the search of cosmological anisotropies.

For now, we use XBO as a proxy for relaxed clusters, as already discussed. We attempt to calibrate the R − T relation and its anisotropy map. We create 10⁵ randomly drawn bootstrap subsamples (same process as in M20), independent of direction. We investigate the correlation between the A_RT and the median XBO for every subsample. The results are shown in Fig. C.3. We then calibrate the A_RT anisotropy map based on the median XBO of every sky region and the observed correlation between the two. Although the induced uncertainties due to this calibration are too large to allow for any meaningful conclusions, this offers a useful example for potential future applications of the R relations. Finally, the behavior of the apparent R − T anisotropies after this correction is applied is very interesting. The previously anisotropic region now becomes milder. The (l, b)∼(260°, −10°) direction starts showing a lower A_RT behavior, which is consistent with the previous results and the scenario of a cosmological anisotropy or a large BF.

C.3. The Y_SZ −  L_BCG relation

Both quantities of this relation are unaffected by absorption effects due to the infrared and submillimeter wavelengths, while they depend on the cosmological parameters in the same manner. Accounting for the overall best-fit slope, a very weak dependence of the normalization on the angular diameter distance remains (). In detail, a 15% variation in H₀ would only lead to a < 3% variation in A_YLLBCG. Similarly, a 1000 km s⁻¹ BF at z = 0.05 would only cause a 0.5% change in A_YLLBCG. It is clear then that only sample-related biases can create observable anisotropies in this scaling relation.

As a result of the large scatter and the limited number of clusters, we consider θ = 90° cones to scan the sky. The anisotropy of the Y_SZ − L_BCG normalization is shown in Fig. C.4. The relation shows a ∼30% variation mostly between the Northern and Southern Galactic Hemispheres, but with a negligible statistical significance of 1.4σ. The relation can be considered statistically isotropic, as expected, with no strong biases toward a direction. The normalization map is shown in Fig. C.4.

Appendix D: Scaling relations and anisotropies as functions of different selection cuts

In this section, we repeat the analysis for the L_X − Y_SZ, Y_SZ − T, and L_BCG − T relations, for different Y_SZ S/N and redshift cuts. The summarized result is that our main conclusions remain unchanged.

D.1. The L_X − Y_SZ relation

For the L_X − Y_SZ relation, a lower Y_SZ threshold of S/N ≥ 4.5 was applied in the default analysis. This resulted in 460 clusters. If we consider two lower thresholds instead, namely S/N ≥ 2 and S/N ≥ 3, we have 1095 and 747 clusters, respectively. The 3σL_X − Y_SZ contours for the three different cases are compared to Fig. D.1. There is a > 5σ shift in the best-fit L_X − Y_SZ relation going from S/N > 4.5 to S/N > 2, while the scatter increases by ∼100%. The S/N > 2 fit, however, is dominated by several systematics. To begin with, the L_X residuals are strongly correlated with z, N_Htot, and Y_SZ, as shown in Fig. D.2. These probably point toward biases in the X-ray selection process and the need for a broken power law to describe L_X − Y_SZ. The same figure show that these systematic behaviors are not present for the S/N > 4.5 case.

Despite these issues, when we scan the sky to detect L_X − Y_SZ anisotropies, we get a similar map as for the default case. The anisotropy significance map (Fig. D.3) illustrates that there are no excess X-ray absorption issues toward (l, b)∼(280°, −15°) where the cosmological anisotropies are found. At the same time, the only anisotropic region again appears to be toward (l, b)∼(128° , + 21°), being ∼14% fainter than the rest of the sky at a 2.5σ level, similar to the default case.

D.2. The Y_SZ − T relation

For the Y_SZ − T relation, we used a lower Y_SZ threshold of S/N ≥ 2 in the default analysis, resulting in 260 clusters. If we increase this lower threshold to S/N ≥ 3 or S/N ≥ 4.5, we are left with 242 and 190 clusters, respectively. The best-fit Y_SZ − T does not change as a function of S/N, as shown in the middle panel of Fig. D.1. The Y_SZ residuals remain independent of the cluster physical properties as well. When we repeat the anisotropy analysis for S/N ≥ 3 and for S/N ≥ 4.5, we obtain almost identical anisotropy results as for the default case. The A_YT variance map is shown in Fig. D.4. Thus, the observed anisotropies are independent of Y_SZ selection effects.

D.3. The Y_SZ − T relation using the Planck values instead

To ensure that the observed Y_SZ − T anisotropies do not emerge because of some unknown directional bias in our own Y_SZ measurements, we repeat the analysis using the Y_SZ values from the PSZ2 catalog. Owing to the higher Y_SZS/N > 4.5 threshold that PSZ2 applied, 206 clusters were matched with our M20 sample. This S/N cut also leads to higher median T and Y_SZ compared to our default analysis. Thus, we adopt C_Y = 85 kpc² and C_X = 6 keV. For the best-fit parameters, we find A_YT = 0.940 ± 0.041, B_YT = 2.074 ± 0.095, σ_int = 0.184 ± 0.019, and σ_int = 0.232 ± 0.024. The scatter is ∼35% larger than in our case (when Y_SZS/N > 4.5) as a result of the different adopted R₅₀₀ between the two independent analyses. This highlights the need to use our own measurement, in which we increase the number of clusters while decreasing the scatter. The slope differs by 2.2σ. The Y_SZ − T scaling relation is plotted in the top panel of Fig. D.5.

Performing the Y_SZ − T sky scanning, we obtain roughly the same anisotropies with the PSZ2 Y_SZ values as we did with our measurements. The A_YT map is shown in the bottom panel of Fig. D.5. The maximum anisotropy is found toward (l, b)=(254°, −22°), where H₀ = 64.6 ± 2.5 km s⁻¹, at a 2.9σ tension with the rest of the sky. This region is only 15° from the most anisotropic region as found in our default analysis, which demonstrates that the Y_SZ − T anisotropies are independent of the adopted Y_SZ catalog.

D.4. The L_BCG − T relation

For the L_BCG scaling relations, we removed clusters at z < 0.03 since they appear systematically brighter than expected. We also excluded all the z > 0.15 clusters, since the redshift evolution of L_BCG remains unknown. We applied these cuts to all L_BCG scaling relations to be consistent. However, these redshift issues are not as important in the L_BCG − T relation. In Fig. D.6, the L_BCG residuals are plotted as a function of z. Only the seven clusters with z < 0.02 are systematically upscattered, while z > 0.15 show the same behavior as the less distant clusters. Thus, we repeat our analysis using all the 259 clusters, regardless of their redshift. The best-fit L_BCG − T relation remains similar to the default analysis, with the parameter uncertainties being naturally smaller, due to the larger number of clusters. The anisotropic behavior of the relation remains the same as before. Interestingly, the statistical significance of the L_BCG − T anisotropy signal rises from 1.9σ to 2.1σ. Both of these results are shown in Fig. D.6.

Appendix E: More details about performed tests

E.1. Zone of avoidance bias

As discussed in Sect. 8.3, the ZoA gap does not introduce any significant bias to our results. In this section we provide some plots related to that discussion. For simplicity, we only provide the relative images for the L_X − T relation since results from the other relations are similar. This is because the samples used are simulated and isotropic and that the cluster positions are the same across scaling relations. Thus, any effects come only from the applied cluster weighting during the sky scanning. We stress again that the observed anisotropies remain unaffected when this weighting is omitted and the ZoA gap is invisible to the used algorithm, highlighting the lack of bias coming from ZoA in the real data.

In the top panel of Fig. E.1, the Galactic latitude of the most anisotropic region for every simulated sample is shown. As explained before, ∼43% of these regions lie within ZoA instead of the ∼33%, which would be the expectation if no bias was present. A nonflat distribution of the galactic latitudes is expected even in the completely bias-free case because the sky area covered by each bin is not the same. Hence, the probability of the most anisotropic region to be located within each bin changes proportionally.

In the bottom panel of Fig. E.1, we compare the relative difference of the observed anisotropy (in terms of σ) for two cases. The first is when the results from the analysis above are considered, without any data in the ZoA. The second is when we fill the ZoA with simulated samples, as described in Sect. 8.3. As discussed there, the average maximum anisotropy signal is increased by 14 ± 8% when the ZoA clusters are excluded. This is mostly due to the different number of available data in the two cases.

Finally, in Fig. E.2 the number of clusters within each cone with θ = 75° is shown. As expected, close to the Galactic center fewer clusters are included in the cones, while the most clusters are found for the Galactic pole cones. The main anisotropic region of our analysis at (l, b)∼(280°, −15°) shows an average number of clusters.

E.2. MCMC fitting

In Sect. 7 we discussed the results obtained by the MCMC fitting. In this appendix, we provide some more details and plots about that test. Owing to the high number of free parameters, we performed 2 × 10⁷ iterations of the chain, with a burn in period of 10⁴. A variable step size was used for every parameter, randomly drawn from the same Gaussian distribution with a zero mean, and a standard deviation of σ = 0.05. Every 10⁴ iterations, σ = 0.5 was used a single time to fully explore the possibility of multiple likelihood maxima within the complicated parameter space. An acceptance rate of 18% was reached. The chain was run multiple times from varying initial positions and step size source distributions to ensure that the same results were reached each time.

In Fig. E.3, we plot the 3σ parameter space for u_z and u_T (redshift and temperature power indexes), which are the only two parameters with a significant impact on the best-fit normalization. As discussed before, their anticorrelation is strong, which generally cancels out any strong effects in the normalization owing to the flux-limited sample and the Eddington bias. The 3σ parameter space for u_NH and u_f (N_Htot and flux) is also plotted as a representative example for the rest of the free parameters. There is no significant effect from these cluster properties (u_NH ∼ u_f ∼ 0) and no correlation between them. This demonstrates that the rest of the cluster parameters do not induce any biases in the observed normalization of a cluster subsample.

E.3. Sky fraction of Chandra and XMM-Newton clusters in our sample

In Fig. E.4, we show the fraction of the clusters for which Chandra or XMM-Newton were used to determine T. Specifically, the fraction is defined as the Chandra clusters minus the XMM-Newton clusters over the sum. The results are discussed in Sect. 8.5.

E.4. Effects of bulk flows to adopted apparent R₅₀₀

If a BF is present, then the distance of a cluster is miscalculated. As a result, its L_X, and its mass M₅₀₀, its physical R₅₀₀ (Mpc), and its apparent θ₅₀₀ (arcmin) are also misinterpreted (see Piffaretti et al. 2011 and Sect. 2.4 of M20). Thus, the region within which we measure T, Y_SZ, and L_X also changes, possibly yielding different results. This is a negligible effect, however, and is not taken into account. In M20 (Appendix B), we showed that for H₀ anisotropies, θ₅₀₀ remains practically unchanged. That is because the opposite changes of R₅₀₀ and D_A cancel each other out. Here, we test if a BF can cause significant changes in the adopted θ₅₀₀ and how this would affect our used parameter values. We combine the relation of Arnaud et al. (2010) with the fact that and that D_L = (1 + z)²D_A. Then, θ₅₀₀ is written as

(E.1)

We typically find BFs of ∼1000 km s⁻¹. For clusters at z = 0.05, this would cause a ≲6% change in θ₅₀₀. For z = 0.1, this change would be ≲2%. For T, the percent changes are roughly the same as in θ₅₀₀ (see Appendix B in M20). Therefore, for the BF case we considered, and for low-z clusters, we expect a ≲6% bias. As an example, let us estimate how this would affect the anisotropies of the L_X − T relation. If clusters appear less luminous (or closer) than expected as a result of a BF (as for the main anisotropic region of our results), this would mean that the distance is underestimated, hence θ₅₀₀ is slightly overestimated. Thus, the measured T is actually slightly underestimated, since it was measured in an annulus further from the center than planned, where clusters are generally cooler. If we measured the correct, higher T, then these clusters would appear even fainter than before compared to the expectations. Therefore, the observed anisotropies would be amplified. Consequently, any small bias that is introduced to θ₅₀₀ owing to an existing BF would eventually suppress the BF signal, which is the opposite of what we would need to explain the anisotropies.

The Y_SZ measurements also suffer the same percent changes as θ₅₀₀ (Fig. E.5). For the same BF case as above, the employed MMF technique leads to a slight overestimation of Y_SZ, which again suppress the BF signal. However, this effect is insignificant for two reasons. Firstly, because of the slope of the Y_SZ − T relation (B_YT = 2.546), the changes of Y_SZ are less important than changes in T (which we already saw that can mildly smooth out the anisotropies). Secondly, the Y_SZ − T anisotropies are much larger than the Y_SZ changes owing to a falsely assumed θ₅₀₀.

Moreover, the measured T value goes into the measurement of Y_SZ when relativistic effects are considered. However, this dependence is very weak since a 1 keV change in T would only lead to a ∼1% change in Y_SZ, and thus a BF would practically not affect Y_SZ in that regard.

For the L_X measurements we cannot make quantitive predictions of the change they suffer for a different θ₅₀₀ input since we did not conduct the measurements. However, we expect L_X to remain almost unchanged, since only the outskirts of the assumed cluster size slightly changes in case of a BF. The vast majority of X-ray cluster emission however comes from well within this area. Even if the change was not negligible, the same effect as for T and Y_SZ is expected, where the miscalculation of L_X actually makes us underestimate the BF signal rather than create it.

It is evident that these changes are only minimal compared to the observed anisotropy amplitudes. It is not trivial to account for these changes because iterative Y_SZ and T measurements and BF estimations would be needed, and these changes would further enhance the observed anisotropies. Thus we ignore these changes.

All Tables

All Figures