A Test of the Cosmological Principle with Quasars

We study the large-scale anisotropy of the Universe by measuring the dipole in the angular distribution of a flux-limited, all-sky sample of 1.3 million quasars observed by the Wide-field Infrared Survey Explorer (WISE). This sample is derived from the new CatWISE2020 catalog, which contains deep photometric measurements at 3.4 and 4.6 um from the cryogenic, post-cryogenic, and reactivation phases of the WISE mission. While the direction of the dipole in the quasar sky is similar to that of the cosmic microwave background (CMB), its amplitude is over twice as large, rejecting the canonical, exclusively kinematic interpretation of the CMB dipole with a p-value of 10^{-4} (3.9 sigma), the highest significance achieved to date in such studies. Our results are in conflict with the cosmological principle, a foundational assumption of the concordance Lambda$CDM model.


INTRODUCTION
The standard Friedmann-Lemaître-Robertson-Walker (FLRW) cosmology is based on the cosmological principle, which posits that the universe is homogeneous and isotropic on large scales. This assumption is supported by the smoothness of the CMB, which has temperature fluctuations of only ∼ 1 part in 100,000 on small angular scales. These higher multipoles of the CMB angular power spectrum are attributed to Gaussian density fluctuations created in the early universe with a nearly scaleinvariant spectrum, which have grown through gravitational instability to create the large-scale structure in the present universe. The dipole anisotropy of the CMB is however much larger, being ∼ 1 part in 1000 as observed in the heliocentric rest frame. This is interpreted as due to our motion with respect to the rest frame in which the CMB is isotropic, and is thus called the kinematic dipole. This motion is usually attributed to the gravitational effect of the inhomogeneous distribution of matter on local scales, originally dubbed the "Great Attractor" (see, e.g., Dressler 1991).
A consistency check would be to measure the concomitant effects on higher multipoles of the CMB anisotropy (Challinor & van Leeuwen 2002); however, even the precise measurements of these by Planck allow up to 40% of the observed dipole to be due to effects other than the Solar System's motion (see discussion in Schwarz et al. 2016). According to galaxy counts in large-scale surveys the universe is sensibly homogeneous when averaged over scales larger than 100 Mpc, as is indeed expected from considerations of structure formation in the concordance ΛCDM model. Hence the reference frame of matter at still greater distances should converge to that of the CMB; i.e. the dipole in the distribution of cosmologically distant sources, induced by our motion via special relativistic aberration and Doppler shifting effects, should align both in direction and in amplitude with the CMB dipole. Independent measurements of the distant matter dipole are therefore an important test of the cosmological principle, and equivalently of the standard model of cosmology. Ellis & Baldwin (1984) proposed that such a test be done using counts of radio sources. These are typi-cally active galactic nuclei (AGNs) at moderate redshift (z ∼ 1), so locally clustered sources (z < 0.1), which can introduce an additional dipole in the distribution of matter (e.g., Tiwari & Nusser 2016), are not a significant contaminant. Consider a population of sources with power-law spectra S ν ∝ ν −α , and integral source counts per unit solid angle dN/dΩ (> S ν ) ∝ S −x ν , above some limiting flux density S ν . If we are moving with velocity v (≪ c) with respect to the frame in which these sources are isotropically distributed, then being "tilted observers" we should see a dipole anisotropy of amplitude (Ellis & Baldwin 1984): The advent of the 1.4 GHz NRAO VLA Sky Survey (NVSS; Condon et al. 1998), which contains ∼ 1.8 million sources, enabled the first estimates of the matter dipole anisotropy (Blake & Wall 2002;Singal 2011;Gibelyou & Huterer 2012;Tiwari & Nusser 2016). To improve sky coverage, data was added from other radio surveys, e.g. the 325 MHz Westerbork Northern Sky Survey (WENSS; Rengelink et al. 1997;Rubart & Schwarz 2013), the 843 MHz Sydney University Molonglo Sky Survey (SUMMS; Mauch et al. 2003;Colin et al. 2017;Tiwari & Aluri 2019) and the 150 MHz TIFR GMRT Sky Survey (TGSS; Bengaly et al. 2018;Singal 2019). However as was first noted by Singal (2011), while the direction of the matter dipole is consistent with that of the CMB, its amplitude is several times larger.
In this Letter, we report the first independent measurement of the dipole in the angular distribution of distant quasars using mid-infrared data from the Wide-field Infrared Survey Explorer (WISE; Wright et al. 2010), which surveyed the sky at 3.4 µm, 4.6 µm, 12µm, and 22 µm (W1, W2, W3, and W4). This provides a measurement of the dipole that is independent of the radio survey-based results, as WISE is a space mission with its own unique scanning pattern, not constrained by the same observational systematics that affect groundbased surveys, such as declination limits or atmospheric effects. While WISE, along with 2MASS, has been used before to set useful constraints on the matter dipole (Gibelyou & Huterer 2012;Yoon et al. 2014;Alonso et al. 2015;Bengaly et al. 2017;Rameez et al. 2018), these studies were of relatively nearby galaxies (z ∼ 0.05 − 0.1) where contamination from local sources can be significant and has to be carefully accounted for. In Section 2, we detail the quasar sample that we use, and we introduce our methodology in Section 3. Our results are presented in Section 4, and we discuss their significance for cosmology in Section 5.

