© ESO 2020

1. Introduction

The isotropy of the Universe on sufficiently large scales is a fundamental pillar of the standard model of cosmology. The most important consequence of isotropy is that the expansion rate of the Universe as well as the physical properties of all astronomical objects must be the same regardless of the direction in the sky. Due to the high significance of this hypothesis, it is necessary that it is robustly scrutinized and tested against different cosmological probes using the latest data samples.

What was initially introduced as a repercussion of general relativity and the Friedmann-Lemaître-Robertson-Walker (FLRW) metric was later supported by observations. The most crucial of them is arguably the cosmic microwave background (CMB) as observed by Cosmic Background Explorer (COBE, Efstathiou et al. 1992), the Wilkinson Microwave Anisotropy Probe (WMAP, Bennett et al. 2013), and the Planck (Planck Collaboration XVI 2016) telescopes. CMB shows a remarkable isotropy at small angular scales (high multipoles), whilst some anisotropies are still present in lower multipoles. The most prominent one is the so-called CMB dipole which, if one assumes its purely kinematic origin, is caused by the Doppler shift due to the motion of our Solar System with respect to the CMB rest frame. This indicates that the Solar System moves toward the Galactic coordinates (l, b)  ∼  (264°, + 48°) with a peculiar velocity of ∼370 km s⁻¹ (Kogut et al. 1993; Fixsen et al. 1996) with respect to the CMB rest frame. Another anisotropic feature present in the CMB is the dipole power asymmetry detected in both WMAP and Planck with a significance of   ∼  2−3.5σ toward (l, b) ∼ (230°, − 20°) (Eriksen et al. 2004; Hanson & Lewis 2009; Bennett et al. 2011; Akrami et al. 2014; Planck Collaboration XXIII 2014; Planck Collaboration XVI 2016). Its nature still remains relatively unclear. The interpretation of the significance of these results differ between papers. Other potential challenges for the isotropy of the Universe found in the CMB is the parity asymmetry, ∼3σ toward (l, b) = (264°, − 17°) and the unexpected quadrupole-octopole alignment and the existence of the Cold Spot at (l, b) = (210°, −57°) (Tegmark et al. 2003; Vielva et al. 2004; Kim & Naselsky 2010; Aluri & Jain 2012; Cai et al. 2014; Planck Collaboration XXIII 2014; Planck Collaboration XVI 2016; Schwarz et al. 2016).

The CMB is a great tool to study the behavior of the early Universe. However, it can be quite challenging to extract information about the directional behavior of the late Universe based on that probe. This problem intensifies when one considers that according to ΛCDM the late Universe is dominated by dark energy whose effects are not directly present in the CMB spectrum. Moreover, since the nature of dark energy is still completely unknown, one can only make assumptions about the isotropic (or not) behavior of dark energy. Consequently, it becomes clear that other cosmological probes, at much lower redshifts than the CMB, are needed in order to search for possible anisotropies in the late Universe. There are indeed many probes that have been used for such tests.

For instance, Type Ia supernovae (SNIa) have been extensively used to test the isotropy of the Hubble expansion with many results reporting no significant deviation from the null hypothesis (Lin et al. 2016; Andrade et al. 2018; Sun & Wang 2018; Wang & Wang 2018). Other studies however, claim to detect mild-significance anisotropies in the SNIa samples that more often than not approximately match the direction of the CMB dipole (Schwarz & Weinhorst 2007; Antoniou & Perivolaropoulos 2010; Colin et al. 2011, 2017; Mariano & Perivolaropoulos 2012; Kalus et al. 2013; Appleby et al. 2015; Bengaly et al. 2015; Javanmardi et al. 2015; Migkas & Plionis 2016). Generally, the reported results from SNIa strongly depend on the used catalog. Moreover, the robustness of SNIa as probes for testing the isotropy of the Universe based on the current status of the relative surveys has been recently challenged (Colin et al. 2017; Beltrán Jiménez et al. 2015; Rameez 2019). This is mainly because of the highly inhomogeneous spatial distribution of the data (most SNIa in the latest catalogs lie close to the CMB dipole direction), their sensitivity to the applied kinematic flow models which readjust their measured heliocentric redshifts, as well as the assumptions that go into the calibration of their light curves.

Other probes that have been used to pinpoint possible anisotropies or inconsistencies with the ΛCDM model are the X-ray background (Shafer & Fabian 1983; Plionis & Georgantopoulos 1999), the distribution of optical (Javanmardi & Kroupa 2017; Sarkar et al. 2019) and infrared galaxies (Yoon et al. 2014; Rameez et al. 2018), the distribution of distant radio sources (Condon 1988; Blake & Wall 2002; Singal 2011; Rubart & Schwarz 2013; Tiwari & Nusser 2016; Bengaly et al. 2018, 2019; Colin et al. 2019), gamma-ray bursts (Řípa & Shafieloo 2017; Andrade et al. 2019), peculiar velocities of galaxy clusters (Kashlinsky et al. 2008, 2010, 2011; Watkins et al. 2009; Atrio-Barandela et al. 2015) and of SNIa (Appleby et al. 2015) etc. While some of them find no statistically significant challenges for the null hypothesis of isotropy, others provide results which are unlikely to occur within the standard cosmological model framework. The use of standard sirens for such tests in the future has also been proposed (e.g., Cai et al. 2018).

Since the outcome of the search for a preferred cosmological direction remains ambiguous, new and independent methods for such tests should be introduced and applied to the latest data samples. In Migkas & Reiprich (2018; hereafter M18), the use of the directional behavior of the galaxy cluster X-ray luminosity-temperature relation is described as a cosmological probe. It is well-known that galaxy clusters are the most massive gravitationally bound systems in the universe, strongly emitting X-ray photons due to the large amounts of hot gas they contain (∼10% of their total mass) in their intra-cluster medium (ICM). Their physical quantities follow tight scaling relations, for which Kaiser (1986) provided mathematical expressions. Specifically, the correlation between the X-ray luminosity (L_X) and the ICM gas temperature (T) of galaxy clusters is of particular interest since it can be used to trace the isotropy of the Universe, which is a new concept for such cosmological studies. The general properties of the L_X–T scaling relation have been extensively scrutinized in the past by several authors (e.g., Vikhlinin et al. 2002; Pacaud et al. 2007; Pratt et al. 2009; Mittal et al. 2011; Reichert et al. 2011; Hilton et al. 2012; Maughan et al. 2012; Bharadwaj et al. 2015; Lovisari et al. 2015; Giles et al. 2016; Zou et al. 2016; Migkas & Reiprich 2018; Ebrahimpour et al. 2018).

In a nutshell, the gas temperature, the flux and the redshift of a galaxy cluster do not require any cosmological assumptions in order to be measured¹. Using the flux and the redshift together with the luminosity distance, through which the cosmological parameters come into play, one can obtain the luminosity of a cluster. The luminosity however can also be predicted (within an uncertainty range) based on the cluster gas temperature. Hence, adjusting the cosmological parameters, one can make the two luminosity estimations match. This can be repeatedly applied to different sky patches in order to test the consistency of the obtained values as a function of the direction. The full detailed physical motivation behind this is discussed in M18. There, it is shown that the directional behavior of the normalization of the L_X–T relation strictly follows the directional behavior of the cosmological parameter values. This newly introduced method to test the Cosmological Principle (CP) could potentially prove very effective due to the very homogeneous sky coverage of many galaxy cluster samples (in contrast to SNIa samples), the plethora of available data as well as large upcoming surveys such as eROSITA (Predehl et al. 2016) which will allows us to measure thousands of cluster temperatures homogeneously (Borm et al. 2014). For studying the isotropy of the Universe with the future surveys, it is of crucial importance that any existing systematic biases that could potentially affect the L_X and T measurements of galaxy clusters would have been identified and taken into account by then.

In this paper, we construct and use a new galaxy cluster sample in order to identify regions that share a significantly different L_X–T relation compared to others. This could lead to pinpointing an anisotropy in the Hubble expansion or discover previously unknown factors which could potentially affect X-ray measurements of any kind. Except for the high quality observations and measurements, another advantage of our sample is the small overlap with the XCS-DR1 (Mehrtens et al. 2012) and ACC (Horner 2001) samples used in M18. There are only three common clusters between our sample and XCS-DR1 (<1%) and only a 30% overlap with ACC. Since we reanalyze the XCS-DR1 and ACC sample in this paper as well, using the same methods we use for our sample in order to perform a consistent comparison, all the common clusters between the different catalogs are excluded from XCS-DR1 and ACC. Our cluster sample does not suffer from any strong archival biases in contrast to ACC. Throughout this paper we use a ΛCDM cosmology with H₀ = 70 km s⁻¹ Mpc⁻¹, Ω_m = 0.3 and Ω_Λ = 0.7 unless stated otherwise.

The paper is organized as follows: In Sect. 2 we describe the construction of the sample and how we derive the properties of the clusters. In Sect. 3 we explain the steps we follow for the data reduction and the spectral analysis of the observations. In Sect. 4 we present the modeling of the L_X–T relation together with the parameter fitting procedure. We also explain how we identify possible spatial anisotropies and assign their statistical significance. In Sect. 5 we present the first results, including the overall L_X–T results and the 1-dimensional (1D) and 2-dimensional (2D) anisotropies of our sample. In Sect. 6 we investigate several X-ray and cluster-related causes that could possibly produce the observed anisotropic signal in our sample. In Sect. 7 we examine the case where the anisotropies in our sample have a cosmological origin instead, assuming there are no systematics associated with X-ray photons or cluster properties. In Sect. 8 we combine our results with those obtained from ACC and XCS-DR1. We express these joint-analysis anisotropies in cosmological terms. In Sects. 9 and 10 we discuss our findings and their implications, compare with other studies and summarize.

2. Sample selection

The cluster sample used in this work is a homogeneously selected one based on the Meta-Catalogue of X-ray detected Clusters of galaxies (MCXC, Piffaretti et al. 2011). Initially it consisted of the 387 galaxy clusters above an unabsorbed flux cut of f_{0.1−2.4 keV} ≥ 5 × 10⁻¹² ergs s⁻¹ cm⁻² excluding the Galactic plane (|b| ≤ 20°), the Magellanic clouds and the Virgo cluster area. The flux for every MCXC cluster is found based on the given X-ray luminosity and redshift, combined with a reversed K-correction (the necessary temperature input was found by the L_X–T relation of Reichert et al. 2011). The MCXC luminosities (and therefore the calculated fluxes on which our sample is based) were corrected for absorption based on HI measurements (see Sects. 2.2 and 2.3).

The parent catalogs of the clusters we use are The ROSAT extended Brightest Cluster Sample (eBCS, Ebeling et al. 2000), The Northern ROSAT All-Sky (NORAS) Galaxy Cluster Survey (Böhringer et al. 2000) and the ROSAT-ESO Flux-Limited X-Ray (REFLEX) Galaxy Cluster Survey Catalog (Böhringer et al. 2004). They are all based on the ROSAT All-Sky Survey (RASS, Voges et al. 1999).

Another selection criterion was for the clusters to have good quality Chandra (Weisskopf et al. 2000) or XMM-Newton (Jansen et al. 2001) public observations (as of July 2019). This criterion is satisfied for 331 clusters. The rest 56 clusters for which such observations were not available have a sparse sky distribution and similar average properties with the 331 clusters (as discussed later in the paper) and therefore their inclusion is not expected to alter the results.

We reduced these 331 available observations, analyzed them and extracted cluster properties as described below. This sample includes a large fraction (∼85%) of the eeHIFLUGCS (extremely expanded HIghest X-ray FLUx Galaxy Cluster Sample Reiprich 2017; Pacaud et al., in prep.) sample clusters. eeHIFLUGCS is a complete, purely X-ray flux-limited sample with similar selection criteria.

The final sample we use for this paper consists of 313 galaxy clusters. The other 18 clusters are not used because of the following reasons. Firstly, we excluded 11 clusters from our analysis that we identified as apparent multiple systems (out of a total of 15). This is due to the fact that clusters in rich environments tend to be systematically fainter than single clusters (M18 and references therein), thus biasing the final results. Moreover, some of these seemingly multiple systems are located at different redshifts but projected in our line of sight as real double and triple systems (see Ramos-Ceja et al. 2019, hereafter R19). When these different components are accounted as one system in MCXC the flux of the “single cluster” is overestimated. As a result, some of these systems falsely overcome the selected flux limit while none of their true individual components has the necessary flux to be included in our catalog. On the other hand, there are cases where one of the individual extended components has enough flux to be kept in our sample and at the same time it is located at a different redshift than the other components of the system (it does not belong to the same rich environment). In this case, these extended sources were kept in the catalog and their L_X values were adjusted correspondingly (see Sect. 2.3) while the other component(s) were excluded. Since we use Chandra and XMM-Newton observations while the MCXC L_X–T values come from ROSAT observations, some minor inconsistencies between these values and our own measurements are expected. In order to account for this, we consider it as an extra source of uncertainty and adjust the confidence levels of the final cluster fluxes and luminosities accordingly as described later in the paper.

Furthermore, we identify nine clusters to be strongly contaminated by point sources, likely Active Galactic Nuclei (AGN). We confirm that by the absence of significant extended emission around the suspected point sources and by fitting a power-law (constraining the power index and the normalization) and an apec model to the spectra of the bright part of the point source using XSPEC (Arnaud 1996). If the pow model returns a better fit than the apec model we mark the source as an AGN. Moreover, we search the literature for known stars and AGNs at these positions. For three of these nine clusters (A2055, A3574E, RXCJ1840.6−7709) the point sources are located close to the (bright) cores of the clusters and cannot be deblended. Thus, we chose to exclude these clusters since their MCXC L_X values would be overestimated and would add extra bias to our analysis. For another cluster (A1735) there is a strong AGN source and a galaxy cluster with extended emission with an angular separation of ∼7′. This system has been identified as one cluster in the MCXC catalog centered at the AGN position, which has the higher contribution to the X-ray flux. Therefore, this system was also excluded since the single extended emission component does not surpass the necessary flux limit. For the rest five clusters, using the Chandra and XMM-Newton images the point sources are easily distinguishable from the cluster emission while the MCXC objects are centered close to the extended emission centers. For these five systems we calculate the flux of the point source using its spectra from one of the two aforementioned telescopes and a pow model, and subtract it from the MCXC flux in order to see if the extended emission alone overcomes the flux limit. This procedure results in the exclusion of three clusters (A0750, A0901, A2351), while the other two (A3392, S0112) stay above the desired flux limit and are considered in our analysis after appropriately decreasing their MCXC luminosity values (Sect. 2.3).

Since ∼50% of the clusters included in our sample have been observed by both Chandra and XMM-Newton, we decide to analyze these common clusters with the former. This is due to the fact that Chandra data are generally less flared than XMM-Newton data. As a result, 237 clusters are analyzed using Chandra observations while 76 clusters are processed using XMM-Newton observations. For both telescopes, we extract and fit the spectra within the energy range of 0.7−7 keV. The cross calibration of the two satellites is discussed in Sect. 2.5. Using Chandra data for ∼75% of our sample offers another advantage as well. As mentioned before, in M18 two samples are used: ACC which consists only of ASCA observations, and XCS-DR1 which consists only of XMM-Newton observations. Subsequently, mostly using a third independent telescope to built our sample and study the anisotropy of the L_X–T eliminates any systematics that might occur in the results because of telescope-specific reasons.

Consequently, the sample with which this analysis is performed consists of 313 single galaxy clusters. For these clusters we have self-consistently measured their gas temperatures T and their uncertainties, as well as their metallicities Z and their X-ray redshifts z. Furthermore, we know their optical spectroscopic z, their fluxes f_X and their luminosities L_X within the 0.1−2.4 keV energy range together with their uncertainties, their Galactic (l, b) and Equatorial (RA, Dec) coordinates and the atomic and molecular hydrogen column density N_Htot in their direction. The exact information for every parameter and where it comes from is described in the following subsections. Their spatial distribution together with the redshift value used for each cluster can be seen in Fig. 1. The vast majority of these 313 galaxy clusters are included in the eeHIFLUGCS sample.

2.1. Redshift

For 264 out of the 313 clusters the given MCXC redshifts are used. We have checked that all these clusters have at least seven galaxies with optical spectroscopic redshifts in the NASA/IPAC Extragalactic Database (NED) and that agree with the assigned MCXC redshift. The median number of galaxies per cluster is 52 for these cases. For seven other clusters we reassigned a redshift based on the already-existing optical spectroscopic data when the offset between the apparently correct redshift value and the MCXC redshift is Δz ≥ 0.007, corresponding to ∼2000 km s⁻¹.

The remaining 42 clusters either do not have enough optical spectroscopic data in order to trust the given redshift or the distribution of the galaxy redshifts of the cluster is inconclusive. In that case, the redshift of the cluster was determined from the available X-ray data. For that, we extracted and fit the spectra within the 0−0.2 R₅₀₀ ² circle and within the 0.2−0.5 R₅₀₀. Using two apec models (one for each cluster region) with the temperature and metallicity parameters free to vary for both, the redshift is also fit simultaneously but linked for the two regions (same z for both regions). In Sect. 3 the technical details of the spectral fitting process are discussed.

