Free Access
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/0004-6361/202140296
Published online 01 June 2021

© ESO 2021

1. Introduction

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

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

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

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

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

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

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

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

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

In this work, alongside LX and T, we also measure and use the total integrated Compton parameter YSZ, the half-light radii of the clusters R, and the infrared luminosity LBCG 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 H0 = 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 best-fit parameters for the ten scaling relations are presented. In Sect. 5, we study the anisotropic behavior of scaling relations, which are sensitive only to X-ray absorption effects and not cosmological factors. In Sect. 6, the focal point of this work is presented, namely, the anisotropy of scaling relations that are sensitive to cosmological anisotropies and BFs. In Sects. 8 and 7, the possible systematic biases are discussed, and the comparison of our results to those from isotropic Monte Carlo simulations is presented. Finally, in Sects. 9 and 10 the discussion and conclusions of this work are given.

2. Sample and measurements

As a general basis, we use the eeHIFLUGCS sample. The only exception is the use of the full Meta-Catalog of X-ray detected Clusters of galaxies (MCXC, Piffaretti et al. 2011) for the LX − YSZ relation. This is done to maximize the available clusters since both LX and YSZ are available beyond eeHIFLUGCS. From eeHIFLUGCS, a slightly different cluster subsample is used for different scaling relations, where all available measurements for both quantities that enter the scaling relation are considered; the vast majority of the subsamples naturally overlap. For the four scaling relations that include T, we use the same sample as in M20. In a nutshell, the sample includes 313 galaxy clusters, the vast majority of which are part of the eeHIFLUGCS sample. It is a relatively low redshift sample with median z = 0.075 (z ∈ [0.004, 0.45]), which includes objects with an X-ray flux of fX, 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 X-ray luminosity LX, temperature T, redshift z, R500, and metallicities of the core Zcore and 0.2−0.5 R500 annulus Zout are obtained as described in M20. The only change we make in the sample is the spectral fit of NGC 5846. When Zout 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 LX − T plane. Thus we repeat the spectral fitting with a fixed value of Zout = 0.400 Z (the median of the sample), which returns T = 0.927 ± 0.013 keV.

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

2.1. Total integrated Compton parameter Y5R500

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

(1)

where σT is the Thompson cross section, me is the electron rest mass, c is the speed of light, kB is the Boltzmann constant, ne the electron number density, Te the electron temperature, and l.o.s. denotes the line of sight toward the cluster. The product of kBne(r)Te(r) represents the electron pressure profile Pe(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 R500 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 Y5R500, where Y5R500 = ∫y dΩ, where Ω is the chosen solid angle.

We estimated Y5R500 from the latest version of the 100, 143, 217, 353, 545, and 857 GHz all-sky maps delivered by the Planck High Frequency Instrument (HFI; Planck Collaboration III 2020), which were published during the final data release in 2018 (R3.00; Planck Collaboration I 2020). The Y5R500 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 Y5R500 obtained using standard MMFs derived from the non-relativistic tSZ spectrum without additional kSZ removal and there were completely negligible changes in the results. Further details of the processes used to extract the Y5R500 values are presented in detail in Erler et al. (2019, hereafter E19) and in Appendix A.1.

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

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

(2)

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

2.2. Half-light radius R

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

To measure Rapp the ROSAT All-Sky Survey (RASS) maps were used. The count-rate growth curves were extracted for all eeHIFLUGCS clusters, and additionally for all MCXC clusters with fX, 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 count-rate growth curves (i.e., the boundaries of cluster X-ray emission) were determined. The radii Rapp were then found, and converted to R using the default cosmology. The exact details of this process will be presented in Pacaud et al. (in prep.).

The point spread function (PSF) of the XRT/PSPC imager of ROSAT varies with the off-axis angle and the photon energies. In general, it is ≤1.5′. To avoid strong biases due to PSF smearing, we excluded all the clusters with Rapp ≤ 2′. This left us with 418 cluster measurements with a median Rapp = 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 Rapp values that are rather insensitive to PSF effects because of its better spatial resolution.

2.3. Near-infrared BCG luminosity LBCG

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

For this work, we use the magnitudes coming from the 2MASS catalog for two reasons. First and foremost, 2MASS returns the largest number of available BCGs for our sample. Out of the 387 clusters of eeHIFLUGCS, we detected the BCG in 2MASS for 331 of them. Secondly, the infrared BCG luminosities are not strongly sensitive to extinction effects, minimizing the potential risk of unaccounted absorption biases. We exclude all BCGs with z < 0.03, as they appear to be systematically overluminous based on all LBCG scaling relations, indicating a flattening of the LBCG 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 LBCG 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 LBCG measurements with a median redshift of 0.069.

2.4. X-ray determined NH, Xray

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

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

For the 1430 MCXC clusters not included in the M20 sample, we started from their LX 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) NHtot 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 YSZ 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 YSZS/N < 2. This results in 168 clusters with X-ray luminosity and temperature values; 113 of these clusters have a YSZ measurement with S/N > 2. Cross-checking the results of a completely independent cluster sample with the results from our sample is crucial to understand the origin of the observed anisotropies (e.g., to exclude that sample selection effects may bias our results). The properties of the ACC sample as we use it are already given in detail by M18 and M20. In a nutshell, the sample mainly consists of massive galaxy clusters, spanning across z ∈ [0.009, 0.839], with a median redshift of 0.226 (for the 113 clusters with YSZ measurements). The temperatures are obtained by a single-thermal model for the whole cluster, while the X-ray luminosities Lbol are given in the bolometric band within the R200 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 Lbol for the total absorption similarly to the LX corrections of our sample. The necessary R500 values to measure YSZ were obtained using the mass-temperature scaling relation of Reichert et al. (2011). The latter mostly uses XMM-Newton-derived temperatures, which might differ from the ASCA temperatures in general. Thus, we compared the ACC temperatures with our temperatures for the common clusters between the two samples and applied the necessary calibration factors before calculating R500.

3. Scaling relations

We study ten cluster scaling relations in total, namely the LX − YSZ, LX − T, YSZ − T, LX − LBCG, LBCG − T, R − YSZ, R − LX, R − LBCG, R − T, and YSZ − LBCG relations. The first nine relations are sensitive to either additionally needed X-ray absorption corrections (1), BFs (2), or possible cosmological anisotropies (3). The YSZ − LBCG 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.

Table 1.

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

(3)

where CY 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, AYX is the normalization of the relation, CX is the calibration term for the X quantity, and BYX is the slope of the relation. The calibration terms CY and CX are taken to be close to the median values of Y and X, respectively. These values are shown in Table 2, together with the assumed (self-similar) values of γYX.

Table 2.

Best-fit 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

(4)

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 best-fit scaling relations (i.e., all the T scaling relations). In this case we assume the normalization of the scaling relation to be known and direction-independent. This is a reasonable assumption since this quantity is associated with intrinsic properties of the clusters and as such there is no obvious reason why it should spatially vary. In this case, the free parameters we wish to constrain are the Hubble constant H0 or the BF amplitude and direction uBF, while the slope is treated as a free, nuisance parameter. To do so, we minimize the following equation:

(5)

where D is either the luminosity distance DL (e.g., for LX − T) or the angular diameter distance DA (e.g., for YSZ − T); Dobs is the observed distance given the measurements Y and X and their exact scaling relation Y − X; Dth is the theoretically expected distance based on the cosmological parameters (e.g., H0); the redshift z and the existing BF uBF; σ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 Dobs enters every scaling relation differently. Generally, it is given by

(6)

In this equation, DH0 = 70 is the distance found for H0 = 70 km s−1 Mpc−1, which enters in the calculation of Y. The use of DH0 = 70 in Eq. (6) cancels out this cosmological contribution to Y. For this, we need for LX − T, YSZ − T, R − T, and LBCG − T, respectively. This leaves us with X-ray flux (in the rest frame of the cluster), apparent size in arcmin, Y5R500, and BCG flux respectively. Moreover, Dth is given by

(7)

where the (1 + z)±1 factor depends on if we are considering DL or DA. 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 zobs is the observed redshift (converted to the CMB rest frame), uBF 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 uBF is changed from “−” to “+” to avoid confusion with negative velocities.

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

3.2. Bulk flow detection

We follow two different methods to estimate the best-fit BF. Firstly, we fit the scaling relations to all the clusters, adding a BF component in our fit, as in Eqs. (5) and (7). The free parameters are the direction and amplitude of the BF along with the intrinsic scatter. The parameters AYX and BYX 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 Dth) is affected by the BF, depending on the angle between the direction of the cluster and the BF direction. Since the spatial distribution of clusters is nearly uniform, any applied BF does not significantly affect the best-fit AYX and BYX, but it does affect the scatter. This procedure essentially minimizes the average Y residuals around each scaling relation; that is, we get for a smaller σint, YX; the procedure is thus labeled the “minimum residuals” (MR) method.

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

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

3.3. Parameter uncertainties

To estimate the uncertainties of the best-fit parameters, we use a bootstrap resampling method with replacement. We randomly draw 5000 resamplings from the studied (sub)sample of clusters with the same size, simultaneously constrain their best-fit parameters, and obtain the final distribution of the latter. The quoted parameter uncertainties refer to the 68.3% credible interval defined by the positive and negative 34.1th percentile of the distributions with respect to the best-fit value of the full (sub)sample. Each parameter value distribution is considered separately from the rest5. 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 best-fit parameters with their uncertainties, as described in the two previous sections. We then create a color-coded full sky map based on the best-fit parameters of every direction.

There are only two differences with the method followed in M20. Firstly, in this work we leave the slope BYX free to vary, instead of fixing it to its best-fit value for the full sample. That way, the parameter of interest is marginalized over BYX. In M20 we demonstrated that this choice does not significantly alter our results compared to the case in which BYX 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 however6, 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 H0 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 best-fit values of the fitted parameters together with their uncertainties for every direction (cone) in the sky. We then identify the region that shares the most statistically significant anisotropy from the rest of the sky (similar to a dipole anisotropy)7. We assess the statistical significance of the deviation between them by

(8)

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

3.6. Monte Carlo simulations

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

3.7. Summary of statistical improvements compared to M20

