© The Authors 2026

1. Introduction

The Λ cold dark matter (ΛCDM) model, validated by decades of observational data, serves as the standard framework of modern cosmology. This model successfully reconciles the observed accelerated expansion of the Universe with the theory of general relativity, resting on two foundational pillars: the theory of gravity itself and the cosmological principle (CP).

The CP, a core postulate of the concordance ΛCDM model, posits that the Universe is statistically homogeneous and isotropic on scales larger than ∼100 Mpc. This assumption forms the basis of the standard Friedmann-Lemaître-Robertson-Walker (FLRW) cosmology.

Direct support for large-scale isotropy comes from the high degree of uniformity observed in the cosmic microwave background (CMB), which appears remarkably isotropic on small angular scales. Complementary evidence is provided by a variety of independent cosmological tracers, including large-scale galaxy surveys (Gibelyou & Huterer 2012; Sarkar et al. 2018; Franco et al. 2023), the X-ray background (Wu et al. 1999), the spatial distribution of galaxy clusters (Bengaly et al. 2016), radio galaxies (Blake & Wall 2002), and Type Ia supernovae (Gupta & Saini 2010; Sorrenti et al. 2023). Collectively, these observations robustly corroborate the statistical isotropy of the Universe on the largest observable scales.

Amidst this backdrop of large-scale isotropy, the CMB itself exhibits a hierarchy of anisotropies. It contains small-scale fluctuations at the level of approximately 1 part in 10⁵, and a far more prominent dipole anisotropy of amplitude ∼1 part in 10³ as measured in the heliocentric frame. Given its distinctive amplitude and dipolar pattern, this dominant signal is almost universally attributed to a kinematic effect–the peculiar motion of the Solar System relative to the CMB rest frame. The latest measurements infer a velocity of v = 369.825 ± 0.11 km/s towards the Galactic coordinates (l, b) = (264.021° ,48.253° ) (Hinshaw et al. 2009; Aghanim et al. 2020).

A critical test of the CP involves verifying the purely kinematic origin of the CMB dipole. This is achieved by measuring a corresponding dipole anisotropy in the distribution of distant matter. If the CP holds, the peculiar velocity inferred from the matter rest frame must be consistent with that derived from the CMB. The Ellis–Baldwin test provides a direct implementation of this idea (Ellis & Baldwin 1984). It posits that if the CMB dipole is indeed due to our motion through a homogeneous and isotropic universe, then a statistically congruent dipole must be present in the number counts of distant sources, such as radio galaxies. The velocity inferred from these counts should agree with the CMB-derived value. Consequently, under the CP, the large-scale distribution of matter should appear isotropic when viewed from the CMB rest frame. Any significant deviation would necessitate a physical explanation and could potentially challenge the validity of the CP itself.

However, this expectation has been challenged by several findings over the past two decades. Singal (2011) reported a significant discrepancy, finding the dipole amplitude from discrete radio sources to be approximately four times larger than that of the CMB, albeit with a possibly consistent direction. A congruent result was obtained by Rubart & Schwarz (2013) using the NVSS and WENSS surveys, who confirmed a dipole direction aligned with the CMB but an amplitude larger by a factor of approximately three. Extending the analysis to Type Ia supernovae, Singal (2022) similarly found a dipole amplitude about four times greater than the CMB value. In a distinct finding from the Pantheon+SH0ES supernova data, Sorrenti et al. (2023) reported a dipole amplitude consistent with the CMB but a direction differing with a statistical significance exceeding 3σ. Earlier, using the Union2.1 supernovae data, Yang et al. (2013) also found a preferred dipole direction via a dipole-fitting approach.

Secrest et al. (2021) used quasars (QSOs) from the CatWISE2020 catalogue and found a dipole direction offset by ∼27° from the CMB dipole, but with an amplitude over twice as large as expected, rejecting the kinematic interpretation at the 4.9σ level. Separately, Dam et al. (2023) also reported a significant tension through a Bayesian analysis of QSO data, measuring a dipole amplitude approximately 2.7 times larger than the kinematic expectation with a statistical significance of 5.7σ. Most recently, an analysis combining NVSS, RACS-low, and LoTSS-DR2 radio surveys reported a dipole amplitude 3.67 ± 0.49 times larger than the kinematic expectation, a 5.4σ discrepancy (Böhme et al. 2025). Secrest et al. (2022) further analysed radio galaxies and QSOs; they suggest that if the CMB dipole is assumed to be fully kinematic, the reference frames of these distant sources may exhibit an intrinsic anisotropy.