QUASAR SAMPLE
Because of the unique power of mid-infrared photometry to pick out AGNs, WISE may be used to create reliable AGN/quasar catalogs based on mid-infrared color alone (e.g., Secrest et al. 2015). We require an AGN sample optimized for cosmological studies, so the objects should preferably be quasars: AGN-dominated and at moderate or high-redshift (z 0.1; cf., Tiwari & Nusser 2016). The sample should cover as much of the celestial sphere as is possible to minimize the impact of missing (or masked) regions, and be as deep as possible to contain the largest number of objects and thus have the greatest statistical power.
We created a custom quasar sample from the new CatWISE2020 data release (Eisenhardt et al. 2020), which contains sources from the combined 4-band cryo, 3-band cryo, post-cryo NEOWISE, and reactivation NEOWISE-R data.
The CatWISE2020 catalog is 0.71 mag and 0.45 mag deeper in W1 and W2 than the previous AllWISE catalog. We select all sources in the CatWISE2020 catalog with valid measurements in W1 and W2, which are the most sensitive to AGN emission (e.g., Stern et al. 2012). We cut out any sources with possible saturation, as well as sources flagged as suffering from possible contaminants. To select AGNs, we impose the color cut W1 − W2 ≥ 0.8 (Stern et al. 2012), which ensures that the spectral energy distribution is AGN-dominated, following a power-law distribution (S ν ∝ ν −α ) that is insensitive to heavy dust reddening at shorter wavelengths. This yields a raw sample of 174,701,084 objects.
We then remove low-redshift AGNs by excluding sources in the 2MASS extended source catalog (XSC; Jarrett et al. 2000), which contains nearly all galaxies not directly behind the Galactic plane out to z ∼ 0.1 (Jarrett 2004). We made the sample uniform (equal depth) across the sky by addressing several known causes of non-uniformity in the WISE data. The first is a decrement of high quality measurements along the Galactic plane where source confusion is prevalent. This can be mitigated by masking the sky below some Galactic latitude; we find |b| > 30 • to be effective in completely removing non-uniformity because of the Galaxy. The second is poor-quality photometry near clumpy and resolved nebulae both in our Galaxy (e.g., planetary nebulae) and in nearby galaxies such as the Magellanic Clouds and Andromeda. We remove these by masking out 6 times the 20 mag arcsec −2 isophotal radii from the 2MASS Large Galaxy Atlas (LGA; Jarrett et al. 2003). The third is a decrement of sources, and the presence of image artifacts, near bright stars, caused by density sup-pression in their vicinity. 1 We find that circular masks with 2MASS K band-dependent radii log 10 (r/deg) = −0.134K − 0.471 effectively remove these. In all, we masked 265 sky regions, plus the Galactic plane. To avoid any directional source count bias we mirror the masks by 180 • on the celestial sphere.
We calculate spectral indices α of our sources in the W1 band by obtaining power-law fits of the form S ν = kν −α , where k is the normalization. We produced a lookup table to determine α based on W1 − W2, by calculating synthetic AB magnitudes following Equation 2 of Bessell & Murphy (2012). The WISE magnitudes are on the Vega magnitude system, so we convert from the AB system using the offsets m AB −m Vega = 2.673, 3.313 for W1 and W2, respectively. These WISE offsets correspond to the constant of −48.60 associated with the definition of the synthetic AB magnitude. The normalisation k is calculated by inverting the equation for the synthetic magnitude and using the observed W1 AB magnitude. Finally, we calculate the isophotal frequency, at which the flux density S ν equals its mean value within the passband, using Equation A19 in Bessell & Murphy (2012). As our sample was constructed with the cut W1 − W2 ≥ 0.8, the distribution peaks at α ∼ 1 and extends to steeper slopes, with a mean value of 1.26. Distributions of spectral indices and fluxes for our final sample of sources are shown in Figure 1. The corresponding mean isophotal frequency is 8.922 × 10 13 Hz, with a dispersion of 0.19%. We select a magnitude cut of 9 > W1 > 16.4 (Vega), equivalent to a flux density cut of 77.77 > S ν > 0.09 mJy, to fix the over-density of fainter sources along overlaps in the WISE scanning pattern, most prevalent at the ecliptic poles where they converge. After removing low-z AGNs, applying the sky masks, and making the flux density cut, our final sample has 1,314,428 AGNs, which we show in Figure 2.
To estimate the distribution of AGN redshifts, we select those within SDSS Stripe 82, a 275 deg 2 region of the sky scanned repeatedly by the SDSS, thus achieving an increase of depth of ∼ 2 mag (Annis et al. 2014). In the specObj table for SDSS DR16, 2 Stripe 82 contains ∼ 4.4 times more objects with spectroscopic r-band magnitudes fainter than 20 (AB) than a comparable sky region in the SDSS main footprint. We use a sub-region of Stripe 82 between −42 • < R.A. < 45 • , which lies outside the |b| < 30 • Galactic plane mask we employ, and which was observed by the Extended Baryon Oscillation Spectroscopic Survey (eBOSS; Dawson et al.   2016), yielding even deeper spectral coverage. There are 14,387 CatWISE AGNs in this region. For photometric information, we cross-match these with the Dark Energy Survey (DES), Data Release 1 (Abbott et al. 2018), which achieved an i-band depth of 23.44 mag (AB). Using a 10 ′′ match for completeness, we find counterparts for 14,343 (99.7%) of the CatWISE AGNs. Matching the DES counterpart coordinates onto specObj table to within 1 ′′ for fiber coverage, we find 8612 matches (60%). The unmatched objects are 0.3 mag fainter in W2 than the matched objects on average, suggesting that they are slightly less luminous or slightly more distant (or both). However, their mean r − W2 value, a measure of AGN obscuration level (e.g., Yan et al. 2013), is ∼ 1.8 mag redder than the mean of the matched sample, implying that the unmatched objects are simply too faint at visual wavelengths for SDSS. Indeed, while about one-third of the full DES-matched sample has r − W2 > 6 mag (Vega), in line with expectations from the literature for the prevalence of type-2 AGNs (Yan et al. 2013), 77% of the unmatched sample have r − W2 > 6. This indicates that the objects in our sample without SDSS spectra are predominantly type-2 systems, an effect of the orientation of the AGN with respect to the line of sight, and so the matched objects may be used to estimate the distribution of redshifts for the full sample. We find a mean redshift of 1.2, with 99% having z > 0.1, i.e. the Cat-WISE AGN sample is not contaminated by low-redshift AGNs. The redshift distribution of our sample is shown in Figure 3.