In this section, we summarize the improvements in the statistical analysis of this work compared to M20: (1) The parameter uncertainties at every stage of this work are found by bootstrap resampling instead of Δχ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 best-fit values. (3) The statistical weighting of clusters during the anisotropy searching is relaxed to avoid creating any artificial anisotropies. (4) Uncertainties of the anisotropy directions are provided. (5) Monte Carlo simulations are carried out to further assess the statistical significance of the anisotropies.

4. General behavior of the ten scaling relations

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

thumbnail Fig. 1.

Best fits of the 10 scaling relations studied in this work. The green area represents the 1σ limits of the best-fit area. From top left panel to the right: LX − YSZ, LX − LBCG, LX − T, YSZ − T, LBCG − T, R − LX, R − YSZ, R − LBCG, R − T, YSZ − LBCG, R − YSZ, Lbol − YSZ (ACC), and YSZ − T (ACC) relations.

4.1. The LX − T relation

The full analysis of the LX − T relation for our sample is presented in detail in M20. The only changes compared to the M20 results are the slightly changed T for NGC 5846 (see Sect. 2), and the use of bootstrap resampling for estimating the parameter uncertainties. Accounting for these changes, the best-fit values for the LX − 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 Lbol − 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 Lbol axis, using bootstrap for estimating the parameter uncertainties. Since the bolometric X-ray luminosity is used for ACC (within R200), both ALT and BLT are larger than the results of our sample. Also, σint is ∼38% larger for ACC.

4.2. The LX − YSZ relation

The LX − YSZ scaling relation has been studied in the past (e.g., Morandi et al. 2007; Planck Collaboration X 2011; Planck Collaboration XI 2011; De Martino & Atrio-Barandela 2016; Ettori et al. 2020; Pratt & Bregman 2020) mostly using YSZ from Planck and LX 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 LX − YSZ having the lowest scatter among all the scaling relations used in this study.

As mentioned in Sect. 2.1, we measured YSZ for 1095 MCXC clusters with S/N > 2. In studying the LX − YSZ 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 NHtot or high z tend to be significantly fainter on average than clusters with high NHtot or low z, respectively. Surprisingly, the same behavior persists even when the original MCXC LX values and the YSZ 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 YSZ 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 LX − YSZ solutions and no systematic behaviors are observed. Additionally, the intrinsic scatter of the LX − YSZ relation decreases drastically compared to cases with lower S/N thresholds, allowing us to put precise constraints on the best-fit parameters and the possibly observed anisotropies of the relation. For these 460 clusters, the best-fit values are in full agreement with past studies within the uncertainties. The scatter we obtain however is lower than most past studies, most probably resulting from the use of the same R500 between LX and YSZ in our analysis. The observed slope BLY ∼ 0.93 is slightly larger than the self-similar prediction of BLY = 0.8.

For ACC, no trends for the YSZ 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 LX − YSZ 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 LX − LBCG relation

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

4.4. The R − LX relation

The relation between R and LX 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 − LX 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 R500 − LX relation (self-similar prediction of γ ∼ 1.3) is not expected since a redshift evolution between the half-light radius R and R500 is also expected (e.g., owing to the time-varying, CC cluster fraction).

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

Finally, the slope is lower than the self-similar prediction for the R500 − LX 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 YSZ − T relation

The YSZ − 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 YSZ measurements with S/N > 2 from the M20 sample. No systematic differences in the YSZ − T relation are observed for different cluster subsamples with different properties. The YSZ residuals remain consistent with zero with increasing T, z, NHtot, and other cluster parameters, while AYT, BYT 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 YSZ − T relation (see Appendix D.2 for more details). The best-fit parameter values that we obtain for these clusters are in line with previous findings. The value of the slope agrees with the self-similar prediction (BYT = 2.5), while the scatter is lower than that for the LX − T relation. Finally, a single power law perfectly describes the relation since a change in AYT, all and BYT, 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 R500).

4.6. The YSZ − LBCG relation

The YSZ − LBCG 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 YSZ 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 YSZ S/N, with high YSZ S/N clusters tending to be up-scattered. This behavior persists even with an increasing S/N threshold. Although the best-fit BYLBCG stays unchanged for different S/N cuts, the AYLBCG, 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 YSZ − LBCG relation offers only limited insight into the anisotropy analysis, we adopt S/N > 2, but suggest caution because of the aforementioned dependence of AYLBCG. The slope lies close to linearity, while the scatter is ∼10% smaller than the LX − LBCG scaling relation.

4.7. The R − YSZ relation

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

4.8. The LBCG − T relation

Similar to the LX − LBCG relation, the LBCG − T relation has not been the focus of many studies in the past. Bharadwaj et al. (2014) studied the scaling of LBCG with the mass of clusters. The latter is however obtained through a cluster mass–temperature scaling relation. For the LBCG − 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 LBCG − T relation and the other simultaneously fit parameters. As a result, from the 259 clusters for which we measured T and LBCG, 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 LBCG − T relation is considerably large, namely ∼70% larger than the LX − 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 LX − T scaling (∼2.8). The slope has also a larger relative uncertainty (∼10%) than the normalization (∼5%).

4.9. The R − LBCG relation

The R − LBCG 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 LBCG redshift cuts. The intrinsic scatter of the relation is similar to the other R relations, while the slope is slightly larger than the R − LX and R − YSZ relations. The same behavior for the R residuals is observed as in the R − LX and R − YSZ scaling relations.

4.10. The R − T relation

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

5. Anisotropies due to unaccounted X-ray absorption effects

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

5.1. LX − YSZ anisotropies

5.1.1. Our sample

The anisotropies of the LX − YSZ relation are the focal point of the search for hidden X-ray effects that we were not aware of in the past. This relation exhibits the lowest scatter and the largest number of clusters among all scaling relations studied in this work. This allows for precise pinpointing of its anisotropies. Even more importantly, LX − YSZ is almost completely insensitive to any spatial H0 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 H0, no significant change in the best-fit ALYBLY would be observed. Analytically, as and , their ratio (considering their best-fit BLY = 0.928) would be

(9)

Thus, a ∼15% spatial variation of H0 would cause a nondetectable ∼2% variation in ALY. Moreover, the LX − YSZ 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 ALY of this region.

Therefore, any statistically significant anisotropies in the LX − YSZ relation should mainly be caused by unaccounted X-ray absorption effects acting on LX. 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 ALY. 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 NHtot needed to explain this discrepancy is ∼3.9 ± 1.1 × 1020 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 ALY and the anisotropy significance maps of LX − YSZ are represented in Fig. 2. We can see that there are no anomalously bright regions. This indicates that there are no regions with significantly lower-than-solar metallicities of the Galactic material. Assuming availability of a much larger number of clusters with LX and YSZ measurements, ideally extending to low Galactic latitudes, this scaling relation could potentially be used as a new probe of the ISM metallicity.

thumbnail Fig. 2.

Normalization anisotropy maps (left) and the respective statistical significance maps of the anisotropies (right), for LX − YSZ (top), joint LX − YSZ (our sample + ACC, middle), and LBCG − 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 LX − 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 LX − T anisotropies are not caused by unmodeled Galactic effects.

5.1.2. Joint analysis of the LX − YSZ relation for our sample and ACC

If the LX − YSZ 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 Lbol, the flux of the ACC clusters was initially measured within the 0.5–2 keV energy range, and thus is sensitive to X-ray absorption effects. Therefore, jointly analyzing the two samples should provide us with better insight into any possible X-ray absorption issues.

To combine the results of the two independent cluster samples, we perform a joint likelihood analysis. The applied method is that followed in M20 (Sect. 8.2), where we determine the overall apparent H0 variation. In this work, the joint parameter is the LX − YSZ normalization for every region over the normalization of the full sample (ALY/ALY, all), marginalized over the slope. We extract the posterior likelihood of ALY/ALY, 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 NHtot ∼ 3.5 ± 1.7 × 1020 cm−2 would be needed, or alternatively a Z ∼ 1.49 ± 0.23 Z for the already-detected Galactic gas and dust.

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

5.2. LX − LBCG anisotropies

Following the same reasoning as for LX − YSZ, the LX − LBCG scaling relation cannot detect cosmological anisotropies or BFs, since both quantities depend on DL and the slope is close to unity. Also, LBCG is not expected to suffer from any excess extinction, since the near-infrared Ks filter of 2MASS shows a nearly negligible sensitivity on extinction (Schlafly & Finkbeiner 2011), and the LBCG values were already corrected for the known extinction effects. Consequently, the only origin of any observed (statistically significant) anisotropies should be a stronger true X-ray absorption than the adopted value, affecting LX. 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 LX − LBCG are shown in the bottom panel of Fig. 2.

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

5.3. Comparison between NH, Xray and NHtot

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

Firstly, we consider the region within 45° around (l, b)=(281°, −16°). This is the most anisotropic and faintest region of the LX − T relation as found in M20, when only our sample is considered. For the 16 clusters lying within this region We find that the best-fiT relation is cm−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 NHtot clusters have a larger, true NH,Xray, then the clusters of the tested region would be less biased, since they have a larger median NHtot, compared to the rest of the sky. Consequently, there is no indication that an untraced, excess X-ray absorption is the cause behind the LX − T anisotropies.

thumbnail Fig. 3.

Comparison between the NHtot from W13 and the NHtot constrained by our X-ray spectral analysis. The equality line is indicated in purple; the 1σ space of the best-fit 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 X-ray 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 X-ray absorption correction.

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

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

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

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

Finally, using the LX − YSZ relation for our sample and ACC, we identified a region, (l, b)∼(122° , + 8°), which might imply that an extra soft X-ray absorption correction is needed for the clusters there. However, the statistical significance of this result is only ∼2σ, while the LX − LBCG and the NH,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 LX values of the clusters lying within this region.

6. Anisotropies due to bulk flows or H0 variations

We have established that strong biases in the X-ray measurements related to previously uncalibrated absorption do not seem to exist, particularly toward (l, b)∼(300°, −20°). The anisotropies observed in M20 thus cannot arise from such issues. Two more possible origins of such anisotropies are BF motions or spatial variations of cosmological parameters. The first adds a systematic, non-negligible, local velocity component on the measured redshift of the objects toward a particular direction. If not accounted for, it leads to a systematic over- or underestimation of the distances of the clusters, and hence of their cosmology-dependent quantities. The same is true if the real underlying values of the cosmological parameters, for example, H0, vary from region to region. This spatial H0 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 H0 variation to explain the observations.