In order to make sure that the obtained X-ray redshifts are trustworthy, we also determined them for the 264 clusters with safe optical redshifts. The comparison of the optical and X-ray redshifts is displayed in the top panel of Fig. 2. By comparing the well-known optical redshifts with the X-ray ones, we see that there is only a very small intrinsic scatter of Δz = 0.0014 (∼420 km s⁻¹) with the agreement being remarkable. It is noteworthy that only 10% of the clusters have a deviation of more than Δz/(1 + z_opt) ≥ 0.01 between the X-ray and the good quality optical redshift, whilst only 3% deviate by more than Δz/(1 + z_opt) ≥ 0.02. Therefore, using the X-ray determined redshifts for the 42 clusters without optically spectroscopic redshifts seems to introduce no bias. The final redshift distribution is shown in the bottom panel of Fig. 2. Moreover, in the appendix is further shown that using or not the clusters with X-ray redshifts one derives very similar results. The redshift distribution of our sample covers the z = (0.003−0.45) range while the median redshift of the sample is z = 0.075 (z = 0.072 for the excluded clusters). All the aforementioned redshifts are heliocentric. The clusters for which we changed their z values compared to the ones from MCXC are displayed in Table C.1 with a star (*) next to their names.

2.2. Hydrogen column density

The value of N_H enters in both the determination of the L_X (as done in the parent catalogs) and in the T determination that is performed in this analysis. Hence, an inaccurate treatment of the input N_H values could potentially bias both parameters mostly in the opposite direction and eventually affect our results.

In the calculation of the L_X of every cluster, the REFLEX and NORAS catalogs used the neutral hydrogen N_HI values coming from Dickey & Lockman (1990; thereafter DL90), while the (e)BCS catalog used again N_HI values but as given in Stark et al. (1992). As shown in Baumgartner & Mushotzky (2006) and Schellenberger et al. (2015; S15 hereafter) the total hydrogen column density N_Htot (neutral + molecular hydrogen) starts to get significantly larger than N_HI for N_HI ≥ 6 × 10²⁰ cm⁻². If this is not taken into account it would result in a misinterpretation of the total X-ray absorption due to Galactic material and hence to underestimating L_X while generally overestimating T. In order to account for this effect, we used the N_Htot values as given by Willingale et al. (2013; hereafter W13) in all the spectra fittings we performed and for correcting the MCXC L_X values as described in the following subsections. The comparison between the N_HI used in the parent catalogs and the N_Htot values we use is displayed in Fig. 3.

One can see that for certain clusters N_Htot <  N_HI, something that seems counter-intuitive. This happens because in the calculation of N_Htot, W13 use the N_HI from the LAB survey (Kalberla et al. 2005) and not from DL90. The LAB survey has a better resolution and tends to give slightly lower N_HI values than DL90 for the same sky positions, something that can create these small inconsistencies for clusters where the molecular hydrogen is not yet high enough. As stated above, we corrected these inconsistencies for all the 313 clusters. The LAB survey covers the velocity range of (−450 km s⁻¹, +400 km s⁻¹), within which all neutral hydrogen is supposed to be detected. This velocity range naturally propagates to the N_Htot values we use and any amount of hydrogen outside of this velocity range is not accounted for. The median N_Htot for the 313 clusters is 3.81 × 10²⁰ cm⁻² (4.35 × 10²⁰ cm⁻² for the excluded clusters).

2.3. Luminosity

We chose to use ROSAT luminosity measurements for the simple reason that the entire R₅₀₀ area of the clusters is observed in the RASS. On the contrary, the field of view (FOV) of XMM-Newton and (especially) Chandra does not cover the full R₅₀₀ for most of our clusters³. It has been shown though that the ROSAT L_X values are fully consistent with the ones from XMM-Newton within the ROSAT L_X uncertainties (Böhringer et al. 2007; Zhang et al. 2011) (see Appendix A.9 for further tests)

Our estimates for the cluster X-ray luminosities used the reported X-ray luminosities in the MCXC catalog as a baseline. These luminosities were homogenized for systematics between the different parent catalogs, and were aperture-corrected to reflect the flux within R₅₀₀ (for more details see Piffaretti et al. 2011, Chap. 3.4.1). For the relative luminosity uncertainty σ_LX we assumed ${\sigma }_{L_{\mathrm{X}}}\sim {\displaystyle \frac{1}{\sqrt{C_{\mathrm{N}}}}}$, where C_N is the RASS counts from the parent catalogs. We applied further corrections to the MCXC cluster luminosities.

Firstly, we calculated K-correction factors to account for the redshifted source spectrum when observed in the observer reference frame. We derived these factors in two iterative steps in XSPEC using the scaling relations by Reichert et al. (2011) for the input temperatures, and the 49 updated redshifts. The changes that occurred are much smaller than σ_LX. Secondly, the L_X values were adjusted accordingly for the 49 clusters for which we used new redshifts. The uncertainties of the fitted X-ray z (both statistical and intrinsic) were propagated to the L_X uncertainties and added in quadrature to the already existing ones. Next, we corrected the MCXC L_X for changes in the soft-band X-ray absorption by using the combined molecular and neutral hydrogen column density values as described in Sect. 2.2. We first derived an absorbed L_X by reversing the absorption correction from Piffaretti et al. (2011), and then derived updated unabsorbed L_X values in XSPEC employing the cluster temperatures and metallicities derived in this work. Finally, the redshift-derived distances of nearby clusters might be biased by peculiar velocities. For the five most nearby clusters (50 h⁻¹ Mpc, z ≤ 0.0116) we used redshift independent distance measurements from NED (published within the last 20 years) to derive the L_X from the unabsorbed, k-corrected flux. The standard deviation of these distance measurements was propagated to the uncertainty of the luminosity. The average change in the distance compared to the redshift distances is ∼7%. For one cluster (S0851), no redshift independent distance was available and thus we adopted the redshift-derived distance but added an uncertainty of 250 km s^{−1 4} due to possible peculiar motions, which propagated in the σ_LX as well.

The comparison between the MCXC L_X values and the values used in this analysis is shown in Fig. 4. As seen there, for 301 clusters (96% of the sample) the change in L_X is ≤25%. Since the intention of this paper is to look for spatial anisotropies of the L_X–T relation in the sky we need to ensure that we do not introduce any directional bias through all of our L_X corrections. To this end, we compare the fraction L_X, ours/L_X, MCXC (the latter is the value given by MCXC) throughout the sky and we find it to be consistent within ±4%. Similar results are obtained if one considers the fraction L_X, MCXC/L_X, parent, where L_X, parent is the value given in the parent catalogs (more details in Appendix A.8).

The final L_X range of the clusters we use is L_X = (1.15 × 10⁴² − 3.51 × 10⁴⁵) erg s⁻¹ while the median value is 1.45 × 10⁴⁴ erg s⁻¹ (1.10 × 10⁴⁴ erg s⁻¹ for the excluded objects). The median σ_LX is 10.6%.

2.4. Cluster radius

The radii of the clusters were used for selecting the region within which the temperature and metallicity were measured. The R₅₀₀ values of the clusters as given in MCXC were determined using the MCXC L_X value and the X-ray luminosity-mass scaling relation L_X − M as given in Arnaud et al. (2010), which results in $L_{\mathrm{X}}\sim R_{500}^{4.92}$. Since we applied certain changes to the MCXC L_X values it is expected that also the respective R₅₀₀ should (slightly) change. Therefore, we used the same scaling relation to calculate the new R₅₀₀ in Mpc units. After that, the appropriate conversion to arcmin units was required, using the angular diameter distance (D_A). Since there is only a weak R₅₀₀ dependency on L_X when the latter changes only because of alternations in the N_H absorption, the new R₅₀₀ value does not significantly differ from the MCXC one (all changes ≤5%, since D_A remains fixed). However, when L_X changes because of a modification in the used z value then also D_A changes, as well as the normalized Hubble parameter E(z) = [Ω_m(1 + z)³ + Ω_Λ]^1/2 and the critical density of the Universe ρ_c(z), which are also included in the L_X − R₅₀₀ relation. This has a stronger impact on the final $R_{500}^{\mathrm{arcmin}}$ (R₅₀₀ in arcmin units) than the absorption case alone. Nonetheless, the new $R_{500}^{\mathrm{arcmin}}$ eventually changes by more than 10% only for 5 out of 313 clusters. At the same time, 294 clusters (>94% of the sample) show a ≤4% relative change in $R_{500}^{\mathrm{arcmin}}$. Therefore, the direct use of the MCXC $R_{500}^{\mathrm{arcmin}}$ values is practically equivalent to our new values.

2.5. Temperature

We determined the temperature of each cluster within the 0.2−0.5 R₅₀₀ annulus of every cluster in order to have self-consistent temperature measurements that reflect the other cluster properties (e.g., L_X) in a similar way. The cores of the clusters are excluded due to the presence of cool-cores, which significantly bias the temperature measurement and potentially increase the scatter of the L_X–T relation (e.g., Hudson et al. 2010). As previously stated, here we used the N_Htot values in order to fit the spectra and obtain T. Since the new R₅₀₀ values do not considerably vary compared to the MCXC ones, we used the latter for the spectra extraction with one exception: when the difference between the two values was >10% (only five clusters as stated above) then we used the redetermined R₅₀₀ value. Generally, the temperature shows only a weak dependance to such small changes in R₅₀₀ since the vast majority of the spectral extraction regions remains unchanged. The relative difference in the obtained temperature for these five clusters when we use both ours and the MCXC R₅₀₀ is by average ∼8%.

It has been shown that the Chandra and XMM-Newton telescopes have systematic differences in the constrained temperature values for the same clusters (S15). Thus, one has to take into account these biases when using temperature measurements from both telescopes. To this end, for the 76 clusters in our sample for which we use XMM-Newton data since we do not have Chandra data, we converted their measured temperatures to Chandra temperatures adopting the conversion relation found in S15. To further check the consistency of this conversion, we applied this test ourselves choosing 15 clusters in an XMM-Newton temperature range of 1.2−8.5 keV (same range as for the 76 XMM-Newton clusters) which have been observed by both telescopes. We constrained their temperatures with both instruments and we find that the best-fit relation provided by S15 still returns satisfactory results (Fig. A.8). The final temperature range of these 313 clusters is T = (0.83−19.24) keV with the median value being T = 4.5 keV, while the median uncertainty is σ_T = 4.9%. Clearly the temperature range is considerably wide which significantly helps the purposes of our study. The suspiciously high temperature of 19.23 keV occurs for the galaxy cluster Abell 2163 (A2163) when the Asplund et al. (2009) abundance table is used. A2163 also lies in a high absorption region with N_Htot = 2 × 10²¹ cm⁻². For other abundance tables, A2163 returns a temperature of ∼13−16 keV. This large difference is mainly driven by the phabs absorption model and the change of the Helium and Oxygen abundances. Generally, the average difference between the obtained temperature values between different abundance tables do not vary by more than ∼3−5% and thus, this cluster is a special case. For consistency reasons we use the 19.23 keV value. Excluding this cluster from the sample does not affect our results significantly since it is not an outlier in the L_X–T plane.

2.6. Metallicity

The metallicity of each cluster was determined simultaneously with the temperature. The two different telescopes are not shown to give systematically different metallicity values for the same clusters (although the intrinsic scatter of the comparison is relatively large), and thus no conversion between the XMM-Newton and Chandra values was needed. The metallicity range of the used sample is Z = (0.04−0.87) Z_⊙ except for two clusters with $Z=1.{53}_{-1.07}^{+1.50}{Z}_{\odot }$ and $Z=3.{86}_{-0.71}^{+0.61}{Z}_{\odot }$, where the metallicity determination of the latter is clearly inaccurate while the former is consistent with typical metallicity values within the 1σ uncertainties. Excluding these two clusters from our analysis has no significant effects on the derived results. Finally, the median value is Z = 0.37 Z_⊙ with the median uncertainty being σ_Z = 23.1%.

3. Data reduction and spectral fitting

3.1. Data reduction

The exact data reduction process slightly differs for the two instruments. For the Chandra analysis, we followed the standard data reduction tasks using the CIAO software package (version 4.8, CALDB 4.7.6). A more detailed description is given in Schellenberger & Reiprich (2017; S17). For the XMM-Newton analysis, we followed the exact same procedure as described in detail by R19. In a nutshell, every observation was treated for solar flares, anomalous state of CCDs (Kuntz & Snowden 2008), instrumental background and exposure correction. For both instruments, the X-ray emission peak was determined and used as the centroid for the spectral analysis⁵, while bright point sources (AGNs and stars), extended structures unrelated to the cluster of interest (e.g., background clusters) and extended substructure sources were masked automatically and later by hand in a visual inspection. For this analysis, the HEASOFT 6.20, XMMSAS v16.0.0 and XSPEC v12.9.1 software packages were used.

3.2. Background modeling

For the Chandra clusters, complementary to the S17 process, the ROSAT All-Sky survey maps in seven bands (Snowden et al. 1997) were used to better constrain the X-ray background components. The background value in each of the seven bands was determined within 1° −2° around the cluster.

For the XMM-Newton clusters, the only difference with the process described in R19 is the background spectra extraction region. The X-ray sky background was obtained when possible, from all the available sky region in the FOV outside of 1.6 R₅₀₀ from the cluster’s center. In this case, no cluster emission residuals were added in the background modeling. This was done mostly for the clusters located at z ≳ 0.1 which have a small apparent angular size in the sky. For most clusters, a partial overlap of the background extraction area with the 1.6 R₅₀₀ circle is inevitable and thus, an extra apec component to account for the cluster emission residuals was added during the spectral fitting, with its temperature and metallicity free to vary. The normalizations of the background model components were also left free to vary during the cluster spectra fitting as described in detail in R19.

3.3. Spectral fitting

For the spectral fitting, the same methods were used as in R19 (apec × phabs + emission and fluorescence lines) with only some small differences which are described here. Firstly, the 0.7−7 keV energy range was used for all spectral fittings for both instruments. This way we managed to exclude the emission lines close to 0.6 keV which originate from the Solar Wind Charge Exchange and cosmic X-ray background (S15, R19 and references therein). Moreover, we avoid the events produced by the fluorescent lines at 7.5, 8 and 8.6 keV which appear in the spectra of the pn detector of XMM-Newton. Furthermore, Chandra has a small effective area for energies higher than 7 keV. For all the spectral fits the Asplund et al. (2009) abundance table was used. Finally, for the 237 Chandra clusters the best-fit parameters of the spectral model were determined from an MCMC chain within XSPEC, while for the 76 XMM-Newton clusters the χ²-statistic was used.

4. The L_X–T scaling relation

For obtaining the best-fit values of the L_X–T relation parameters and comparing them for clusters located in different directions in the sky, we use a similar approach to M18. Here the strong dependance of the L_X on the cosmological parameters should be stressed again, combined with the fact that T can be measured without any cosmological assumptions (see appendix for the exact R₅₀₀ and T dependance on the chosen cosmology).

4.1. Form of the L_X–T scaling relation

We adopt a standard power-law form of the L_X–T scaling relation as shown below:

$$\begin{array}{c}\hfill \frac{L_{\mathrm{X}}}{{10}^{44}\mathrm{erg}{\mathrm{s}}^{-1}}E{(z)}^{-1}=A\times {\left(\frac{T}{4\mathrm{keV}}\right)}^B,\end{array}$$(1)

where the term E(z) = [Ω_m(1 + z)³ + Ω_Λ]^1/2 scales L_X accordingly to account for the redshift evolution of the L_X–T scaling relation. The scaling of the temperature term was chosen to be close to the median T = 4.5 keV. The exact constant scaling of the L_X values (10⁴⁴ erg s⁻¹) is not important since it is only a multiplication factor of the normalization. The exact scaling correction for the redshift evolution of the L_X–T relation [E(z)] is also not particularly significant for consistent redshift distributions and low-z samples like our own, as discussed later in the paper. In order to constrain the best-fit parameters the χ²-minimization method is used and applied to the logarithmic form of the L_X–T relation,

$$\begin{array}{c}\hfill log{L}_{\mathrm{X}}^{\prime }=logA+Blog{T}^{\prime }.\end{array}$$(2)

Here, $L_{\text{X}}^{\prime }$ and T′ are defined as

$$\begin{array}{c}\hfill L_{\mathrm{X}}^{\prime }={\displaystyle \frac{L_{\mathrm{X}}}{{10}^{44}\mathrm{erg}{\mathrm{s}}^{-1}}}E{(z)}^{-1}\mathrm{and}{T}^{\prime }={\displaystyle \frac{T}{4\mathrm{keV}}}·\end{array}$$(3)

4.2. Linear regression

The exact form of the χ²-statistic used to find the best-fit A, B, σ_int and H₀ values is given by

$$\begin{array}{c}\hfill {\chi }_{\mathrm{L}}^2={\displaystyle \sum _{i=1}^N}\frac{{(log\left[L_{\mathrm{X},\mathrm{obs}}^{\prime }\right]-log\left[L_{\mathrm{X},\mathrm{th}}^{\prime }(T^{\prime },\mathit{p})\right])}^2}{{{\sigma }_{logL,i}}^2+{B^2\times {\sigma }_{logT,i}}^2+{{\sigma }_{\mathrm{int}}}^2},\end{array}$$(4)

where N is the number of clusters used for the fit, $L_{\text{X},\text{obs}}^{\prime }$ and T′ are the measured luminosity and temperature values respectively (scaled as explained above), $L_{\text{X},\text{th}}^{\prime }$ is the theoretically expected value for the luminosity based on the measured temperature in addition to the fitted parameters p (A and B, or H₀). Furthermore, σ_{logL, T, i} are the Gaussian logarithmic uncertainties which are derived in the same way as in M18⁶, while σ_int (which was not included in M18) accounts for the intrinsic scatter of the relation.