Other studies also hint at potential violations of statistical isotropy. For instance, isotropy tests performed on the SDSS Luminous Red Galaxy sample indicate a possible breakdown on scales as large as ∼300 h⁻¹ Mpc (Park et al. 2017). In the CMB itself, a persistent hemispherical power asymmetry remains significant at an approximate the 3σ level in the latest data (Eriksen et al. 2004; Planck Collaboration VII 2020; Planck Collaboration XVI 2016). Notably, several of these anomalous signals, including the dipole in radio polarization and the alignment of the CMB low-order multipoles, indicate a preferred direction closely aligned with the CMB dipole (Ralston & Jain 2004). Furthermore, the discovery of enormous cosmic structures–such as large QSO groups (Clowes et al. 2013), correlated QSO orientations (Friday et al. 2022), and giant cosmic voids (Haslbauer et al. 2020; Keenan et al. 2013)–challenges the assumption of homogeneity at the largest scales.

The cumulative evidence from these diverse studies, which utilize different tracers and methodologies and sometimes yield conflicting results, calls into question the foundational assumption of large-scale isotropy and homogeneity inherent in the CP. In this context, we employed an independent cosmological probe to estimate the Solar System’s peculiar motion by analysing the dipole anisotropy in the sky distribution of QSO redshifts. The redshift dipole offers a direct measurement of the peculiar velocity relative to the matter reference frame traced by the QSOs. This method circumvents the need to account for local over- or under-densities (e.g. – a local void or specific clustering) and other large-scale structures (LSSs) known to potentially bias traditional number-count dipole measurements (Oayda et al. 2024; Rameez et al. 2018; Tiwari & Nusser 2016; Rubart et al. 2014).

This study aims to employ the high-precision, wide-coverage QSO data from the Dark Energy Spectroscopic Instrument (DESI) Data Release 1 (DR1) to perform a tomographic estimation of the dipole modulation across redshift bins, following the methodology of da Silveira Ferreira & Marra (2024). We used the DESI DR1 QSO sample to infer the Solar System’s peculiar velocity via the aberration effect, thereby performing a direct test of the CP through the Ellis–Baldwin test. This test confronts the kinematic dipole measured in the matter frame (traced by distant QSOs) with that of the cosmic rest frame (defined by the CMB).

This paper is structured as follows. Section 2 details the DESI and eBOSS QSO catalogues and their sky footprints used in the analysis. Section 3 presents the formalism for the redshift dipole, outlines the least-squares estimator used to measure the dipole amplitude, justifies the separate analysis of the northern Galactic caps (NGC) and the southern Galactic caps (SGC) to account for systematic depth variations, and introduces Fisher’s method for evaluating the overall statistical significance. Finally, Section 4 presents the main results and discusses their implications, concluding that the QSO dipole inferred from DESI DR1 is consistent with the kinematic interpretation of the CMB dipole, thereby finding no compelling evidence against the CP.

2. Data

We aim to test the CP using LSS tracers from recent spectroscopic surveys. To ensure robust results, we selected catalogues based on data quality, minimal known systematics, substantial sky coverage, and the ability to perform survey comparisons. Specifically, we utilized the final QSO clustering catalogues from the SDSS Extended Baryon Oscillation Spectroscopic Survey (eBOSS) and the early data release of the Dark Energy Spectroscopic Instrument (DESI). These catalogues provide observational weights, redshifts, and angular positions, and their derived cosmological parameters are consistent with the ΛCDM framework.

We analysed three distinct QSO samples:

(i) The eBOSS DR14 QSO Catalogue (Abolfathi et al. 2018): This catalogue contains 194 801 objects in the redshift range 0.04 < z < 6.95. After removing extreme outliers, 90% of the data lie within 0.42 < z < 3.28, leaving 192 853 objects. To mitigate potential boundary effects, we further restricted our analysis to the well-populated range 0.4 < z < 2.8, resulting in a final sample of 187 954 QSOs. The spectra cover the wavelength region 3610−10 140 Å at a spectral resolution in the range 1300 < R < 2500 (Abolfathi et al. 2018).