6.1. The LX − T relation

Our inference of LX 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 LX values of the clusters across the sky based on their T and the globally calibrated LX − T and attribute any directionally systematic deviations to BFs or H0 anisotropies.

6.1.1. Our sample

The anisotropies of the LX − 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 ALT/ALT, 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.

thumbnail Fig. 4.

Top: same as in Fig. 2, for LX − T. Bottom: H0 anisotropy map derived from the joint LX − T (our sample+ACC).

Cosmological anisotropies and BFs. We assume that the cause of the observed tension is an anisotropic H0 value in a Universe without any BFs. We would need H0 = 66.0 ± 1.7 km s−1 Mpc−1 toward (l, b)∼(274°, −9°), and H0 = 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.

thumbnail Fig. 5.

Hubble diagram of galaxy clusters as derived by the LX − T (top) and the YSZ − 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 best-fit lines are indicated with the same color.

We now assume that a BF motion is the sole origin of the observed anisotropies, where H0 is isotropic. Based on the MR method, we find a BF of uBF = 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 uBF = 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 uBF = 1170 ± 400 km s−1 toward (l, b)=(262°±52° , + 2°±26°), and uBF = 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 uBF = 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 LX − 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 uBF = 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 uBF = 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 LX − T anisotropies are not subject to a specific redshift bin but consistently extend throughout the z range.

6.1.2. Joint analysis of the LX − 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 H0, similar to Fig. 23 of M20. The results are plotted in Fig. 4. The free parameter H0 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 H0 = 66.2 ± 1.6 km s−1 Mpc−1, while for the rest of the sky we gets H0 = 72.7 ± 1.5 km s−1 Mpc−1. While the posterior H0 range is identical to that found in M20, the significance of the anisotropies do not exceed 3σ, owing to the more conservative methodology followed, along with the non-use of XCS-DR1, in this work. The peak anisotropy is close to the result from our sample since it dominates the joint fit, which is also revealed because the σ level compared to our sample alone increases only slightly. There is a second, weaker peak on the maximum anisotropy position of ACC.

In terms of a BF motion, there is no meaningful way to combine the two independent data sets in an analytical way similar to the H0 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 YSZ − T relation

The anisotropies of the YSZ − T relation are presented in this work for the first time. The quantity YSZ 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 LX − T relation applies. There are three advantages to the YSZ − T relation. Firstly, the scatter is clearly smaller than the LX − T relation allowing for more precise constraints. Secondly, YSZ 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 LX. Thus, in practice no YSZ − T anisotropies can occur from unaccounted Galactic absorption issues. Finally, owing to the applied redshift evolution of YSZ − T, the latter is more sensitive to BFs than LX − T.

6.2.1. Our sample

As a consequence of the low observed scatter and the 260 clusters with quality YSZ 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 LX − 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 YSZ − T scatter and the narrower scanning cones.

This result strongly demonstrates that the observed LX − T anisotropies are not due to X-ray absorption issues. The AYT/AYT, all and the sigma maps are shown in Fig. 6. We should note that there is a mild correlation between the LX and YSZ scatter compared to T due to the physical state of galaxy clusters. If the origin of the anisotropies was sample-related (e.g., a surprisingly strong archival bias), we would expect to see similar anisotropies in LX − T and YSZ − T. However, this would still not explain the large statistical significance of the YSZ − 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 YSZ values from PSZ2 instead, we obtain similar anisotropic results (Appendix D.3).

thumbnail Fig. 6.

Top: same as in Fig. 2 for YSZ − T. Bottom: H0 anisotropy map (left) and the respective significance map (right) derived from the joint YSZ − T (our sample+ACC, bottom).

Cosmological anisotropies and BFs. We now investigate the necessary H0 variations to fully explain the observed anisotropies. We would need H0 = 62.3 ± 2.2 km s−1 Mpc−1 toward (l, b)∼(268°, −16°) and H0 = 72.1 ± 0.9 km s−1 Mpc−1 for the rest of the sky. The result is consistent within 1.4σ with the H0 value obtained from the joint LX − 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 YSZ − 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 YSZ value, and in the redshift evolution of the relation. However, the anisotropies appear to be even stronger now compared to the LX − T relation.