The latter is fitted iteratively, starting from 0 and increasing step-by-step until there is a combination of p that gives ${\chi }_{\text{red}}^2\sim 1$, as in Maughan (2007), Maughan et al. (2012), Zou et al. (2016) etc. Under certain conditions, this procedure might return slightly underestimated σ_int values. However, this should not be a concern since the exact values of σ_int are not of particular importance for this analysis and they are only used to derive trustworthy parameter uncertainties from our χ² model.

Additionally, the 1σ uncertainties of the fitted parameters are based on the standard $\Delta {\chi }^2={\chi }^2-{\chi }_{\text{min}}^2$ limits (Δχ² ≤ 1 or 2.3 for one or two fitted parameters respectively). In the case of the slope being free to vary, the projection of the x-axis uncertainties to the y-axis also varies. This fitting method is comparable to the BCES Y|X fitting method described by Akritas & Bershady (1996).

Finally, we should stress that H₀ and A cannot be simultaneously constrained since they are degenerate. One can put absolute constraints only on the product $A\times H_0^2$. Therefore, one needs to fix one of the parameters to investigate the behavior of the other. In Sect. 5 we use a fixed H₀ = 70 km s⁻¹ Mpc⁻¹ to investigate the behavior of A. In Sect. 7 we fix A to its best-fit value and study the directional behavior of H₀ through the χ²-minimization procedure described above⁷.

4.3. Pinpointing anisotropies via sky scanning

With the purpose of studying the consistency of the fitted parameters throughout the sky and identifying specific sky patches that seem to show a significantly different behavior than the rest, we follow the method described below. We consider a cone of a given radius θ (we use θ = 45°, 60°, 75° and 90°) and we only consider the clusters that lie within this cone. For instance, if we choose a θ = 60° cone centered at (l, b) = (150°, 30°) then the subsample of clusters consists of all the clusters with an angular separation of ≤60° from these specific coordinates. By fitting the L_X–T scaling relation to these clusters, we obtain the normalization (or H₀), slope and intrinsic and total scatter for these clusters. The extracted best-fit value for the fitted parameter is assigned at these coordinates.

Shifting this cone throughout the full sky in steps of Δl = 1° and Δb = 1° in Galactic coordinates⁸, we can obtain the desired parameter values for every region of the sky. We additionally apply a statistical weighting on the clusters based on their angular separation from the center of the cone. This is given by simply dividing their uncertainties by $cos({\displaystyle \frac{{\theta }_1}{\theta }}\times 90°)$, where θ₁ is the above-mentioned angular separation. Hence, the weighting cos term is calibrated in such way that it shifts from 1 to 0 as we move from the center of the cone to its boundaries, independently of the angular size of the cone. This enlargement of the uncertainties results in an artificial decrease of the σ_int which is not of relevance here since, as explained before, σ_int mostly acts as a nuisance parameter. Nevertheless, we perform tests to ensure this does not bias our results, as explained in Sect. 9.2.

All the A maps are plotted based on the A/A_all value, where A_all is the best-fit A when all the clusters are used independently of the direction. Finally, the maps have the same color scale for easier comparison, except for the θ = 45° cone maps for which the color scale is enlarged for better visualization.

4.4. Statistical significance and sigma maps

With the desired best-fit values and their uncertainties for every sky region at hand, it is easy to identify the direction that shows the most extreme behavior and assess the statistical significance of their deviation. For quantifying the latter in terms of number of sigma for two different subsamples we use:

$$\begin{array}{c}\hfill \mathrm{No}.\mathrm{of}\sigma ={\displaystyle \frac{{\mathit{p}}_1-{\mathit{p}}_2}{\sqrt{{\sigma }_{{\mathit{p}}_1}^2+{\sigma }_{{\mathit{p}}_2}^2}}},\end{array}$$(5)

where p_1, 2 are the best-fit values for the two different subsamples and ${\sigma }_{{\mathit{p}}_{1,2}}^2$ are their uncertainties⁹.

Each time we constrain the anisotropic amplitude of the most extreme dipole in the sky, while we also compare the two most extreme regions in terms of the fitted parameter, regardless of their angular separation. This is done by calculating the statistical deviation (in terms of σ) between all the different cone subsamples. The two sky regions for which the largest deviation (highest no. of σ) is found between them, are the ones reported in the following sections as “the most extreme regions”. In addition, a percentage value (%) is displayed next to each σ deviation. This value comes from the difference of these two extreme regions over the best-fit value for the full sample.

In order to create the significance maps, we use Eq. (5) to compare the best-fit result p₁ of every cone with the best-fit result p₂ of the rest of the sky. The obtained sigma value is assigned to the direction at the center of the cone. All the significance maps have the same color scale for easier comparison. For the majority of cases, the two most extreme regions as defined above match the highest σ regions in the significance maps.

Finally, as an extra test we also create 10⁵ realizations using the bootstrap resampling method in order to check the probability of the extreme results to randomly occur independently of the sky direction. The followed procedure is described in detail in Sect. 9.2. Using all these estimates, one can determine the consistency with what one would expect in an isotropic universe.

5. Results

5.1. The L_X–T scaling relation for the full sky

Before we search for apparent anisotropies in the sky we constrain the behavior of the L_X–T scaling relation for the full sample. We do not account for any selection biases in our analysis since we believe that their effects are not important for this work. This is because we wish to study the relative L_X–T differences between different sky regions (or from the overall best-fit line). If we indeed corrected for selection effects we would constrain the “true” underlying L_X–T relation which would not represent our data (but the true distribution). This might cause wrong estimates for the relative L_X–T differences. Therefore, we need to constrain the relation that describes our 313 clusters best. Nevertheless, in Sect. 6.4 we discuss the possible effects of selection systematics and find that there is no indication that they compromise our results.

We use the aforementioned 313 clusters and fit Eq. (1) obtaining the best-fit normalization and slope of the L_X–T relation as well as its intrinsic scatter. The results are:

$$\begin{array}{c}\hfill A_{\mathrm{all}}=1.{114}_{-0.040}^{+0.044},B_{\mathrm{all}}=2.102\pm 0.064,\mathrm{and}{\sigma }_{\mathrm{int}}=0.242\mathrm{dex}.\end{array}$$(6)

The statistical uncertainties for A and B are limited to ∼3−4% which highlights the precision of our results based on the number and the quality of the data, combined with the large covered temperature range of the clusters. Moreover, the total scatter (statistical + intrinsic) of the L_X–T is σ_total = 0.262 dex, which means that the statistical uncertainties of the clusters contribute to only 7% of the total scatter. The L_X–T fit of our sample is displayed in Fig. 5 (top panel).

The best-fit slope is slightly steeper than the expected value from the self-similar model. Naively, the slope best-fit value (∼2.1) might seem surprising since most studies find a L_X–T slope of B ∼ 2.3−3.6 (see references in Sect. 1). However, the exact value depends on multiple aspects such as the energy range for which L_X was measured, the instrument used, the sample selection (when no bias correction is applied), the temperature distribution of the used sample, the cluster radius within which parameters were measured etc. Generally, it is expected that bolometric L_X values return a steeper slope than soft band L_X values, such as the 0.1−2.4 keV band we use. This happens due to the fact that the bolometric emissivity ϵ of the ICM for thermal bremsstrahlung (which is the dominant emission process for T ≳ 3 keV) is $ϵ\propto n_{\text{e}}^2{T}^{0.5}$ (where n_e is the electron density), while ϵ in the soft band (0.1−2.4 keV as used here) is rather independent of T for T ≳ 3 keV ($ϵ\propto n_{\text{e}}^2$). Therefore, one very roughly expects the slope of the L_X–T relation to be smaller by ∼0.5 in the 0.1−2.4 keV band; that is L_X ∼ T^1.5 in the self-similar case. In general, the best-fit L_X–T relation tends to change slightly when one corrects for selection biases (see references in Sect. 1 about the L_X–T relation).

If we now convert all the measured temperatures to XMM-Newton temperatures¹⁰ using the relation given in S15, we obtain a slope of B ∼ 2.38, shifting by 2.8σ compared to our main result, while the normalization remains the same. This result is consistent with the previously reported values that used the 0.1−2.4 keV luminosities (e.g., Chen et al. 2007; Eckmiller et al. 2011; Lovisari et al. 2015)¹¹. In the bottom panel of Fig. 5 the 68.3% and 99.7% confidence levels (1 and 3σ respectively) of the fitted parameters are shown for Chandra-converted temperatures and XMM-Newton-converted temperatures. It should be clear that this is done just for the sake of comparison and that the Chandra temperatures are used for the rest of the paper.

A comparison between our results and the derived L_X–T scaling relation from other works is also shown in the top panel of Fig. 5. We note that L_X corresponds to the 0.1−2.4 keV energy band for all the compared studies, while the results from the BCES (Y|X) fitting method were used when available. Additionally, the L_X–T results for the full samples were used without any bias corrections. From this comparison, it is clear that all the derived results agree in the L_X–T normalization value. In terms of the slope, our Chandra fit is more consistent with the high-T part of the distribution, while our XMM-Newton fit is quite similar to the results of Eckmiller et al. (2011) and Lovisari et al. (2015).

5.2. 1-dimensional anisotropies

As a first test for the potentially anisotropic behavior of our galaxy cluster sample, we recreate the normalization against the Galactic longitude plot as presented in M18 (Fig. 3 in that paper). For this test, we consider regions centered at l with a width of Δl = 90°. At the same time, the whole Galactic latitude range b ∈ (− 90°, + 90°) is covered by every region.

Firstly, we allow both A and B to vary simultaneously. The behavior of these two parameters as functions of the Galactic longitude are displayed in Fig. 6. One sees that the slope remains relatively constant throughout the sky, varying only by 18% from its lowest to highest value, and with a relatively low dispersion. Also, the largest deviation between any two independent sky regions is limited to 1.8σ. No obvious systematic trend in the slope as a function of the galactic longitude can be seen since all the regions return slope values consistent with the full sample at ≤1.2σ. At the same time, this variation for the normalization reaches 31% with a higher dispersion and a clear trend with galactic longitude, while the strongest tension between two independent sky regions appears to be 3.2σ.

Based on these results, the slope is kept fixed at the best-fit value for the whole sample and only the normalization of the L_X–T relation is free to vary. In the top left panel of Fig. 7 the best-fit normalization value A for every region is displayed with respect to the best-fit A_all for the full sky (all 313 clusters). The same is also done for ACC and XCS-DR1 with the results displayed in the top and bottom left panels of Fig. 7 respectively. The only difference with the M18 results for these two samples is that here the intrinsic scatter term is taken into account as well during the fitting as shown in Eq. (4).

Surprisingly enough, the pattern in the behavior of the L_X–T normalization for our sample strongly resembles the results of both ACC and XCS-DR1, despite being almost independent with XCS-DR1, sharing only ∼30% of the clusters with ACC and following different analysis strategies. Specifically, the region with the most anisotropic behavior compared to the rest of the sky (2.9σ significance) is the one with the lowest A lying within l ∈ [210°, 300°]. This region exactly matches the findings of M18 for XCS-DR1, while the lowest A region for ACC is separated by 40°. Here we should remind the reader that the 313 clusters we use share only three common clusters with XCS-DR1 and 104 with ACC as these samples were used in M18. The opposite most extreme behavior (highest A) is detected in l ∈ [−20°, 70°] (same brightest region in ACC as well, 25° away from XCS-DR1’s brightest region) with a deviation of 2.5σ compared to the rest of the sky and 3.4σ compared to the lowest-A region, which is similar to the two other samples.

In order to verify that the observed behavior is not caused by the few common clusters between our sample and ACC or XCS-DR1, we exclude all 104 of them from our sample and repeat the analysis. The result is shown in the bottom left panel of Fig. 7. One can see that this systematic trend persists and does not significantly depend on the common clusters between the two samples. The region with the largest deviation from the rest of the sky remains the same as for the full sample with an even higher significance of 3.1σ. This striking similarity between the three different samples in the 1D search for anisotropies should be investigated in more depth in order for its exact reason to be identified.

5.3. 2-dimensional investigation

In order to identify the exact regions with the highest degree of anisotropy, we should consider every possible direction in the sky. Different size regions should be considered as well, thus systematic behaviors can be detected. To this end, we use scanning cones (solid angles), as described in Sects. 4.3 and 4.4. The slope is fixed to the best-fit value since the variations of the normalization are much stronger. This choice does not bias our results, as shown in Sect. 6.5.

5.3.1. θ = 90° cone

To begin with, we choose a scanning cone with θ = 90°, meaning we divide the sky in all the possible hemisphere¹² combinations. The lowest number of clusters in any hemisphere is 109 toward the (l, b) = (150°, −2°) direction, with 204 clusters located in the opposite hemisphere. Constraining A for every hemisphere, one obtains the A and significance color maps displayed in the top left panels of Figs. 8 and 9 respectively.

As shown in the plots, there is mainly one low A region within ∼20° from (l, b) ∼ (270°, −5°) with a rather strong behavior. There is also one high A peak. It is noteworthy that the two most extreme regions in the map (deep purple and bright yellow) are located close to the Galactic plane (within ∼20°), where there are no observed clusters. The clusters toward the Galactic center in particular seem to be overluminous compared to other sky regions. Of course the color differences are not visually strong here since the color scale was chosen based on the largest deviations, appearing in later maps.

In detail, the most extreme hemispheres are found at (l, b) = (272°, −8°) with A = 1.062 ± 0.048 and at (l, b) = (47°, + 22°) with A = 1.236 ± 0.047. The angular separation between them is 135° while they deviate by 2.59σ (16 ± 6%) from each other. Although they are not completely independent, the contribution of the common clusters is not the same to both subsamples due to the applied statistical weighting based on the distance of every cluster from the center of the cone. The most extreme dipole (2 independent subsamples separated by 180° in the sky) appears at (l, b) = (230°, −20°) with a significance of 1.90σ. This dipole is separated by 75° from the CMB dipole, although a dipole interpretation is obviously not reflecting the maximum apparent anisotropies in that case.

5.3.2. θ = 75° cone

If an anisotropy toward one direction exists, the clusters lying close to that direction would be the most affected ones and as we move further away from that direction, the anisotropic effect on clusters would fade. Therefore, such anisotropic behaviors are better studied if one uses smaller solid angles in the sky. To this end, we decrease the radius of the scanning cone first to θ = 75°. Indeed, the fluctuations of A as well as the significance of the anisotropies increase, while the general behavior of the directional anisotropies in the θ = 75° map however remains relatively unchanged compared to the previous map with a larger cone. The results are displayed in the top right panels of Figs. 8 and 9.

For the θ = 75° cones, A varies from A = 0.999 ± 0.050 at (l, b) = (274°, −22°) to A = 1.288 ± 0.061 toward (l, b) = (17°, −9°). These two regions are separated by 99° and deviate from each other by 3.64σ (26 ± 7%). Furthermore, the most significant dipole that appears is the one centered at (l, b) = (263°, −21°) at 3.22σ, 68° away from the CMB dipole.

5.3.3. θ = 60° cone

Further decreasing the size of the solid angles, we use θ = 60° cones. We see that the behavior of the sky regions suffers some changes whilst staying generally consistent with the previous results. The most prominent change is the existence of a low A region close to (l, b) ∼ (120°, + 20°), although its statistical significance (as displayed in Fig. 9) is lower than the other, main low A region since it only contains ∼45 clusters. Another change in the 60° map is that the brightest part of the sky is shifted toward (l, b) ∼ (170°, −10°). However, as one can clearly see in Fig. 9, the most statistically significant region with a high normalization remains in the same area as in the previous cases, namely toward (l, b) = (34°, + 4°) with 78 clusters and A = 1.346 ± 0.069.

At the same time, the lowest normalization value A = 0.940 ± 0.051 is located at (l, b) = (281°, −16°) (84 clusters). The most extreme regions deviate from each other by 4.73σ (36 ± 8%), which constitutes a considerably strong tension. The most extreme dipole in this case is found toward (l, b) = (260°, −36°) with a statistical significance of 3.77σ.

If we now exclude these two most extreme low and high A regions and their 159 individual clusters from the rest of the sky, we are left with 154 clusters. Performing the fit on these clusters, we obtain A = 1.138 ± 0.048. We see that the rest of the sky is at a 2.49σ tension with the bright region toward (l, b) = (34°, + 4°), and at a 2.85σ tension with the faint region toward (l, b) = (281°, −16°). Thus, the anisotropic behavior of the faint region is somewhat more statistically significant than the behavior of the bright region.

5.3.4. θ = 45° cone

The last cone we use has θ = 45°. Since there are many regions mostly close to the Galactic plane, with fewer clusters than needed in order to obtain a trustworthy result, we enforce an extra criterion. We only consider regions with ≥35 clusters¹³. The A and σ maps are shown in the bottom right panels of Figs. 8 and 9 respectively. The white regions show the regions without enough clusters for a reliable fit. The most extreme regions are found toward (l, b) = (280°, + 1°) (42 clusters) and (l, b) = (32°, + 14°) (40 clusters) with A = 0.822 ± 0.067 and A = 1.413 ± 0.095. The statistical discrepancy between the rises to 5.08σ (53 ± 10%) being the most statistically significant result up to now.