(ii) The eBOSS DR16 QSO Catalogue (Ross et al. 2020): This catalogue contains 343 708 QSOs with a tight and uniform redshift distribution in the range 0.8 < z < 2.2, selected from the final SDSS-IV/eBOSS quasar sample (Lyke et al. 2020). We used the entire specified range without further cuts. The spectra cover 3600−10 400 Å with a spectral resolution of R ∼ 2000 (Lyke et al. 2020). At the largest QSO sample at the time, eBOSS DR16 provided important data for cosmological analysis.

(iii) The DESI DR1 QSO Catalogue (DESI Collaboration 2025): The Dark Energy Spectroscopic Instrument (DESI) is a Stage IV spectroscopic survey designed to obtain tens of millions of galaxy and QSO redshifts to constrain dark energy and cosmological models (Adame et al. 2025). Its first data release (DR1) encompasses the first year of main survey observations. The public large-scale structure catalogue for QSOs contains 1 223 170 objects with confident redshift measurements in the range 0.8 < z < 3.5. The spectra cover the wavelength region 3600−9824 Å at a spectral resolution in the range 2000 < R < 5200 (DESI Collaboration 2025).

The construction of this LSS catalogue, detailed in Adame et al. (2025), involves critical steps to correct for systematic variations in the observed number density and to ensure its suitability for cosmological clustering analyses, as follows.

Systematic weights: The catalogue incorporates weights to correct for two main sources of spurious density fluctuations: (1) imaging systematics (e.g. seeing, depth, and stellar density) affecting the input target selection, and (2) variations in the ability to confidently measure redshifts from DESI spectra, which depend on observing conditions and instrumental performance.

Reference randoms: A matched set of synthetic ‘randoms’ is provided to accurately model the survey selection function, geometry, and fibre assignment incompleteness.

Validation: The two-point clustering measurements from the final LSS catalogue have been compared to high-fidelity simulations of the DESI survey. They are found to be in general statistical agreement, typically within a 2% factor in the inferred galaxy bias, confirming that observational systematic uncertainties are sub-dominant to statistical uncertainties for large-scale analyses.

For our dipole analysis, we used the total systematic weight w_i provided in the catalogue, which combines the corrections mentioned above. To avoid potential biases from sparse edges of the redshift distribution, we restricted our analysis to the well-characterized range 0.8 < z < 3.0, yielding a final sample of 1 176 570 QSOs.

The redshift range and the total number of sources for each final sample are summarized in Table 1. The redshift distributions are visualized in Figure 1 as scatter plots with a bin width of Δz = 0.001

Fig. 1. Redshift distributions of the QSO catalogues. Left to right: eBOSS DR14 (orange), eBOSS DR16 (green), and DESI DR1 (blue). The DR14 distribution is broad (0.4 < z < 8) but sparse at the extremes; we use the well-populated range 0.4 < z < 2.8. The DESI DR1 distribution spans 0.8 < z < 3.5; we use 0.8 < z < 3.0 to avoid edge effects. All histograms use a bin width of Δz = 0.001. Note that the y-axis scales differ between panels. The primary overlap region for the eBOSS samples is 0.8 < z < 2.2, while DESI significantly overlaps in the range 0.8 < z < 2.3.

The sky footprints of the three catalogues are shown in Figure 2. To address potential hemispheric systematic biases from instrumentation or calibration, we performed a detailed comparison of the redshift distributions between the NGC and the SGC. As shown in Figure 3, we used bootstrap resampling to estimate the mean redshift in each cap. The NGC mean is taken as a reference (centred at zero). While the shapes of the NGC and SGC distributions are similar, they do not perfectly overlap. The systematic offset is approximately 0.008 for the eBOSS samples and 0.001 for the DESI samples, both of which are of the same order as the expected kinematic signal β ≡ v/c ∼ 0.001. This significant offset could bias a combined NGC+SGC dipole estimate. Therefore, we implemented an unbiased estimation method, described in Section 3.3, which analyses the two hemispheres separately.