Assuming the apparent anisotropies are due to a BF affecting the entire sample, using the MR method we find uBF = 950 ± 340 km s−1 toward (l, b)=(263°±39°, −22°±20°). For the redshift bin z ∈ [0, 0.07], we obtain uBF = 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 best-fit results remain within the uncertainties without an amplitude decay. For the z ∈ [0.07, 0.12] bin, we obtain uBF = 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 uBF = 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 uBF = 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 uBF = 1200 ± 350 km s−1 toward (l, b)=(254°±22°, −18°±16°). For z > 0.09, the amplitude is reduced (uBF = 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 YSZ − T anisotropies reveal similar BFs than the LX − 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 YSZ − T relation for our sample and ACC

The ACC sample shows a very similar anisotropic YSZ − T behavior to our sample (Appendix B.3). We combine the two independent samples and their 376 different clusters to jointly constrain the apparent H0 anisotropies with the YSZ − T relation. The H0 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 H0 = 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 LX − T results. Finally, a dipole form of the anisotropy is apparent in the sigma maps.

6.3. The LBCG − T relation

The final scaling relation that can potentially trace cosmological anisotropies and BFs is the LBCG − 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 LX − T and YSZ − T constitute the main disadvantages of LBCG − 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 LBCG − 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 LX − T and YSZ − T results offers additional confirmation of the existence of the physical phenomenon causing this.

Cosmological anisotropies and BFs. In terms of H0, we obtain H0 = 66.2 ± 2.5 km s−1 Mpc−1 for that direction and H0 = 72.6 ± 2.3 km s−1 Mpc−1 for the opposite hemisphere. The H0 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 uBF = 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 LX − T and YSZ − T.

6.4. Combined anisotropies of LX − T, YSZ − T, and LBCG − T

We now combine the information from the LX − T, YSZ − T, and LBCG − T anisotropies, for our sample and ACC, into one single H0 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 H0 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 LX, YSZ, and LBCG. This is of course true for cosmological anisotropies and BFs. This also requires that there are no unaccounted X-ray absorption effects that bias the T measurement, which should be the case according to the analysis presented in Sect. 5. If a strong temperature outlier exists due to chance, it would not significantly affect the results by propagating to all scaling relations. The effect of the outlier would be buried under the average behavior of the rest of the clusters in its region, if no reason for an anisotropy exists. The posterior H0 likelihood from the entire region is then combined with the same results from the other scaling relations.

Secondly, a moderately correlated scatter exists between LX and YSZ, which relates to the physical state of galaxy clusters. We discuss this in detail in Sect. 8.6. When the average properties of clusters are similar from sky region to sky region (and any anisotropies arise purely from cosmological effects), this correlated scatter is not expected to bias our combined results. If, however, a strongly inhomogeneous spatial distribution of cluster populations (e.g., CC and non-CC clusters) exist, then the correlation between LX and YSZ 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 H0 spatial variation.

6.4.2. Apparent H0 anisotropy from joint analysis

When LX − T, YSZ − T, and LBCG − T results for both our sample and ACC are combined, we obtain H0 = 66.5 ± 1.0 km s−1 Mpc−1 toward (l, b)=(273°±40°, −11°±27°), and H0 = 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 H0 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 H0 uncertainty from which probes that attempt to put absolute constraints on H0 suffer. However, it is not needed for the relative spatial differences to which we are interested in. The joint H0 map together with the statistical significance map are shown in Fig. 7. The observed anisotropy forms a dipole.

thumbnail Fig. 7.

H0 anisotropy map as derived from the joint analysis of LX − T, YSZ − T, and LBCG − T relations for both samples.

7. Comparison with isotropic Monte Carlo simulations

The well-established statistical methods used up to now provide reliable parameter uncertainties and a robust estimation of the rareness of the observed anisotropies. Nonetheless, biases we are not aware of, can still be present in our analysis and lead to overestimated anisotropy signals. The same is true for cosmic variance. To further investigate this, we need to apply our analysis to simulated isotropic MC simulations. For the LX − T, YSZ − T, and LBCG − 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 LX − T relation. We keep the coordinates, redshifts, and T (together with σlog T) of the simulated clusters fixed to the real values. This is done to incorporate all the possible effects that could create anisotropies, including the spatial distribution of the real clusters. Based on the best-fit scaling relation for the full, real sample, we calculate the predicted LX value. We then add a random offset to the latter, drawn from a log-normal distribution with a standard deviation of . Since low-mass systems might exhibit a larger intrinsic scatter, we consider three different values of . We divide the sample in three equally sized subsamples according to their T value, and for every subsample we constrain σint, LT 9 given the best-fit ALT and BLT 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 LX value. We ensure that the posterior distribution of the best-fit values of the simulated samples follow the input values for every parameter.

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

A simulated, isotropic LBCG value is also drawn for the 196 clusters with 0.03 < z < 0.15 in the same way as for LX − T. The infrared BCG luminosity does not show any correlation in its scatter with LX and YSZ, and therefore no extra procedures are needed. We finally repeat the Lbol − T and YSZ − 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.

thumbnail Fig. 8.

Histograms of the statistical significance of the maximum anisotropy as detected in 10 000 isotropic MC simulations for the LX − T, YSZ − T, LBCG − T, Lbol − T (ACC), and YSZ − T (ACC) scaling relations. The vertical black line represents the results from the real data; the p-value and the direction are also shown.

7.2.1. Our sample

For the LX − 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 YSZ − 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 YSZ − T than for LX − 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 LX − T and YSZ − 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 LX and YSZ values has been taken into account as explained before.

For the LBCG − 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 LX − T and YSZ − 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 LX, YSZ, and LBCG.

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 Lbol − 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 YSZ − T relation yields p = 0.055 alone, which is more statistically significant than Lbol − 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 X-ray and sample-related issues were tested and discussed for the LX − T relation in M20. Here we explore some additional effects that might undermine the significance of our results.

8.1. Cool-core and morphologically relaxed clusters

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

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

thumbnail 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 LX − T, YSZ − T, and LBCG − T

Relaxed and disturbed clusters do not show any meaningful difference in their YSZ − T and LBCG − T normalization, as shown in Fig. 10. Therefore, the possible CC bias is irrelevant for these two scaling relations. For LX − 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 LX − T anisotropy region, hence no considerable bias is expected.

thumbnail Fig. 10.

3σ (99.7%) parameter space of the normalization and slope of the YSZ − T (left), LBCG − T (middle), and LX − T (right) relations for relaxed (purple) and disturbed (green) clusters.

This was shown in M18, where using core-excised LX did not have any effect on the LX − 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 LX − T anisotropies, we create 105 randomly drawn bootstrap subsamples (same process as in M20), independent of direction. We investigate the correlation between the ALT 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 ALT (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 ALT.

thumbnail Fig. 11.

Correlation between the best-fit ALT (over the full sample’s best-fit value) and the median XBO for every of the 105 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

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

The MB correction to be applied in the normalization of the LX − T for low-z, flux-limited samples such as our own is (Vikhlinin et al. 2009). Consequently, for the (l, b)∼(274°, −9°) region we have MBcorr ≈ 0.594 ± 0.048, and for the rest of the sky MBcorr ≈ 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 MBcorr uncertainties, there is only a 0.007% (4σ) probability that the full amplitude of the LX − T anisotropy is a result of the MB.

Similar to the effects of the scatter, we might expect that larger LX uncertainties could lead to larger normalizations and vice versa. However, the faintest region again has a slightly larger average LX 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 LX − 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 YSZ − 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 LX, YSZ, and T (i.e., mass) distributions between low and high z. Low- and high-mass systems might share slightly different scaling laws. If the z distributions of the compared cluster subsamples were significantly different, this could lead to a biased comparison and artificial anisotropies. However, in M20 we showed that the z and T distributions of the (l, b)∼(280°, −15°) are very similar to the rest of the sky, as expected for a homogeneously selected sample. This is further confirmed by our physical properties test in Sect. 8.4. Subsequently, the Eddington bias is rather irrelevant to our analysis.

8.3. Zone of avoidance bias

The fact that the direction of the observed anisotropies mostly lies close to the Galactic plane gap, or the zone of avoidance (ZoA; |b ≤ 20° |) is another issue that might raise concerns about the underlying origin of these anisotropies. In Sect. 5, we demonstrated that uncalibrated X-ray absorption issues are not the reason behind LX − T anisotropies and definitely cannot explain the YSZ − 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 LX − T, YSZ − T, and LBCG − T relations of our sample, this time without applying any weighting on the cluster uncertainties; these weights were based on each the distance of each cluster from the center of each cone. As a result, the scanning algorithm does not take into account the spatial distribution of clusters and does not “see” the ZoA gap. The only information the scanning algorithm reads is the number of clusters in each cone. The cones close to the most anisotropic region happen to have average numbers of objects (Fig. E.2). Even if fewer clusters were included in the cones, in M20 we showed that the best-fit normalization is completely independent of the number of clusters.

When repeating the analysis for LX − 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 YSZ − 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 LBCG − 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 LX − T and YSZ − 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 ALT. In this work, we wish to find out if such a behavior of ALT (or AYT) could also result from a certain combination of average cluster properties. To do so, we construct 106 random bootstrap cluster subsamples of random size (10−60% of the total size of the sample), as done in M20. Except for the best-fit ALT, BLT and σint, LT, 12 more average properties of each subsample are also derived12. We express ALT as a function of all the parameters pi (normalized by their sample mean pmean), and their power-law index vi, as follows:

(10)

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

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

To determine if the observed LX − T anisotropies are caused by the physical cluster properties of the (l, b)∼(274°, −9°), we predict its ALT using Eq. (10) and its mean cluster properties. We obtain ALT ∼ 1.12 and ALT ∼ 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 YSZ − T returns similar results.

8.5. Temperature calibration issues

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

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

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

Past studies (i.e., Nagarajan et al. 2019, and references therein) have reported a positive correlation between the scatter of LX and YSZ 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 LBCG. Considering the 260 clusters with both YSZ and LX measurements, we confirm that their residuals with respect to T are correlated and have a Pearson’s correlation coefficient of rcorr = 0.668 ± 0.089. This is shown in Fig. 12. The best-fit line has a slope of ∼0.63, while the total scatter of the relation is ∼0.16 dex.

thumbnail Fig. 12.

Correlation between the YSZ scatter from the YSZ − T relation and the LX scatter from the LX − 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 YSZ − T and LX − 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 ALT and AYT 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 LBCG with respect to T does not show any meaningful correlation with the scatter of LX or YSZ.

8.7. Selection biases

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

9. Discussion and conclusions

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

9.1. Known systematics do not explain the observed anisotropies

The LX − T anisotropy toward (l, b)∼(280°, −15°) that was originally observed in M20 could in principle be due to uncalibrated X-ray effects, such as excess X-ray absorption toward this sky direction. With our new data and analysis strategy, especially from the LX − YSZ relation, which shows no signs of extra X-ray absorption toward the region in question, we can now decisively exclude this possibility. The lack of excess X-ray absorption effects toward that region is further confirmed by the analysis of ACC and by the LX − LBCG relation of our sample. The comparison between the X-ray-determined NHtot and the W13-based NHtot also does not indicate such an effect. Furthermore, the same anisotropies persist in scaling relations that are insensitive to unaccounted X-ray absorption, such as YSZ − T and LBCG − 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 YSZ − T, LBCG − T, and the ACC sample, which is not flux-limited. It was also shown that the fraction of morphologically relaxed and disturbed clusters does not spatially vary strongly enough to have any effect in the apparent anisotropy of the LX − T relation. Even if it did, this would have no strong effect on the YSZ − 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 H0 anisotropies and BFs, will be employed to better understand the precision and accuracy with which our methodology can detect such phenomena.

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

9.2. Cosmological phenomena behind the tension

For now, there are just two obvious remaining explanations for the persisting cluster scaling relation anisotropies. The first is a spatial variation of H0 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 H0 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.

Table 3.

Maximum anisotropy direction for every scaling relation, together with the needed H0 relative variation to fully explain the anisotropy.

Table 4.

Best-fit BFs for every scaling relation, method, and redshift bin considered in this work.

thumbnail 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 LX − T, YSZ − T, and LBCG − T relations, respectively. The blue diamonds and triangles correspond to the ACC Lbol − T and YSZ − 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 high-z clusters are needed. The upcoming eRASS catalogs will provide us with thousands of such clusters. These distant objects will be of crucial importance when applying the techniques presented in this work, since we will be able to tell if, and at which scale, the cluster behavior converges to isotropy. The kSZ effect, which is caused by the peculiar motion of clusters in the CMB frame, can also be utilized to differentiate cosmological anisotropies from BFs. The limited number of clusters in our sample, together with the low redshifts and the large angular cluster sizes, did not allow us to extract useful information for this effect by stacking filtered kSZ maps of objects that presumably move toward a common direction. The reason for this lies in the large power of the primary CMB anisotropies at the scales that are relevant to our sample as well as the limited angular resolution and sensitivity provided by Planck. Next-generation SZ instruments will significantly improve over Planck’s sensitivity and angular resolution and allow for direct measurements of BFs using the kSZ effect.

9.2.1. Anisotropic Hubble expansion

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

Similar large-scale anisotropies in the distribution of high-z radio and infrared sources have also been found (e.g., Tiwari et al. 2015; Colin et al. 2017; Bengaly et al. 2018; Rameez et al. 2018; Siewert et al. 2020), usually implying a ∼600−1500 km s−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 X-ray point sources, tracing the matter distribution out to very large scales.

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

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

9.2.2. Large-scale bulk flows

The consistent observation of cluster anisotropies across several scaling relations could hint at the existence of large BFs that by far exceed the BF scales predicted by ΛCDM. In Fig. 13, a ≳600 km s−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 high-z clusters are necessary to derive safe conclusions at these scales.

The ΛCDM model predicts negligible BF amplitudes at scales of ≳250 Mpc. Flow motions of similar, large amplitudes have been reported in the past for galaxy clusters (e.g., Lauer & Postman 1994; Hudson et al. 1999; Kashlinsky et al. 2008, 2010; Atrio-Barandela et al. 2015), while Osborne et al. (2011), Mody & Hajian (2012), and Planck Collaboration Int. XIII (2014) fail to see any large-scale BFs in the CMB kSZ data. Studies using thousands of galaxies consistently detect a BF toward the same direction as us, but with a ∼3−5 times smaller amplitude. However, it also seems to extend to scales larger than ∼300 Mpc, up to where the galaxy surveys are usable. For instance, Watkins et al. (2009) and Magoulas et al. (2016) find a ∼400 km s−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 larger-than-expected scales. However, they again stress the limitations of current SNIa samples for such studies and that conclusive answers cannot be given. Other studies have stressed that the usual Newtonian treatment of peculiar velocity fields can significantly underestimate the inferred BFs (e.g., Filippou & Tsagas 2021; Tsaprazi & Tsagas 2020). From all that, it is evident that it is still not very clear if the detected BF motions are entirely in agreement with the standard model and to which scale the BF persists.

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

10. Summary

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

If the observed anisotropies were the result of an anisotropic expansion rate, we would need a ∼9% spatial variation of H0 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.


1

Measured by ROSAT.

3

X does not depend on z and H0 when X = temperature T.

4

Our analysis and conclusions are rather insensitive to a (small) systematic over- or underestimation of σint, YX. This was tested by repeating our anisotropy analysis using 20% smaller or larger σint, YX.

5

Thus, this process is equivalent to marginalizing one parameter over the rest.

6

Only when there are strongly up- or down-scattered clusters close to the center of the cone with small measurement uncertainties, and the number of clusters within the cone is relatively small.

7

In M20 we identified the two regions with the most extreme, opposite behaviors. The currently adopted method leads to slightly reduced anisotropy signals and is part of the more conservative approach we follow in this work.

8

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.

9

We find σint, LT = (0.29, 0.23, 0.16) for the T < 3.5 keV, 3.5 keV < T < 5.8 keV, and T > 5.8 keV clusters respectively.

10

Even if we completely ignore the directions of the anisotropies and only consider the amplitudes, the probability still is p = 2.5 × 10−5.

11

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 best-fit 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 best-fit normalization). The statistical significance of this relative difference (i.e., anisotropy) accounts of course for σint through the normalization uncertainties.

12

These are the average redshift, temperature, < 0.2 R500 (core) temperature, flux, core metallicity, 0.2−0.5 R500 metallicity, NHtot, RASS time exposure, X-ray peak-BCG offset, number of clusters, original source catalog (REFLEX/BCS = 1, NORAS = 2), and instrument used for T measurement (XMM-Newton = 1, Chandra = 2).

13

The existing literature is too large to be fully included in this work. Therefore we focus only on the most recent results or on results that were not already mentioned in M20. The latter paper gives many additional studies that find consistent results with this work.

14

The used SZ center is the position of the brightest pixel within a 15 × 15 arcmin2 box centered around the X-ray peak.

15

In M20 we considered θ = 75° because the slope was kept fixed and fewer cluster per regions were needed to sufficiently determine the normalization. In this work the slope is left free to vary.

Acknowledgments

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

References

  1. Abbott, T., Abdalla, F., Allam, S., et al. 2018, ApJS, 239, 18 [Google Scholar]
  2. Allen, S. W., Evrard, A. E., & Mantz, A. B. 2011, ARA&A, 49, 409 [NASA ADS] [CrossRef] [Google Scholar]
  3. Andrade, U., Bengaly, C. A. P., Alcaniz, J. S., & Capozziello, S. 2019, MNRAS, 490, 4481 [NASA ADS] [CrossRef] [Google Scholar]
  4. Appleby, S., Shafieloo, A., & Johnson, A. 2015, ApJ, 801, 76 [NASA ADS] [CrossRef] [Google Scholar]
  5. 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]
  6. Arnaud, M., Pratt, G. W., Piffaretti, R., et al. 2010, A&A, 517, A92 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  7. Asplund, M., Grevesse, N., Sauval, A. J., & Scott, P. 2009, ARA&A, 47, 481 [NASA ADS] [CrossRef] [Google Scholar]
  8. Atrio-Barandela, F., Kashlinsky, A., Ebeling, H., Fixsen, D. J., & Kocevski, D. 2015, ApJ, 810, 143 [NASA ADS] [CrossRef] [Google Scholar]
  9. Bender, A. N., Kennedy, J., Ade, P. A. R., et al. 2016, MNRAS, 460, 3432 [Google Scholar]
  10. Bengaly, C. A. P., Maartens, R., & Santos, M. G. 2018, J. Cosmol. Astropart. Phys., 2018, 031 [Google Scholar]
  11. Bengaly, C. A. P., Maartens, R., Randriamiarinarivo, N., & Baloyi, A. 2019, J. Cosmol. Astropart., 2019, 025 [Google Scholar]
  12. Betoule, M., Kessler, R., Guy, J., et al. 2014, A&A, 568, A22 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  13. Bharadwaj, V., Reiprich, T. H., Schellenberger, G., et al. 2014, A&A, 572, A46 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  14. Bharadwaj, V., Reiprich, T. H., Lovisari, L., & Eckmiller, H. J. 2015, A&A, 573, A75 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  15. Blanton, M. R., & Roweis, S. 2007, AJ, 133, 734 [Google Scholar]
  16. Bleem, L. E., Stalder, B., de Haan, T., et al. 2015, ApJS, 216, 27 [NASA ADS] [CrossRef] [Google Scholar]
  17. Bleem, L. E., Bocquet, S., Stalder, B., et al. 2020, ApJS, 247, 25 [CrossRef] [Google Scholar]
  18. Böhringer, H., Chon, G., & Collins, C. A. 2020, A&A, 633, A19 [EDP Sciences] [Google Scholar]
  19. Bolejko, K., Nazer, M. A., & Wiltshire, D. L. 2016, J. Cosmol. Astropart., 2016, 035 [Google Scholar]
  20. Boruah, S. S., Hudson, M. J., & Lavaux, G. 2020, MNRAS, 498, 2703 [Google Scholar]
  21. Carrick, J., Turnbull, S. J., Lavaux, G., & Hudson, M. J. 2015, MNRAS, 450, 317 [Google Scholar]
  22. Chang, Z., & Lin, H.-N. 2015, MNRAS, 446, 2952 [CrossRef] [Google Scholar]
  23. Chilingarian, I. V., & Zolotukhin, I. Y. 2011, MNRAS, 419, 1727 [Google Scholar]
  24. Chluba, J., Nagai, D., Sazonov, S., & Nelson, K. 2012, MNRAS, 426, 510 [NASA ADS] [CrossRef] [Google Scholar]
  25. Colin, J., Mohayaee, R., Sarkar, S., & Shafieloo, A. 2011, MNRAS, 414, 264 [NASA ADS] [CrossRef] [Google Scholar]
  26. Colin, J., Mohayaee, R., Rameez, M., & Sarkar, S. 2017, MNRAS, 471, 1045 [NASA ADS] [CrossRef] [Google Scholar]
  27. Colin, J., Mohayaee, R., Rameez, M., & Sarkar, S. 2019, A&A, 631, L13 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  28. Dai, D.-C., Kinney, W. H., & Stojkovic, D. 2011, J. Cosmol. Astropart., 2011, 015 [Google Scholar]
  29. Das, K. K., Sankharva, K., & Jain, P. 2021, J. Cosmol. Astropart., submited [arXiv:2101.11016] [Google Scholar]
  30. De Martino, I., & Atrio-Barandela, F. 2016, MNRAS, 461, 3222 [Google Scholar]
  31. Deng, H.-K., & Wei, H. 2018, Eur. Phys. J. C, 78, 755 [CrossRef] [EDP Sciences] [Google Scholar]
  32. Dąbrowski, M. P., & Wagner, F. 2020, Eur. Phys. J. C, 80, 676 [EDP Sciences] [Google Scholar]
  33. Eckert, D., Molendi, S., & Paltani, S. 2011, A&A, 526, A79 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  34. Erler, J., Basu, K., Chluba, J., & Bertoldi, F. 2018, MNRAS, 476, 3360 [NASA ADS] [CrossRef] [Google Scholar]
  35. Erler, J., Ramos-Ceja, M. E., Basu, K., & Bertoldi, F. 2019, MNRAS, 484, 1988 [Google Scholar]
  36. Ettori, S., Lovisari, L., & Sereno, M. 2020, A&A, 644, A111 [EDP Sciences] [Google Scholar]
  37. Feindt, U., Kerschhaggl, M., Kowalski, M., et al. 2013, A&A, 560, A90 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  38. Filippou, K., & Tsagas, C. G. 2021, Ap&SS, 366, 4 [Google Scholar]
  39. Fitzpatrick, E. L. 1999, PASP, 111, 63 [NASA ADS] [CrossRef] [Google Scholar]
  40. Fosalba, P., & Gaztañaga, E. 2021, MNRAS, 504, 5840 [Google Scholar]
  41. Furnell, K. E., Collins, C. A., Kelvin, L. S., et al. 2018, MNRAS, 478, 4952 [NASA ADS] [CrossRef] [Google Scholar]
  42. Harrison, E. R. 1974, ApJ, 191, L51 [NASA ADS] [CrossRef] [Google Scholar]
  43. Haslbauer, M., Banik, I., & Kroupa, P. 2020, MNRAS, 499, 2845 [Google Scholar]
  44. Hasselfield, M., Hilton, M., Marriage, T. A., et al. 2013, J. Cosmol. Astropart., 2013, 008 [Google Scholar]
  45. Hilton, M., Sifón, C., Naess, S., et al. 2021, ApJS, 253, 3 [Google Scholar]
  46. Hoffman, Y., Courtois, H. M., & Tully, R. B. 2015, MNRAS, 449, 4494 [Google Scholar]
  47. Horner, D. J. 2001, PhD Thesis, University of Maryland College Park, USA [Google Scholar]
  48. Horvath, I., Szécsi, D., Hakkila, J., et al. 2020, MNRAS, 498, 2544 [Google Scholar]
  49. Hu, J. P., Wang, Y. Y., & Wang, F. Y. 2020, A&A, 643, A93 [EDP Sciences] [Google Scholar]
  50. Hudson, M. J., Smith, R. J., Lucey, J. R., Schlegel, D. J., & Davies, R. L. 1999, ApJ, 512, L79 [NASA ADS] [CrossRef] [Google Scholar]
  51. Hudson, M. J., Smith, R. J., Lucey, J. R., & Branchini, E. 2004, MNRAS, 352, 61 [Google Scholar]
  52. Hudson, D. S., Mittal, R., Reiprich, T. H., et al. 2010, A&A, 513, A37 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  53. Huterer, D. 2020, ApJ, 904, L28 [Google Scholar]
  54. Itoh, N., Kohyama, Y., & Nozawa, S. 1998, ApJ, 502, 7 [Google Scholar]
  55. Kaiser, N. 1986, MNRAS, 222, 323 [Google Scholar]
  56. 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]
  57. Kaiser, N., Burgett, W., Chambers, K., et al. 2010, in Ground-based and Airborne Telescopes III, Int. Soc. Opt. Photon., 7733, 77330E [Google Scholar]
  58. Kashlinsky, A., Atrio-Barandela, F., Kocevski, D., & Ebeling, H. 2008, ApJ, 686, L49 [Google Scholar]
  59. Kashlinsky, A., Atrio-Barandela, F., Ebeling, H., Edge, A., & Kocevski, D. 2010, ApJ, 712, L81 [Google Scholar]
  60. Kazantzidis, L., & Perivolaropoulos, L. 2020, Phys. Rev. D, 102, 023520 [CrossRef] [Google Scholar]
  61. Keenan, R. C., Barger, A. J., Cowie, L. L., et al. 2012, ApJ, 754, 131 [NASA ADS] [CrossRef] [Google Scholar]
  62. Korkidis, G., Pavlidou, V., Tassis, K., et al. 2020, A&A, 639, A122 [EDP Sciences] [Google Scholar]
  63. Lauer, T. R., & Postman, M. 1994, ApJ, 425, 418 [NASA ADS] [CrossRef] [Google Scholar]
  64. Lavaux, G., Afshordi, N., & Hudson, M. J. 2013, MNRAS, 430, 1617 [NASA ADS] [CrossRef] [Google Scholar]
  65. Li, M., Pan, J., Gao, L., et al. 2012, ApJ, 761, 151 [Google Scholar]
  66. Lopes, P. A. A., Trevisan, M., Laganá, T. F., et al. 2018, MNRAS, 478, 5473 [NASA ADS] [CrossRef] [Google Scholar]
  67. Lovisari, L., Schellenberger, G., Sereno, M., et al. 2020, ApJ, 892, 102 [CrossRef] [Google Scholar]
  68. 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]
  69. Manolopoulou, M., Hoyle, B., Mann, R. G., Sahlén, M., & Nadathur, S. 2021, MNRAS, 500, 1953 [Google Scholar]
  70. Mathews, G. J., Rose, B. M., Garnavich, P. M., Yamazaki, D. G., & Kajino, T. 2016, ApJ, 827, 60 [Google Scholar]
  71. Maughan, B. J. 2007, ApJ, 668, 772 [NASA ADS] [CrossRef] [Google Scholar]
  72. Maughan, B. J., Giles, P. A., Randall, S. W., Jones, C., & Forman, W. R. 2012, MNRAS, 421, 1583 [NASA ADS] [CrossRef] [Google Scholar]
  73. Migkas, K., & Reiprich, T. H. 2018, A&A, 611, A50 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  74. Migkas, K., Schellenberger, G., Reiprich, T. H., et al. 2020, A&A, 636, A15 [CrossRef] [EDP Sciences] [Google Scholar]
  75. Mittal, R., Hudson, D. S., Reiprich, T. H., & Clarke, T. 2009, A&A, 501, 835 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  76. Mittal, R., Hicks, A., Reiprich, T. H., & Jaritz, V. 2011, A&A, 532, A133 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  77. Mody, K., & Hajian, A. 2012, ApJ, 758, 4 [NASA ADS] [CrossRef] [Google Scholar]
  78. Mohayaee, R., Rameez, M., & Sarkar, S. 2020, ArXiv e-prints [arXiv:2003.10420] [Google Scholar]
  79. Mohr, J. J., & Evrard, A. E. 1997, ApJ, 491, 38 [NASA ADS] [CrossRef] [Google Scholar]
  80. Mohr, J. J., Reese, E. D., Ellingson, E., Lewis, A. D., & Evrard, A. E. 2000, ApJ, 544, 109 [NASA ADS] [CrossRef] [Google Scholar]
  81. Morandi, A., Ettori, S., & Moscardini, L. 2007, MNRAS, 379, 518 [NASA ADS] [CrossRef] [Google Scholar]
  82. Nagai, D., Kravtsov, A. V., & Vikhlinin, A. 2007, ApJ, 668, 1 [Google Scholar]
  83. Nagarajan, A., Pacaud, F., Sommer, M., et al. 2019, MNRAS, 488, 1728 [CrossRef] [Google Scholar]
  84. Osborne, S. J., Mak, D. S. Y., Church, S. E., & Pierpaoli, E. 2011, ApJ, 737, 98 [NASA ADS] [CrossRef] [Google Scholar]
  85. Paliathanasis, A., & Leon, G. 2020, Eur. Phys. J. C, 80, 589 [EDP Sciences] [Google Scholar]
  86. Pavlidou, V., Korkidis, G., Tomaras, T. N., & Tanoglidis, D. 2020, A&A, 638, L8 [CrossRef] [EDP Sciences] [Google Scholar]
  87. Peery, S., Watkins, R., & Feldman, H. A. 2018, MNRAS, 481, 1368 [Google Scholar]
  88. Piffaretti, R., Arnaud, M., Pratt, G. W., Pointecouteau, E., & Melin, J.-B. 2011, A&A, 534, A109 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  89. Planck Collaboration X. 2011, A&A, 536, A10 [Google Scholar]
  90. Planck Collaboration XI. 2011, A&A, 536, A11 [Google Scholar]
  91. Planck Collaboration XXIX. 2014, A&A, 571, A29 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  92. Planck Collaboration XXVII. 2016, A&A, 594, A27 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  93. Planck Collaboration I. 2020, A&A, 641, A1 [CrossRef] [EDP Sciences] [Google Scholar]
  94. Planck Collaboration III. 2020, A&A, 641, A3 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  95. Planck Collaboration Int. XIII. 2014, A&A, 561, A97 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  96. Pratt, C. T., & Bregman, J. N. 2020, ApJ, 890, 156 [Google Scholar]
  97. Pratt, G. W., Croston, J. H., Arnaud, M., & Böhringer, H. 2009, A&A, 498, 361 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  98. Qin, F., Howlett, C., Staveley-Smith, L., & Hong, T. 2019, MNRAS, 482, 1920 [Google Scholar]
  99. Rameez, M., Mohayaee, R., Sarkar, S., & Colin, J. 2018, MNRAS, 477, 1772 [NASA ADS] [CrossRef] [Google Scholar]
  100. Reichert, A., Böhringer, H., Fassbender, R., & Mühlegger, M. 2011, A&A, 535, A4 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  101. Reiprich, T. H. 2017, Astron. Nachr., 338, 349 [NASA ADS] [CrossRef] [Google Scholar]
  102. Reiprich, T. H., & Böhringer, H. 2002, ApJ, 567, 716 [NASA ADS] [CrossRef] [Google Scholar]
  103. Rossetti, M., Gastaldello, F., Ferioli, G., et al. 2016, MNRAS, 457, 4515 [Google Scholar]
  104. Rubart, M., Bacon, D., & Schwarz, D. J. 2014, A&A, 565, A111 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  105. Salehi, A., Farajollahi, H., Motahari, M., et al. 2020, Eur. Phys. J. C, 80, 753 [EDP Sciences] [Google Scholar]
  106. Salehi, A., Yarahmadi, M., & Fathi, S. 2021, MNRAS, 504, 1304 [Google Scholar]
  107. Schellenberger, G., Reiprich, T. H., Lovisari, L., Nevalainen, J., & David, L. 2015, A&A, 575, A30 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  108. Schlafly, E. F., & Finkbeiner, D. P. 2011, ApJ, 737, 103 [NASA ADS] [CrossRef] [Google Scholar]
  109. Schlegel, D. J., Finkbeiner, D. P., & Davis, M. 1998, ApJ, 500, 525 [NASA ADS] [CrossRef] [Google Scholar]
  110. Scrimgeour, M. I., Davis, T. M., Blake, C., et al. 2016, MNRAS, 455, 386 [CrossRef] [Google Scholar]
  111. Secrest, N. J., von Hausegger, S., Rameez, M., et al. 2021, ApJ, 908, L51 [Google Scholar]
  112. Shamir, L. 2020, PASA, 37, e053 [Google Scholar]
  113. Shanks, T., Metcalfe, N., Chehade, B., et al. 2015, MNRAS, 451, 4238 [Google Scholar]
  114. Shanks, T., Hogarth, L. M., & Metcalfe, N. 2019a, MNRAS, 484, L64 [NASA ADS] [CrossRef] [Google Scholar]
  115. Shanks, T., Hogarth, L. M., Metcalfe, N., & Whitbourn, J. 2019b, MNRAS, 490, 4715 [NASA ADS] [CrossRef] [Google Scholar]
  116. Siewert, T. M., Schmidt-Rubart, M., & Schwarz, D. J. 2020, ArXiv e-prints [arXiv:2010.08366] [Google Scholar]
  117. Skrutskie, M., Cutri, R., Stiening, R., et al. 2006, AJ, 131, 1163 [NASA ADS] [CrossRef] [Google Scholar]
  118. Soltis, J., Farahi, A., Huterer, D., & Liberato, C. M. 2019, Phys. Rev. Lett., 122, 091301 [Google Scholar]
  119. Spallicci, A. D. A. M., Helayël-Neto, J. A., López-Corredoira, M., & Capozziello, S. 2021, Eur. Phys. J. C, 81, 4 [EDP Sciences] [Google Scholar]
  120. Springob, C. M., Magoulas, C., Colless, M., et al. 2014, MNRAS, 445, 2677 [Google Scholar]
  121. Sun, Z. Q., & Wang, F. Y. 2019, Eur. Phys. J. C, 79, 783 [CrossRef] [EDP Sciences] [Google Scholar]
  122. Tiwari, P., Kothari, R., Naskar, A., Nadkarni-Ghosh, S., & Jain, P. 2015, Astropart. Phys., 61, 1 [NASA ADS] [CrossRef] [Google Scholar]
  123. Tsaprazi, E., & Tsagas, C. G. 2020, Eur. Phys. J. C, 80, 757 [EDP Sciences] [Google Scholar]
  124. Tully, R. B., Pomarède, D., Graziani, R., et al. 2019, ApJ, 880, 24 [NASA ADS] [CrossRef] [Google Scholar]
  125. Verde, L., Kamionkowski, M., Mohr, J. J., & Benson, A. J. 2001, MNRAS, 321, L7 [NASA ADS] [CrossRef] [Google Scholar]
  126. Vikhlinin, A., Burenin, R. A., Ebeling, H., et al. 2009, ApJ, 692, 1033 [NASA ADS] [CrossRef] [Google Scholar]
  127. Watkins, R., & Feldman, H. A. 2015, MNRAS, 447, 132 [NASA ADS] [CrossRef] [Google Scholar]
  128. Watkins, R., Feldman, H. A., & Hudson, M. J. 2009, MNRAS, 392, 743 [NASA ADS] [CrossRef] [Google Scholar]
  129. Whitbourn, J. R., & Shanks, T. 2014, MNRAS, 437, 2146 [NASA ADS] [CrossRef] [Google Scholar]
  130. 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]
  131. Wright, E. L. 1979, ApJ, 232, 348 [NASA ADS] [CrossRef] [Google Scholar]
  132. Wright, E. L., Eisenhardt, P. R., Mainzer, A. K., et al. 2010, AJ, 140, 1868 [NASA ADS] [CrossRef] [Google Scholar]
  133. York, D. G., Adelman, J., Anderson, J. E., Jr, et al. 2000, AJ, 120, 1579 [Google Scholar]
  134. Zhang, Y.-Y., Reiprich, T. H., Schneider, P., et al. 2017, A&A, 599, A138 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  135. Zhao, D., Zhou, Y., & Chang, Z. 2019, MNRAS, 486, 5679 [CrossRef] [Google Scholar]
  136. Zitrin, A., Bartelmann, M., Umetsu, K., Oguri, M., & Broadhurst, T. 2012, MNRAS, 426, 2944 [NASA ADS] [CrossRef] [Google Scholar]

