Issue 
A&A
Volume 649, May 2021



Article Number  A151  
Number of page(s)  38  
Section  Cosmology (including clusters of galaxies)  
DOI  https://doi.org/10.1051/00046361/202140296  
Published online  01 June 2021 
Cosmological implications of the anisotropy of ten galaxy cluster scaling relations
^{1}
ArgelanderInstitut für Astronomie, Universität Bonn, Auf dem Hügel 71, 53121 Bonn, Germany
email: kmigkas@astro.unibonn.de
^{2}
Center for Astrophysics  Harvard & Smithsonian, 60 Garden Street, Cambridge, MA 02138, USA
^{3}
Deutsches Zentrum für Luft und Raumfahrt e.V. (DLR) Projektträger, JosephBeuysAllee 4, 53113 Bonn, Germany
^{4}
Institut für Astronomie und Astrophysik, Sand 1, 72076 Tübingen, Germany
^{5}
Max Planck Institute for Extraterrestrial Physics, Gießenbachstraße 1, 85748 Garching bei München, Germany
^{6}
INAF – Osservatorio Astronomico di Bologna, Via Piero Gobetti, 93/3, 40129 Bologna, BO, Italy
Received:
7
January
2021
Accepted:
22
March
2021
The hypothesis that the late Universe is isotropic and homogeneous is adopted by most cosmological studies, including studies of galaxy clusters. The cosmic expansion rate H_{0} is thought to be spatially constant, while bulk flows are often presumed to be negligible compared to the Hubble expansion, even at local scales. The effects of bulk flows on the redshift–distance conversion are hence usually ignored. Any deviation from this consensus can strongly bias the results of such studies, and thus the importance of testing these assumptions cannot be understated. Scaling relations of galaxy clusters can be effectively used for this testing. In previous works, we observed strong anisotropies in cluster scaling relations, whose origins remain ambiguous. By measuring many different cluster properties, several scaling relations with different sensitivities can be built. Nearly independent tests of cosmic isotropy and large bulk flows are then feasible. In this work, we make use of up to 570 clusters with measured properties at Xray, microwave, and infrared wavelengths to construct ten different cluster scaling relations and test the isotropy of the local Universe; to our knowedge, we present five of these scaling relations for the first time. Through rigorous and robust tests, we ensure that our analysis is not prone to generally known systematic biases and Xray absorption issues. By combining all available information, we detect an apparent 9% spatial variation in the local H_{0} between (l, b)∼(280°_{−35°}^{+35°}, −15°_{−20°}^{+20°}) and the rest of the sky. The observed anisotropy has a nearly dipole form. Using isotropic Monte Carlo simulations, we assess the statistical significance of the anisotropy to be > 5σ. This result could also be attributed to a ∼900 km s^{−1} bulk flow, which seems to extend out to at least ∼500 Mpc. These two effects will be indistinguishable until more highz clusters are observed by future allsky surveys such as eROSITA.
Key words: cosmology: observations / largescale structure of Universe / galaxies: clusters: general / Xrays: galaxies: clusters / methods: statistical
© 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 socalled 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 powerlaw forms. Some of the measured cluster properties depend on the assumed values of the cosmological parameters (e.g., Xray 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 Xray 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 Xray 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 Xray absorption, bulk flows (BFs), and Hubble expansion anisotropies. Firstly, the existence of yet undiscovered excess Xray 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 redshiftbased 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 largescale 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 indepth 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_{0}, 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 lowz 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 Xray 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 halflight 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_{0} = 70 km s^{−1} Mpc^{−1}, Ω_{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 bestfit parameters for the ten scaling relations are presented. In Sect. 5, we study the anisotropic behavior of scaling relations, which are sensitive only to Xray 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 MetaCatalog of Xray 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 Xray flux of f_{X, 0.1−2.4 keV} ≥ 5 × 10^{−12} erg s^{−1} 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 Xray luminosity L_{X}, temperature T, redshift z, R_{500}, and metallicities of the core Z_{core} and 0.2−0.5 R_{500} 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 Xray 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
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_{B}n_{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_{500} 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 allsky 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 nonrelativistic 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 signaltonoise 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 Xray 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
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. Halflight radius R
The apparent angular size R_{app} in the sky within which half of the total Xray 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 halflight 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^{−12} erg s^{−1} cm^{−2}. This led to 438 measurements.
To measure R_{app} the ROSAT AllSky Survey (RASS) maps were used. The countrate 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^{−12} erg s^{−1} cm^{−2}. By a combination of applied iterative algorithms and visual validation by six different astronomers, the plateau of the countrate growth curves (i.e., the boundaries of cluster Xray 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 offaxis 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. Nearinfrared 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 nearinfrared (NIR) data from the SDSS (York et al. 2000), PanSTARRS (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)^{2}. All the galaxy magnitudes were corrected for Galactic extinction and the proper kcorrection 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. Xray determined N_{H, Xray}
We determined the total hydrogen column density N_{H,Xray} using the Xray 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. Xray 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_{SZ}S/N < 2. This results in 168 clusters with Xray luminosity and temperature values; 113 of these clusters have a Y_{SZ} measurement with S/N > 2. Crosschecking 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 singlethermal model for the whole cluster, while the Xray luminosities L_{bol} are given in the bolometric band within the R_{200} 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_{500} values to measure Y_{SZ} were obtained using the masstemperature scaling relation of Reichert et al. (2011). The latter mostly uses XMMNewtonderived 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_{500}.
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 Xray 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.
Possible anisotropy causes that can be traced by each cluster scaling relation.
The form of the studied scaling relations between some measured cluster quantities Y and X is given by
where C_{Y} is the calibration term for the Y quantity, the term E(z)=[Ω_{m}(1 + z)^{3} + Ω_{Λ}]^{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 (selfsimilar) values of γ_{YX}.
Bestfit parameters of the 10 scaling relations.
3.1. Linear regression
Similar to M20, to constrain the scaling relation parameters we perform a linear regression in the logarithmic space using a χ^{2} 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
where N is the number of clusters, , ^{3}, σ_{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 ^{4}. σ_{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 bestfit scaling relations (i.e., all the T scaling relations). In this case we assume the normalization of the scaling relation to be known and directionindependent. 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_{0} 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:
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_{0}); 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
In this equation, D_{H0 = 70} is the distance found for H_{0} = 70 km s^{−1} Mpc^{−1}, 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 Xray flux (in the rest frame of the cluster), apparent size in arcmin, Y_{5R500}, and BCG flux respectively. Moreover, D_{th} is given by
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_{0} variations express relative differences between regions. The absolute H_{0} values cannot be constrained by cluster scaling relations only since a fiducial H_{0} value was initially assumed to calibrate the relation. Hence, the H_{0} anisotropy range always extends around the initial H_{0} choice. Apparent H_{0} 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 bestfit 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 bestfit 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 bestfit 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 fluxlimited 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 lowz objects and do not clearly reflect the largerscale 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 lowz 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 bestfit 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 bestfit 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 bestfit value of the full (sub)sample. Each parameter value distribution is considered separately from the rest^{5}. 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., Δχ^{2} 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 bestfit parameters with their uncertainties, as described in the two previous sections. We then create a colorcoded full sky map based on the bestfit 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 bestfit 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^{6}, 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θ_{1} term, where θ_{1} 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_{0} anisotropy to the measured cluster distance scale with cosθ_{1}. 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 bestfit 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)^{7}. We assess the statistical significance of the deviation between them by
where p_{1, 2} are the bestfit 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 Δχ^{2} 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 bestfit 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 bestfit results of all the scaling relations is given in Table 2, while the scaling relations themselves are plotted in Fig. 1.
Fig. 1. Best fits of the 10 scaling relations studied in this work. The green area represents the 1σ limits of the bestfit area. From top left panel to the right: L_{X} − Y_{SZ}, L_{X} − L_{BCG}, L_{X} − T, Y_{SZ} − T, L_{BCG} − T, R − L_{X}, R − Y_{SZ}, R − L_{BCG}, R − T, Y_{SZ} − L_{BCG}, R − Y_{SZ}, L_{bol} − Y_{SZ} (ACC), and Y_{SZ} − T (ACC) relations. 
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 bestfit 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 χ^{2}minimization in the L_{bol} axis, using bootstrap for estimating the parameter uncertainties. Since the bolometric Xray luminosity is used for ACC (within R_{200}), 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 & AtrioBarandela 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 bestfit parameters and the possibly observed anisotropies of the relation. For these 460 clusters, the bestfit 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_{500} between L_{X} and Y_{SZ} in our analysis. The observed slope B_{LY} ∼ 0.93 is slightly larger than the selfsimilar 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 lowz clusters, respectively, to compare the L_{X} − L_{BCG} relation between cool and noncoolcore (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 bestfit slope of our analysis is less steep than the bestfit slope found by Bharadwaj et al. (2014); in contrast to this work the latter authors used the bolometric Xray 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_{500} − L_{X} relation (selfsimilar prediction of γ ∼ 1.3) is not expected since a redshift evolution between the halflight radius R and R_{500} is also expected (e.g., owing to the timevarying, 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 downscattered 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 nonnegligible correlation is observed between the residuals and the apparent halflight radius R_{app}. The clusters with the lowest R_{app} appear to be downscattered in the R − L_{X} plane, and vice versa. A mild correlation of the R residuals is also observed with the offset between the Xray 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). Coolcore clusters are expected to strongly bias the scaling relations involving R, since they emit most of their Xray photons from near their centers. As a result, the halflight radius is lower compared with nonCC clusters at a fixed mass. More details on that can be found in Appendix C.1.
Finally, the slope is lower than the selfsimilar prediction for the R_{500} − 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 bestfit parameter values that we obtain for these clusters are in line with previous findings. The value of the slope agrees with the selfsimilar 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_{500}).
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 upscattered. This behavior persists even with an increasing S/N threshold. Although the bestfit 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 Xray 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 bestfit 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 XrayBCG 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 bestfit 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 Xray absorption effects
In this section we study the scaling relations, whose observed anisotropies could be caused purely by previously unknown soft Xray effects, such as extra absorption from “dark”, metalrich, 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 Xray 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 Xray 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_{0} 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_{0}, no significant change in the bestfit A_{LY}B_{LY} would be observed. Analytically, as and , their ratio (considering their bestfit B_{LY} = 0.928) would be
Thus, a ∼15% spatial variation of H_{0} 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^{−1} 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 Xray 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^{20} cm^{−2} (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 lowerthansolar 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.
Fig. 2. Normalization anisotropy maps (left) and the respective statistical significance maps of the anisotropies (right), for L_{X} − Y_{SZ} (top), joint L_{X} − Y_{SZ} (our sample + ACC, middle), and L_{BCG} − T (bottom). All the maps in this work are shown in a Hammer projection. 
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 Xray absorption effects. Therefore, jointly analyzing the two samples should provide us with better insight into any possible Xray 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_{0} 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^{20} cm^{−2} would be needed, or alternatively a Z ∼ 1.49 ± 0.23 Z_{⊙} for the alreadydetected Galactic gas and dust.
As such, the tension could be attributed to chance and not necessarily to an unaccounted Xray 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 nearinfrared 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 Xray 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 anisotropyinducing 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^{20} cm^{−2}. Alternatively, a metal abundance of Z ∼ 2.4 ± 0.8 Z_{⊙} of the alreadydetected hydrogen cloud would also alleviate this small tension. Thus, the L_{X} − L_{BCG} relation does not show any indications of previously unknown Xray 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 Xray 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 Xray 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 bestfiT relation is cm^{−2}. The Xray 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 upscatter of the remaining N_{H,Xray} values, does not necessarily mean that the true Xray 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 directiondependent. It should be borne in mind that this is an approximate, complementary test as a result of its limitations and not a standalone check for excess Xray 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 bestfiT relation is cm^{−2}. 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 Xray absorption is the cause behind the L_{X} − T anisotropies.
Fig. 3. Comparison between the N_{Htot} from W13 and the N_{Htot} constrained by our Xray spectral analysis. The equality line is indicated in purple; the 1σ space of the bestfit function for the region of interest (green, black points) and for the rest of the sky (cyan, red points) are also shown. The reason for the systematic difference is given in the main text. The region of interest is within 45° from (l, b)=(281°, −16°) (top) and from (l, b)=(122° , + 8°) (bottom). 
The second region we consider is within 45° around (l, b)=(122° , + 8°). This region shows mild indications of uncalibrated true Xray absorption that differs from the W13 values. Using its 11 clusters, we find cm^{−2}. The results (Fig. 3) are again consistent with the rest of the sky and do not reveal any signs of a biased applied Xray absorption correction.
5.4. Overall conclusion on possible Xray 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 Xray 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 AllSky 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 Xray absorption correction.
If a mysterious hard Xray absorbing material existed that did not affect soft Xrays, 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 Xray 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_{0} variations
We have established that strong biases in the Xray 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, nonnegligible, 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 cosmologydependent quantities. The same is true if the real underlying values of the cosmological parameters, for example, H_{0}, vary from region to region. This spatial H_{0} 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_{0} 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_{0} 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.
Fig. 4. Top: same as in Fig. 2, for L_{X} − T. Bottom: H_{0} anisotropy map derived from the joint L_{X} − T (our sample+ACC). 
Cosmological anisotropies and BFs. We assume that the cause of the observed tension is an anisotropic H_{0} value in a Universe without any BFs. We would need H_{0} = 66.0 ± 1.7 km s^{−1} Mpc^{−1} toward (l, b)∼(274°, −9°), and H_{0} = 72.1 ± 1.4 km s^{−1} Mpc^{−1} for the rest of the sky. The respective Hubble diagrams are compared in the top panel of Fig. 5.
Fig. 5. Hubble diagram of galaxy clusters as derived by the L_{X} − T (top) and the Y_{SZ} − T (bottom) relations. The clusters from the most anisotropic region of each scaling relation are shown (blue) together with the clusters from the rest of the sky (red). The bestfit lines are indicated with the same color. 
We now assume that a BF motion is the sole origin of the observed anisotropies, where H_{0} is isotropic. Based on the MR method, we find a BF of u_{BF} = 980 ± 300 km s^{−1} 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^{−1} 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^{−1} toward (l, b)=(262°±52° , + 2°±26°), and u_{BF} = 1040 ± 570 km s^{−1} 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^{−1} 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^{−1} 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)^{8} and only consider the 170 clusters at larger distances (z > 0.067). For these, we obtain u_{BF} = 620 ± 310 km s^{−1} 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_{0}, similar to Fig. 23 of M20. The results are plotted in Fig. 4. The free parameter H_{0} seems to vary within ∼65−76 km s^{−1} Mpc^{−1}. The most anisotropic region is found toward (l, b)=(284°, −4°). The best fit of this region is H_{0} = 66.2 ± 1.6 km s^{−1} Mpc^{−1}, while for the rest of the sky we gets H_{0} = 72.7 ± 1.5 km s^{−1} Mpc^{−1}. While the posterior H_{0} 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 nonuse of XCSDR1, 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_{0} 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 Xray 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 samplerelated (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).
Fig. 6. Top: same as in Fig. 2 for Y_{SZ} − T. Bottom: H_{0} anisotropy map (left) and the respective significance map (right) derived from the joint Y_{SZ} − T (our sample+ACC, bottom). 
Cosmological anisotropies and BFs. We now investigate the necessary H_{0} variations to fully explain the observed anisotropies. We would need H_{0} = 62.3 ± 2.2 km s^{−1} Mpc^{−1} toward (l, b)∼(268°, −16°) and H_{0} = 72.1 ± 0.9 km s^{−1} Mpc^{−1} for the rest of the sky. The result is consistent within 1.4σ with the H_{0} 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^{−1} toward (l, b)=(263°±39°, −22°±20°). For the redshift bin z ∈ [0, 0.07], we obtain u_{BF} = 1060 ± 390 km s^{−1} toward (l, b)=(254°±42°, −17°±19°), similar to the results from the full sample. By gradually expanding the redshift range, the bestfit 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^{−1} 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^{−1} 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^{−1} toward (l, b)=(254°±22°, −18°±16°). For z > 0.09, the amplitude is reduced (u_{BF} = 720 ± 380 km s^{−1}), 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_{0} anisotropies with the Y_{SZ} − T relation. The H_{0} 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_{0} = 63.4 ± 2.5 km s^{−1} Mpc^{−1}, 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_{0}, we obtain H_{0} = 66.2 ± 2.5 km s^{−1} Mpc^{−1} for that direction and H_{0} = 72.6 ± 2.3 km s^{−1} Mpc^{−1} for the opposite hemisphere. The H_{0} 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^{−1} 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_{0} 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_{0} 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 Xray 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_{0} 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 nonCC 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_{0} spatial variation.
6.4.2. Apparent H_{0} 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_{0} = 66.5 ± 1.0 km s^{−1} Mpc^{−1} toward (l, b)=(273°±40°, −11°±27°), and H_{0} = 72.8 ± 0.6 km s^{−1} Mpc^{−1} 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_{0} 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_{0} uncertainty from which probes that attempt to put absolute constraints on H_{0} suffer. However, it is not needed for the relative spatial differences to which we are interested in. The joint H_{0} map together with the statistical significance map are shown in Fig. 7. The observed anisotropy forms a dipole.
Fig. 7. H_{0} anisotropy map as derived from the joint analysis of L_{X} − T, Y_{SZ} − T, and L_{BCG} − T relations for both samples. 
7. Comparison with isotropic Monte Carlo simulations
The wellestablished 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 bestfit 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 lognormal distribution with a standard deviation of . Since lowmass 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} ^{9} given the bestfit 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 bestfit 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 bestfit Y_{SZ} − T relation and T. Then, we further add the expected scatter of Y_{SZ} based on the observed bestfit 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 lognormal distribution combining the statistical uncertainties of the observed L_{X} and Y_{SZ} residuals and the intrinsic scatter. We constrain the bestfit 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.
Fig. 8. Histograms of the statistical significance of the maximum anisotropy as detected in 10 000 isotropic MC simulations for the L_{X} − T, Y_{SZ} − T, L_{BCG} − T, L_{bol} − T (ACC), and Y_{SZ} − T (ACC) scaling relations. The vertical black line represents the results from the real data; the pvalue and the direction are also shown. 
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^{−4}). 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^{−4}, 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^{−4}). Thus, we obtain p = 3.5 × 10^{−6}. 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^{−8}. 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 Xray and samplerelated 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. Coolcore and morphologically relaxed clusters
Coolcore clusters have a strong central peak in their surface brightness profile, and they are known to be intrinsically brighter in Xrays 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 Xray 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_{500} as morphologically relaxed and possibly CC. We also consider those with XBO > 0.08 R_{500} 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).
Fig. 9. Fractional difference between disturbed and relaxed clusters over all the clusters for every sky patch of the extragalactic sky. 
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.
Fig. 10. 3σ (99.7%) parameter space of the normalization and slope of the Y_{SZ} − T (left), L_{BCG} − T (middle), and L_{X} − T (right) relations for relaxed (purple) and disturbed (green) clusters. 
This was shown in M18, where using coreexcised 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^{5} 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}.
Fig. 11. Correlation between the bestfit A_{LT} (over the full sample’s bestfit value) and the median XBO for every of the 10^{5} bootstrap subsamples. 
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
Fluxlimited samples such as eeHIFLUGCS suffer from the socalled 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 Xray 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%^{11}. 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 Xray fluxlimited samples.
The MB correction to be applied in the normalization of the L_{X} − T for lowz, fluxlimited 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 highmass 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 Xray 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 bestfit 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^{6} random bootstrap cluster subsamples of random size (10−60% of the total size of the sample), as done in M20. Except for the bestfit A_{LT}, B_{LT} and σ_{int, LT}, 12 more average properties of each subsample are also derived^{12}. We express A_{LT} as a function of all the parameters p_{i} (normalized by their sample mean p_{mean}), and their powerlaw index v_{i}, as follows:
The 1.132 term corresponds to the bestfit 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 bestfit 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^{6} 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 bestfit 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 fluxlimited samples and as a result of the Eddington bias, highz clusters are also highT 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 Xray 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 XMMNewton (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, lowscatter relation however, which we used to convert XMMNewton temperatures to the equivalent Chandra temperatures. Therefore, we should not expect any dependence of the results on the fraction of Chandra and XMMNewton clusters. This is confirmed by the analysis in Sect. 8.4, where this fraction does not strongly correlate with the bestfit 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 XMMNewton 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 bestfit line has a slope of ∼0.63, while the total scatter of the relation is ∼0.16 dex.
Fig. 12. Correlation between the Y_{SZ} scatter from the Y_{SZ} − T relation and the L_{X} scatter from the L_{X} − T relation. 
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 bestfit scaling relations as a frame of reference. The biascorrected 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 biascorrected 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 Xray effects, such as excess Xray 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 Xray absorption toward the region in question, we can now decisively exclude this possibility. The lack of excess Xray 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 Xraydetermined N_{Htot} and the W13based N_{Htot} also does not indicate such an effect. Furthermore, the same anisotropies persist in scaling relations that are insensitive to unaccounted Xray 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 fluxlimited. 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_{0} 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, directiondependent, 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_{0} 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^{−1} extending to ≳500 Mpc scales (assuming H_{0} 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.
Maximum anisotropy direction for every scaling relation, together with the needed H_{0} relative variation to fully explain the anisotropy.
Bestfit BFs for every scaling relation, method, and redshift bin considered in this work.
Fig. 13. Top left: bulk flow amplitude and its 68.3% uncertainty as a function of the redshift radius of the used spherical volumes (or the comoving distance radius). Circles and squares correspond to the MR and MA methods, respectively. The black, red, and green points correspond to the L_{X} − T, Y_{SZ} − T, and L_{BCG} − T relations, respectively. The blue diamonds and triangles correspond to the ACC L_{bol} − T and Y_{SZ} − T relations, respectively. The results above z ≳ 0.16 are expected to be dominated by more local clusters and thus are overestimated. Top right: bulk flow amplitude and its 68.3% uncertainty as a function of the median redshift of each shell (or the comoving distance), together with the standard deviation of the redshift distribution. The color coding is the same as before. The data points with the same color and shape are independent to each other. Bottom: examples of BF directions as found for different redshift bins and methods. The color coding is the same as before. All directions agree with each other within 1σ. 
Unfortunately, the above explanations are not distinguishable with the current data sets and more highz 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. Nextgeneration 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^{13}. 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 largescale 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^{−1} velocity, similar to our BF results.
Similar largescale anisotropies in the distribution of highz 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^{−1} 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 Xray 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 maximumlikelihood analysis. In contrast, Soltis et al. (2019) find no evidence of departure from H_{0} 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. Largescale 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^{−1} 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 highz 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; AtrioBarandela et al. 2015), while Osborne et al. (2011), Mody & Hajian (2012), and Planck Collaboration Int. XIII (2014) fail to see any largescale 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^{−1} BF at ∼100−150 h^{−1} 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^{−1} 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^{−1} for 200 h^{−1} 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^{−1} 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 largerthanexpected 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 matterdominated 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_{0} 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_{0} variation due to these voids is slightly lower than the H_{0} 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 Xrays, 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 Xray 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_{0} to reconcile with the observations. Alternatively, we would need a BF motion of ∼900 km s^{−1}, 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.
This is the scale that many studies consider for studying BFs or local voids (e.g., Betoule et al. 2014; Carrick et al. 2015); see later discussion.
We represent σ_{int} in terms of percentile and not dex to make the following calculations clearer: σ_{int} mirrors the standard deviation of the distribution of clusters around the bestfit line. The relative difference of the two subsamples though corresponds to the standard error difference, which is the uncertainty of the mean of this distribution (i.e., the bestfit normalization). The statistical significance of this relative difference (i.e., anisotropy) accounts of course for σ_{int} through the normalization uncertainties.
These are the average redshift, temperature, < 0.2 R_{500} (core) temperature, flux, core metallicity, 0.2−0.5 R_{500} metallicity, N_{Htot}, RASS time exposure, Xray peakBCG offset, number of clusters, original source catalog (REFLEX/BCS = 1, NORAS = 2), and instrument used for T measurement (XMMNewton = 1, Chandra = 2).
Acknowledgments
We thank the anonymous referee for their constructive comments that helped us improve our manuscript. KM is a member of the MaxPlanck International School for Astronomy and Astrophysics (IMPRS) and of the BonnCologne Graduate School for Physics and Astronomy (BCGS), and thanks for their support. GS acknowledges support through NASA Chandra grant GO516137X. LL acknowledges financial contribution from the contracts ASIINAF Athena 201927HH.0, “Attività di Studio per la comunità scientifica di Astrofisica delle Alte Energie e Fisica Astroparticellare” (Accordo Attuativo ASIINAF n. 201714H.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
 Abbott, T., Abdalla, F., Allam, S., et al. 2018, ApJS, 239, 18 [Google Scholar]
 Allen, S. W., Evrard, A. E., & Mantz, A. B. 2011, ARA&A, 49, 409 [NASA ADS] [CrossRef] [Google Scholar]
 Andrade, U., Bengaly, C. A. P., Alcaniz, J. S., & Capozziello, S. 2019, MNRAS, 490, 4481 [NASA ADS] [CrossRef] [Google Scholar]
 Appleby, S., Shafieloo, A., & Johnson, A. 2015, ApJ, 801, 76 [NASA ADS] [CrossRef] [Google Scholar]
 Arnaud, K. A. 1996, in Astronomical Data Analysis Software and Systems V, eds. G. H. Jacoby, & J. Barnes, ASP Conf. Ser., 101, 17 [Google Scholar]
 Arnaud, M., Pratt, G. W., Piffaretti, R., et al. 2010, A&A, 517, A92 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Asplund, M., Grevesse, N., Sauval, A. J., & Scott, P. 2009, ARA&A, 47, 481 [NASA ADS] [CrossRef] [Google Scholar]
 AtrioBarandela, F., Kashlinsky, A., Ebeling, H., Fixsen, D. J., & Kocevski, D. 2015, ApJ, 810, 143 [NASA ADS] [CrossRef] [Google Scholar]
 Bender, A. N., Kennedy, J., Ade, P. A. R., et al. 2016, MNRAS, 460, 3432 [Google Scholar]
 Bengaly, C. A. P., Maartens, R., & Santos, M. G. 2018, J. Cosmol. Astropart. Phys., 2018, 031 [Google Scholar]
 Bengaly, C. A. P., Maartens, R., Randriamiarinarivo, N., & Baloyi, A. 2019, J. Cosmol. Astropart., 2019, 025 [Google Scholar]
 Betoule, M., Kessler, R., Guy, J., et al. 2014, A&A, 568, A22 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bharadwaj, V., Reiprich, T. H., Schellenberger, G., et al. 2014, A&A, 572, A46 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bharadwaj, V., Reiprich, T. H., Lovisari, L., & Eckmiller, H. J. 2015, A&A, 573, A75 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Blanton, M. R., & Roweis, S. 2007, AJ, 133, 734 [Google Scholar]
 Bleem, L. E., Stalder, B., de Haan, T., et al. 2015, ApJS, 216, 27 [NASA ADS] [CrossRef] [Google Scholar]
 Bleem, L. E., Bocquet, S., Stalder, B., et al. 2020, ApJS, 247, 25 [CrossRef] [Google Scholar]
 Böhringer, H., Chon, G., & Collins, C. A. 2020, A&A, 633, A19 [EDP Sciences] [Google Scholar]
 Bolejko, K., Nazer, M. A., & Wiltshire, D. L. 2016, J. Cosmol. Astropart., 2016, 035 [Google Scholar]
 Boruah, S. S., Hudson, M. J., & Lavaux, G. 2020, MNRAS, 498, 2703 [Google Scholar]
 Carrick, J., Turnbull, S. J., Lavaux, G., & Hudson, M. J. 2015, MNRAS, 450, 317 [Google Scholar]
 Chang, Z., & Lin, H.N. 2015, MNRAS, 446, 2952 [CrossRef] [Google Scholar]
 Chilingarian, I. V., & Zolotukhin, I. Y. 2011, MNRAS, 419, 1727 [Google Scholar]
 Chluba, J., Nagai, D., Sazonov, S., & Nelson, K. 2012, MNRAS, 426, 510 [NASA ADS] [CrossRef] [Google Scholar]
 Colin, J., Mohayaee, R., Sarkar, S., & Shafieloo, A. 2011, MNRAS, 414, 264 [NASA ADS] [CrossRef] [Google Scholar]
 Colin, J., Mohayaee, R., Rameez, M., & Sarkar, S. 2017, MNRAS, 471, 1045 [NASA ADS] [CrossRef] [Google Scholar]
 Colin, J., Mohayaee, R., Rameez, M., & Sarkar, S. 2019, A&A, 631, L13 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Dai, D.C., Kinney, W. H., & Stojkovic, D. 2011, J. Cosmol. Astropart., 2011, 015 [Google Scholar]
 Das, K. K., Sankharva, K., & Jain, P. 2021, J. Cosmol. Astropart., submited [arXiv:2101.11016] [Google Scholar]
 De Martino, I., & AtrioBarandela, F. 2016, MNRAS, 461, 3222 [Google Scholar]
 Deng, H.K., & Wei, H. 2018, Eur. Phys. J. C, 78, 755 [CrossRef] [EDP Sciences] [Google Scholar]
 Dąbrowski, M. P., & Wagner, F. 2020, Eur. Phys. J. C, 80, 676 [EDP Sciences] [Google Scholar]
 Eckert, D., Molendi, S., & Paltani, S. 2011, A&A, 526, A79 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Erler, J., Basu, K., Chluba, J., & Bertoldi, F. 2018, MNRAS, 476, 3360 [NASA ADS] [CrossRef] [Google Scholar]
 Erler, J., RamosCeja, M. E., Basu, K., & Bertoldi, F. 2019, MNRAS, 484, 1988 [Google Scholar]
 Ettori, S., Lovisari, L., & Sereno, M. 2020, A&A, 644, A111 [EDP Sciences] [Google Scholar]
 Feindt, U., Kerschhaggl, M., Kowalski, M., et al. 2013, A&A, 560, A90 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Filippou, K., & Tsagas, C. G. 2021, Ap&SS, 366, 4 [Google Scholar]
 Fitzpatrick, E. L. 1999, PASP, 111, 63 [NASA ADS] [CrossRef] [Google Scholar]
 Fosalba, P., & Gaztañaga, E. 2021, MNRAS, 504, 5840 [Google Scholar]
 Furnell, K. E., Collins, C. A., Kelvin, L. S., et al. 2018, MNRAS, 478, 4952 [NASA ADS] [CrossRef] [Google Scholar]
 Harrison, E. R. 1974, ApJ, 191, L51 [NASA ADS] [CrossRef] [Google Scholar]
 Haslbauer, M., Banik, I., & Kroupa, P. 2020, MNRAS, 499, 2845 [Google Scholar]
 Hasselfield, M., Hilton, M., Marriage, T. A., et al. 2013, J. Cosmol. Astropart., 2013, 008 [Google Scholar]
 Hilton, M., Sifón, C., Naess, S., et al. 2021, ApJS, 253, 3 [Google Scholar]
 Hoffman, Y., Courtois, H. M., & Tully, R. B. 2015, MNRAS, 449, 4494 [Google Scholar]
 Horner, D. J. 2001, PhD Thesis, University of Maryland College Park, USA [Google Scholar]
 Horvath, I., Szécsi, D., Hakkila, J., et al. 2020, MNRAS, 498, 2544 [Google Scholar]
 Hu, J. P., Wang, Y. Y., & Wang, F. Y. 2020, A&A, 643, A93 [EDP Sciences] [Google Scholar]
 Hudson, M. J., Smith, R. J., Lucey, J. R., Schlegel, D. J., & Davies, R. L. 1999, ApJ, 512, L79 [NASA ADS] [CrossRef] [Google Scholar]
 Hudson, M. J., Smith, R. J., Lucey, J. R., & Branchini, E. 2004, MNRAS, 352, 61 [Google Scholar]
 Hudson, D. S., Mittal, R., Reiprich, T. H., et al. 2010, A&A, 513, A37 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Huterer, D. 2020, ApJ, 904, L28 [Google Scholar]
 Itoh, N., Kohyama, Y., & Nozawa, S. 1998, ApJ, 502, 7 [Google Scholar]
 Kaiser, N. 1986, MNRAS, 222, 323 [Google Scholar]
 Kaiser, N., Aussel, H., Burke, B. E., et al. 2002, in Survey and Other Telescope Technologies and Discoveries, Int. Soc. Opt. Photon., 4836, 154 [Google Scholar]
 Kaiser, N., Burgett, W., Chambers, K., et al. 2010, in Groundbased and Airborne Telescopes III, Int. Soc. Opt. Photon., 7733, 77330E [Google Scholar]
 Kashlinsky, A., AtrioBarandela, F., Kocevski, D., & Ebeling, H. 2008, ApJ, 686, L49 [Google Scholar]
 Kashlinsky, A., AtrioBarandela, F., Ebeling, H., Edge, A., & Kocevski, D. 2010, ApJ, 712, L81 [Google Scholar]
 Kazantzidis, L., & Perivolaropoulos, L. 2020, Phys. Rev. D, 102, 023520 [CrossRef] [Google Scholar]
 Keenan, R. C., Barger, A. J., Cowie, L. L., et al. 2012, ApJ, 754, 131 [NASA ADS] [CrossRef] [Google Scholar]
 Korkidis, G., Pavlidou, V., Tassis, K., et al. 2020, A&A, 639, A122 [EDP Sciences] [Google Scholar]
 Lauer, T. R., & Postman, M. 1994, ApJ, 425, 418 [NASA ADS] [CrossRef] [Google Scholar]
 Lavaux, G., Afshordi, N., & Hudson, M. J. 2013, MNRAS, 430, 1617 [NASA ADS] [CrossRef] [Google Scholar]
 Li, M., Pan, J., Gao, L., et al. 2012, ApJ, 761, 151 [Google Scholar]
 Lopes, P. A. A., Trevisan, M., Laganá, T. F., et al. 2018, MNRAS, 478, 5473 [NASA ADS] [CrossRef] [Google Scholar]
 Lovisari, L., Schellenberger, G., Sereno, M., et al. 2020, ApJ, 892, 102 [CrossRef] [Google Scholar]
 Magoulas, C., Springob, C., Colless, M., et al. 2016, in The Zeldovich Universe: Genesis and Growth of the Cosmic Web, eds. R. van de Weygaert, S. Shandarin, E. Saar, & J. Einasto, 308, 336 [Google Scholar]
 Manolopoulou, M., Hoyle, B., Mann, R. G., Sahlén, M., & Nadathur, S. 2021, MNRAS, 500, 1953 [Google Scholar]
 Mathews, G. J., Rose, B. M., Garnavich, P. M., Yamazaki, D. G., & Kajino, T. 2016, ApJ, 827, 60 [Google Scholar]
 Maughan, B. J. 2007, ApJ, 668, 772 [NASA ADS] [CrossRef] [Google Scholar]
 Maughan, B. J., Giles, P. A., Randall, S. W., Jones, C., & Forman, W. R. 2012, MNRAS, 421, 1583 [NASA ADS] [CrossRef] [Google Scholar]
 Migkas, K., & Reiprich, T. H. 2018, A&A, 611, A50 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Migkas, K., Schellenberger, G., Reiprich, T. H., et al. 2020, A&A, 636, A15 [CrossRef] [EDP Sciences] [Google Scholar]
 Mittal, R., Hudson, D. S., Reiprich, T. H., & Clarke, T. 2009, A&A, 501, 835 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mittal, R., Hicks, A., Reiprich, T. H., & Jaritz, V. 2011, A&A, 532, A133 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mody, K., & Hajian, A. 2012, ApJ, 758, 4 [NASA ADS] [CrossRef] [Google Scholar]
 Mohayaee, R., Rameez, M., & Sarkar, S. 2020, ArXiv eprints [arXiv:2003.10420] [Google Scholar]
 Mohr, J. J., & Evrard, A. E. 1997, ApJ, 491, 38 [NASA ADS] [CrossRef] [Google Scholar]
 Mohr, J. J., Reese, E. D., Ellingson, E., Lewis, A. D., & Evrard, A. E. 2000, ApJ, 544, 109 [NASA ADS] [CrossRef] [Google Scholar]
 Morandi, A., Ettori, S., & Moscardini, L. 2007, MNRAS, 379, 518 [NASA ADS] [CrossRef] [Google Scholar]
 Nagai, D., Kravtsov, A. V., & Vikhlinin, A. 2007, ApJ, 668, 1 [Google Scholar]
 Nagarajan, A., Pacaud, F., Sommer, M., et al. 2019, MNRAS, 488, 1728 [CrossRef] [Google Scholar]
 Osborne, S. J., Mak, D. S. Y., Church, S. E., & Pierpaoli, E. 2011, ApJ, 737, 98 [NASA ADS] [CrossRef] [Google Scholar]
 Paliathanasis, A., & Leon, G. 2020, Eur. Phys. J. C, 80, 589 [EDP Sciences] [Google Scholar]
 Pavlidou, V., Korkidis, G., Tomaras, T. N., & Tanoglidis, D. 2020, A&A, 638, L8 [CrossRef] [EDP Sciences] [Google Scholar]
 Peery, S., Watkins, R., & Feldman, H. A. 2018, MNRAS, 481, 1368 [Google Scholar]
 Piffaretti, R., Arnaud, M., Pratt, G. W., Pointecouteau, E., & Melin, J.B. 2011, A&A, 534, A109 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration X. 2011, A&A, 536, A10 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XI. 2011, A&A, 536, A11 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XXIX. 2014, A&A, 571, A29 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XXVII. 2016, A&A, 594, A27 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration I. 2020, A&A, 641, A1 [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration III. 2020, A&A, 641, A3 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration Int. XIII. 2014, A&A, 561, A97 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Pratt, C. T., & Bregman, J. N. 2020, ApJ, 890, 156 [Google Scholar]
 Pratt, G. W., Croston, J. H., Arnaud, M., & Böhringer, H. 2009, A&A, 498, 361 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Qin, F., Howlett, C., StaveleySmith, L., & Hong, T. 2019, MNRAS, 482, 1920 [Google Scholar]
 Rameez, M., Mohayaee, R., Sarkar, S., & Colin, J. 2018, MNRAS, 477, 1772 [NASA ADS] [CrossRef] [Google Scholar]
 Reichert, A., Böhringer, H., Fassbender, R., & Mühlegger, M. 2011, A&A, 535, A4 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Reiprich, T. H. 2017, Astron. Nachr., 338, 349 [NASA ADS] [CrossRef] [Google Scholar]
 Reiprich, T. H., & Böhringer, H. 2002, ApJ, 567, 716 [NASA ADS] [CrossRef] [Google Scholar]
 Rossetti, M., Gastaldello, F., Ferioli, G., et al. 2016, MNRAS, 457, 4515 [Google Scholar]
 Rubart, M., Bacon, D., & Schwarz, D. J. 2014, A&A, 565, A111 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Salehi, A., Farajollahi, H., Motahari, M., et al. 2020, Eur. Phys. J. C, 80, 753 [EDP Sciences] [Google Scholar]
 Salehi, A., Yarahmadi, M., & Fathi, S. 2021, MNRAS, 504, 1304 [Google Scholar]
 Schellenberger, G., Reiprich, T. H., Lovisari, L., Nevalainen, J., & David, L. 2015, A&A, 575, A30 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Schlafly, E. F., & Finkbeiner, D. P. 2011, ApJ, 737, 103 [NASA ADS] [CrossRef] [Google Scholar]
 Schlegel, D. J., Finkbeiner, D. P., & Davis, M. 1998, ApJ, 500, 525 [NASA ADS] [CrossRef] [Google Scholar]
 Scrimgeour, M. I., Davis, T. M., Blake, C., et al. 2016, MNRAS, 455, 386 [CrossRef] [Google Scholar]
 Secrest, N. J., von Hausegger, S., Rameez, M., et al. 2021, ApJ, 908, L51 [Google Scholar]
 Shamir, L. 2020, PASA, 37, e053 [Google Scholar]
 Shanks, T., Metcalfe, N., Chehade, B., et al. 2015, MNRAS, 451, 4238 [Google Scholar]
 Shanks, T., Hogarth, L. M., & Metcalfe, N. 2019a, MNRAS, 484, L64 [NASA ADS] [CrossRef] [Google Scholar]
 Shanks, T., Hogarth, L. M., Metcalfe, N., & Whitbourn, J. 2019b, MNRAS, 490, 4715 [NASA ADS] [CrossRef] [Google Scholar]
 Siewert, T. M., SchmidtRubart, M., & Schwarz, D. J. 2020, ArXiv eprints [arXiv:2010.08366] [Google Scholar]
 Skrutskie, M., Cutri, R., Stiening, R., et al. 2006, AJ, 131, 1163 [NASA ADS] [CrossRef] [Google Scholar]
 Soltis, J., Farahi, A., Huterer, D., & Liberato, C. M. 2019, Phys. Rev. Lett., 122, 091301 [Google Scholar]
 Spallicci, A. D. A. M., HelayëlNeto, J. A., LópezCorredoira, M., & Capozziello, S. 2021, Eur. Phys. J. C, 81, 4 [EDP Sciences] [Google Scholar]
 Springob, C. M., Magoulas, C., Colless, M., et al. 2014, MNRAS, 445, 2677 [Google Scholar]
 Sun, Z. Q., & Wang, F. Y. 2019, Eur. Phys. J. C, 79, 783 [CrossRef] [EDP Sciences] [Google Scholar]
 Tiwari, P., Kothari, R., Naskar, A., NadkarniGhosh, S., & Jain, P. 2015, Astropart. Phys., 61, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Tsaprazi, E., & Tsagas, C. G. 2020, Eur. Phys. J. C, 80, 757 [EDP Sciences] [Google Scholar]
 Tully, R. B., Pomarède, D., Graziani, R., et al. 2019, ApJ, 880, 24 [NASA ADS] [CrossRef] [Google Scholar]
 Verde, L., Kamionkowski, M., Mohr, J. J., & Benson, A. J. 2001, MNRAS, 321, L7 [NASA ADS] [CrossRef] [Google Scholar]
 Vikhlinin, A., Burenin, R. A., Ebeling, H., et al. 2009, ApJ, 692, 1033 [NASA ADS] [CrossRef] [Google Scholar]
 Watkins, R., & Feldman, H. A. 2015, MNRAS, 447, 132 [NASA ADS] [CrossRef] [Google Scholar]
 Watkins, R., Feldman, H. A., & Hudson, M. J. 2009, MNRAS, 392, 743 [NASA ADS] [CrossRef] [Google Scholar]
 Whitbourn, J. R., & Shanks, T. 2014, MNRAS, 437, 2146 [NASA ADS] [CrossRef] [Google Scholar]
 Willingale, R., Starling, R. L. C., Beardmore, A. P., Tanvir, N. R., & O’Brien, P. T. 2013, MNRAS, 431, 394 [NASA ADS] [CrossRef] [Google Scholar]
 Wright, E. L. 1979, ApJ, 232, 348 [NASA ADS] [CrossRef] [Google Scholar]
 Wright, E. L., Eisenhardt, P. R., Mainzer, A. K., et al. 2010, AJ, 140, 1868 [NASA ADS] [CrossRef] [Google Scholar]
 York, D. G., Adelman, J., Anderson, J. E., Jr, et al. 2000, AJ, 120, 1579 [Google Scholar]
 Zhang, Y.Y., Reiprich, T. H., Schneider, P., et al. 2017, A&A, 599, A138 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Zhao, D., Zhou, Y., & Chang, Z. 2019, MNRAS, 486, 5679 [CrossRef] [Google Scholar]
 Zitrin, A., Bartelmann, M., Umetsu, K., Oguri, M., & Broadhurst, T. 2012, MNRAS, 426, 2944 [NASA ADS] [CrossRef] [Google Scholar]
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 Xray 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 twodimensional 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 nonrelativistic 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 bandpasscorrected 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 tSZmeasured center of each cluster by the position of the brightest pixel within a 15′ × 15′ box around the Xray coordinates. Using the normalized cluster template, the value of is converted to Y_{5R500} by
where y(θ) is the yparameter profile of the cluster that has been normalized to unity, θ is the radial angular coordinate in the plane of the sky, and θ_{500} is the apparent R_{500} in arcmin. The uncertainties of θ_{500} (which are not provided by MCXC) were ignored since their effect is negligible. The θ_{500} values were calculated through the Xray luminosity–mass relation of Arnaud et al. (2010) and the angular diameter distance to each cluster. The only meaningful contribution to the θ_{500} uncertainties comes from the scatter of the scaling relation. However, owing to the shallow dependence of R_{500} to L_{X} (), a fluctuation of L_{X} within the scatter would only cause a minor shift in θ_{500} (see M20). Additionally, as shown in Fig. E.5, even shifts of 6% in R_{500} would not cause a major change in Y_{5R500}, certainly less than the statistical uncertainties. Thus, it is reasonable to ignore any uncertainties of θ_{500} 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 nonrelativistic 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 Xray 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 nonrelativistic 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} = θ_{500}/c_{500}, where c_{500} 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 XMMNewton 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^{14}. 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.
Fig. A.1. Comparison between the default Y_{5R500} values in this work based on the method of E19 and the Y_{5R500} values as given in the PSZ2 catalog for the 566 clusters in common. The equality line is shown in green. 
A.2. Details on obtaining the nearinfrared 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_{500} from the Xray 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 kcorrection. The kcorrection 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 kcorrection code from Blanton (Blanton & Roweis 2007).
A.3. Details on the Xray determined N_{H, Xray}
To determine N_{H,Xray}, we used the exact same spectral fitting procedure as described in M20 when the Xray 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_{500} region and the 0.2 to 0.5 R_{500} 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 XMMNewton clusters. When we performed the N_{H,Xray} − N_{Htot} comparison, we saw that the data residuals compared to the bestfit 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 Xray absorption taking place toward that direction, it would correspond to an extra N_{Htot} ∼ 7.9 ± 3.4 × 10^{20} cm^{−2}. 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 Xray 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 Xray 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^{15}. 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_{0} = 61.3 ± 4.2 km s^{−1} Mpc^{−1} toward (l, b)∼(318°, −9°) and H_{0} = 72.2 ± 2.4 km s^{−1} Mpc^{−1} for the rest of the sky. The obtained H_{0} 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^{−1} 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^{−1} 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.
Fig. B.1. Normalization anisotropy maps and the respective statistical significance maps of the anisotropies for the L_{bol} − Y_{SZ} (top), the L_{bol} − T (middle), and the Y_{SZ} − T (bottom) scaling relations for the ACC sample. The L_{bol} − Y_{SZ} traces only unaccounted for Xray absorption effects, while the L_{X} − T and Y_{SZ} − T behavior mirrors cosmological anisotropies and BFs. 
B.3. The Y_{SZ} − T relation
We repeat the Y_{SZ} − T anisotropy analysis with only the 113 ACC clusters with Y_{SZ}S/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_{0}, we obtain H_{0} = 60.6 ± 3.6 km s^{−1} Mpc^{−1} toward the most anisotropic region and H_{0} = 73.4 ± 1.9 km s^{−1} Mpc^{−1} 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^{−1} 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^{−1} 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_{0} 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 Xray absorption effects. The latter can only affect the determination of R_{500}, which is based on the absorptionaffected flux. However, , so a 20% bias in L_{X} would only cause a < 4% bias in R_{500}, 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 halflight 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.
Fig. C.1. 3σ (99.7%) parameter space of the normalization and slope of the R − T relation for relaxed (purple) and disturbed (green) 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.
Fig. C.2. Normalization anisotropy map for R − L_{X} (left) and R − Y_{SZ} (right). 
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^{5} 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.
Fig. C.3. Top: correlation between the bestfit normalization A_{RT} (over the full bestfit value of the full sample) and the median XBO for every of the 10^{5} bootstrap subsamples. Bottom: A_{RT} anisotropy map after calibrating for the existing correlation with the XBO. 
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 bestfit slope, a very weak dependence of the normalization on the angular diameter distance remains (). In detail, a 15% variation in H_{0} would only lead to a < 3% variation in A_{YLLBCG}. Similarly, a 1000 km s^{−1} BF at z = 0.05 would only cause a 0.5% change in A_{YLLBCG}. It is clear then that only samplerelated 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.
Fig. C.4. Normalization anisotropy map for the Y_{SZ} − L_{BCG} relation. The statistical significance of the observed anisotropies is ≤1.4σ, and thus the relation is statistically isotropic. 
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 bestfit 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 Xray 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.
Fig. D.1. 3σ (99.7%) parameter space of the normalization and slope of the L_{X} − Y_{SZ} (top), Y_{SZ} − T (middle), and Y_{SZ} − L_{BCG} (bottom) relations for S/N > 2 (purple), S/N > 3 (green), and S/N > 4.5 (black). 
Fig. D.2. L_{X} residuals for the L_{X} − Y_{SZ} relation for S/N > 2 (top) and S/N > 4.5 (bottom) as a function of z (left), N_{Htot} (middle), and Y_{SZ} (right). The green strips indicate the bestfit line with their 1σ uncertainties. It is evident that for S/N > 4.5, the slope is always consistent with zero, indicating no dependence of the residuals on the cluster physical properties. 
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 Xray 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.
Fig. D.3. Statistical significance map of the L_{X} − Y_{SZ} anisotropies for S/N > 2. 
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 bestfit 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.
Fig. D.4. Normalization anisotropy map of the Y_{SZ} − T relation for S/N > 3. 
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_{SZ}S/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^{2} and C_{X} = 6 keV. For the bestfit 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_{SZ}S/N > 4.5) as a result of the different adopted R_{500} 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.
Fig. D.5. Top: Y_{SZ} − T relation when using the PSZ2 Y_{SZ} values, together with its 1σ bestfit function (green). Bottom: normalization anisotropy map of the Y_{SZ} − T relation when the PSZ2 values for Y_{SZ} are used. 
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_{0} = 64.6 ± 2.5 km s^{−1}, 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 bestfit 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.
Fig. D.6. Top: L_{BCG} residuals of the L_{BCG} − T fit, as a function of the BCG redshift. The green stripe corresponds to the bestfit function within 1σ. Middle: 3σ (99.7%) parameter space of the normalization and slope of the L_{BCG} − T relation, for all clusters (green), and for 0.03 < z < 0.15 clusters (purple). Bottom: normalization anisotropy map of the L_{BCG} − T relation when all clusters are considered independent of their redshift. 
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 biasfree 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.
Fig. E.1. Top: distribution of the Galactic latitude of the most anisotropic regions as detected in the 10 000 isotropic simulated samples for the L_{X} − T relation. The distance weighting during the L_{X} − T fitting was used. The bins close to the ZoA cover a larger portion of the sky and naturally more anisotropies are expected to be detected there even if no bias existed. 
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.
Fig. E.2. Number of clusters per θ = 75° cone for the 313 clusters used in the L_{X} − T relation. 
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^{7} iterations of the chain, with a burn in period of 10^{4}. 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^{4} 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 bestfit normalization. As discussed before, their anticorrelation is strong, which generally cancels out any strong effects in the normalization owing to the fluxlimited 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.
Fig. E.3. Top: distribution of the galactic latitude of the most anisotropic regions as detected in the 10 000 isotropic simulated samples for the L_{X} − T relation. The distance weighting during the L_{X} − TT fitting was used. The bins close to the ZoA cover a larger portion of the sky and naturally more anisotropies are expected to be detected there even if no bias existed. 
E.3. Sky fraction of Chandra and XMMNewton clusters in our sample
In Fig. E.4, we show the fraction of the clusters for which Chandra or XMMNewton were used to determine T. Specifically, the fraction is defined as the Chandra clusters minus the XMMNewton clusters over the sum. The results are discussed in Sect. 8.5.
Fig. E.4. Spatial variation of fraction of Chandra clusters minus XMMNewton clusters over the sum. 
E.4. Effects of bulk flows to adopted apparent R_{500}
If a BF is present, then the distance of a cluster is miscalculated. As a result, its L_{X}, and its mass M_{500}, its physical R_{500} (Mpc), and its apparent θ_{500} (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_{0} anisotropies, θ_{500} remains practically unchanged. That is because the opposite changes of R_{500} and D_{A} cancel each other out. Here, we test if a BF can cause significant changes in the adopted θ_{500} 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)^{2}D_{A}. Then, θ_{500} is written as
We typically find BFs of ∼1000 km s^{−1}. For clusters at z = 0.05, this would cause a ≲6% change in θ_{500}. For z = 0.1, this change would be ≲2%. For T, the percent changes are roughly the same as in θ_{500} (see Appendix B in M20). Therefore, for the BF case we considered, and for lowz 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 θ_{500} 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 θ_{500} 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 θ_{500} (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 θ_{500}.
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 θ_{500} 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 Xray 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.
Fig. E.5. Relative change of the measured Y_{5R500} when the input R_{500} (and θ_{500}) is increased (red) or decreased (blue) by 6%. 
All Tables
Maximum anisotropy direction for every scaling relation, together with the needed H_{0} relative variation to fully explain the anisotropy.
Bestfit BFs for every scaling relation, method, and redshift bin considered in this work.
All Figures
Fig. 1. Best fits of the 10 scaling relations studied in this work. The green area represents the 1σ limits of the bestfit area. From top left panel to the right: L_{X} − Y_{SZ}, L_{X} − L_{BCG}, L_{X} − T, Y_{SZ} − T, L_{BCG} − T, R − L_{X}, R − Y_{SZ}, R − L_{BCG}, R − T, Y_{SZ} − L_{BCG}, R − Y_{SZ}, L_{bol} − Y_{SZ} (ACC), and Y_{SZ} − T (ACC) relations. 

In the text 
Fig. 2. Normalization anisotropy maps (left) and the respective statistical significance maps of the anisotropies (right), for L_{X} − Y_{SZ} (top), joint L_{X} − Y_{SZ} (our sample + ACC, middle), and L_{BCG} − T (bottom). All the maps in this work are shown in a Hammer projection. 

In the text 
Fig. 3. Comparison between the N_{Htot} from W13 and the N_{Htot} constrained by our Xray spectral analysis. The equality line is indicated in purple; the 1σ space of the bestfit function for the region of interest (green, black points) and for the rest of the sky (cyan, red points) are also shown. The reason for the systematic difference is given in the main text. The region of interest is within 45° from (l, b)=(281°, −16°) (top) and from (l, b)=(122° , + 8°) (bottom). 

In the text 
Fig. 4. Top: same as in Fig. 2, for L_{X} − T. Bottom: H_{0} anisotropy map derived from the joint L_{X} − T (our sample+ACC). 

In the text 
Fig. 5. Hubble diagram of galaxy clusters as derived by the L_{X} − T (top) and the Y_{SZ} − T (bottom) relations. The clusters from the most anisotropic region of each scaling relation are shown (blue) together with the clusters from the rest of the sky (red). The bestfit lines are indicated with the same color. 

In the text 
Fig. 6. Top: same as in Fig. 2 for Y_{SZ} − T. Bottom: H_{0} anisotropy map (left) and the respective significance map (right) derived from the joint Y_{SZ} − T (our sample+ACC, bottom). 

In the text 
Fig. 7. H_{0} anisotropy map as derived from the joint analysis of L_{X} − T, Y_{SZ} − T, and L_{BCG} − T relations for both samples. 

In the text 
Fig. 8. Histograms of the statistical significance of the maximum anisotropy as detected in 10 000 isotropic MC simulations for the L_{X} − T, Y_{SZ} − T, L_{BCG} − T, L_{bol} − T (ACC), and Y_{SZ} − T (ACC) scaling relations. The vertical black line represents the results from the real data; the pvalue and the direction are also shown. 

In the text 
Fig. 9. Fractional difference between disturbed and relaxed clusters over all the clusters for every sky patch of the extragalactic sky. 

In the text 
Fig. 10. 3σ (99.7%) parameter space of the normalization and slope of the Y_{SZ} − T (left), L_{BCG} − T (middle), and L_{X} − T (right) relations for relaxed (purple) and disturbed (green) clusters. 

In the text 
Fig. 11. Correlation between the bestfit A_{LT} (over the full sample’s bestfit value) and the median XBO for every of the 10^{5} bootstrap subsamples. 

In the text 
Fig. 12. Correlation between the Y_{SZ} scatter from the Y_{SZ} − T relation and the L_{X} scatter from the L_{X} − T relation. 

In the text 
Fig. 13. Top left: bulk flow amplitude and its 68.3% uncertainty as a function of the redshift radius of the used spherical volumes (or the comoving distance radius). Circles and squares correspond to the MR and MA methods, respectively. The black, red, and green points correspond to the L_{X} − T, Y_{SZ} − T, and L_{BCG} − T relations, respectively. The blue diamonds and triangles correspond to the ACC L_{bol} − T and Y_{SZ} − T relations, respectively. The results above z ≳ 0.16 are expected to be dominated by more local clusters and thus are overestimated. Top right: bulk flow amplitude and its 68.3% uncertainty as a function of the median redshift of each shell (or the comoving distance), together with the standard deviation of the redshift distribution. The color coding is the same as before. The data points with the same color and shape are independent to each other. Bottom: examples of BF directions as found for different redshift bins and methods. The color coding is the same as before. All directions agree with each other within 1σ. 

In the text 
Fig. A.1. Comparison between the default Y_{5R500} values in this work based on the method of E19 and the Y_{5R500} values as given in the PSZ2 catalog for the 566 clusters in common. The equality line is shown in green. 

In the text 
Fig. B.1. Normalization anisotropy maps and the respective statistical significance maps of the anisotropies for the L_{bol} − Y_{SZ} (top), the L_{bol} − T (middle), and the Y_{SZ} − T (bottom) scaling relations for the ACC sample. The L_{bol} − Y_{SZ} traces only unaccounted for Xray absorption effects, while the L_{X} − T and Y_{SZ} − T behavior mirrors cosmological anisotropies and BFs. 

In the text 
Fig. C.1. 3σ (99.7%) parameter space of the normalization and slope of the R − T relation for relaxed (purple) and disturbed (green) clusters. 

In the text 
Fig. C.2. Normalization anisotropy map for R − L_{X} (left) and R − Y_{SZ} (right). 

In the text 
Fig. C.3. Top: correlation between the bestfit normalization A_{RT} (over the full bestfit value of the full sample) and the median XBO for every of the 10^{5} bootstrap subsamples. Bottom: A_{RT} anisotropy map after calibrating for the existing correlation with the XBO. 

In the text 
Fig. C.4. Normalization anisotropy map for the Y_{SZ} − L_{BCG} relation. The statistical significance of the observed anisotropies is ≤1.4σ, and thus the relation is statistically isotropic. 

In the text 
Fig. D.1. 3σ (99.7%) parameter space of the normalization and slope of the L_{X} − Y_{SZ} (top), Y_{SZ} − T (middle), and Y_{SZ} − L_{BCG} (bottom) relations for S/N > 2 (purple), S/N > 3 (green), and S/N > 4.5 (black). 

In the text 
Fig. D.2. L_{X} residuals for the L_{X} − Y_{SZ} relation for S/N > 2 (top) and S/N > 4.5 (bottom) as a function of z (left), N_{Htot} (middle), and Y_{SZ} (right). The green strips indicate the bestfit line with their 1σ uncertainties. It is evident that for S/N > 4.5, the slope is always consistent with zero, indicating no dependence of the residuals on the cluster physical properties. 

In the text 
Fig. D.3. Statistical significance map of the L_{X} − Y_{SZ} anisotropies for S/N > 2. 

In the text 
Fig. D.4. Normalization anisotropy map of the Y_{SZ} − T relation for S/N > 3. 

In the text 
Fig. D.5. Top: Y_{SZ} − T relation when using the PSZ2 Y_{SZ} values, together with its 1σ bestfit function (green). Bottom: normalization anisotropy map of the Y_{SZ} − T relation when the PSZ2 values for Y_{SZ} are used. 

In the text 
Fig. D.6. Top: L_{BCG} residuals of the L_{BCG} − T fit, as a function of the BCG redshift. The green stripe corresponds to the bestfit function within 1σ. Middle: 3σ (99.7%) parameter space of the normalization and slope of the L_{BCG} − T relation, for all clusters (green), and for 0.03 < z < 0.15 clusters (purple). Bottom: normalization anisotropy map of the L_{BCG} − T relation when all clusters are considered independent of their redshift. 

In the text 
Fig. E.1. Top: distribution of the Galactic latitude of the most anisotropic regions as detected in the 10 000 isotropic simulated samples for the L_{X} − T relation. The distance weighting during the L_{X} − T fitting was used. The bins close to the ZoA cover a larger portion of the sky and naturally more anisotropies are expected to be detected there even if no bias existed. 

In the text 
Fig. E.2. Number of clusters per θ = 75° cone for the 313 clusters used in the L_{X} − T relation. 

In the text 
Fig. E.3. Top: distribution of the galactic latitude of the most anisotropic regions as detected in the 10 000 isotropic simulated samples for the L_{X} − T relation. The distance weighting during the L_{X} − TT fitting was used. The bins close to the ZoA cover a larger portion of the sky and naturally more anisotropies are expected to be detected there even if no bias existed. 

In the text 
Fig. E.4. Spatial variation of fraction of Chandra clusters minus XMMNewton clusters over the sum. 

In the text 
Fig. E.5. Relative change of the measured Y_{5R500} when the input R_{500} (and θ_{500}) is increased (red) or decreased (blue) by 6%. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.