Fig. 2. Sky footprints of the QSO samples in Mollweide projection (HEALPix format with Nside = 64). Left to right: eBOSS DR14, eBOSS DR16, and DESI DR1. The maps show the number count per pixel for each catalogue.

3. Method

3.1. Redshift dipole

The kinematic redshift dipole arises from the Doppler shift due to our peculiar motion relative to the cosmic rest frame. In an isotropic universe, the observed redshift distribution of distant sources should exhibit a dipolar anisotropy modulated by this motion. The dominant contribution to the observed redshift dipole Δ is our peculiar velocity v, expressed as β = v/c, where c is the speed of light. The theoretical relationship is given by Nadolny et al. (2021):

$$\begin{array}{c}\hfill \mathbf{\Delta }=-\mathit{\beta }+{\mathbf{\Delta }}_{\mathrm{cd}}+\mathcal{O}({\beta }^2),\end{array}$$(1)

where Δ_cd is the clustering dipole from large-scale structure, and 𝒪(β²) represents higher-order terms

Due to our motion, the observed frequency (ν) of light from a source in direction $\widehat{\mathit{n}}$ is related to its rest-frame frequency ν_rest by

$$\begin{array}{c}\hfill \nu ={\nu }_{\mathrm{rest}}\delta (\widehat{\mathit{n}},\mathit{\beta }),\mathrm{where}\delta (\widehat{\mathit{n}},\mathit{\beta })=\frac{1+\widehat{\mathit{n}}·\mathit{\beta }}{\sqrt{1-{\beta }^2}}·\end{array}$$(2)

Using c = λν and the definition of redshift z = (λ_observe − λ_emit)/λ_emit, we derived the Lorentz boost relation between the observed redshift z′ (in the observer’s frame) and the intrinsic redshift z (in the CMB rest frame):

$$\begin{array}{c}\hfill 1+z\prime =(1+z)\delta {(\widehat{\mathit{n}},\mathit{\beta })}^{-1}.\end{array}$$(3)

The term z represents the expected redshift in an isotropic universe.

For high-redshift sources (z ≳ 1), such as the QSOs used here, the clustering dipole contribution Δ_cd is significantly suppressed compared to the kinematic term (Tiwari & Nusser 2016). Therefore, to first approximation, we assumed Δ ≈ −β when translating a measured dipole into a peculiar velocity.

3.2. Estimator

Following the tomographic methodology of da Silveira Ferreira & Marra (2024), we estimated the kinematic dipole by analysing the redshift distribution in fine bins. This approach mitigates the impact of large-scale clustering noise and allows us to probe potential redshift evolution of the signal.

We divided the QSO sample into narrow redshift bins of width Δz = 0.001. For each bin, we modelled the observed redshift z′_i as a Doppler-shifted version of an intrinsic, isotropic redshift field. The goal is to determine the dipole vector Δ = (Δ_x, Δ_y, Δ_z) defined in the Cartesian coordinate system transformed from Galactic coordinates (x towards the Galactic centre, y in the direction of Galactic rotation, and z towards the North Galactic Pole), and a monopole offset m that best match this model. This is achieved by minimizing a per-bin χ² statistic using a least-squares estimator. In practice, we employed the fit_dipole function from the HEALPix package (Gorski et al. 2005), which performs this fit on a spherical harmonic representation of the data within the bin.

The core of the estimator compares the observed redshifts with a simulated isotropic field. For a trial dipole Δ, we first computed the Lorentz boost factor $\delta ({\widehat{\mathit{n}}}_i,-\mathbf{\Delta })$ for each QSO i in direction ${\widehat{\mathit{n}}}_i$. We then constructed a simulated intrinsic redshift for the bin, $z_{\mathrm{bin}}^{\mathrm{simul}}$, by removing the estimated Doppler shift from all objects and computing the weighted mean:

$$\begin{array}{c}\hfill z_{\mathrm{bin}}^{\mathrm{simul}}=\frac{{\sum }_i^N{w}_i(1+z{\prime }_i)\delta ({\widehat{\mathit{n}}}_{\mathit{i}},-\mathbf{\Delta })}{{\sum }_i^N{w}_i}-1.\end{array}$$(4)