Appendix A: Measurements of YSZ, LBCG, R, and NH,Xray

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

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

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

(A.1)

where y(θ) is the y-parameter profile of the cluster that has been normalized to unity, θ is the radial angular coordinate in the plane of the sky, and θ500 is the apparent R500 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 X-ray 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 R500 to LX (), a fluctuation of LX 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 R500 would not cause a major change in Y5R500, certainly less than the statistical uncertainties. Thus, it is reasonable to ignore any uncertainties of θ500 for measuring Y5R500.

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

The four Y5R500 values do not differ significantly to each other ( < 7.2%), and any selection for the default Y5R500 used in this work results in the same conclusions with only minimal numerical fluctuations. We choose as the default Y5R500 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 Y5R500 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 YSZ comes from DA and the (biased) redshift–distance conversion.

We compare the values of our default Y5R500 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 Y5R500 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 Y5R500 as a function of MMF filter templates of different cluster size, parameterized through the scale radius θs = θ500/c500, where c500 is the concentration parameter of the cluster. However, for the vast majority of clusters, the θs − Y5R500 plane is left largely unconstrained, for which reason the Planck Collaboration employs an XMM-Newton derived θ − Y5R500 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 well14. 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 Y5R500 can therefore still appear. The comparison does not significantly change with the use of Y5R500 values based on the other three E19 approaches we described in Appendix A.1.