Dipole Estimator
We determine the dipole of our sample with the 3dimensional linear estimator: wherer i is the unit vector pointing to source i, and N is the sample size. This estimator simply calculates the mean resultant length and direction of the N unit vectors and is agnostic with regard to the true underlying signal (e.g., Fisher et al. 1987), as opposed to other estimators (e.g., Blake & Wall 2002;Bengaly et al. 2019) which explicitly seek a dipolar pattern. However, if the signal has a dipolar form then Equation 2 generally has a bias in both amplitude and direction (Rubart & Schwarz 2013) induced by Poisson noise and masking. We account for amplitude bias in our results as well as in the estimates of their significance using Monte Carlo methods, correcting for directional bias as discussed in Appendix A. We further confirm our results by employing the quadratic estimator D q which does not suffer from bias and is evaluated by minimising the quantity (e.g., Bengaly et al. 2019): where n p denotes the number of sources in sky pixel p (r p being the unit vector to the pixel) and the sum is to be taken over all unmasked pixels (in which n is the average number of sources). Due to significantly higher computational expense for the quadratic estimator, we run simulations only with the linear estimator.

Mock data and statistical significance
We generate mock samples of N init vectors drawn from a statistically isotropic distribution, whose directions are subsequently modified by special relativistic aberration according to an observer boosted with velocity β. Each sample is then masked with the same mask that was applied to the data (Figure 2). In order to respect the exact distribution of flux values and spectral indices in the data, we assign to each simulated source a flux density S ν and a spectral index α drawn at random from their empirical distributions before applying the flux density cut (Figure 1). The sampled fluxes are now modulated depending on source position, velocity β, and α. Lastly, only sources with S ν > S ν,cut are retained, and the number of remaining sources is finally reduced to that of the true sample, N , through random selection.
Under the null hypothesis that the measured dipole D l is a consequence of our motion with respect to a frame shared by both quasars and the CMB, we generate a set of mock skies according to the above recipe. For each random choice we record D sim l , and correct for its directional bias using Equations A3 and A4. The fraction of mock skies with amplitude | D l | larger than our empirical sample, and with angular distance from the CMB dipole closer than our sample, gives the p-value with which the null hypothesis is rejected. Note that the effect on our results of the distributions of flux and spectral index (Figure 1) is automatically included via the bootstrap approach employed for our simulations.