This $z_{\mathrm{bin}}^{\mathrm{simul}}$ represents our best estimate of what the redshift would be in the CMB rest frame for that bin, given the current Δ. The N denotes the total number of objects. The total weights w_i = w_i, sys ⋅ w_i, fail are derived from the LSS catalogues, where w_i, sys represents the systematic weight and w_i, fail the redshift failure weight. We emphasize that the weight definition is consistent for both the eBOSS and DESI catalogues used in this analysis. Next, we applied the boost to this intrinsic value to predict what we should observe if Δ were correct: $z{″}_{i,\mathrm{bin}}=(1+z_{\mathrm{bin}}^{\mathrm{simul}})\delta {(\widehat{{\mathit{n}}_{\mathit{i}}},-\mathbf{\Delta })}^{-1}-1$. The χ² for the bin measures the discrepancy between the actual observations z′_i and these predictions z″_i, bin:

$$\begin{array}{c}\hfill {\chi }_{\mathrm{bin}}^2(\mathbf{\Delta })=\frac{{\sum }_i^N{w}_i{[1+z{\prime }_i-(1+z{″}_{i,\mathrm{bin}})]}^2}{{\sum }_i^N{w}_i}·\end{array}$$(5)

The estimator iteratively adjusts Δ to minimize χ_bin². The minimization process aims to fit a reference frame that is as stationary as possible relative to the CMB. This was achieved by removing our peculiar velocity. It should be noted that, since the peculiar velocity of the observed objects is unknown, we relied solely on the Earth’s peculiar velocity and the statistical characteristics of the material reference frame for the fit. As mentioned in Section 3.4, the distant material reference frame should remain relatively stationary with respect to the CMB as a whole.

3.3. Unbiased estimation

As established in Section 2 and Figure 3, significant systematic differences exist between the NGC and SGC. A joint fit to the combined data, assuming a single global monopole, would be biased because the mean redshift (the monopole) differs between hemispheres. To obtain an unbiased estimate of the dipole – which should be a global vector – we must account for these hemispheric offsets.

Fig. 3. Assessment of NGC and SGC systematic differences via bootstrap resampling. For each QSO sample (left: eBOSS DR14; centre: eBOSS DR16; right: DESI DR1), we repeatedly (1000 times) drew bootstrap resamples of 28 000 data points from each hemisphere’s dataset, computing the mean redshift for each resample. The distributions of these resampled means are plotted as histograms with a bin width of 0.000333. The y-axis represents the count of resamples whose mean redshift falls within the corresponding bin. All distributions are plotted relative to the average redshift of the total NGC dataset (vertical line at zero). The x-axis shows the redshift difference Δz (dimensionless) relative to the total NGC mean. The offset between the NGC (centred) and SGC distributions indicates a hemisphere-level systematic discrepancy. The amplitude of this offset (∼0.008 for eBOSS, ∼0.001 for DESI) is comparable to the expected kinematic dipole signal (β = v/c ∼ 10−3), necessitating separate hemispheric analysis (Section 3.3).

The key to our unbiased estimation is to perform the dipole fit separately on the NGC and SGC datasets while enforcing consistency in the resulting dipole vector. If the two Galactic caps are analysed together, the dipole calculation would be adjusted relative to the combined monopole of both hemispheres, introducing systematic differences of the same order–or larger–than the expected dipole signal (cΔ ∼ 300 km/s). We therefore analysed the NGC and the SGC independently.

In practice, within each redshift bin, we minimized χ_bin, NGC²(Δ) and χ_bin, SGC²(Δ), allowing the fit in each hemisphere to account for its own redshift weighting. Crucially, we then required that the dipole vector Δ that minimizes the combined χ_total² in Equation (6) be identical for both hemispheres. This approach effectively avoids the bias introduced by merging the data; it extracts a single, coherent kinematic dipole signal consistent across the full sky.

We implemented this by performing the minimization separately for the NGC and the SGC samples within each redshift bin. This yields two χ² functions: χ_bin, NGC²(Δ) and χ_bin, SGC²(Δ), each already marginalized over their respective hemispheric monopoles. The total χ² for the entire dataset is then the sum over all bins and both hemispheres:

$$\begin{array}{c}\hfill {\chi }_{\mathrm{total}}^2(\mathbf{\Delta })={\displaystyle \sum [{\chi }_{\mathrm{bin},\mathrm{NGC}}^2(\mathbf{\Delta })+{\chi }_{\mathrm{bin},\mathrm{SGC}}^2(\mathbf{\Delta })]}.\end{array}$$(6)

Minimizing χ_total²(Δ) yields the dipole estimate Δ that best fits the data while being insensitive to the NGC-SGC mean redshift offset.

3.4. Fisher’s statistical method

To quantitatively assess the tension between our measured dipole and the prediction from the CMB kinematic dipole (Δ_CMB), we needed a single statistical measure that combines the evidence from all three Cartesian components. Fisher’s combined probability test (Fisher 1925) provides a rigorous framework for this.

For each Cartesian component k ∈ x, y, z, we defined a null hypothesis H₀^(k): the measured dipole component Δ_k is consistent with the CMB-predicted value Δ_CMB, k. We computed the corresponding p-value, P_k, which represents the probability of obtaining a value as extreme as or more extreme than our measurement Δ_k if H₀^(k) were true. Assuming Gaussian errors, P_k is derived from the two-tailed test statistic |Δ_k − Δ_CMB, k|/σ, where $\sigma =\sqrt{{\sigma }_{{\mathrm{\Delta }}_k}^2+{\sigma }_{{\mathrm{\Delta }}_{\mathrm{CMB},k}}^2}\approx {\sigma }_{{\mathrm{\Delta }}_k}$, since the CMB dipole is measured with high precision.

Following the CP and the work of Ellis & Baldwin (1984), the matter reference frame defined by distant objects is expected to be at rest relative to the CMB on large scales, meaning that an isotropic and homogeneous reference frame should not systematically bias any particular dipole component. Based on this, we assumed that the tests for the three orthogonal components are independent, Fisher’s method combines the p-values into a single test statistic:

$$\begin{array}{c}\hfill X^2=-2{\displaystyle \sum _{k=1}^3}ln(P_k).\end{array}$$(7)

If all individual null hypotheses are true and the test statistics |Δ_k − Δ_CMB, k| are independent, the test statistic X² follows a χ² distribution with 2 × 3 = 6 degrees of freedom.

The final combined p-value, P_combined, was calculated as the survival function of the χ₆² distribution at the observed X². A very small P_combined indicates a low probability that all three CMB dipole component predictions are simultaneously consistent with our measurements. This combined p-value was then converted into a significance level in units of Gaussian standard deviations (σ). This metric provides a comprehensive and conservative measure of the overall discrepancy between the large-scale structure dipole and the CMB dipole.

4. Discussion and conclusion

In this study, we presented a test of the CP by measuring the large-scale structure dipole anisotropy using QSO samples from eBOSS DR14 (0.4 < z < 2.8), DR16 (0.6 < z < 2.2), and DESI DR1 (0.8 < z < 3.0). Our analysis, employing the methodology of da Silveira Ferreira & Marra (2024), allows for a direct comparison of the kinematic dipole inferred from the distribution of QSOs with the canonical CMB dipole. The key results, summarized in Table 2 and Fig. 4, reveal a complex picture that warrants careful scientific discussion.

Fig. 4. Empirical distribution plots of the DESI DR1 redshift dipole in its Cartesian components (m, cΔx, cΔy, cΔz). The red curve represents the dipole expected from the CMB dipole. The x-axis represents velocity (in kilometres per second), the y-axis represents distribution density, and the shaded regions denote the 1σ confidence intervals. Notably, the cΔy component shows small deviation from the CMB prediction, while the 1σ ranges of both cΔx and cΔz fully encompass the CMB dipole value.

From the DESI DR1 data, we infer a bulk flow velocity of v = 443.8 ± 204.1 km/s towards (l, b) = (107.4° ±86.8° ,28.4° ±45.2° ). While the velocity magnitude is consistent with the CMB kinematic value of 370 km/s within 1.56σ, the best-fit direction shows a substantial offset from the CMB dipole (l = 264.02° ,b = 48.25° ). However, a component-by-component statistical assessment using Fisher’s method yields an overall tension of 3.45σ for DESI, primarily driven by a deviation in the Δ_y component. To interpret these seemingly contradictory signals – a consistent velocity magnitude but offset Δ_y direction and component level tension – we considered the evidence in a broader context.