thumbnail Fig. A.1.

Comparison between the default Y5R500 values in this work based on the method of E19 and the Y5R500 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 near-infrared BCG luminosity LBCG

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 R500 from the X-ray center of the cluster. Both of these values are taken by the MCXC catalog. Thirdly, the redshift of the BCG must be within a certain range Δz of the cluster redshift. If the BCG redshift was estimated via photometric measurements, then Δz ≤ 0.02, whilst the allowed redshift difference for spectroscopic redshifts was Δz ≤ 0.01. If multiple BCG candidates were suggested for one cluster, the following criteria were considered. The BCG should possess a rather extended envelope and it should be selected by optical surveys over NIRs.

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

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

To determine NH,Xray, we used the exact same spectral fitting procedure as described in M20 when the X-ray redshift was fitted. This time the NH,Xray is left free to vary instead. Its value was linked between the spectra from the 0 to 0.2 R500 region and the 0.2 to 0.5 R500 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 NH,Xray for 213 out of 237 Chandra clusters and 27 out of 76 XMM-Newton clusters. When we performed the NH,Xray − NHtot comparison, we saw that the data residuals compared to the best-fit model behaved strongly as a function of the NH,Xray measurement uncertainty. Clusters with large statistical uncertainties were significantly downscattered, increasing the scatter of the relation. These were mostly low NH,Xray clusters. The 0.7 keV cut removes most of the power for detecting NHtot (which shows up more strongly at lower energies). So when the NHtot 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 NH,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 Lbol − YSZ, Lbol − T, and YSZ − 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 LX − YSZ relation