RESULTS
Our sample of 1,314,428 quasars exhibits a dipole with amplitude: D l = 0.0173. Correcting for the directional bias induced by the mask employed, we find that it points in the direction: (l, b) = (234. • 1, 29. • 2). This is 29. • 8 from the direction of the CMB dipole (l, b = 264. • 021, 48. • 253;Planck Collaboration et al. 2018). However, when the expected dipole is simulated assuming the kinematic interpretation of the CMB dipole, only 4 out of 40,000 such simulations give D sim l with an amplitude larger than the observed value (left panel, Figure 4) and within 29. • 8 of the CMB dipole direction as for our sample. We can therefore reject the null hypothesis with a p-value of 10 −4 corresponding to a significance of 3.9σ.
If we assume that the anomalous quasar dipole is still of kinematic origin, albeit with a velocity different from that inferred from the CMB, we can estimate its directional uncertainty. To avoid bias, we first compute the dipole with the quadratic estimator D q , which gives D q = 0.01629 towards (l, b) = (234. • 0, 27. • 4). The corresponding velocity from Equation 1, with (median) α = 1.17 and index x = 1.7 at the flux density cut, is v = 861 km s −1 . A set of 15,000 simulations with this input velocity is then performed to find the directional uncertainty. The right panel of Figure 4 shows this as a patch around the (consistent) dipole direction obtained with both estimators.

DISCUSSION
The CatWISE AGN sample exhibits an anomalous dipole, oriented similarly to the CMB dipole but over twice as large. Whereas a "clustering dipole" is expected from correlations in the spatial distribution of the sources, this can be estimated knowing their autocorrelation function (or power spectrum) and distribution in redshift (see Appendix B). It is smaller by a factor of ∼ 60 than the dipole we observe in these higher redshift quasars.
The unique statistical power of our study has allowed us to confirm the anomalously large matter dipole suggested in previous work, which used objects selected at a different wavelength (radio), using surveys completely independent of WISE, viz. NVSS, WENNS, SUMMS, and TGSS. The ecliptic scanning pattern of WISE has no relationship with the CMB dipole, so there is no rea-son to suspect that the dipole we measure in the Cat-WISE AGN catalog is an artifact of the survey.
After Ellis & Baldwin (1984) proposed this important observational test of the cosmological principle, agreement was initially claimed between the dipole anisotropy of the CMB and that of radio sources (Blake & Wall 2002). If the rest frame of distant AGNs is indeed that of the CMB, it would support the consensus that there exists a cosmological standard of rest, related to quantities measured in our heliocentric frame via a local special relativistic boost. This underpins modern cosmology: for example, the observed redshifts of Type Ia supernovae are routinely transformed to the "CMB frame". From this it is deduced that the Hubble expansion rate is accelerating (isotropically), indicating dominance of a cosmological constant, and this has led to today's concordance ΛCDM model. If the purely kinematic interpretation of the CMB dipole that underpins the above procedure is in fact suspect, then so are the important conclusions that follow from adopting it. In fact, as observed in our heliocentric frame, the inferred acceleration is essentially a dipole aligned approximately with the local bulk flow of galaxies and towards the CMB dipole (Colin et al. 2019), so cannot be due to a cosmological constant.
If it is established that the distribution of distant matter in the large-scale universe does not share the same reference frame as the CMB, then it will become imperative to ask whether the differential expansion of space produced by nearby nonlinear structures of voids and walls and filaments can indeed be reduced to just a local boost (Wiltshire et al. 2013). Alternatively the CMB dipole may need to be interpreted in terms of new physics, e.g. as a remnant of the pre-inflationary universe (Turner 1991). Gunn (1988) had noted that this issue is closely related to the bulk flow observed in the local universe, which in fact extends out much further than is expected in the concordance ΛCDM model (e.g., Colin et al. 2011;Feindt et al. 2013). Further work is needed to clarify these important issues.
As Ellis & Baldwin (1984) emphasized, a serious disagreement between the standards of rest defined by distant quasars and the CMB may require abandoning the standard FLRW cosmology itself. The importance of the test we have carried out can thus not be overstated.