Internal consistency and data evolution: The progression from the highly anomalous eBOSS DR14 result to the CMB-consistent eBOSS DR16 result underscores the critical impact of data quality. The DR14 discrepancy is likely attributable to its smaller sample size and higher redshift dispersion, limitations mitigated in DR16. This evolution cautions against over-interpreting significant tensions from any single data release, including the 3.45σ value from DESI DR1.

Methodological validation and external comparison: Our successful replication of the da Silveira Ferreira & Marra (2024) analysis on the eBOSS DR16 QSO sample yields a consistent result. We derive a bulk flow of v = 477.6 ± 224.4 km/s towards (l, b) = (261.4° ±117.3° , − 0.4° ±34.6° ), in agreement with the CMB kinematic dipole at a 2.1σ level (see Fig. A.1). This is in close agreement with their original finding of $v={196}_{-79}^{+92}$ km/s towards $(l,b)=({298}_{-52}^{°+29},{50}_{-62}^{°+8})$, which they reported as consistent with the CMB kinematic expectation at a 2σ level. This independent verification validates our analytical framework and, crucially, confirms that the methodology itself does not inherently produce high-tension results. The discrepancy, therefore, arises specifically with the DESI DR1 sample, suggesting sample-specific factors may be at play.

The nature of the 3.45σ tension in DESI: The significant Fisher combined tension for DESI warrants careful dissection. As Fig. 4 shows, it stems largely from the Δ_y component, whose associated unipolar amplitude differs from the expected value by nearly 300 km/s. This deviation may originate from the data itself, potentially challenging the CP. This pattern–a component-isolated deviation–is more indicative of residual systematic effects or unusual sample variance within DESI’s specific scan pattern and calibration, rather than a coherent, all-sky anisotropy expected from a genuine breakdown of isotropy.

The central question is whether the measured dipole reflects our motion relative to the Hubble flow (kinematic origin, supporting CP) or an intrinsic anisotropy (challenging CP). The amplitude of the DESI-inferred bulk flow is consistent with the kinematic expectation, which is a necessary condition for the kinematic interpretation. However, the offset in the best-fit direction, though not statistically significant given the large uncertainties, introduces ambiguity. Claims of a CP violation require not only high statistical significance but also robustness across multiple independent probes, elimination of plausible systematics, and a consistent directional signal across different redshift slices and tracers. Our DESI result, characterized by a component-level tension and a direction that is not robustly constrained, currently falls short of these stringent criteria.

In conclusion, our analysis employed the inferred bulk flow velocity–derived from the dipole–as the critical test. The amplitude of this velocity, at |v| = 443.8 ± 204.1 km/s, is consistent with the kinematic interpretation of the CMB dipole within 1.56σ. However, its direction (l, b) = (107.4° ±86.8° ,28.4° ±45.2° ) shows a notable offset from the CMB dipole, though with large uncertainties. The significant statistical tension (3.45σ) is primarily localized to the Δ_y component, which exhibits a large deviation. Crucially, the Δ_x and Δ_z components are in good agreement with expectations (see Fig. 4). This pattern–a coherent amplitude, a directionally ambiguous vector with two well-behaved components–suggests that the tension is more likely driven by residual systematics or sample variance in DESI DR1, rather than by cosmological anisotropy. Therefore, the current measurements do not provide compelling evidence against the CP. Future, larger datasets will be essential to reduce uncertainties and clarify the origin of the Δ_y anomaly.

Data availability

The eBOSS data are available at https://data.sdss.org/sas/dr16/eboss. The DESI data are available at https://data.desi.lbl.gov/public/dr1.

Acknowledgments

This work was supported by the National SKA Program of China (Grants Nos. 2022SKA0110200 and 2022SKA0110203). This work was also supported by Xiaofeng Yang’s ZHISHAN Distinguished Professor startup funding of Henan University.

References

Appendix A: Methodological validation and additional analyses

This appendix provides supporting details for the dipole analysis presented in the main text, focusing on methodological consistency and the exploration of potential systematics.