We study the Lbol − YSZ 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 Lbol − YSZ relation. If however there was some extra X-ray absorption taking place toward that direction, it would correspond to an extra NHtot ∼ 7.9 ± 3.4 × 1020 cm−2. This value is consistent with the LX − YSZ result of our sample within 1.1σ. The maximum anisotropy region however is located 90° away from the LX − YSZ direction. This could indicate that the normalization variation is due to statistical noise and does not reveal any unknown absorption. Alternatively, it could be attributed to the incapability of the θ = 90° cones to accurately pinpoint the direction where the hidden X-ray absorption takes place, especially since the direction uncertainties are considerably large. Based on all the above, it cannot be excluded that the mild anisotropies seen in the ACC sample are due to chance.

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

B.2. The Lbol − T relation

For the Lbol − 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 sky15. 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 Lbol − 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 H0 = 61.3 ± 4.2 km s−1 Mpc−1 toward (l, b)∼(318°, −9°) and H0 = 72.2 ± 2.4 km s−1 Mpc−1 for the rest of the sky. The obtained H0 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 uBF = 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 uBF = 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.

thumbnail Fig. B.1.

Normalization anisotropy maps and the respective statistical significance maps of the anisotropies for the Lbol − YSZ (top), the Lbol − T (middle), and the YSZ − T (bottom) scaling relations for the ACC sample. The Lbol − YSZ traces only unaccounted for X-ray absorption effects, while the LX − T and YSZ − T behavior mirrors cosmological anisotropies and BFs.

B.3. The YSZ − T relation

We repeat the YSZ − T anisotropy analysis with only the 113 ACC clusters with YSZS/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 AYT of this region compared to the rest of the sky is 35 ± 11%. The direction is also identical to that obtained with the Lbol − T relation with a larger statistical significance.

In terms of H0, we obtain H0 = 60.6 ± 3.6 km s−1 Mpc−1 toward the most anisotropic region and H0 = 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 LX − T analysis results.

For the BF scenario and the MR method, we obtain uBF = 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 uBF = 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 − LX, R − T, R − YSZ, R − LBCG, and YSZ − LBCG 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 −  LX, R − T, R − YSZ, and R −  LBCG 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 − LX, R − YSZ, and R − LBCG 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 H0 would result in a 7−10% variation of ARX. 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 ART ∝ DA dependency. For the same underlying effect, R − T shows half the normalization variation than LX − T, YSZ − T, and LBCG − T. Despite the weaker signal, this test can still be a valuable to detect cosmological anisotropies and BFs when applied to larger cluster samples, and when a more precise modeling of the systematics is feasible.

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

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

thumbnail 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 − LX and R − YSZ are shown in Fig. C.2. The map for R − LBCG 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 ARX. 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.

thumbnail Fig. C.2.

Normalization anisotropy map for R − LX (left) and R − YSZ (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 105 randomly drawn bootstrap subsamples (same process as in M20), independent of direction. We investigate the correlation between the ART and the median XBO for every subsample. The results are shown in Fig. C.3. We then calibrate the ART 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 ART behavior, which is consistent with the previous results and the scenario of a cosmological anisotropy or a large BF.

thumbnail Fig. C.3.

Top: correlation between the best-fit normalization ART (over the full best-fit value of the full sample) and the median XBO for every of the 105 bootstrap subsamples. Bottom: ART anisotropy map after calibrating for the existing correlation with the XBO.

C.3. The YSZ −  LBCG relation

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

As a result of the large scatter and the limited number of clusters, we consider θ = 90° cones to scan the sky. The anisotropy of the YSZ − LBCG 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.

thumbnail Fig. C.4.

Normalization anisotropy map for the YSZ − LBCG 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 LX − YSZ, YSZ − T, and LBCG − T relations, for different YSZ S/N and redshift cuts. The summarized result is that our main conclusions remain unchanged.

D.1. The LX − YSZ relation

For the LX − YSZ relation, a lower YSZ 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σLX − YSZ contours for the three different cases are compared to Fig. D.1. There is a > 5σ shift in the best-fit LX − YSZ 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 LX residuals are strongly correlated with z, NHtot, and YSZ, as shown in Fig. D.2. These probably point toward biases in the X-ray selection process and the need for a broken power law to describe LX − YSZ. The same figure show that these systematic behaviors are not present for the S/N > 4.5 case.

thumbnail Fig. D.1.

3σ (99.7%) parameter space of the normalization and slope of the LX − YSZ (top), YSZ − T (middle), and YSZ − LBCG (bottom) relations for S/N > 2 (purple), S/N > 3 (green), and S/N > 4.5 (black).

thumbnail Fig. D.2.

LX residuals for the LX − YSZ relation for S/N > 2 (top) and S/N > 4.5 (bottom) as a function of z (left), NHtot (middle), and YSZ (right). The green strips indicate the best-fit 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 LX − YSZ anisotropies, we get a similar map as for the default case. The anisotropy significance map (Fig. D.3) illustrates that there are no excess X-ray absorption issues toward (l, b)∼(280°, −15°) where the cosmological anisotropies are found. At the same time, the only anisotropic region again appears to be toward (l, b)∼(128° , + 21°), being ∼14% fainter than the rest of the sky at a 2.5σ level, similar to the default case.

thumbnail Fig. D.3.

Statistical significance map of the LX − YSZ anisotropies for S/N > 2.

D.2. The YSZ − T relation

For the YSZ − T relation, we used a lower YSZ threshold of S/N ≥ 2 in the default analysis, resulting in 260 clusters. If we increase this lower threshold to S/N ≥ 3 or S/N ≥ 4.5, we are left with 242 and 190 clusters, respectively. The best-fit YSZ − T does not change as a function of S/N, as shown in the middle panel of Fig. D.1. The YSZ 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 AYT variance map is shown in Fig. D.4. Thus, the observed anisotropies are independent of YSZ selection effects.

thumbnail Fig. D.4.

Normalization anisotropy map of the YSZ − T relation for S/N > 3.

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

To ensure that the observed YSZ − T anisotropies do not emerge because of some unknown directional bias in our own YSZ measurements, we repeat the analysis using the YSZ values from the PSZ2 catalog. Owing to the higher YSZS/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 YSZ compared to our default analysis. Thus, we adopt CY = 85 kpc2 and CX = 6 keV. For the best-fit parameters, we find AYT = 0.940 ± 0.041, BYT = 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 YSZS/N > 4.5) as a result of the different adopted R500 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 YSZ − T scaling relation is plotted in the top panel of Fig. D.5.

thumbnail Fig. D.5.

Top: YSZ − T relation when using the PSZ2 YSZ values, together with its 1σ best-fit function (green). Bottom: normalization anisotropy map of the YSZ − T relation when the PSZ2 values for YSZ are used.