Additionally, the most extreme dipole is centered at (l, b) = (255°, −53°) with a significance of 4.22σ. However, it should be beared in mind that many regions that appeared to have the maximum dipoles for other cones are excluded now due to low number of clusters. This could lead the maximum dipole to shift toward lower Galactic latitudes on the low normalization side. Moreover, due to the low number of clusters in these regions the results are more sensitive to outliers, especially when these outliers are located close to the center of the regions where they have more statistical weight than other clusters. Nevertheless, the large statistical tension cannot be neglected.

Excluding once again the two extreme regions from the rest of the sample we are left with 234 clusters which have a best-fit of A = 1.107 ± 0.041. Thus, they are in a 2.96σ tension with the brightest region and in a 3.63σ with the faintest region. Once again, the region toward (l, b) = (280°, + 1°) seems to be more anisotropic than the one toward (l, b) = (24°, + 16°).

5.3.5. Overview of results

As a summary of the above, we identify the clear existence of a region with galaxy clusters appearing systematically fainter than expected based on their temperature measurements. This region is roughly located at (l, b) ∼ (277° ±5°, −11° ±12°). On the contrary, the systematically brightest region is found toward (l, b) ∼ (32° ±15°, + 8° ±17°). Their angular separation in the sky is ∼115°. The statistical tension between these two regions rises significantly while narrower cones are considered, reaching ∼5σ for the smaller cones. The same is true for the dipole anisotropies, going up to ∼4σ. Interestingly enough, the same behavior for this sky patch is also detected for ACC and XCS-DR1 (see Sect. 8). Another interesting trend is the systematically bright region at (l, b) ∼ (175° ±15°, + 5° ±20°, which appears to have the same behavior in all three maps with the larger scanning cones. Unfortunately, not enough available clusters lie there for the 45° map to return reliable results.

The most statistically significant dipole anisotropy is consistently found toward (l, b) ∼ (253° ±13°, −32° ±15°), lying ∼30° ±25° away from the systematically fainter sky region. Finally, the correlation of these results with the CMB dipole is not strong since the faintest regions of our analysis are found ∼55° −75° away from the CMB’s corresponding dipole end, while the strongest anisotropic dipoles of the L_X–T relation are located ∼70° −85° away from the CMB one.

6. Possible X-ray and cluster-related causes and consistency of anisotropies

Galaxy clusters are complex systems where many aspects of physics come into play when one wishes to analyze them. Thus, we have to investigate if the apparent anisotropies are caused by any systematic effects. With a purpose of trying to identify the reason behind these strong L_X–T anisotropies, we perform an in-depth analysis using different subsamples of the 313 clusters which are chosen based on their physical properties. If the best-fit L_X–T relation of galaxy clusters significantly differs for clusters with different physical parameters (e.g., low and high T or Z clusters, different N_Htot values etc.), a nonuniform sky distribution of such clusters could create artificial anisotropies.

6.1. Excluding galaxy groups and low-z clusters

It has been shown that the low-T clusters (mainly galaxy groups) can sometimes exhibit a slightly different L_X–T behavior compared to the most massive and hotter systems (e.g., Lovisari et al. 2015, and references therein). We wish to test if this possibly different behavior has any effects on the apparent anisotropies.

Hence, we first excluded all the systems below T ≤ 2.5 keV. Moreover, all the clusters within ∼130 Mpc (z ≤ 0.03) were excluded in order to avoid the peculiar velocity effects on the measured redshift (the vast majority of these clusters are already excluded based on the T ≤ 2.5 keV limit). This resulted in the exclusion of 67 objects.

We applied the necessary correction to convert our heliocentric redshifts to “CMB frame” redshifts. This conversion is not expected to cause any significant changes in our results for two reasons. Firstly, the spatial distribution of our sample is rather uniform, therefore only ∼25% of this subsample’s clusters are located within 30° from the CMB dipole for which this correction might have a notable impact. Secondly, due to the low-z cut we apply here, the CMB frame redshift correction is much smaller than the cosmological recession velocity (≲4%). Hence, the final propagated correction to the L_X values is far less than the observed anisotropies. Nevertheless, we transformed the redshifts for the sake of completeness.

When we fit the L_X–T relation to all the 246 clusters with T >  2.5 keV and z >  0.03, we obtain the following best-fit values:

$$\begin{array}{c}\hfill A=1.{114}_{-0.041}^{+0.047},B=2.096\pm 0.078,\mathrm{and}{\sigma }_{\mathrm{int}}=0.218\mathrm{dex}.\end{array}$$(7)

It is noteworthy that A and B remain unchanged compared to the case where all the 313 clusters are considered. This indicates that a single power law model can be an efficient option for fitting our sample. The most clear difference of this subsample fitting is the decrease of the intrinsic scatter by 10%. The total scatter also goes down by the same factor (σ_tot = 0.236 dex). The 3σ solution spaces for the entire sample and for these 246 clusters are displayed in Fig. A.6, being entirely consistent.

Performing the 2D scanning of the full sky using θ = 75°, the A map shown in the bottom panel of Fig. 10 is produced. A similar pattern with the previous maps persist, although there are some changes. The main differences are that the –statistically insignificant– bright region toward (l, b) ∼ (170°, −10°) vanishes whilst the behavior of the faint region toward (l, b) ∼ (120°, + 10°) seems to be amplified. Despite of that, the most statistically significant low-A regions approximately remains in the sky patch that was found before, toward (l, b) = (288°, −35°) with A = 1.016 ± 0.045 (110 clusters). The most extreme high-A sky region is again consistent with our previous findings, lying at (l, b) = (10°, + 16°) with A = 1.371 ± 0.061 (113 clusters). The statistical discrepancy between these two results is 4.68σ (32 ± 7%), not being alleviated by the exclusion of these groups and local clusters. Their angular separation in the sky is 93°. The most extreme dipole for this map is found toward (l, b) = (196°, −34°) with 3.27σ, shifted compared to the previously found most extreme dipole regions by ∼45° ±28°.

To further scrutinize the effects of low-T systems on our results, as well as the effects of local clusters and their peculiar velocities, we wish to restrict our sample even more by expanding the lower limits of T and z. To this end, we excluded all the 115 objects with T <  3 keV or z <  0.05 (∼210 Mpc). This left us with 198 clusters. The best-fit results are:

$$\begin{array}{c}\hfill A=1.{172}_{-0.046}^{+0.053},B=2.049\pm 0.077,\mathrm{and}{\sigma }_{\mathrm{int}}=0.205\mathrm{dex}.\end{array}$$(8)

The best-fit L_X–T relation slightly changes compared to the full sample results, but remains consistent within 1.1σ. At the same time, σ_int further decreases, being 15% lower than the full sample’s σ_int. In Fig. A.6, the comparison between the 3σ solution spaces for the full sample and for these 198 clusters is displayed. In the panel of Fig. 11 the A map is displayed for this subsample of clusters, with a θ = 75° cone. The significance map is shown in the bottom panel of the same figure.

The behavior of A throughout the sky remains consistent with the previous results, even after excluding more low-T clusters and using only clusters with z >  0.05 with CMB-frame z values. The lowest A = 1.081 ± 0.054 is found toward (l, b) = (286°, −36°) (85 clusters) while the highest A = 1.445 ± 0.070 is located toward (l, b) = (9°, + 15°) (91 clusters). The statistical tension between these two results is 4.12σ (31 ± 8%). The most extreme dipole on the other hand is centered toward (l, b) = (223°, −47°) with a relatively low significance of 2.27σ. This highlights the fact that the most extreme behavior in the sky is not found in a dipole form, and this becomes more obvious as we go to higher redshifts. Consequently, it is quite safe to conclude that this anisotropic behavior is caused neither by the galaxy groups or the local clusters nor by the use of heliocentric redshifts.

6.2. Different cluster metallicities

A slightly nonsimilar behavior of the L_X–T relation for varying metallicities of clusters can be expected mainly due to two factors. Firstly, in the parent catalogs from which our cluster sample has been constructed, the conversion of the count-rate to flux was done by using a fixed metal abundance of 0.3 Z_⊙. When the true metallicity of a cluster deviates from this fixed value, small biases can propagate in the flux and luminosity determination. In general, the measured luminosity of clusters with Z >  0.3 Z_⊙ might be eventually slightly underestimated. However, this overestimation is only minimal. For instance, for ΔZ ∼ 0.4 Z_⊙ between fixed and true Z, the final flux changes by ∼0−2%, where the exact change depends on the other cluster parameters, such as the temperature.

The second and most important factor is that clusters with higher Z values tend to be intrinsically brighter when the rest of the physical parameters are kept constant. This can be shown through an apec model simulation in XSPEC. Even for a small deviation of ΔZ ∼ 0.1 Z_⊙ the flux of a cluster can fluctuate by ≳17% for a cluster with T ≲ 1 keV, while this fluctuation becomes only ≲1% for a cluster with T ≳ 8 keV. Therefore, a randomly different metallicity distribution between different sky regions could in principle cause small anisotropies. However, in order for the observed anisotropies to be purely caused by that, strong inhomogeneities in the metallicity distribution should exist, which, if detected, would be a riddle of its own.

6.2.1. Core metallicities within 0−0.2 × R₅₀₀

Galaxy clusters do not show a single metallicity component. Since we wish to focus first on the effects that a varying metallicity could have on the luminosity, we consider the metallicity of the core of the cluster (Z_core, where by “core” we mean 0−0.2 × R₅₀₀) from where the bulk of the X-ray luminosity comes from. It is also expected that the clusters with the higher Z_core values would have a higher fraction of cool-core members, which are generally more luminous than non cool-core clusters for the same T (e.g., Mittal et al. 2011).

In order to investigate the behavior of the L_X–T relation as a function of the metallicity of the galaxy clusters, we divided our sample into three subsamples based on their Z_core value. Our only criterion for this division was the equal number of clusters in each subsample. These subsamples are 105 clusters with Z_core ≤ 0.452 Z_⊙, 104 clusters with 0.452 Z_⊙ <  Z_core ≤ 0.590 Z_⊙ and 104 clusters with Z_core >  0.590 Z_⊙. For each subsample, we perform the fitting letting A and B to vary. The following results are not particularly sensitive to the exact Z_core limits.

The 1σ solution spaces for each subsample are shown in Fig. 12. One can see that all the three subsamples share a very similar L_X–T solution. The maximum statistical deviation of ∼1.12σ is found between the two subsamples with the lowest and highest Z_core, with the latter being slightly more luminous on average. Furthermore, the intrinsic scatter for the two subsamples with the lower Z_core is σ_int ∼ 0.260 dex while for the high-Z_core subsample is σ_int ∼ 0.197 dex.

Although it is not expected that the high Z_core subsample would cause any apparent A anisotropies with a possibly nonhomogeneous spatial sky coverage, for the sake of completeness we excluded all the 104 clusters with Z_core >  0.590 Z_⊙ and scanned the sky again with a 75° radius cone. The produced A and significance maps are illustrated Fig. 13. The obtained directional behavior of A completely matches the results of the full sample. The lowest A = 0.927 ± 0.064 and highest A = 1.274   ±   0.071 are found toward (l, b) = (264°, −18°) (83 clusters) and (l, b) = (30°, + 23°) (88 clusters) respectively. Their deviation is 3.63σ (32 ± 9%), staying unchanged despite the smaller number of available clusters. The most extreme dipole is found toward (l, b) = (261°, −20°) with 3.41σ significance.

6.2.2. Outer metallicities within 0.2–0.5 × R₅₀₀

The metallicity Z_out of the 0.2−0.5 × R₅₀₀ annulus might not affect the final L_X as strongly as the core metallicity. However, it could in principle correlate with the measured temperature of a galaxy cluster since these two quantities were fitted simultaneously. To check if there is an inconsistent L_X–T behavior based on Z_out we follow the same procedure as for Z_core, dividing the full sample into three subsamples similarly with before. These subsamples are 105 clusters with Z_out ≤ 0.320 Z_⊙, 104 clusters with 0.320 Z_⊙ <  Z_out ≤ 0.426 Z_⊙ and 104 clusters with Z_out >  0.426 Z_⊙.

In the top panel of Fig. 14 the 99.7% (3σ) solution spaces for the three subsamples are shown. It is obvious that the 104 clusters with the highest Z_out share a significantly different L_X–T solution than the 105 clusters with the lowest Z_out. The statistical deviation between these two subsamples is ∼4.3σ. However, as shown in the bottom panel of Fig. 14, the main source of deviation are the local, low temperature groups. Excluding objects with T <  2 keV and z <  0.02, the deviation between low and high Z clusters drops to 2.5σ, which is still a nonnegligible tension. At the same time, the medium Z_out subsample seems to be consistent with the low Z subsamples while also being in tension with the high Z clusters. Furthermore, the intrinsic scatter remains similar for all three Z_out subsamples (∼0.230−0.245 dex).

With the purpose of examining whether the strong A anisotropies are affected in any way by these Z-dependent different L_X–T behaviors, we once again excluded all the 104 clusters with Z_out >  0.426 Z_⊙ and performed the usual sky scanning with a 75° radius cone. In Fig. 15 the results are displayed.

The similarity with the full sample θ = 75° result is striking. The lowest and highest A directions are (l, b) = (270°, −14°) (80 clusters) and (l, b) = (24°, + 15°) (93 clusters). The direction (l, b) = (174°, −12°) (58 clusters) is actually brighter by ∼3% but its statistical significance is lower, which is similar to the results of previous maps. The A values of the most extreme regions are A = 1.035 ± 0.069 and A = 1.390 ± 0.066 respectively with a tension of 3.72σ (30 ± 8%), not relieved despite the exclusion of the Z_out subsample with the significantly different L_X–T behavior. The most extreme dipole is located toward (l, b) = (265°, −16°) but with a lower statistical significance of 2.26σ. Consequently, the derived anisotropies persist when the high Z clusters are excluded, while the significance of the dipolar anisotropy drops by ∼1σ compared to the full sample results.

6.3. Absorption correction

Another possible systematic effect resulting in the observed anisotropies could be the inaccurate treatment of the N_Htot column density correction in our apec model. This could lead to systematic differences in the L_X–T values of clusters in directions with different N_H. Since also the most extreme regions always lie within 35° from the Galactic plane, we have to ensure that the apparent anisotropies are not caused by such effects. There are two main cases for which a systematic bias could be introduced through the absorption correction and they are described in the following subsections.

6.3.1. Consistency throughout N_Htot range

The first case is that the N_Htot value does not trace the true absorption consistently throughout the full N_Htot range. Thus, clusters in regions with different amounts of hydrogen get a systematically different boost in their L_X–T values after the applied correction.

This can be easily checked by comparing the L_X–T scaling relation for the clusters with low and high N_Htot. To this end, we divided our sample into three subsamples of equal size based on their N_Htot values. These samples are the 105 clusters with N_Htot ≤ 2.53 × 10²⁰ cm⁻², the next 104 clusters with N_Htot ≤ 5.16 × 10²⁰ cm⁻² and finally the 104 clusters with N_Htot >  5.16 × 10²⁰ cm⁻².

We fit the full L_X–T relation for these three independent subsamples. As shown in the top panel of Fig. 16, clusters in high and low N_H regions show completely consistent L_X–T behaviors with each other. The only noticeable difference between the three subsamples is their intrinsic scatter. Going from the low to the high N_Htot subsamples, the intrinsic scatter is σ_intr = 0.214, 0.242 and 0.258 dex respectively. This is not surprising since the high N_Htot clusters undergo stronger corrections based on the molecular hydrogen column densities of W13. However, one should not forget that these molecular hydrogen values are approximations and thus some random scatter around the true values is expected, which then propagates to the L_X values. In any case, this does not constitute any source of N_Htot-related bias since the overall A and B behavior is similar for different N_Htot values (see also Sects. 6.4.1 and 6.4.2).

As an extra test, we excluded the 104 clusters with the highest absorption (N_Htot >  5.16 × 10²⁰ cm⁻²) and repeated the 2D sky scanning with θ = 75° in order to see if we observe the same anisotropies. The results are shown in Fig. 17. The previously detected anisotropic behavior persists, with the lower A = 0.879 ± 0.059 (65 clusters) found toward (l, b) = (242°, −27°), which is consistent within ∼30° from the previous findings. The brightest part of the sky remains unchanged compared with the full sample case, namely toward (l, b) = (35°, −15°) (93 clusters) with A = 1.368   ±   0.069. These two regions share a statistical tension of 5.39σ (45 ± 8%), the most statistically significant result we found up to now. The most extreme dipole anisotropy on the other hand is found toward (l, b) = (221°, −33°) with 3.55σ.

Subsequently, the detected > 3.5−4σ apparent anisotropies not only do not result due to the different amounts of absorbing material throughout the sky and its effects on X-ray photons, but they significantly increase to a > 5σ level when the 104 clusters with the highest absorption are excluded. This is mostly due to the decrease of the intrinsic scatter of the clusters left, which leads to a decrease in the final A uncertainties.

6.3.2. Extra absorption from undetected material or varying metallicity of the Galactic material

The second case is that the exact amount of X-ray absorbing material is not accurately known and a higher or lower absorption correction is needed than the one applied. Such problems could occur for example if not all the absorbing material in the line of sight of a galaxy cluster has been detected by the radio surveys such as LAB, either because it is outside of the velocity range of the radio survey or for other unknown reasons (e.g., more than expected hydrogen in ionized or molecular form).

Another possible reason could be the varying metal abundance of the ISM throughout the Galaxy. The applied X-ray absorption correction is mostly applied as this: the amount of hydrogen detected is used as a proxy for the total amount of absorbing material that exists toward a given direction. The elements of this material that contribute the most in the absorption of the X-ray photons are helium¹⁴ and metals such as oxygen, neon, silicon etc. Based on the detected N_Htot value, a Solar metal abundance is assumed for the Galactic interstellar medium (ISM) in every direction in order to quantify the number of metals absorbing X-ray radiation. However, throughout the Galaxy the true metal abundance might diverge from this approximation since there are metal-rich and metal-poor regions. Consequently, the same amount of detected hydrogen could correspond to different amounts of X-ray absorption from metals, which is not taken into account by our current absorption correction models. It needs to be checked if the apparent anisotropies could in principle be caused by such effects.

In order to test this, one can estimate the needed absorption using two ways. Firstly, one can calculate the necessary “true” N_Htot in order to fully explain the observed anisotropies. Secondly, one can fit the extracted X-ray cluster spectra and leave N_Htot to vary. Then, the obtained best-fit N_Htot can be compared with the ones we use, which come from W13.

Necessary N_Htot to fully explain L_X–T anisotropies. Any existence of X-ray absorbing “dark clouds” in certain parts of the sky could potentially explain the observed anisotropies. Approximately quantifying how much extra (or less) hydrogen column density would cause such an effect, one could see if this value can realistically be missed by radio surveys. Since any such clouds is unlikely to cover significantly large portions of the sky (more than 90° width) it is more appropriate to first look for them in the smallest radius cones.

To this end, we consider the region with the lowest normalization A for the θ = 45° cone at (l, b) = (280°, + 1°) which includes 42 clusters. Its A value is 29 ± 7% lower than the rest of the sky (the fitting for the rest of the sky is performed without any distance weight) and therefore its clusters would need to be more luminous by the same degree in order to be consistent with an isotropic behavior. To quantify how much extra N_Htot is required to make these clusters more luminous by ∼29% and explain the apparent anisotropy, we performed the following:

We selected 25 clusters from the region of interest with varying temperatures and metallicities. For each cluster we used an apec model in XSPEC, reproducing the current absorbed L_X value. For several N_Htot values we found the new unabsorbed L_X. For the same clusters, we refit their X-ray spectra and constrained the new T for the same N_Htot values as above. For every cluster and for every N_Htot change compared to the W13 values, we thus knew the relative change of L_X and T compared to the standard values. Next, we were able to find the average relative change of L_X and T for every tested N_Htot. Of course this change is not identical for every cluster since it depends on the exact T and Z. However, the actual average value of this change (∼20%) is much larger than its variation between clusters with different properties (±4%). Assuming the slope to be B = 2.102, we could obtain the relative change of A for every sky region based on the tested N_Htot values. Consequently, we found how much extra or less N_Htot one would actually need toward the apparently anisotropic sky regions in order to explain their behavior.

We find that an extra N_Htot = 3.3  ±  0.9  ×  10²⁰ cm⁻² is required in order to make the above-mentioned low A region consistent with the rest of the sky. For its 37 clusters, the average N_Htot is ∼7 × 10²⁰ cm⁻². Thus, the final N_Htot which would explain the low A value of this region is ∼48% larger than its current value. If we express this difference in terms of metal abundance of the existing Galactic ISM (and not just larger amounts of ISM material), the absorbing elements toward that direction should have a metallicity of Z ∼ 1.5 Z_⊙ to create such apparent anisotropies due to extra absorption.

The same analysis for the bright region toward (l, b) = (24°, + 16°) (which deviates by 31 ± 10% from the rest of the sky) yields that 3.5 ± 1.1 × 10²⁰ cm⁻² less hydrogen would be needed toward that direction. This would mean that we falsely applied a higher absorption correction, systematically increasing the unabsorbed luminosities of the clusters lying in that part of the sky. For these 42 clusters included in that cone, the average total hydrogen column density is ∼6 × 10²⁰ cm⁻² (the actual individual values vary significantly within 0.9−20.6 × 10²⁰ cm⁻²). Therefore, ∼60% less absorbing material should exist toward this direction, in order for the A value to match the rest of the sky. This seems considerably unlikely. In terms of varying metallicity, the hydrogen cloud that was detected there should be metal-poor (Z ∼ 0.4 Z_⊙) to explain the obtained discrepancy.

This seems rather unlikely since it has been shown that toward the central bulges of spiral galaxies, and the Milky Way specifically, the metallicity of the ISM is expected to be higher (Boissier & Prantzos 1999; Schönrich & McMillan 2017; Spina et al. 2017 than the metallicity in regions further away from the Galactic center. While these studies focus more on the Galactic plane and we do not use any clusters within 20° from the latter, they indicate that the high-A regions are expected to be more metal-rich than the low-A regions, instead of ∼3−4 times more metal-poor (which could potentially explain the anisotropies). Since X-ray absorption models do not account for these effects that could potentially bias the extracted cluster properties, further testing will be needed in the future. Following the same reasoning for the θ = 60° cones results in quantitatively very similar results¹⁵.

Another possible explanation for the behavior of the low-A regions would be the existence of nearby dwarf galaxies (e.g., McConnachie 2012), that contain sufficient amounts of X-ray absorbing material to cause such dimming to the clusters, for which we do not account for. However, these systems would need to fulfill some conditions such as having a large apparent size in the sky and containing absorbing material not detected by LAB. For the latter to happen, the absorbing material would need to either have a line-of-sight velocity outside of the LAB range or its hydrogen content to be limited compared to the existing metals (as explained before). Even though this “hidden” absorption by nearby galactic systems would still not explain the behavior of the bright, high-A regions, the statistical significance of the latter would drop since clusters from other parts of the sky would see an increase in their L_X.

As an overview of this analysis, we see that such large differences between the detected and the true amount of N_Htot (if this is the only reason behind the apparent anisotropies), are relatively difficult to occur, but definitely worth further checking. The necessary metallicities of the ISM to explain the behavior of the anisotropic regions seem quite unlikely as well, since one would expect oversolar metallicities close to the Galactic center, and not undersolar ones.

Free to vary N_H results from literature. A direct way to check if any of the above cases seems possible to explain our results is to try to estimate the absorption using only the X-ray spectra independently of the N_Htot measurements that were used above. This can be done by leaving the N_Htot parameter free to vary when fitting the cluster spectra. Comparing these estimations to the W13 N_Htot values, one can see if there is a systematic difference for sky regions that show extreme A behavior, indicating lower or higher absorption than the one previously adopted. These potential differences can reflect either differences in the actual amount of the ISM material as calculated before or differences of the true metallicity of that material, compared to the universally-assumed Solar one. Of course if different instruments are used for this estimation (e.g., Chandra and XMM-Newton), calibration issues must be taken into account.

In our case, the spectral fitting was performed only with a fixed N_Htot as described in previous sections. The results for a varying N_Htot will be presented in future work. For now, we use the N_Htot measurements as obtained by Lovisari & Reiprich (2019) who fit the X-ray spectra of 207 nearby galaxy groups and clusters, determining their metallicity radial profiles using only XMM-Newton observations. We assume that the determination of N_Htot does not strongly depend on the physical properties of the fitted cluster spectra, and thus only the sky coordinates of each object is of interest for our test. There are 142 overlapping clusters between these 207 clusters and our 313 clusters.

We plot the difference D = N_H, free − N_Htot as a function of N_Htot for two regions: a cone with 45° radius centered at (l, b) ∼ (273°, −19°) and the same cone centered at (l, b) ∼ (26°, + 9°). The selected coordinates are the average values of all the results from the analysis up to now. The results are illustrated in Fig. 18.

As one can see, there are no significant differences of D between the two extreme regions¹⁶. Moreover, the existing deviations between the free absorption and the previously assumed one are not large enough to explain the apparent anisotropies. For the low-A region, the median value of D is D_med = −0.28 × 10²⁰ cm⁻², while it would need to be D ≳ +3   ×   10²⁰ cm⁻² to alleviate the existing statistical tension. For the high-A region, the result is D_med = −0.22 × 10²⁰ cm⁻², while the anisotropies could be explained if D ≲ −4 × 10²⁰ cm².

All in all, using the X-ray cluster spectra of Lovisari & Reiprich (2019) as a first indication we see that the true, total absorption, does not seem to significantly deviate from the adopted absorption from W13. If a significant deviation between true and “currently-measured” absorption is detected in the future, it could potentially explain the observed anisotropies of the L_X–T relation in the sky. For now though, such an explanation seems unlikely.

6.3.3. Absorption from the Magellanic system

The Magellanic system is comprised of the Large Magellanic Cloud (LMC) galaxy, the Small Magellanic Cloud (SMC) galaxy, the Magellanic Stream (MS), the Magellanic Bridge (MB) and finally the Leading Arm (LA). All these objects are known to contain sufficient amounts of neutral hydrogen which could potentially interfere with our measurements if not taken into account. However, the very vast majority of the hydrogen of the Magellanic system is well within the velocity range for which the N_Htot values are extracted (velocity range of LAB). Thus, it should already be included in the results of the Sect. 6.3 and taken into account during the correction of the L_X–T values for the absorption.

In addition, it has been shown (e.g., D’Onghia & Fox 2016; Choudhury et al. 2016, and references therein) that the Magellanic system is metal-poorer (∼10%−50%) than the Solar metallicity assumed in the absorption correction models. As a result, the latter will overestimate the absorption effects caused by the Magellanic system¹⁷, eventually overestimating the unabsorbed L_X values of the clusters in these regions. Combining with the fact that the Magellanic system mostly covers sky regions in which clusters appear to be fainter than expected (low A anisotropies), one sees that is unlikely that the Magellanic system has any effects on our anisotropic results. Nevertheless, we examine all the above-mentioned components to see where they lie in the sky and if they correlate with the anisotropic behavior we observe.

The LMC is located at (l, b) ∼ (281°, −33°), within the low normalization regions, and moving away from us with +262 km s⁻¹ (McConnachie 2012). This velocity implies that the LAB survey would have not detected only the neutral hydrogen with a peculiar velocity of ≥ + 140 km s⁻¹ compared to the LMC center, and toward our line of sight. Moreover, its N_HI distribution is peaked close to the stellar population, covering a “circle” with a ∼3° −4° radius in the sky, centered at the above coordinates D’Onghia & Fox (2016; thereafter D16). From the 313 galaxy clusters we use, only two are within 15° from LMC, but only one is dimmer than expected based on its temperature (however within the intrinsic scatter limits). The N_Htot value toward the LMC as given by W13 is 3 × 10²¹ cm⁻². Therefore, based on all the above, it is safe to conclude that the LMC system does not bias our analysis since we would need multiple systems to be affected by that and appear underluminous.

The SMC is located at (l, b) ∼ (303°, −44°) (where N_Htot = 3 × 10²¹ cm⁻²), further away than LMC but still moderately close the low normalization regions. Its line of sight velocity compared to us is 145.6 km s⁻¹, therefore the greatest parts of its hydrogen components are expected to have been accounted for from the LAB survey. Its angular size is ∼50% of LMC and the N_HI distribution still seems to be mostly concentrated within its optical counterpart (D16). From our sample, only two clusters are within 10° of SMC and five are within 15°. From these five clusters, three have minimal random residuals from the best-fit L_X–T relation for the whole sample, one is up-scattered and the last one is low-scattered. Thus, once again we can safely assume (mainly because of its low relative velocity and the normal L_X–T behavior of the few clusters) that SMC does not cause any significant bias to our results.

The MS extends over 100° on the sky, starting from LMC and spreading toward the south Galactic pole. Then it moves up to Galactic latitudes of b ∼ −40°, for l ∼ 100° (Fig. 1 in D16) covering 2700 deg² in total. However, its N_HI density linearly decreases more than 20°away from LMC. Generally, MS is not expected to affect our results for two reasons. Firstly, the fraction of the sky it covers does not seem to correlate with low L_X–T normalization regions, since it lies at fairly low Galactic latitudes and to positions where high normalization regions are also located. Secondly, according to D16, its velocity range varies within −450 km s⁻¹ to +180 km s⁻¹, so practically all of its hydrogen component is expected to be accounted for in the LAB survey.

The MB connects the LMC and SMC systems and unlike MS, it contains a stellar component as well. Its central coordinates are (l, b) ∼ (294°, −37°) and there are only three clusters within 15° of these coordinates. However, all these three clusters are already included in the 10° circles of LMC and SMC we considered before. The main velocity of MB with respect to us is ∼225 km s⁻¹, therefore it should be included in the LAB results.

Finally, the LA extends in the opposite direction compared to the MS, starting from LMC and extending up to the northern Galactic hemisphere (for l ∼ 250° −280°). Its angular size is ∼60° (D16). It is the least massive component of the Magellanic System while its velocity range is relatively constant, ∼180−270 km s⁻¹. The fact that a large part of the LA lies within 20° of the Galactic plane, where no galaxy clusters exist in our sample, combined with its velocity range, indicates that no bias can occur for our results.

6.4. Systematics, selection effects, and correlation of results with cluster properties

Some cluster properties are usually associated with potential systematic effects. For instance, one might expect that clusters with a lower RASS exposure time might be generally up-scattered and vice versa. This is due to the fact that brighter clusters are more likely to be detected than fainter ones for the same exposure time. On the other hand, when the RASS exposure time is large enough, fainter clusters should also meet the detection thresholds set by the parent catalogs. If such a systematic indeed exists it might translate to higher and lower A values for the two cases respectively, creating artificial anisotropies. As shown in the next sections, this has no impact on our results.

Similar systematics might occur near the X-ray flux limit of our sample. For a similar temperature, intrinsically brighter clusters are more likely than other clusters to overpass the flux limit and be included in the final sample (Malmquist bias). Thus, if this applies to our sample and an excess of such clusters exist within a sky region, this will possibly result to higher A. As again shown in the following sections, our sample and analysis do not suffer of such effects.

A third possible systematic is the detection of clusters in high N_Htot regions since the difference between N_Htot and N_HI can have an impact on the selection of every X-ray flux-limited sample (including the parent catalogs). These selections are based on the unabsorbed flux, corrected only for N_HI. This flux will be underestimated for clusters lying in high N_Htot and therefore they might not overcome the flux threshold set by each sample. A ≳10% underestimation is expected for regions with N_Htot − N_HI >  4 × 10²⁰ cm⁻². The Galactic plane is usually excluded from such cluster selection processes. Nevertheless, there are still many sky regions with high N_Htot within which clusters could be missed. This could cause unaccounted selection biases and affect the completeness of the samples (which is not important to our study as explained before).

On the other hand, the high N_Htot clusters that overcome the flux limit of a sample might be intrinsically brighter in average (in order to be detected even though their estimated flux is biased low). This is “revealed” only when we correct their L_X for the molecular absorption. As a result, these clusters might be upscattered in the L_X–T plane and potentially jeopardize our results. In Sect. 6.3 we showed however that there is no such bias in our analysis. In the next sections we provide further evidence for this.

6.4.1. Correlations between A and subsample average parameters

As a generalization of the above, the apparent anisotropies in the behavior of the L_X–T relation could be in principle caused by different cluster subpopulations in the different regions. In practice, this would mean that these subpopulations might have different (average) physical properties, leading to the derived directional behavior of A. To investigate this, we perform a bootstrap resampling analysis. We drew 10⁵ random subsamples of 65 clusters (typical number for the θ = 60° cones) independently of the direction. For every subsample, we find the best-fit A and the weighted mean of the temperature, redshift, core and outer metallicity, flux, luminosity, N_Htot, intrinsic scatter and RASS exposure time. Thus, we can study if any correlation between the average values and A exists. Such a correlation, combined with a different parameter distribution in the most extreme regions could (at least partially) explain the observed anisotropies. Additionally, we created another 10⁵ subsamples with a random number of clusters (between 35 and 170) in order to test if a correlation between the number of data and A exists.

In order to check for any possible correlations, we plot the A value against all these parameters parameters for every one of the 10⁵ subsamples. The most characteristic of these plots are displayed in Fig. 19, while the rest can be found in (Fig. A.1). We also calculate the Pearson’s correlation coefficient r_corr given by Eq. (9).

$$\begin{array}{c}\hfill r_{\mathrm{corr}}=\frac{n\sum x_i{y}_i-\sum x_i\sum y_i}{\sqrt{n\sum x_i^2-{(\sum x_i)}^2}\sqrt{n\sum y_i^2-{(\sum y_i)}^2}},\end{array}$$(9)

where x and y are the two quantities for which we wish to study their correlation, and n = 10⁵ is the number of subsamples used. Depending on the obtained value of r_corr one can assess if there is indeed some correlation between the value of A and some average property of the different cluster subsamples.

The only physical parameters that seem to mildly correlate with the behavior of the L_X–T normalization A are the average subsample luminosity (r_xy = +0.239 ± 0.006) and the average 0.2−0.5 × R₅₀₀ metallicity (Z_out, (r_corr = −0.248 ± 0.005)). As shown, the best-fit A tends to slightly increase for clusters with a higher average L_X, which is expected. Subsamples with randomly more up-scattered clusters in the L_X–T plane will naturally return a higher A. Hence, pure randomness can produce similar small anisotropies, which are highly unlikely to explain the observed spatial L_X–T anisotropies (see additionally Sect. 9.2).

If randomly up- and down-scattered L_X values were not the reason behind the weak A − L_X correlation, then a similar (positive) correlation would exist between A and average T. However, there is no correlation between these parameters as shown in Fig. A.1 (r_corr = +0.005 ± 0.006). Here we should also note that the weighted T average of the highest and lowest A regions is similar, shifting between ∼5.1 keV and ∼5.7 keV for both of them, depending on the cone radius. Hence, we can safely conclude that the observed anisotropies do not arise due to any different T distribution.

On the other hand, a slight decrease in A can be seen with an increasing average Z_out. One can clearly conclude though that the obtained anisotropies of A cannot be attributed to this mild correlation since the highest and lowest A regions have similar average Z_out values (Z_out ∼ 0.40 Z_⊙).

A rather weak correlation is observed between A and average redshift (r_xy = +0.198 ± 0.007), although the brightest and faintest regions have again a very similar average z (∼0.09−0.1 for both). This indicates that the exact choice of the redshift evolution parametrization in the L_X–T relation (e.g., E(z)⁻¹) is not particularly important, since any parametrization would approximately have the same effect in the two most anisotropic regions. Indeed, if one tries different x priors for the E(z)^x term, the significance of the final anisotropies fluctuates only by ±0.1σ compared to the used self-similar case (Appendix A.1).

Another weak correlation of A is observed with σ_int (r_corr = −0.194   ±   0.006). However, this trend cannot explain the apparent anisotropies since σ_int does not strongly differ for the highest and lowest A regions (0.10 dex and 0.13 dex¹⁸ respectively for the θ = 60° cones). No correlation is observed between the A value and N_Htot, no. of clusters, RASS exposure time, Z_core and flux.

As an overview, no correlation of A with an average parameter (including systematics) can explain the apparent anisotropies. In the future, the correlation of A with combinations of these average parameters will be explored as a possible explanation behind the discrepancies, even if this seems unlikely based on the results up to now.

6.4.2. L_X–T fitting residuals as functions of cluster properties

As a further, secondary check, we tested the correlation between the cluster properties and their log L_X residuals compared to the overall best-fit L_X–T model¹⁹. As expected, similar results with the bootstrap analysis were obtained. Thus, to avoid repetition we do not go into a detailed presentation of all the results, but instead focus on the ones usually related to systematic biases. The full discussion is found in Appendix A.2.

In a nutshell, no strong systematic behavior of the log L_X residuals is observed for varying RASS exp. time, N_Htot and z (Fig. A.2). A mild systematic behavior of the residuals exists in terms of the flux and the statistical uncertainties (σ_stat) of the clusters. However, this has no effect in the derived anisotropies since the strongly anisotropic sky regions have similar flux and σ_stat distributions. Finally, the residuals versus the outer cluster metallicity are also displayed in the same figure since there is a mild systematic behavior between these quantities. We already showed that this does not significantly affect our results in previous sections.

6.5. Fixed slope vs free slope

In our analysis until now we fixed the slope to its best-fit value for every subcategory of clusters, before we study the spatial anisotropies of A. This choice is motivated by the fact that B does not significantly fluctuate throughout the sky for the 1D analysis (Fig. 6, similar results obtained for a 2D scanning). Moreover, a significant correlation between A and B is not expected, due to the pivot point of the L_X–T relation being close to the median T. To investigate the possible biases that a fixed B introduces to our analysis, we perform the following:

Case 1: we scan the sky using θ = 75° while we treat B as a nuisance parameter. We allow B to vary simultaneously with A, within its 2σ limits from its overall best-fit value. We then marginalize over B to study the spatial behavior of A. The 1σ uncertainties of A are again extracted based on the Δχ² ≤ 1 limits since there is only one parameter of interest.

Case 2: We repeat the procedure but this time we allow B to vary freely. We again study the A anisotropies and quantify the statistical significance using the Δχ² ≤ 2.3 limits (2 parameters of interest) for the 1σ parameter uncertainties.

For Case 1, the maximum anisotropy is found between the regions (l, b) ∼ (272°, −21°) (A = 0.977 ± 0.050) and (l, b) ∼ (26°, −13°) (A = 1.274 ± 0.062). The statistical significance of the tension is 3.74σ (27 ± 7%), slightly larger than before despite the marginalization over B. One sees that the results are entirely equivalent to the case where B is kept fixed, in terms of both statistical significance and direction. This strongly demonstrates the robustness of our method and the independence of the A constraints from B.

For Case 2, the A map is portrayed in Fig. 20 (top panel). The spatial fluctuations of A slightly intensify (∼34% between the most extreme values) and its directional pattern remains the same as when B is kept fixed. This once more illustrates that the derived A sky pattern does not depend on the true B values of the different sky regions.

As already shown for the 1D analysis, the fluctuations of B are smaller (∼19%) than the ones of A. Every sky region is consistent within < 2σ with the rest of the sky, making the behavior of B fairly consistent throughout the sky. The largest A anisotropy is found toward the same regions as for Case 1 (drifting by < 6°) and now slightly drops to 2.78σ (30 ± 11%). This is due to the enlarged uncertainties obtained from the 2-parameter Δχ² ≤ 2.3 limits. These two regions return very similar slope values (B = 2.256 ± 0.110 and B = 2.109 ± 0.119 for lowest and highest A respectively). Their L_X–T plot is displayed in the bottom panel of Fig. 20.

These results confirm that the choice of keeping the slope fixed for the bulk of our analysis does not introduce any biases in the directional behavior of A or the statistical significance of the observed anisotropies.

7. Cosmological constraints

Many reasons that could potentially lead to a biased anisotropic behavior of the L_X–T relation were tested until now. These tests explored the possibility that the apparent anisotropies could appear due to systematic differences of the subsamples in different patches of the sky or that unknown effects could influence the observed X-ray photons coming from specific extragalactic regions. Furthermore, we tested if possible systematics, such as RASS exposure time, Malmquist bias close to the flux limit etc., could bias our results. The observed anisotropies seem to be consistent and are not significantly alleviated by such tests.

During our analysis up to now, we assumed fixed cosmological parameters toward all the directions in the sky when deriving the normalization and slope of the L_X–T relation. On the contrary, one can reasonably assume that the physics within the ICM of galaxy clusters that determine the correlation between L_X and T should be the same regardless of the direction. As a result, the true normalization and slope of the L_X–T relation should not depend on the coordinates and should be fixed to their best-fit values.

Consequently, the last thing to be checked is if any apparent anisotropies could occur because of an anisotropic Hubble expansion. In practice, this would mean that the luminosity distance would differ toward varying directions for a fixed z. These differences can be expressed in terms of the cosmological parameter H₀ which enters in the luminosity distance through the conversion of the X-ray flux to luminosity.

To explain the behavior of faint L_X–T regions we would need a higher D_L for the same z (thus higher L_X). For z ≲ 0.3 it is known that $D_{\mathrm{L}}\approx {\displaystyle \frac{1}{H_0}}[z+{\displaystyle \frac{z^2}{2}}(1-q_0)+O(z^3)]$, where q₀ is the deceleration parameter. Hence, a lower H₀, implying a lower current expansion rate, would return a higher D_L for a fixed z. The same would be true for a more negative q₀, implying a higher acceleration rate²⁰.

One could also study the directional behavior of Ω_m, but for low redshift objects (like the clusters we use), D_L is not very sensitive to this parameter. Therefore, large deviations from region to region would be needed in order to explain the anisotropies (see results of M18). Moreover, Ω_m variations would have a different effect to higher and lower redshift (and thus temperature) clusters, changing both the normalization and the slope of the L_X–T relation. Other effects that have to be taken into consideration in such a case is the higher (lower) matter density of the Universe toward different directions, leading to more (less) structures. Structured environment can alter the behavior of the L_X–T relation as we have shown in M18. Based on that, a robust directional study of Ω_m is not ideal for our sample and method. On the other hand, the effect of H₀ on D_L does not depend on z and hence variations of smaller amplitude than Ω_m could result in the observed anisotropies. Also, since H₀ variations will have the same effect on the L_X of every cluster independently of its z (and thus T), the slope of the L_X–T will remain unchanged.

Therefore, we fix A = 1.114 and B = 2.102 and fit H₀ as the only free parameter (together with σ_intr) as described in Sect. 4.2. It should be noted that here we investigate the relative change of H₀ due to spatial anisotropies of the cosmic expansion, and absolute values of H₀ are arbitrary.

The H₀ map as produced using θ = 75° is portrayed in Fig. 21. As one can see the H₀ and A maps show exactly the same behavior for the reasons explained above. We do not plot the significance map in this case since it is identical to the significance map of A using θ = 75°, as shown in Fig. 9.

The apparently maximum acceleration direction is found toward (l, b) = (274°, −22°) with H₀ = 66.20 ± 1.72 km s⁻¹ Mpc⁻¹ while the most extreme opposite behavior is found at (l, b) = (17°, −9°) with H₀ = 75.17 ± 1.81 km s⁻¹ Mpc⁻¹. Their deviation from each other is at 3.59σ (13 ± 4%). The most extreme dipole is centered at (l, b) = (263°, −21°) with a 3.15σ significance. One sees that these three directions completely match the directions for the normalization analysis using the 75° radius cones, highlighting the ability of the normalization of the L_X–T relation to trace possible cosmological anisotropies. Moreover, the sigma values also match the ones from the normalization map as expected.

8. Combination with ACC and XCS-DR1

8.1. H₀ results for each sample

For the ACC and XCS-DR1 samples only the 1D analysis is presented in M18. In order to see if the behavior of the L_X–T relation for these two samples is comparable to the one of our sample in the 2D space, we repeat the analysis described in this paper using these two samples. Prior to the analysis, the 104 common clusters between our sample and ACC are excluded from the latter. Since the focus of the paper is on the cluster sample we build and use, here we only present the results for θ = 75°. Nevertheless, we cannot use more narrow cones to either sample. This is due to the fact that ACC does not have enough clusters (168 in total, after the exclusion of the 104 common clusters) for such small cones, while the spatial distribution of the XCS-DR1 clusters is not entirely uniform. This results in the number of clusters falling below 30 for many regions when 60° cones are used. The results for both samples can be seen in Fig. 22.

For ACC the highest and lowest H₀ (brightest and faintest respectively in terms of A) regions are at (l, b) = (77°, + 15°) and (l, b) = (317°, −14°) respectively. These two directions are relatively consistent with the general behavior of our sample, with a ∼40° separation compared to the results of our sample. We remind the reader that the two samples are completely independent. The H₀ values of these extreme regions are H₀ = 78.76 ± 4.15 km s⁻¹ Mpc⁻¹ and H₀ = 58.12 ± 2.68 km s⁻¹ Mpc⁻¹ deviating by 4.18σ (30 ± 7%). Their angular separation is 122°, which is similar to the ones for the extreme regions of our sample. The most extreme dipole is found toward (l, b) = (327°, −21°) with a 3.68σ. It is noteworthy that the low H₀ is much more statistically significant in this sample than the high H₀ region, indicating a monopole anisotropy. One obtains similar results for ACC when B is left free to vary within its 2σ limits as a nuisance parameter, but with a decreased statistical significance. In that case, the statistical significance of the anisotropy slightly drops to 3.12σ (from 4.18σ for a fixed B) toward similar sky directions.

For XCS-DR1 the most extreme regions are located at (l, b) = (31°, + 25°) (brightest) and (l, b) = (281°, + 24°) (faintest) separated by 117°. Their respective H₀ values are H₀ = 77.91 ± 2.20 km s⁻¹ Mpc⁻¹ and H₀ = 63.56 ± 2.32 km s⁻¹ Mpc⁻¹ deviating by 4.52σ (21 ± 5%). One can see that this discrepancy is larger than the one in our sample or the one obtained from ACC. However, XCS-DR1 has some properties that might lead to overestimating the anisotropies between different sky region. For instance, overluminous clusters tend to have smaller statistical uncertainties (M18), and when these clusters are in the center of the cones (higher statistical weight), this can lead to artificially high H₀ (or A). Thus, one has to be conservative when interprenting the statistical significance of the anisotropies found in the XCS-DR1 sample.

Interestingly, the direction for the lowest H₀ (which corresponds to the maximum cosmic expansion rate) is separated only by 28° from the CMB dipole. Since XCS-DR1 is a high redshift sample (median z ∼ 0.35), naively one would not expect any effects on the XCS-DR1 results due to the peculiar velocity of the Solar System compared to the CMB frame and therefore there is no obvious reason why these two directions should be close. The most extreme dipole for XCS-DR1 is located toward (l, b) = (211°, + 14°) with a 2.75σ significance.

Now we allow B to vary within its best-fit 2σ limits. Some changes are observed, although the general directional behavior of H₀ remains relatively consistent. The statistical significance of the maximum anisotropies significantly decreases from 4.52σ to 2.82σ. This is due to the fact that the median T = 2.7 keV of the XCS-DR1 sample is smaller than the pivot point (4 keV) of the L_X–T relation, and thus A and B values are more correlated than for our sample (or ACC). These small differences between the results of the two cases can be avoided if one chooses the pivot point to be ∼2.7 keV for the XCS-DR1 modeling. Despite of these small alternations, the most extreme region is still found toward (l, b) ∼ (292°, + 23°), only ∼10° away from the previously found direction, and still with a ∼3σ significance.

8.2. Combining the H₀ results for the three samples

Remarkably, ACC and XCS-DR1 roughly agree with our sample on their L_X–T anisotropic behavior despite the fact that they do not share any common clusters. While at first sight it might seem that the H₀ maps of ACC and XCS-DR1 look different, the location of their most extreme regions is still consistent within ∼40° −55°.

In total, they contain 842 different galaxy clusters. Consequently, any constraints on the fitted parameters would be much stronger if we combined them. While, the normalization values of the three samples are quite different (cluster populations, used energy range for L_X, T constrain method etc. vary significantly), H₀ is a global parameter that should not depend on specific samples or even cosmological probes. The normalization and slope values of the three different samples can be set in such way so the best-fit H₀ value considering the entire sample is H₀ = 70 km s⁻¹ Mpc⁻¹. Nevertheless, we see that the three samples return a different H₀ range. As shown before ACC and XCS-DR1 show a larger variation of H₀ (± ∼ 20%) than our sample (± ∼ 9%). This correlates with the larger scatter of the other two samples and it can be attributed to randomness (since the H₀ uncertainties of ACC and XCS-DR1 are ∼2−3 times larger than the ones of our sample), reasons that we have not yet identified or a combination of the above (the significance however remains similar for the three samples).

By performing the H₀ scanning analysis, one obtains three different and independent estimations of the likelihood of the H₀ parameter for every region. Multiplying these three likelihoods gives us the combined most likely H₀ value for every region in the sky. In order to consistently use the three samples, we use the smallest possible cone radius (75°) for which we have enough data for all three catalogs in any cone, and we use the same parameter fitting range (H₀ ∈ [50, 90] km s⁻¹ Mpc⁻¹) as well. Therefore, the H₀ map displayed in the left panel of Fig. 23 is obtained, while the significance map is shown in the bottom panel of the same figure (we also overplot the results of other studies, as discussed in Sect. 9.1 and Table 3).

From the combined H₀ results, the lowest value H₀ = 65.20 ± 1.48 km s⁻¹ Mpc⁻¹ occurs toward (l, b) = (303°, −27°) (237 clusters) while the highest value H₀ = 76.64 ± 1.41 km s⁻¹ Mpc⁻¹ is found at (l, b) = (34°, + 26°) (302 clusters). Therefore, the null isotropy hypothesis between these two regions is rejected with a remarkable significance of 5.59σ (16 ± 3%). The angular separation of these two regions is 103°. On the other hand, the strongest dipole occurs toward (l, b) = (265°, −20°) (57° away from the CMB dipole) with a significance of 4.06σ.

We repeat the joint analysis considering the obtained H₀ results from every sample when B was left free to vary as a nuisance parameter. As expected, the overall behavior of H₀ persists with some limited changes. The statistical significance of the maximum anisotropy drops to 4.55σ (from 5.59σ), and is found between (l, b) ∼ (312°, −21°) and (l, b) ∼ (45°, + 21°). Consequently, the choice of keeping B fixed slightly overestimates the exact statistical significance of our findings but does not affect the general conclusion.

All these results demonstrate clearly that the similar anisotropies in all three independent samples are extremely unlikely to be random and that there is an underlying reason causing the L_X–T relation to show a strong directionally depended behavior.

9. Discussion

The significance of cosmic isotropy for the standard cosmological paradigm is undisputed. Designing scrutinizing methods to test this hypothesis is vital since much new information about the Universe can be revealed through such tests.

One can assume that the isotropic expansion of the Universe holds, but a cosmological probe could still consistently show a significantly anisotropic behavior. This could result in the identification of yet unknown factors with a surprisingly strong impact on the data collection, analysis, or both. Since these factors are not accounted for in previous studies using similar wavelengths (e.g., X-rays) or the same astrophysical objects, these biases could in principle extrapolate to many aspects of relative research fields.

For instance, the anisotropy of the L_X–T scaling relation found in this paper could have multiple implications for other studies using X-ray galaxy clusters or other X-ray objects. Since the strong anisotropies do not strongly depend on the specific sample, X-ray satellite etc. the vastly more probable scenario is that the underlying reason is not a sample-specific systematic.

Moreover, the amplitude and direction of the anisotropies are preserved even after excluding several cluster subcategories, such as low-T systems, local clusters, clusters with high absorption, metal-rich clusters, high flux ones etc. Thus, different subpopulations of clusters toward different directions is not a likely explanation as well. Additionally, possible biases due to selection effects do not explain the findings. Therefore, if this is eventually proven to be caused by an unknown (extra)Galactic effect acting on X-ray photons, previously published results would need modifications correcting for this effect.

Such an example would be the galaxy cluster masses obtained through the X-ray luminosity-mass scaling relation L_X − M. If we assume a typical scaling relation slope of L_X ∼ M^1.5, then the masses of the clusters toward the faint regions of our analysis would be underestimated by ∼10−20% while the clusters in the bright regions would end up with masses overestimated by the same amount. As a result, the cosmological parameters obtained via the halo mass function could be biased if the sky coverage of the clusters is not uniform. Even in the latter case, the scatter of the final results would increase. It is characteristic that without applying any statistical weighting in the clusters and performing the sky scanning using θ = 60°, ∼72% of the subsamples in the different directions show a lower σ_int than the full sample results. This indicates the potential increase of the scatter in X-ray scaling relations when the full sky is used as if galaxy clusters were showing the same behavior everywhere. Possibly biased results when the used samples do not cover the full sky homogeneously can clearly occur to any other studies as well, if these use measured X-ray luminosities (or temperatures) of galaxy clusters.

Another useful test would be to study the dependance of these anisotropies on the exact energy range. This will be particularly helpful in order to check if the observed anisotropies could be the result of absorption effects, such as strong variations in the galactic ISM metallicity, metal-rich nearby dwarf galaxies etc. However, in Sect. 6.3.2 we showed that this is unlikely, but further testing is needed. Nevertheless, checking if these L_X–T anisotropies also appear in the hard X-ray band alone, where the absorbing effects are minimal, would provide us with valuable information about their exact nature. This will be feasible with the upcoming eROSITA all-sky survey. Here we should remind the reader that while we only use the 0.1−2.4 keV energy range for L_X, the ACC and XCS-DR1 samples (which also show similar anisotropies), use the bolometric energy range.

In the L_X measurements used for this study the cluster cores are not excluded, since this is very difficult to do with ROSAT data due to its large PSF. It has been shown however that core-excised luminosities scatter less in their scaling with temperature (e.g., Markevitch 1998; Pratt et al. 2009; Maughan et al. 2012, etc.). Such values would be optimal for our analysis since a lower scatter in the L_X–T relation would decrease the uncertainties of the derived A values. This could eventually allow the detection of spatial anisotropies with an even higher statistical significance and strengthen our results. This will be possible with eROSITA data and with possible future XMM-Newton and Chandra-based samples that provide core-excised luminosity values.

The summary of the best-fit A, B and σ_int values is shown in Table 1. The directions of the most extreme regions for every subsample, together with the statistical significance of the anisotropic signal between these two regions and the direction and significance of the most anisotropic dipole are shown in Table 2.

The consistent value of the slope throughout the different subsamples is noteworthy. The largest difference (∼1σ) is found between the Z_out ≤ 0.426 Z_⊙ and the N_Htot ≤ 7.37 × 10²⁰ cm⁻² subsamples. On the other hand, A deviates by ∼2σ between the T >  3 keV, z >  0.05 and the N_Htot ≤ 7.37 × 10²⁰ cm⁻² subsamples, while it is quite consistent between the rest. As expected, the lowest scatter is found for the subsamples with the highest T and z. On the contrary, the subsample where the high Z_core clusters were excluded returns the largest scatter, still consistent though with the other subsamples. We should also note here the significantly lower total scatter of our sample against XCS-DR1 and ACC (after converted to L_X–T scatter).

Generally, as θ decreases, the statistical significance of the results increases, as we are able to pinpoint the anisotropies more effectively. However, the amount of available data is not yet enough to use even narrower angles. This will change with future surveys such as the upcoming all-sky eROSITA survey which will provide us with a larger number of observed clusters with temperature measurements (Borm et al. 2014). Nevertheless, the existence and consistency of these apparent L_X–T anisotropies are already on solid ground granting these results, especially when one combines all three independent samples with a > 5σ anisotropy emerging from this. This holds even when the slope is left free to vary (within a limited range) from region to region, and then marginalized over. It is also quite interesting that the maximum anisotropic directions in almost every tested case seem to prefer an angular separation of ∼80° −120° instead of a dipole form. The most extreme observed dipole anisotropies have a statistical significance of ∼4σ, with an angular distance of ∼50° −100° from the CMB dipole direction. At the same time, the faintest parts of the maps are slightly closer (∼35° −90°) to the corresponding end of the CMB dipole.

Here we should discuss some possible reasons for caution when one interprets the large statistical significance of the observed anisotropies. Firstly, while we have tested a large number of potential X-ray and cluster-related reasons and systematics that might cause such a spatially inconsistent behavior, we only tested them one by one. If one takes into account two or more such reasons simultaneously the statistical tension might decrease. Although it seems improbable that the observed anisotropies can be attributed purely to such effects (since three independent samples show similar behavior), one cannot discard the possibility of an overestimation of the anisotropies due to the (unchecked) combination of systematics. Secondly, the derived statistical significance of the results is based on the Δχ² limits of the fit. While the applied bootstrap method returns similar results (see Sect. 9.2), one still has to consider the so-called cosmic variance. To do so, one can use Monte Carlo simulations to draw similar samples from an inputted isotropic universe and, following the same method as in this paper, check how often such large anisotropies appear. This will be done in future work.

9.1. Comparison with other studies

Except for identifying previously unknown factors that can significantly affect the determination of physical parameters of astrophysical objects as discussed above, testing the isotropy of the Universe has of course another aspect as well. If many independent cosmological probes agree on a similar anisotropic direction and amplitude, while all known biases have been accounted and corrected for, then the hypothesis of cosmic isotropy should be reconsidered. This could eventually lead to a major shift in the standard cosmological model.

The direction we identify as the one with the maximum acceleration (or minimum expansion rate as explained before) if the anisotropies were indeed only of cosmological origin, agrees well with many other studies that used SNIa and other probes to look for possible anisotropies in the Hubble expansion. Several examples of such studies are shown in Table 3, together with their the most anisotropic directions and their significance.

Generally, it is usual that the anisotropies found in SNIa come mostly from z ≲ 0.1 and they are attributed to local bulk flows, arising due to the Shapley supercluster at (l, b) ∼ (306°, + 30°) with z ∼ 0.04−0.05. We should note however that the anisotropic results of SNIa strongly depend on the used sample since studies that have been performed with the latest SNIa compilations tend to find consistency with isotropy as discussed in Sect. 1. Moreover, the rather inhomogeneous SNIa coverage of the sky can create problems in the search of a preferred cosmological axis.

Within an isotropic FLRW background the directions of peculiar velocities are expected to be randomly distributed. However, a coherent bulk flow toward a massive structure due to gravitational attraction, it would affect the redshifts of local objects in a systematic way. If not taken into account, the luminosity distance (calculated through z) of clusters would be over or underestimated depending on their position in the sky. This could inevitably lead to apparent anisotropies arising from local probes.

Although these local flow motions are not expected to extent beyond ∼200 h⁻¹ Mpc, the studies shown in Table 3 (among others) detect bulk flows (or anisotropies) further away than this scale and with amplitudes which are hard to explain within ΛCDM. This detection is performed by different independent probes. The statistical significance however decreases compared to local probes due to the limited number of data in certain sky patches. An example of studying the scale of bulk flows is given in Carrick et al. (2015) who find a 5σ bulk flow of ∼160 km s⁻¹ extending over 200 h⁻¹ Mpc toward (l, b) ∼ (304°, + 6°). The structures that could fully explain such a bulk flow motion have not been identified yet. Moreover, the direction of the anisotropies of more distant probes tends to converge with the one from the CMB dipole, but often with a slightly larger amplitude.

The consistency of the apparent anisotropies beyond ∼210 h⁻¹ Mpc (z >  0.05) can be also seen in our results, where the tension with the null hypothesis of isotropy does not decrease. Another effective test could be to perform our L_X–T anisotropy analysis with clusters at z >  0.2, beyond the effects of the recent large-scale bulk flow detections. Currently there are not enough data for such a test though, but this is expected to change with the upcoming all-sky eROSITA survey. Finally, if the only reason behind the anisotropies we observe in the L_X–T behavior was local or cosmic coherent flow motions, one would expect to retrieve mostly dipole anisotropies, whether we have shown that anisotropies separated by ∼90° −120° are more significant in our analysis. However, a more in-depth testing is needed to draw safe conclusions about this scenario.

9.2. Statistical significance validation by bootstrapping

In order to further investigate the statistical significance of our results and if they could be attributed to pure chance we perform a bootstrap resampling analysis. We consider two cases:

In the first case, we used all the 313 clusters covering the whole sky. We drew 10⁵ random subsamples of the same size as the region we want to test its significance. We assigned random statistical weights in the drawn clusters²¹ to simulate the method we use during the sky scanning analysis. There, the weights were assigned based on the distance of every cluster from the center of the scanning cone. This test demonstrates how often our cluster sample can reproduce such low or high A values randomly and independently of the direction, when having the same number of clusters as in the extreme regions.

In the second case, we excluded the subsample of interest and performed the 10⁵ resamplings based on the rest of the clusters. This way, we can estimate how many times the extreme result of the excluded subsample can occur randomly from data in other directions. The random statistical weighting is used here as well.

Both cases also offer a direct comparison with the deviations occurring from the Δχ² limits, from which the reported statistical significance for every result comes from. In order to have minimal overlapping between the 10⁵ realizations, we choose to perform this analysis for the results occurring for θ = 60°. The number of clusters in the extreme regions is small enough so there is no significant overlapping, while it is large enough to be relatively insensitive to strong outliers.

Drawing and analyzing 10⁵ subsamples of 84 clusters from the full sample, we find that only 0.68% of the results have a lower A than the one found for the (l, b) = (281°, −16°) direction (A = 0.940 ± 0.051). This corresponds to a p−value of p = 0.007 for the null hypothesis, or in a Gaussian significance of 2.71σ. Now we repeat the analysis with a subsample size of 78 clusters, same as the brightest region for a θ = 60° cone toward (l, b) = (34°, + 4°). We find that 10% (p = 0.010, 1.65σ) of the results have a higher A ≤ 1.346 compared to the aforementioned bright region. Therefore, the statistical significance of the fainter region toward (l, b) = (281°, −16°) is much higher in that case. The statistical deviation of these two extreme regions based on the based on the Δχ² limits as shown in Eq. (5) is 4.73σ (Table 2). As found from the bootstrap resampling method however is ∼3.9σ, slightly decreased but still significant. Finally, the most probable value for these 10⁵ realizations is A ∼ 1.118 ± 0.115, which is consistent with the results of the full sample fitting.

For the second case, we first excluded the 84 clusters within 60° from (l, b) = (281°, −16°) and we only considered the rest 229 clusters. Following the same procedure as before, we find that only 0.32% of the subsamples have A ≤ 0.940 (p = 0.003, 2.95σ), the same result as the one we obtained from Δχ² limits in Sect. 6.3.2. Doing the same for the 78 clusters (full sample except these 78) toward (l, b) = (34°, + 4°) we see that an A ≥ 1.346 value is reproduced only for 0.45% of the 100 000 subsamples (p = 0.005, 2.85σ), again consistent with our previous findings. In this case where only clusters from the rest of the sky are considered, the probability of the high A result to occur randomly drops significantly compared to the case where the full sample is used. This indicates that these 78 clusters strongly affect the bootstrap results when all 313 clusters are used.

Repeating the analysis for the results when θ = 45°, first we used the full sample with all the 313 clusters. For the results described in Sect. 5.3, we obtain a probability of 0.45% for the lowest A value to occur randomly from the whole sample, and a probability of 8.71% for the highest A value. The deviation between these two results is ∼3.9σ. This is displayed in the left panel of Fig. 24.

Excluding these two extreme subsamples one at a time as we did in the second case before, the probability decreases to 0.09% (p = 9  ×  10⁻⁴, 3.33σ) for the faintest region and to 0.67% for the brightest region. The results are once again consistent with the deviations obtained from the Δχ² in Sect. 6.3.2.

If one excludes both subsamples simultaneously and only considers the 234 left, it results in a probability of 0.21% and 0.69% of the null hypothesis to reproduce the A results of the faintest and brightest regions respectively. The result can be seen in the right panel of Fig. 24. The deviation of the two extreme results is at ∼5.4σ, slightly higher from the 5.08σ value predicted from the Δχ² limits.

Finally, as discussed in Sect. 4.3, in order to ensure that the artificially low values of σ_int do not affect our results, we repeated the bootstrap analysis without including the σ_int term in our model (Eq. (4)). For that we also found the A value of the extreme regions without accounting for σ_int. The significance of the bootstrap results remains high, decreasing only by ∼5−12% for every case (e.g., from ∼3.9σ to ∼3.7σ for the θ = 60^θ case and from ∼5.4σ to ∼4.8σ for the θ = 45°) since outliers are now more likely to cause extreme A behaviors.

This analysis strongly demonstrates the high statistical significance of the results. Also, it is shown that the apparent anisotropies are very unlikely to be attributed to randomness, as well as that the statistical deviations obtained through Δχ² limits match the ones from bootstrapping.

10. Conclusions

In this work, we constructed and analyzed a new, large homogeneously selected X-ray galaxy cluster sample of 313 objects, with the purpose of probing the anisotropic behavior of the L_X–T scaling relation as first found in Migkas & Reiprich (2018; M18). Through the strong correlation between the X-ray luminosity and temperature and the null hypothesis that the L_X–T behavior must be similar throughout the sky, one can probe the existence of up-to-now unknown factors affecting the behavior of X-ray photons, galaxy clusters or both, for different sky directions. Furthermore, one can estimate how isotropic the Hubble expansion seems to be by constraining the cosmological parameters for different sky patches. This can be done due to the inclusion of the cosmological parameters in the X-ray flux-luminosity conversion, where we take advantage of the fact that the determination of the temperature is cosmology-independent. A necessary requirement however is to verify that no underlying unknown systematics exist, affecting the X-ray observations and the galaxy cluster scaling relations in particular.

We tested the consistency of the L_X–T relation for different directions by scanning the full sky using cones of different sizes, and quantify deviations in terms of the normalization parameter A, or the Hubble constant H₀. A consistent and strong directional behavior of these parameters emerged. Dividing the sky into hemispheres, we first found that the hemisphere with its pole located at (l, b) = (272°, −8°) seems to be fainter (lower A or lower H₀) compared to the opposite hemisphere at a 2.58σ level. With our cluster sample having a quite uniform spatial distribution we could pinpoint apparent anisotropies more effectively with narrower cones. Using cones with 75° down to 45° radius we found that the sky region toward (l, b) ∼ (277°, −11°) systematically returns a lower A/H₀ compared to the sky region toward (l, b) ∼ (32°, + 15°) with a significance of ∼3.6−5σ (99.97−99.9999%). The main bulk of the deviations though come from the faint region rather than being balanced between the two extreme regions.

Surprisingly, the maximum dipole form anisotropies are systematically weaker by ∼0.4−0.9σ compared to these ∼110° anisotropies, although still significant. Moreover, the region close to (l, b) ∼ (170°, + 15°) is also systematically brighter with values comparable to the (l, b) ∼ (25°, + 4°) region but with lower significance due to fewer clusters.

We examined multiple reasons, mostly related to galaxy cluster physics, X-ray analysis and systematic biases, that could provide us with an explanation about the derived anisotropies. For instance, the L_X–T behavior of different cluster population was studied. We found that clusters in low absorption regions show the same behavior with clusters in high absorption regions after the proper corrections have been applied. Excluding the latter subcategory, the anisotropies remain. Moreover, excluding galaxy groups and clusters with T >  3 keV and redshifts of z >  0.05 do not significantly affect our results, as can be seen in Table 2. Dividing our cluster sample according to the metallicity values of our clusters (both core and outer regions) and performing the sky scanning also does not seem to explain our findings. The same is true if one allows the slope to vary within limits during the sky scanning, and then marginalizes over the slope values. We also checked if our analysis is biased by selection effects related to the RASS exposure times of the clusters, the applied flux limit and high molecular hydrogen regions, not finding any indication for such effects. However, all these tests were done one at a time. One can argue that a combination of such effects may partially decrease the high statistical significance of the anisotropies. Of course this is still a presumption since the full magnitude of the anisotropies seems unaffected by the different tests, but it is worth checking in future work.

Furthermore, we discussed the possibility of extragalactic, metal-rich systems causing X-ray absorption that is not accounted for in the LAB survey. The case where the true metal abundance in the Galaxy’s ISM shows strong spatial variations, possibly biasing the applied absorption correction and causing these anisotropies, was also discussed. Even though we showed that the last two cases are unlikely to be the reason behind the apparent L_X–T anisotropies, it is worth checking if this behavior persists also in the case where L_X is only measured in the hard X-ray band. In these photon energies the X-ray absorption is not significant and one would not expect any anisotropies caused by such effects. The eROSITA all-sky survey would be a great tool that will allow us to test that.

As a final test, we created 10⁵ random bootstrap subsamples and investigated the correlation of the average properties of their clusters with the best-fit A value. No strong correlation was found, while the most extreme regions tended to have similar average properties. This bootstrapping method we used further verified the statistical significance of our results, while it hints to the faint sky region as the most statistically unique one. In future work, simulated isotropic samples similar to the one in this work will be used to test the frequency with which such strong anisotropies appear.

Some useful by-products of our analysis have to do with the general L_X–T scaling relation behavior, such as the decrease in the scatter for higher T and z clusters, the slightly larger scatter of low core metallicity clusters compared to the rest, the ∼3σ discrepancy in the L_X–T slope when Chandra or XMM-Newton were used, the excellent agreement between X-ray and optical redshifts, as well as the strong L_X–T inconsistency between clusters with low and high metallicities within the 0.2−0.5  ×  R₅₀₀ annulus.

When our sample is combined with the ACC and XCS-DR1 samples as used in M18, we see that their sky behavior agree well with each other even without having even one common cluster among them. Moreover, the fact that the observations of the three samples come mostly from three different telescopes and the sample have been compiled by different authors and sharing different properties (such as the z distribution) should be kept in mind. Creating a full-sky H₀ map using the 842 individual clusters included in these three catalogs, a ∼5.5σ anisotropy was obtained between the sky regions toward (l, b) ∼ (303°, −27°) (H₀ ∼ 65 km s⁻¹ Mpc⁻¹) and (l, b) ∼ (34°, + 26°) (H₀  ∼  77 km s⁻¹ Mpc⁻¹). These values were obtained keeping the slope fixed. When the slope is free to vary one obtains similar results at a ∼4.5σ level. This could either mean that indeed the explanation of the anisotropies might be of cosmological origin (including strong bulk flows) or that there is a hidden (extra)Galactic factor that affects X-ray cluster measurements independently of the used sample. The direction of the anisotropies strongly correlates with results from other independent probes as shown in Table 3.

The assumption of the isotropic nature of X-ray galaxy cluster scaling relation is common, even though this had not been observationally tested and confirmed before. The possible discovery of systematics which X-ray cluster studies do not account for until now, could considerably alter the way X-ray scaling relations are used and interpreted. If this anisotropic behavior persists in other X-ray wavelengths as well, it could indicate that also other X-ray astronomy studies might need readjustments.

On the other hand, the cosmic isotropy still remains an ambiguous topic since several independent cosmological probes have been found to have an anisotropic behavior recently. While there are results not reporting any significant anisotropies, others claim to detect ∼2−3σ anisotropies either in the local Universe (z ≲ 0.1) or to larger distances. To assess this question, independent methods such as the L_X–T test are needed to be applied and their results to be compared. If no biases are identified as the reason behind the anisotropies we observe and other probes seem to consistently agree, then the explanation might indeed be of cosmological origin. Such examples would be an anisotropic dark energy nature leading to different expansion rates for different directions in the late Universe, coherent bulk flow motions up to certain cosmic scales affecting the cosmological redshift measurements etc. Irrelevantly of the actual reason, studies dealing with X-ray cluster measurements are potentially affected from our findings.

---

Acknowledgments

We would like to thank the anonymous referee for their valuable feedback and their constructive comments that helped solidify the robustness of our results. 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 grants GO4-15129X, GO5-16137X, and AR9-20013X. THR, FP and MERC acknowledge support from the German Aerospace Agency (DLR) with funds from the Ministry of Economy and Technology (BMWi) through grant 50 OR 1514. LL acknowledges support from NASA through contracts 80NSSCK0582 and 80NSSC19K0116. This research has made use of the NASA/IPAC Extragalactic Database (NED) which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration.

References

Appendix A: Extra tests

A.1. Effect of redshift evolution parametrization

During this analysis, we choose a fixed prior of E(z)⁻¹ for the redshift scaling of the L_X–T relation as shown in Eq. (2). If the true scaling is not self-similar and if two separate subsamples have a different redshift distribution, artificial normalization anisotropies might be induced. As already discussed in Sect. 6.4.1 though, the two most anisotropic sky regions (brightest and faintest) share a similar redshift distributions and are not expected to be affected by such possible biases. Also, the cluster redshifts are relatively low and thus E(z) does not rise to high values in order to significantly affect our results.

Nevertheless, we wish to test the dependance of our results on the exact E(z) prior selection. To this end, we repeat the θ = 75° analysis for the full sample, for four different cases, E(z)⁻², E(z)^−1.5, E(z)^−0.5 and for no redshift evolution (keeping B fixed to the best-fit value obtained for every case separately). The location of the most extreme regions, both faintest and brightest, fluctuate only by <9°. The statistical significance for the anisotropy between these regions maximizes for the E(z)⁻² case (3.75σ) and decreases gradually to 3.57σ for the case without any redshift evolution. The A and the significance maps for all four E(z) scenarios do not practically differ from the maps shown in the upper left panels of Figs. 8 and 9 respectively and thus we do not display them.

A.2. Correlations and systematics

As discussed in Sect. 6.4.1 we examine the possible correlations that might exist between A and the average physical parameters of every region subsample. In Fig. A.1 the rest of the correlations are shown. None of the physical parameters seems to have a significantly enough correlation with A to explain the observed anisotropies.

As also explained in Sect. 6.4.1, we look for any possible systematic behavior between the properties of our 313 clusters and their logarithmic luminosity (log L_X) residuals from the overall best-fit model. To this end, we fit the behavior of the residuals against every cluster property. Thus, we can quantify the significance of the possible deviation from the case of no systematic behavior. This method cannot be applied to the bootstrap realizations since the significance depends on the number of data points, and one creates as many realizations (data points) as one wishes.

The log L_X residuals as a function of various cluster properties are displayed in Fig. A.2.

No systematic behavior of the cluster L_X–T residuals arises for different RASS exposure times. This also holds true for the absorption correction measure N_Htot. The only noticeable feature there is toward large N_Htot values where the clusters seem to be weakly upscattered compared to the overall L_X–T best fit model. The effect of this in our observed anisotropies has been quantified in Sect. 6.3.1. There we found that if one excludes these clusters, the significance of the anisotropies actually increases, since these clusters lie in relatively low-A sky regions. The residuals seem to be consistent throughout the z range as well, with the exception of z ≲ 0.025. In Sect. 6.1 we extensively show that excluding these clusters neither alleviate the statistical tension between the oppositely anisotropic regions nor changes their sky direction.

A mild systematic behavior can be seen for high flux clusters, being upscattered in average. This is the opposite behavior than the one expected due to selection biases. In low fluxes, one can see that the residuals are randomly distributed. This limited number of upscatered high flux clusters would only affect our anisotropy results if they were not randomly distributed in the sky (which they are). Despite of that, we excluded the 37 clusters with f >  2.6  ×  10⁻¹¹ erg s⁻¹ cm⁻² (after which this systematic behavior becomes clear) and repeated the sky scanning process with θ = 75°. The maximum anisotropy actually increased from 3.64σ to 3.91σ (for this cone size), and is found between the regions (l, b) = (272°, −26°) and (l, b) = (39°, −7°).

Another mild systematic behavior is observed for the clusters with low statistical uncertainties (${\sigma }_{\mathrm{stat}}=\sqrt{{\sigma }_{logL}^2+{B^2\times {\sigma }_{logT}}^2}$, with B = 2.102), as they tend to be intrinsically brighter than average. Based on the distance-weighing method we follow during the sky scanning process, when such a low σ_stat bright cluster is close to the center of a cone, the best-fit A value of that cone can be biased high to roughly match the behavior of this particular cluster. Consequently larger anisotropies might be obtained. However, this effect is limited in this work due to the inclusion of the intrinsic scatter term σ_intr in our model.

To test this, we excluded the 39 clusters with σ_stat <  0.035 dex and repeated the analysis for θ = 75°. While the most anisotropic low-A region was found again toward (l, b) = (272°,−18°), the brightest region shifted toward (l, b) = (75°, + 22°). The statistical significance of their in-between anisotropy slightly decreased to 2.85σ (from 3.64σ). This small change is expected since for this test we discard from our sample the clusters with the best-quality measurements, marginally increasing the uncertainties of the derived A. Despite of that one sees that the significance of the anisotropies remains high. If we repeat the test for the θ = 60° cones, the maximum anisotropy found is 3.76σ (from 4.73σ initially).

Finally, the clusters with high metallicities in the 0.2−0.5 R₅₀₀ annulus appear to be systematically fainter. This result and its effects on the apparent anisotropies have been extensively discussed in Sect. 6.2.2.

A.3. Chandra-only clusters

In order to make sure that the anisotropic behavior of the L_X–T relation is not the result of a systematic bias coming between Chandra and XMM-Newton clusters (even if we calibrate the temperatures properly as described in the paper), we reproduce some of the A color maps using only the 237 Chandra clusters.

In Figs. A.3 and A.4 is shown that both the 1D and the 2D analysis yield similar results to the full sample. In the 2D map the faint regions tend to shift to lower Galactic latitude. However, the lowest A is found toward (l, b) = (294°, −34°) which is only 10° away from the combined lowest result found when all three independent samples were used. The location of the brightest regions is at (l, b) = (25°, −11°). The statistically deviation between the two most extreme regions is 2.82σ, somewhat decreased compared to the full sample (3.64σ) but not relieved.

A.4. Optical redshifts only

Finally, we check if the use of X-ray redshifts affect our results somehow. This is not expected to happen since there is an excellent agreement between the two types of redshifts as discussed in Sect. 2.1.

As shown in Fig. A.5, the A and B solution space remains identical when we considered all the 313 clusters or only the 271 clusters with optical z. Consequently, clusters with X-ray z do not bias the results and agree well with the rest.

A.5. Solution space excluding low T systems

If one excludes the low T and low z clusters, one can see that they do not strongly affect the overall L_X–T solution of the sample. This is shown in Fig. A.6 where the 3σ solution spaces are shown for the cases where we exclude all the clusters with T <  2.5 keV and z <  0.03 (left) and T <  3 keV and z <  0.05 (right), compared to the solution of the full sample.

While in the first case the two solutions are entirely consistent, in the second case they are consistent within ∼1σ. Therefore, adding to all the tests done in the main sections of the paper, we can safely conclude that these systems do not affect our anisotropic findings.

A.6. REFLEX vs NORAS L_X–T behavior

Our sample consists mainly of clusters from the REFLEX (185 clusters, 59%) and NORAS (105, 34%) clusters. A possible systematic difference between the two catalogs in their cluster population or in the flux measurements could artificially create apparent anisotropies and bias our findings. Here we should note that the two catalogs were constructed by the same team and analyzed in a similar way, so naively significant discrepancies should not be expected.

In order to test the consistency of the L_X–T behavior of the clusters coming from these two catalogs, we fit A and B for both subsamples, and compare the results. The 3σ contour plots for the two subsamples are shown in Fig. A.7.

The best-fit results are consistent for the two subsamples at a 1.4σ level. It is clear that this discrepancy is not the reason behind the observed anisotropies when the full sample is used. The NORAS clusters seem to be slightly more luminous than the REFLEX clusters, but this seems to be due to the existence of the strongest low-A anisotropic region in the REFLEX part of the sky (south ecliptic hemisphere). If one excludes the clusters within 25° from the lowest A sky direction as found for the θ = 75° cones (Table 2), then the discrepancy between the two subsamples drops to 0.8σ, which is negligible.

A.7. Systematic temperature differences between Chandra and XMM-Newton

As explained in Sect. 2.5, the Chandra and XMM-Newton telescopes show systematic differences in the temperature determination. In order to consistently use the measurements from both telescopes, one has to take this into account. To this end, we converted all the temperatures measured with XMM-Newton into “Chandra” temperatures, using the relation found in S15. To verify that this relation sufficiently describes the needed conversion for our sample as well, we measured the temperature of 15 clusters with both instruments and compare the results, which are shown in Fig. A.8.

As shown, the conversion relation found by S15 using 64 clusters is consistent with our results and thus used for the necessary temperature conversions. The statistical uncertainties of the S15 best-fit relation as well as the given scatter are taken into account in the final converted temperature values we use.

A.8. Isotropic L_X processing throughout catalogs

The L_X values have gone through several steps of processing (RASS to REFLEX/NORAS/eBCS to MCXC to our values). If the values suffered an anisotropically biased analysis during this multiprocessing, this would propagate to our results.

Firstly, we need to ensure that the L_X corrections we applied to the respective MCXC values did not introduce any artificial anisotropy. For this purpose, we check the directional behavior of the fraction between our luminosity estimated L_X,  ours and the MCXC L_X,  MCXC. This is done with the same methodology as the A scanning of the sky, for θ = 75° cones. Each cluster is assigned a statistical weight based on its distance from the center of each cone and the average L_X,  ours/L_X,  MCXC is obtained. The produced map is displayed in the upper panel of Fig. A.9. In order to directly compare with the observed anisotropies of the L_X–T, the same color scale is used.

One can see that the corrections we applied to the L_X,  MCXC values did not introduce any spatial anisotropies. The lowest fraction L_X,  ours/L_X,  MCXC = 1.005 is found toward (l, b) ∼ (320°, −46°) while the highest fraction L_X,  ours/L_X,  MCXC = 1.071 is found toward (l, b) ∼ (147°, −15°).

Next, we test the isotropy of the processing step from the parent catalogs to MCXC. We follow the same procedure as before, using the 313 clusters of our sample. As shown in the bottom panel of Fig. A.9, the MCXC homogenization of the original L_{X,  parent} is greatly isotropic. The lowest fraction L_X,  MCXC/L_{X,  parent} = 0.964 is found toward (l, b) ∼ (148°, −17°) while the largest L_X,  MCXC/L_{X,  parent} = 0.999 is found toward (l, b) ∼ (156°, + 54°).

Thus, no anisotropic bias was introduced going from the original catalogs to our sample. The last step to be tested is the original L_X measurement from REFLEX, NORAS and eBCS using the RASS data. Such a procedure clearly cannot be checked unless we remeasure the cluster fluxes from the RASS data ourselves. However, there is no obvious reason why such a directional behavior would exist in the original analysis, especially since for the vast majority of clusters (∼88% of the sample) the analysis was conducted in a self-consistent way by the same authors (REFLEX/NORAS).

A.9. ROSAT vs XMM-Newton L_X measurements

It has already been discussed that ROSAT and XMM-Newton return consistent L_X values for the same clusters. As an additional test, we compare our L_X,  ours values with the one derived by Pratt et al. (2009) (L_{X,  Pratt09}, XMM-Newton values) for the 19 common clusters between the two samples. For that, we calibrated our values using the same z as for L_{X,  Pratt09}. The comparison is portrayed in Fig. A.10.

The weighted mean for the L_X,  ours/L_{X,  Pratt09} fraction is 1.001 ± 0.150 and highlights that the luminosity measurements values based on the two different telescopes and studies agree with each other. There is no cluster more than 2σ away from the 1:1 line. Finally, the five clusters that are located in low-A regions in our analysis do not show a different behavior (weighted mean L_X,  ours/L_{X,  Pratt09} = 0.954 ± 0.116) compared to the rest of the clusters.

Appendix B: R₅₀₀ and temperature dependence on cosmology

The most useful feature of the X-ray galaxy cluster L_X–T scaling relation for cosmological isotropy studies is that the determination of the temperature is insensitive to cosmology. The only way that T can be affected by cosmological parameters is through the angular diameter distance D_A and the apparent size of R₅₀₀. The latter is used to select the area from which the spectrum is extracted and the X-ray cluster parameters are constrained. Below we show that the way R₅₀₀ is determined and used in our work is almost independent of changes in cosmological parameters (in particular H₀ which we fit), which propagates to the T determination.

The apparent size of R₅₀₀ which we use, is in arcmin. Thus, it is equal to the physical size of R₅₀₀ in Mpc over D_A. Moreover, the $R_{500}^{\mathrm{Mpc}}$ is derived based on the $L_{\mathrm{X}}\sim E{(z)}^{7/3}{M}_{500}^{1.64}$ relation of Arnaud et al. (2010) and the fact that $M_{500}~R_{500}^3{H}_0^{2}E{(}^z$. With the measured redshift z of the clusters remaining unchanged and the cosmological parameters Ω_m and Ω_Λ fixed to global values (as done in Sect. 7), this can be written as a function of H₀ as shown in Eq. (B.1).

$$\begin{array}{c}\hfill R_{500}^{\mathrm{arcmin}}={\displaystyle \frac{R_{500}^{\mathrm{Mpc}}}{D_{\mathrm{A}}}}\sim {\left({\displaystyle \frac{M_{500}}{H_0^2}}\right)}^{1/3}\frac{1}{D_{\mathrm{A}}}\sim {\displaystyle \frac{L_{\mathrm{X}}^{0.203}}{H_0^{0.667}}}\frac{1}{D_{\mathrm{A}}}·\end{array}$$(B.1)

The luminosity L_X depends on H₀ only through the luminosity distance D_L ∼ 1/H₀. This dependance writes as $L_{\text{X}}\sim D_{\text{L}}^2\sim 1/H_0^2$. Moreover, it also holds that D_A ∼ 1/H₀. Plugging these two relation in Eq. (B.1) results in:

$$\begin{array}{c}\hfill R_{500}^{\mathrm{arcmin}}\sim {\displaystyle \frac{H_0^{-0.406}}{H_0^{0.667}}}{H}_0\Rightarrow R_{500}^{\mathrm{arcmin}}\sim H_0^{-0.073}.\end{array}$$(B.2)

Consequently, a 20% in H₀, which is similar to the H₀ deviations we obtain in Sect. 7, it would only cause a ∼1% change in $R_{500}^{\mathrm{arcmin}}$ with a similar change in the measured T. At the same time, it would cause a ∼45% change in L_X. Additionally, due to the above, the angular radius within which L_X, 500 is measured does not significantly change as well. Thus, one can safely neglect the impact of H₀ anisotropies on the measured flux through the selection of the apparent radius.

All these strongly demonstrate the usefulness of the L_X–T relation for cosmic anisotropies studies.

Appendix C: Table of galaxy cluster data

All Tables

All Figures