As a primary validation step of our entire analysis pipeline—including data processing, dipole estimation, and uncertainty quantification—we applied it to the eBOSS DR16 QSO sample. This is the same dataset used in the benchmark study of (da Silveira Ferreira & Marra 2024).

We followed the identical procedure outlined in Section 3 of the main text. The sample selection, coordinate transformations, and likelihood analysis for the dipole components were performed without modification.

The resulting probability distributions for the monopole and the three Cartesian dipole components are shown in Fig. A.1 in the main text. The key outcome is that the expected value from the CMB kinematic dipole (red line) falls entirely within the 1σ posterior credible interval for every single component.

Fig. A.1. Empirical distribution plots of the eBOSS DR16 redshift dipole in its Cartesian components (m, cΔx, cΔy, cΔz). The red curve represents the dipole expected from the CMB dipole. The x-axis represents velocity (in km/s), the y-axis represents distribution density, and the shaded regions denote the 1σ confidence intervals. Notably, all components fully encompass the CMB dipole value at the 1σ ranges.

This successful replication serves two crucial purposes: 1. Methodological Robustness: It confirms that our implementation of the (da Silveira Ferreira & Marra 2024) method is correct and yields consistent results on the same data. 2. Baseline Establishment: It provides a clear benchmark for a “consistent” result within the ΛCDM framework. The contrast between this clean consistency (eBOSS DR16) and the component-level tension found in DESI DR1 (Fig. 4) sharpens the discussion on the origin of the latter.

All Tables

All Figures

Fig. 1. Redshift distributions of the QSO catalogues. Left to right: eBOSS DR14 (orange), eBOSS DR16 (green), and DESI DR1 (blue). The DR14 distribution is broad (0.4 < z < 8) but sparse at the extremes; we use the well-populated range 0.4 < z < 2.8. The DESI DR1 distribution spans 0.8 < z < 3.5; we use 0.8 < z < 3.0 to avoid edge effects. All histograms use a bin width of Δz = 0.001. Note that the y-axis scales differ between panels. The primary overlap region for the eBOSS samples is 0.8 < z < 2.2, while DESI significantly overlaps in the range 0.8 < z < 2.3.

In the text

	Fig. 2. Sky footprints of the QSO samples in Mollweide projection (HEALPix format with Nside = 64). Left to right: eBOSS DR14, eBOSS DR16, and DESI DR1. The maps show the number count per pixel for each catalogue.
In the text

Fig. 3. Assessment of NGC and SGC systematic differences via bootstrap resampling. For each QSO sample (left: eBOSS DR14; centre: eBOSS DR16; right: DESI DR1), we repeatedly (1000 times) drew bootstrap resamples of 28 000 data points from each hemisphere’s dataset, computing the mean redshift for each resample. The distributions of these resampled means are plotted as histograms with a bin width of 0.000333. The y-axis represents the count of resamples whose mean redshift falls within the corresponding bin. All distributions are plotted relative to the average redshift of the total NGC dataset (vertical line at zero). The x-axis shows the redshift difference Δz (dimensionless) relative to the total NGC mean. The offset between the NGC (centred) and SGC distributions indicates a hemisphere-level systematic discrepancy. The amplitude of this offset (∼0.008 for eBOSS, ∼0.001 for DESI) is comparable to the expected kinematic dipole signal (β = v/c ∼ 10−3), necessitating separate hemispheric analysis (Section 3.3).

In the text

Fig. 4. Empirical distribution plots of the DESI DR1 redshift dipole in its Cartesian components (m, cΔx, cΔy, cΔz). The red curve represents the dipole expected from the CMB dipole. The x-axis represents velocity (in kilometres per second), the y-axis represents distribution density, and the shaded regions denote the 1σ confidence intervals. Notably, the cΔy component shows small deviation from the CMB prediction, while the 1σ ranges of both cΔx and cΔz fully encompass the CMB dipole value.

In the text

Fig. A.1. Empirical distribution plots of the eBOSS DR16 redshift dipole in its Cartesian components (m, cΔx, cΔy, cΔz). The red curve represents the dipole expected from the CMB dipole. The x-axis represents velocity (in km/s), the y-axis represents distribution density, and the shaded regions denote the 1σ confidence intervals. Notably, all components fully encompass the CMB dipole value at the 1σ ranges.

In the text