Performing the YSZ − T sky scanning, we obtain roughly the same anisotropies with the PSZ2 YSZ values as we did with our measurements. The AYT map is shown in the bottom panel of Fig. D.5. The maximum anisotropy is found toward (l, b)=(254°, −22°), where H0 = 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 YSZ − T anisotropies are independent of the adopted YSZ catalog.

D.4. The LBCG − T relation

For the LBCG 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 LBCG remains unknown. We applied these cuts to all LBCG scaling relations to be consistent. However, these redshift issues are not as important in the LBCG − T relation. In Fig. D.6, the LBCG residuals are plotted as a function of z. Only the seven clusters with z < 0.02 are systematically upscattered, while z > 0.15 show the same behavior as the less distant clusters. Thus, we repeat our analysis using all the 259 clusters, regardless of their redshift. The best-fit LBCG − 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 LBCG − T anisotropy signal rises from 1.9σ to 2.1σ. Both of these results are shown in Fig. D.6.

thumbnail Fig. D.6.

Top: LBCG residuals of the LBCG − T fit, as a function of the BCG redshift. The green stripe corresponds to the best-fit function within 1σ. Middle: 3σ (99.7%) parameter space of the normalization and slope of the LBCG − T relation, for all clusters (green), and for 0.03 < z < 0.15 clusters (purple). Bottom: normalization anisotropy map of the LBCG − 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 LX − T relation since results from the other relations are similar. This is because the samples used are simulated and isotropic and that the cluster positions are the same across scaling relations. Thus, any effects come only from the applied cluster weighting during the sky scanning. We stress again that the observed anisotropies remain unaffected when this weighting is omitted and the ZoA gap is invisible to the used algorithm, highlighting the lack of bias coming from ZoA in the real data.

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

thumbnail 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 LX − T relation. The distance weighting during the LX − 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.

thumbnail Fig. E.2.

Number of clusters per θ = 75° cone for the 313 clusters used in the LX − 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 × 107 iterations of the chain, with a burn in period of 104. 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 104 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 uz and uT (redshift and temperature power indexes), which are the only two parameters with a significant impact on the best-fit normalization. As discussed before, their anticorrelation is strong, which generally cancels out any strong effects in the normalization owing to the flux-limited sample and the Eddington bias. The 3σ parameter space for uNH and uf (NHtot 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 (uNH ∼ uf ∼ 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.

thumbnail 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 LX − T relation. The distance weighting during the LX − 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 XMM-Newton clusters in our sample

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

thumbnail Fig. E.4.

Spatial variation of fraction of Chandra clusters minus XMM-Newton clusters over the sum.

E.4. Effects of bulk flows to adopted apparent R500

If a BF is present, then the distance of a cluster is miscalculated. As a result, its LX, and its mass M500, its physical R500 (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, YSZ, and LX 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 H0 anisotropies, θ500 remains practically unchanged. That is because the opposite changes of R500 and DA 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 DL = (1 + z)2DA. Then, θ500 is written as

(E.1)

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 low-z clusters, we expect a ≲6% bias. As an example, let us estimate how this would affect the anisotropies of the LX − 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 YSZ 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 YSZ, which again suppress the BF signal. However, this effect is insignificant for two reasons. Firstly, because of the slope of the YSZ − T relation (BYT = 2.546), the changes of YSZ are less important than changes in T (which we already saw that can mildly smooth out the anisotropies). Secondly, the YSZ − T anisotropies are much larger than the YSZ changes owing to a falsely assumed θ500.

Moreover, the measured T value goes into the measurement of YSZ 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 YSZ, and thus a BF would practically not affect YSZ in that regard.

For the LX 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 LX to remain almost unchanged, since only the outskirts of the assumed cluster size slightly changes in case of a BF. The vast majority of X-ray cluster emission however comes from well within this area. Even if the change was not negligible, the same effect as for T and YSZ is expected, where the miscalculation of LX 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 YSZ and T measurements and BF estimations would be needed, and these changes would further enhance the observed anisotropies. Thus we ignore these changes.

thumbnail Fig. E.5.

Relative change of the measured Y5R500 when the input R500 (and θ500) is increased (red) or decreased (blue) by 6%.

All Tables

Table 1.

Possible anisotropy causes that can be traced by each cluster scaling relation.

Table 2.

Best-fit parameters of the 10 scaling relations.

Table 3.

Maximum anisotropy direction for every scaling relation, together with the needed H0 relative variation to fully explain the anisotropy.

Table 4.

Best-fit BFs for every scaling relation, method, and redshift bin considered in this work.

All Figures

thumbnail Fig. 1.

Best fits of the 10 scaling relations studied in this work. The green area represents the 1σ limits of the best-fit area. From top left panel to the right: LX − YSZ, LX − LBCG, LX − T, YSZ − T, LBCG − T, R − LX, R − YSZ, R − LBCG, R − T, YSZ − LBCG, R − YSZ, Lbol − YSZ (ACC), and YSZ − T (ACC) relations.

In the text
thumbnail Fig. 2.

Normalization anisotropy maps (left) and the respective statistical significance maps of the anisotropies (right), for LX − YSZ (top), joint LX − YSZ (our sample + ACC, middle), and LBCG − T (bottom). All the maps in this work are shown in a Hammer projection.

In the text
thumbnail Fig. 3.

Comparison between the NHtot from W13 and the NHtot constrained by our X-ray spectral analysis. The equality line is indicated in purple; the 1σ space of the best-fit 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
thumbnail Fig. 4.

Top: same as in Fig. 2, for LX − T. Bottom: H0 anisotropy map derived from the joint LX − T (our sample+ACC).

In the text
thumbnail Fig. 5.

Hubble diagram of galaxy clusters as derived by the LX − T (top) and the YSZ − 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 best-fit lines are indicated with the same color.

In the text
thumbnail Fig. 6.

Top: same as in Fig. 2 for YSZ − T. Bottom: H0 anisotropy map (left) and the respective significance map (right) derived from the joint YSZ − T (our sample+ACC, bottom).

In the text
thumbnail Fig. 7.

H0 anisotropy map as derived from the joint analysis of LX − T, YSZ − T, and LBCG − T relations for both samples.

In the text
thumbnail Fig. 8.

Histograms of the statistical significance of the maximum anisotropy as detected in 10 000 isotropic MC simulations for the LX − T, YSZ − T, LBCG − T, Lbol − T (ACC), and YSZ − T (ACC) scaling relations. The vertical black line represents the results from the real data; the p-value and the direction are also shown.

In the text
thumbnail Fig. 9.

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

In the text
thumbnail Fig. 10.

3σ (99.7%) parameter space of the normalization and slope of the YSZ − T (left), LBCG − T (middle), and LX − T (right) relations for relaxed (purple) and disturbed (green) clusters.

In the text
thumbnail Fig. 11.

Correlation between the best-fit ALT (over the full sample’s best-fit value) and the median XBO for every of the 105 bootstrap subsamples.

In the text
thumbnail Fig. 12.

Correlation between the YSZ scatter from the YSZ − T relation and the LX scatter from the LX − T relation.

In the text
thumbnail 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 LX − T, YSZ − T, and LBCG − T relations, respectively. The blue diamonds and triangles correspond to the ACC Lbol − T and YSZ − 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
thumbnail Fig. A.1.

Comparison between the default Y5R500 values in this work based on the method of E19 and the Y5R500 values as given in the PSZ2 catalog for the 566 clusters in common. The equality line is shown in green.

In the text
thumbnail Fig. B.1.

Normalization anisotropy maps and the respective statistical significance maps of the anisotropies for the Lbol − YSZ (top), the Lbol − T (middle), and the YSZ − T (bottom) scaling relations for the ACC sample. The Lbol − YSZ traces only unaccounted for X-ray absorption effects, while the LX − T and YSZ − T behavior mirrors cosmological anisotropies and BFs.

In the text
thumbnail 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
thumbnail Fig. C.2.

Normalization anisotropy map for R − LX (left) and R − YSZ (right).

In the text
thumbnail Fig. C.3.

Top: correlation between the best-fit normalization ART (over the full best-fit value of the full sample) and the median XBO for every of the 105 bootstrap subsamples. Bottom: ART anisotropy map after calibrating for the existing correlation with the XBO.

In the text
thumbnail Fig. C.4.

Normalization anisotropy map for the YSZ − LBCG relation. The statistical significance of the observed anisotropies is ≤1.4σ, and thus the relation is statistically isotropic.

In the text
thumbnail Fig. D.1.

3σ (99.7%) parameter space of the normalization and slope of the LX − YSZ (top), YSZ − T (middle), and YSZ − LBCG (bottom) relations for S/N > 2 (purple), S/N > 3 (green), and S/N > 4.5 (black).

In the text
thumbnail Fig. D.2.

LX residuals for the LX − YSZ relation for S/N > 2 (top) and S/N > 4.5 (bottom) as a function of z (left), NHtot (middle), and YSZ (right). The green strips indicate the best-fit 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
thumbnail Fig. D.3.

Statistical significance map of the LX − YSZ anisotropies for S/N > 2.

In the text
thumbnail Fig. D.4.

Normalization anisotropy map of the YSZ − T relation for S/N > 3.

In the text
thumbnail Fig. D.5.

Top: YSZ − T relation when using the PSZ2 YSZ values, together with its 1σ best-fit function (green). Bottom: normalization anisotropy map of the YSZ − T relation when the PSZ2 values for YSZ are used.

In the text
thumbnail Fig. D.6.

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

In the text
thumbnail 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 LX − T relation. The distance weighting during the LX − 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
thumbnail Fig. E.2.

Number of clusters per θ = 75° cone for the 313 clusters used in the LX − T relation.

In the text
thumbnail 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 LX − T relation. The distance weighting during the LX − 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
thumbnail Fig. E.4.

Spatial variation of fraction of Chandra clusters minus XMM-Newton clusters over the sum.

In the text
thumbnail Fig. E.5.

Relative change of the measured Y5R500 when the input R500 (and θ500) is increased (red) or decreased (blue) by 6%.

In the text

Current usage metrics show cumulative count of Article Views (full-text 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 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.