Statistical effects of the observer’s peculiar velocity on source number counts

Charles Dalang,<sup>1,2,*</sup> Ruth Durrer,<sup>1</sup> and Fabien Lacasa<sup>1,3,4,5</sup>

<sup>1</sup>Université de Genève, Département de Physique Théorique and Center for Astroparticle Physics, 24 quai Ernest-Ansermet, CH-1211 Genève 4, Switzerland

<sup>2</sup>Queen Mary University of London, Mile End Road, London E1 4NS, United Kingdom

<sup>3</sup>Dipartimento di Fisica e Astronomia “Galileo Galilei,” Università degli Studi di Padova, via Marzolo 8, I-35131 Padova, Italy

<sup>4</sup>INFN, Sezione di Padova, via Marzolo 8, I-35131 Padova, Italy

<sup>5</sup>Université Paris-Saclay, CNRS, Institut d'Astrophysique Spatiale, 91405 Orsay, France

<figure><img src="image_1.png" /></figure>

(Received 7 October 2022; accepted 27 March 2023; published 12 May 2023)

The velocity of the Sun with respect to the cosmic microwave background (CMB) can be extracted from the CMB dipole, provided its intrinsic dipole is assumed to be small in comparison. This interpretation is consistent, within fairly large error bars, with the measurement of the correlations between neighboring CMB multipoles induced by the velocity of the observer, which effectively breaks isotropy. In contrast, the source number count dipole was reported to privilege a velocity of the observer with an amplitude that is about twice as large as the one extracted from the entirely kinematic interpretation of the CMB dipole, with error bars that indicate a more and more significant tension. In this work, we study the effect of the peculiar velocity of the observer on correlations of nearby multipoles in the source number counts. We provide an unbiased estimator for the kinematic dipole amplitude, which is proportional to the peculiar velocity of the observer and we compute the expected signal-to-noise ratio. Assuming full sky coverage, near future experiments can achieve better than 5% constraints on the velocity of the Sun with our estimator.

DOI: 10.1103/PhysRevD.107.103514

# I. INTRODUCTION

Watching the stars on a clear night sky may easily make one wonder about one's place in the Universe. From such a point of view, it seems hard not to think that one is located in a very special place. Be that as it may, the Copernican principle states that humans are not privileged observers of the Universe. Combined with evidence of statistical isotropy in the temperature associated with a black body spectrum from the cosmic microwave background (CMB) up to temperature fluctuations of the order of $10^{-5}$ [1], the Copernican principle points towards the cosmological principle, which is a corner stone of modern cosmology: The Universe is statistically homogeneous and isotropic. The implementation of these strong assumptions allows for the use of the highly symmetric Friedmann-Lemaître-Robertson-Walker (FLRW) metric and the application of perturbation theory to describe the Universe, on sufficiently large scales. This drastically simplifies calculations and has

allowed us to constrain with percent precision the six-parameter standard model of cosmology, the Λ-cold-dark-matter (ΛCDM) model [1].

In this reasoning, one important step was swept under the carpet. To recover statistical isotropy of CMB temperature fluctuations of the order of $10^{-5}$, one must boost the observer to the so-called CMB frame in which the apparently large CMB dipole, which is of order $10^{-3}$, i.e., a hundred times larger than other multipoles, vanishes [1-4]. This is called the “entirely kinematic interpretation” of the CMB dipole and is motivated by standard single field inflation. This canonical model predicts a nearly scale invariant power spectrum of primordial fluctuations of the inflaton field, at the end of a period of quasi-de Sitter expansion. In this context, there is no reason to expect the dipole to be a hundred times larger than other multipoles and one rather expects a primordial intrinsic dipole of the order of $10^{-5}$. On the other hand, one also expects a CMB dipole to be generated by the peculiar velocity of the observer, as pointed out by D. W. Sciama [5] and calculated by P. J. E. Peebles and D. T. Wilkinson [6]. One can easily be tempted to attribute the large dipole to the velocity of the observer with respect to the CMB, potentially absorbing a small intrinsic dipole, which is expected to yield a 1% correction. Absorbing the entire CMB dipole in the velocity leads to an observer velocity $||v_o|| = 369.82 \pm 0.11 \text{ km s}^{-1}$,

<footnote>*c.dalang@qmul.ac.uk</footnote>

<footnote>Published by the American Physical Society under the terms of the Creative Commons Attribution 4.0 International license. Further distribution of this work must maintain attribution to the author(s) and the published article's title, journal citation, and DOI.</footnote>
corresponding to $\beta_{\text{dip}} = \|\mathbf{v}_o\|/c = 0.00123 \pm 0.00036$, pointing towards $\hat{\beta}_{\text{dip}} = (264.021^\circ \pm 0.011^\circ, 48.253^\circ \pm 0.005^\circ)$ in galactic coordinates, where the yearly modulation of ~30 km s⁻<fnref n="1" /> of the Satellite orbiting around the Sun has been removed [1]. While this order of magnitude for the peculiar velocity of the Sun is expected, it is important to recall that an intrinsic CMB dipole is degenerate with a boost to linear order in $\beta$. This implies that a boost to the CMB frame, may not correspond to the frame in which matter fields live at rest on average, if there exists a significant intrinsic cosmic dipole [7,8]. This may potentially have dramatic consequences on the interpretation of the CMB or of supernova data. (See for example [9–12].)

There are several alternative ways to measure the velocity of the observer with respect to the rest frame of matter fields. One can check that the distribution of far enough sources<fnref n="1" /> leads to a kinematic dipole $D_{\text{kin}}$, consistent with the CMB dipole. This was pioneered by G. Ellis and J. Baldwin [15], who determined the dipole in the source number counts per unit solid angle for radio sources with a flux following a power law spectrum in frequency. Different teams reported a source number count dipole with a direction that agrees fairly well with the CMB dipole, but with an amplitude that is about twice as large as expected from the CMB dipole [16–19]. It was suggested by [14] that evolution bias may, at least partially, explain the reported tension. This was further studied by the authors of [20], which also find significant variations in the number count kinematic dipole in the presence of parameter evolution when using different quasar luminosity function models. Combining midinfrared quasars and radio sources, a $5.2\sigma$ tension between these two dipoles was reported in [17]. The authors of [21] reanalyzed the same midinfrared quasars and concluded that neither masking nor parameter evolution could fully explain the reported tension, even if this is subject to further assumptions. Let us, however, also mention that a source number count dipole performed with data from the Very Large Array Sky Survey and the Rapid Australian Square Kilometer Array Pathfinder Continuum Survey was reported to be consistent with the CMB dipole [22], although with much larger error bars. Note that a dipole measurement with 1048 supernovae type 1a from the Pantheon sample was found to point in a direction well aligned with the CMB dipole with an amplitude which is $2.4\sigma$ smaller than the CMB dipole [23]. Carefully analyzing the peculiar velocity field with different catalogs, the authors of [24] found a bulk flow velocity of about 400 km s⁻<fnref n="1" /> extending to ≲150 h⁻<fnref n="1" /> Mpc, which is difficult to accommodate within the standard cosmological paradigm. The cosmic infrared background can also be used to measure the velocity of far away galaxies relative to the observer and a signal-to-noise ratio of 50–100

is forecasted for the Euclid and Roman surveys [25]. Gravitational waves also offer promising ways to measure the kinematic dipole [26,27]. Additionally, we also point out that the degeneracy between intrinsic and kinematic dipole may be broken using the redshift dependent dipoles in future source number count experiments [28]. Fluctuations of the number counts from radio surveys such as the Square Kilometer Array may also be useful to constrain the observer’s peculiar velocity [29].

The CMB, however, appears to be consistent with itself. Indeed, the velocity of the observer which, effectively breaks statistical isotropy, induces correlation of the $l$ and $(l \pm 1)$ multipoles in the CMB [30–32]. This leads to an independent measurement of $\beta$, which was found to be $\beta = 0.00128 \pm 0.00026(\text{stat}) \pm 0.00038(\text{syst})$ in [4], consistent with $\beta_{\text{dip}}$. Slightly tighter constraints, consistent with $\beta_{\text{dip}} = 0.00123$ were obtained in the analysis presented in [33,34]. This shows that the entirely (or at least dominantly) kinematic interpretation of the CMB dipole is consistent with the correlation of neighboring multipoles, although there is still room for an intrinsic CMB dipole, which can make up a significant portion of the observed CMB dipole, without contradicting the observed $l$ and $(l \pm 1)$ correlations [35]. However, a significant bulk flow of galaxy clusters extending up to $300 \ h^{-1} \ \text{Mpc}$ was found in WMAP via the kinematic Sunyaev-Zeldovich effect [36]. Let us also mention that while a boost and an intrinsic dipole are degenerate in the CMB dipole at linear order in $\beta$, second order corrections in $\beta$ give distinct spectral distortions in the CMB monopole and quadrupole, which may allow for a measurement of the intrinsic CMB dipole in futuristic CMB spectral distortions experiments, as pointed out in [37].

In this work, we study the correlation of neighboring multipoles in the source number counts, which yield an independent crosscheck of the validity of the measurement of the kinematic dipole with source number counts. We find that the peculiar velocity of the observer induces correlations between the $l$ and $(l \pm 1)$ multipoles in the source number counts, which are absent for comoving observers in a statistically isotropic Universe. Assuming full sky coverage, we derive an unbiased estimator $\hat{D}_{\text{kin}}$ of the kinematic dipole amplitude, proportional itself to $\beta$, and compute its variance and signal-to-noise ratio. We also comment on the determination of the velocity direction and on limitations from a partial sky survey.

The paper is structured as follows: In Sec. II, we present the setup and the transformation rules of the ingredient quantities under boost. In Sec. III, we detail the computation of the boosted spherical harmonic coefficients of the source number counts for redshift surveys and redshift independent ones. In Sec. IV, we calculate the two-point correlation function for an observer boosted in an otherwise isotropic Universe. In Sec. V, we write down an estimator for the amplitude of the dipole and for the velocity of the observer with respect to the rest frame of distance sources.

<footnote n="1">¹Far enough sources means sources with observed redshift $z \ge 0.1$, such that the intrinsic dipole, predicted in $\Lambda$CDM is small in comparison to the kinematic dipole [13,14].</footnote>
We show that this estimator is unbiased and impose an upper bound for its variance in terms of measured quantities. In Sec. VI, we briefly outline the determination of the orientation of the peculiar velocity. We discuss our findings, the potential for a measurement and its limitations in Sec. VII.

Units are such that $c = 1$, bold symbols indicate three dimensional vectors, hats may indicate unit vectors when used on bold symbols or statistical estimators.

## II. THE SETUP

We consider two observers O and O', which are related by a Lorentz boost of velocity $\beta = ||v_o||/c$. Here $v_o$ is the velocity of O' with respect to O, which is aligned with their respective $\hat{z} = \hat{z}'$ axes. This choice allows for the azimuthal angles of the two observers to coincide, $\varphi = \varphi'$, despite the Lorentz transformation. In Sec. VI, we consider arbitrary directions of the peculiar velocity of the observer. Primed quantities relate to O' and quantities without primes relate to O. We assume that both of these observers live in a FLRW universe described by the line element $ds^2 = a^2(\eta)(-d\eta^2 + dr^2 + r^2(d\theta^2 + \sin^2\theta d\phi^2))$, where $a(\eta)$ is the scale factor, which only depends on conformal time $\eta$ and $r, \theta$ and $\varphi$ are comoving coordinates. Observer O is assumed to be following the Hubble flow, i.e., having zero peculiar velocity and may be called a comoving observer. The motion of observer O' with respect to O affects their measurements of time intervals, cosines of polar angles, polar angles, solid angles, frequencies, and redshift, which transform, respectively, in the following way:

$$dt' = dt\sqrt{1 - \beta^2}, \quad (1)$$

$$\cos \theta' = (\cos \theta + \beta)[1 + \beta \cos \theta]^{-1}, \quad (2)$$

$$\theta' = \theta + \delta\theta, \delta\theta = -\beta \sin\theta + \mathcal{O}(\beta^2), \quad (3)$$

$$d\Omega' = d^2\hat{n}' = \sin(\theta')d\theta'd\varphi' = d\Omega(1 - \beta^2)[1 + \beta \cos \theta]^{-2}, \quad (4)$$

$$\nu' = \nu(1 + \beta \cos \theta)\gamma, \quad (5)$$

$$1 + z' = (1 + z)(1 + \beta \cos \theta)^{-1}\gamma^{-1}, \quad (6)$$

where $z$ and $z'$ denote the redshifts of a photon observed at an angle $\theta$ (respectively, $\theta'$) with respect to $v_o$, and $\gamma = 1/\sqrt{1-\beta^2}$ is the Lorentz factor. Observer O' calls $\hat{n}'$ the direction of an incoming photon which corresponds to $\hat{n}$ for observer O. Those are related by [38]

$$\hat{n}' = \left( \frac{\hat{n} \cdot \hat{\beta} + \beta}{1 + \hat{n} \cdot \hat{\beta}} \right) \hat{\beta} + \frac{\hat{n} - (\hat{n} \cdot \hat{\beta})\hat{\beta}}{\gamma(1 + \hat{n} \cdot \hat{\beta})}. \quad (7)$$

In the next section, we apply these transformation rules to the source number count.

## III. SOURCE NUMBER COUNTS

We express the number of sources $dN$ per unit solid angle $d\Omega$ and per redshift bin $dz$ in the direction $\hat{n} = (\theta, \varphi)$ and at redshift $z$ for sources with flux $S$ (in W m$^{-2}$ Hz$^{-1}$) above a certain threshold $S_*(v_o)$ in some frequency band $[v_o, v_o + d\nu_o]$ as

$$\frac{dN}{d\Omega dz}[\hat{n}, z, S > S_*(v_o)]. \quad (8)$$

For fixed $S_*(v_o)$ and fixed redshift $z$, the real-valued function $\frac{dN}{d\Omega dz} : \mathbb{S}_2 \to \mathbb{R}$ may be expanded in complex spherical harmonics

$$\frac{dN}{d\Omega dz}[\hat{n}, z, S > S_*(v_o)] = \sum_{l=0}^{+\infty} \sum_{m=-l}^{l} a_{lm}(z) Y_{lm}(\hat{n}) \quad (9)$$

with complex valued functions $a_{lm}(z) : \mathbb{R}_+ \to \mathbb{C}$, with $l \in \mathbb{N}$ and $m \in [-l, -l+1, \dots, l-1, l]$. In case the survey lacks redshift information, we also consider the same quantity integrated over redshift, which we denote by

$$\begin{aligned} \frac{dN}{d\Omega}[\hat{n}, S > S_*(v_o)] &\equiv \int_0^{+\infty} dz \frac{dN}{d\Omega dz}[\hat{n}, z, S > S_*(v_o)], \\ &= \sum_{l=0}^{+\infty} \sum_{m=-l}^{l} a_{lm} Y_{lm}(\hat{n}), \end{aligned} \quad (10)$$

and for which the expansion coefficients are redshift independent. Using the orthogonality relations of the spherical harmonics

$$\int_{\mathbb{S}_2} d^2\hat{n} Y_{lm}(\hat{n}) Y_{l'm'}^*(\hat{n}) = \delta_{ll'}\delta_{mm'}, \quad (11)$$

one can obtain the expansion coefficients

$$a_{lm}(z) = \int_{\mathbb{S}_2} d^2\hat{n} \frac{dN}{d\Omega dz}[\hat{n}, z, S > S_*(v_o)] Y_{lm}^*(\hat{n}). \quad (12)$$

The redshift independent coefficients are obtained in a similar way:

$$a_{lm} = \int_{\mathbb{S}_2} d^2\hat{n} \frac{dN}{d\Omega}[\hat{n}, S > S_*(v_o)] Y_{lm}^*(\hat{n}). \quad (13)$$

A masked sky breaks these orthogonality relations, of course. In the following, we distinguish between an isotropic universe and a statistically isotropic universe. In an isotropic Universe, only the monopole, i.e., $a_{00}$ contributes to the sum in Eq. (9). This corresponds to the case
where the number density of sources is exactly the same in every direction. In practice, isotropy is only true in a statistical sense. The expectation value of the number of galaxies is the same in every direction. Within the standard model, gravity acts to cluster sources, such that, in a statistically isotropic Universe, there are in principle an infinite series of $a_{lm}$s that contribute. The $a_{lm}$s of higher $l$ correspond to perturbations on smaller and smaller angular scale. So far, number counts have been analyzed in angular space mainly in photometric surveys which make a 3 × 2 pt analysis of shear, number counts and their cross-correlations, see, e.g., [39,40]. A harmonic space analysis of the number counts from the Dark Energy Survey data is found, e.g., in [41,42]. Importantly, the motion of the observer in an otherwise isotropic Universe contributes a dipole to the source number count such that the number of sources per unit solid angle and redshift observed by $O'$ at redshift $z'$ in direction $\hat{n}'$ reads [14,43]

$$ \frac{dN'}{d\Omega' dz'}[\hat{n}', z', S > S_*(\nu_o)] = (1 + (\hat{\beta} \cdot \hat{n})D_{\text{kin}}(z)) \quad (14) $$

$$ \times \left\langle \frac{dN}{d\Omega dz} [z, \hat{n}, S > S_*(\nu_o)] \right\rangle + \mathcal{O}(\beta^2), \quad (15) $$

where the amplitude of the dipole is given by

$$ D_{\text{kin}}(z) = \left[ 2 + \frac{2[1-x(z)]}{r(z)\mathcal{H}(z)} + \frac{\dot{\mathcal{H}}(z)}{\mathcal{H}^2(z)} - f_{\text{evol}}(z) \right] \beta. \quad (16) $$

Here $\mathcal{H}(z) = \dot{a}(\eta)/a(\eta)$ indicates the conformal Hubble rate, a dot indicates a derivative with respect to conformal time and $r(z)$ is the background comoving distance,

$$ r(z) = \int_0^z \frac{dz}{(1+z)\mathcal{H}(z)}. \quad (17) $$

The magnification bias $x(z)$ [sometimes noted $s(z) = 2x(z)/5$] is defined as

$$ x(z) = - \frac{\partial \ln \langle \frac{dN}{d\Omega dz} [\hat{n}, z, S > S_*] \rangle}{\partial \ln S_*}, \quad (18) $$

and is sometimes defined in its integrated form $\langle \frac{dN}{d\Omega dz} [\hat{n}, z, S > S_*] \rangle \propto S_*^{-x(z)}$. It tracks the number density of sources above a given flux density threshold $S_*$. Intuitively, for positive $x(z)$, the number density of objects above the threshold $S_*$ goes to zero for large enough $S_*$ and diverges for $S_*$ going to zero. A constant $x(z)$ means that this power law does not change with redshift or, alternatively that the population distribution of fluxes is constant

in cosmic time. The evolution bias traces the time evolution of the number of sources per unit comoving volume $dV$

$$ \begin{aligned} f_{\text{evol}}(z) &\equiv \frac{1}{\mathcal{H}} \left\langle \frac{dN}{dV} [r(z), L > L_*] \right\rangle^{-1} \frac{\partial}{\partial \eta} \left\langle \frac{dN}{dV} [r(z), L > L_*] \right\rangle, \\ &= - \frac{\partial \ln \langle \frac{dN}{dV} [r(z), L > L_*] \rangle}{\partial \ln(1+z)}. \end{aligned} \quad (19) $$

Here $L$ and $L_*$ are luminosity densities (in W Hz⁻¹) corresponding to the flux densities $S$ and $S_*$, respectively. The angular brackets<fnref n="2" /> indicate an average over the two-sphere $S_2$

$$ \langle \dots \rangle = \frac{1}{4\pi} \int_{S_2} \dots d^2\hat{n}. \quad (20) $$

The difference between $\hat{n}'$ and $\hat{n}$ results in second order corrections [i.e., $\mathcal{O}(\beta^2)$] of $dN'/(d\Omega' dz')[\hat{n}', S > S_*(\nu_o)]$ in a Universe where the sources are isotropically distributed. Instead, in a Universe that is statistically isotropic, the difference between $\hat{n}$ and $\hat{n}'$ becomes first order in $\beta$. We note $\hat{n}' \hat{n} = (\theta', \phi') = (\theta + \delta\theta, \phi)$. Recall that since we have assumed that $\beta$ is aligned with $\hat{z}$, the azimuthal angle $\phi$ is unaffected by the boost. In a Universe with intrinsic anisotropies, an additional term in the number count plays a role. More precisely, Eq. (25) in [14] becomes

$$ \begin{aligned} & \frac{dN'}{d\Omega' dz'}[\hat{n}', z', S > S_*] \\ &= \frac{dN'}{d\Omega' dz'}(\theta + \delta\theta, \phi, r[z', \hat{n}], L > L'_*[z', \hat{n}, \nu_s]), \end{aligned} \quad (21) $$

$$ \begin{aligned} & \approx \frac{dN}{d\Omega dz}[\hat{n}, r[z], L > L_*] \left( \frac{d\Omega}{d\Omega'} + \frac{dz}{dz'} \right) \\ & \quad + \left. \frac{\partial}{\partial r'} \left( \frac{dN}{d\Omega dz}[\hat{n}, r', L > L_*] \right) \right|_{r'=r[z]} \cdot \delta r[z, \hat{n}] \\ & \quad + \left. \frac{\partial}{\partial L_*'} \left( \frac{dN}{d\Omega dz}[\hat{n}, r, L > L_*'] \right) \right|_{L_*'=L_*} \cdot \delta L_*[z, \hat{n}, \nu_s] \\ & \quad + \left. \frac{\partial}{\partial \theta'} \left( \frac{dN}{d\Omega dz}[\theta', \phi, r, L > L_*'] \right) \right|_{\theta'=\theta} \cdot \delta\theta[\hat{n}], \end{aligned} \quad (22) $$

where in the first line, we have rewritten $\hat{n}' \hat{n} = (\theta', \phi') = (\theta + \delta\theta, \phi)$. We associated a direction dependent comoving distance $r[z', \hat{n}] = r[z'] + \delta r[z', \hat{n}]$ to sources located at fixed observed redshift (with $\delta r[z, \hat{n}] = -\hat{n} \cdot \beta/\mathcal{H}$). We also associated a direction dependent luminosity density

<footnote n="2">²Later in the paper, angular brackets denote expectation values. The context should allow the experienced reader to break this degeneracy.</footnote>
threshold $L'_*[z', \hat{n}, \nu'_s] = L_*[z', \nu'_s] + \delta L_*[z', \hat{n}, \nu'_s]$, which corresponds to a fixed observed flux density threshold of the detector, which is independent of its motion. In the second line, we have changed the “per observed” unit solid angle $d\Omega'$ and redshift interval $dz'$ to “per background” unit solid angle $d\Omega$ and redshift $dz$, according to Eqs. (1)–(6). In the following three lines, we Taylor expand around background quantities the three variables that are affected by the boost, namely, $r$, $L_*$, and $\theta$. The only term that was not accounted for in [14] is the $\partial_{\theta'}$ derivative. The $\partial_{\theta'}$ derivative is equal to zero to first order in $\beta$ in case the number counts are independent of angular direction for observer $O$. This is equivalent to assuming that the number counts are isotropic for observer $O$ and consist strictly of a monopole, which is an assumption of [14] and is sufficient if one is only interested in the dipole generated by the monopole due to the motion of the observer. Here, however we want to study the modification of all multipoles due to the motion of the observer. This is because $\delta\theta(\hat{n}) = -\beta \sin\theta$ is linear in $\beta$ and the $\partial_{\theta'}$ derivative acting on a monopole, which, by definition, is independent of angles, vanishes. However, if there are intrinsic anisotropies, meaning if $a_{lm} \neq 0$ for $l \ge 1$,

then this derivative is nonzero. We have $\delta\theta = -\beta \sin\theta = -\tan(\theta)\beta \cdot \hat{n}$. This $\partial_{\theta'}$ derivative results in an additional dipolar term that comes from relating the number density from direction $\hat{n}'$ for the boosted observer to direction $\hat{n}$ for the observer at rest with respect to the source's rest frame. Therefore, instead of Eq. (14), in a statistically isotropic Universe, we get

$$ \frac{dN'}{d\Omega' dz'}[\hat{n}', z', S > S_*] = (1 + [\cos\theta \mathcal{D}_{\text{kin}}(z) - \beta \sin(\theta)\partial_{\theta}]) \\ \times \frac{dN}{d\Omega dz}[\hat{n}, z, S > S_*]. \quad (23) $$

Assuming that the sources have a luminosity density that follows a frequency power law $L \propto \nu_s^{-\alpha(z)}$ with spectral index $\alpha(z)$, one can integrate over $dz$, write the partial derive of $f_{\text{evol}}(z)$ in terms of a total redshift derivative, relate the luminosity density to a flux density, and integrate by parts (see Sec. II of [14] or Appendix A of [28]) to find

$$ \frac{dN'}{d\Omega'}[\hat{n}', S > S_*] = \int_0^{+\infty} dz \left( 1 + \hat{n} \cdot \beta \left[ 3 + x(z)[1 + \alpha(z)] - \tan(\theta)\partial_{\theta} + (1+z)\frac{d}{dz} \right] \right) \frac{dN}{d\Omega dz}[\hat{n}, z, S > S_*]. \quad (24) $$

Assuming for simplicity that $x(z)$ and $\alpha(z)$ are constant, we find, after integrating by parts and neglecting boundary terms<fnref n="3" /> that

$$ \frac{dN'}{d\Omega'}[\hat{n}', S > S_*] = (1 + [\cos\theta\mathcal{D}_{\text{kin}} - \beta\sin(\theta)\partial_{\theta}]) \\ \times \frac{dN}{d\Omega}[\hat{n}, S > S_*], \quad (25) $$

where the kinematic dipole boils down to the Ellis and Baldwin formula<fnref n="4" /> [15]

$$ \mathcal{D}_{\text{kin}} = [2 + x(1 + \alpha)]\beta. \quad (26) $$

Recall that a constant $x$ means that the flux distribution of sources does not depend on redshift, as we have discussed below Eq. (18). Typical values of these parameters

for radio galaxies are $x \sim 1$ and $\alpha \sim 1$, such that $\mathcal{D}_{\text{kin}} \sim 4\beta$. From now on, we focus on the redshift integrated surveys, which measure $dN'/d\Omega'$, but comparing (23) and (25), one sees that the only change in redshift dependent surveys is the change of $\mathcal{D}_{\text{kin}}$ in (25) to $\mathcal{D}_{\text{kin}}(z)$, defined in (16).

For observer $O'$, the number of sources per unit solid angle may also be expanded in spherical harmonics, although the coefficients will in general be different

$$ \frac{dN'}{d\Omega'}[\hat{n}', S > S_*(\nu_o)] = \sum_{l=0}^{+\infty} \sum_{m=-l}^{l} a'_{lm} Y_{lm}(\hat{n}') . \quad (27) $$

One computes the boosted $a'_{lm}$s by computing the following integrals

$$ a'_{lm} = \int_{\mathbb{S}_2} d^2\hat{n}' \frac{dN'}{d\Omega'}[\hat{n}', S > S_*(\nu_o)] Y_{lm}^*(\hat{n}') . \quad (28) $$

One can express these integrals as follows

<footnote n="3">³Boundary terms may not necessarily vanish. For example, if one works with redshift bins, one may have to include these boundary terms, which are straightforward to compute from Eq. (24).</footnote>

<footnote n="4">⁴Strictly speaking, it is sufficient that $\alpha(x)$ and $x(z)$ are uncorrelated to recover the Ellis and Baldwin formula.</footnote>
$$
\begin{align}
a'_{l'm'} &= \int_{S_2} d^2\hat{n}' \frac{dN'}{d\Omega'} [\hat{n}', S > S_*(\nu_o)] Y_{l'm'}^*(\hat{n}'), \nonumber \\
&= \int_{S_2} d^2\hat{n}' \left[ (1 + [\cos \theta \mathcal{D}_{\text{kin}} - \beta \sin(\theta) \partial_\theta]) \frac{dN}{d\Omega} [\hat{n}, S > S_*(\nu_o)] \right] \cdot Y_{l'm'}^*(\hat{n}'), \nonumber \\
&= \int_{S_2} d^2\hat{n}' \left[ (1 + [\cos \theta \mathcal{D}_{\text{kin}} - \beta \sin(\theta) \partial_\theta]) \sum_{lm} a_{lm} Y_{lm}(\hat{n}) \right] \cdot (-1)^{m'} Y_{l'-m'}(\hat{n}'), \tag{29}
\end{align}
$$

where we have used that $Y_{lm}^*(\hat{n}') = (-1)^m Y_{l(-m)}(\hat{n}')$ and
expanded the number counts in spherical harmonics
according to Eq. (10). We express $Y_{lm}(\hat{n})$ in terms of
the variable $\hat{n}'$ by a Taylor expansion

$$
Y_{lm}(\hat{n}) = Y_{lm}(\hat{n}') + \beta \sin \theta \partial_{\theta} Y_{lm}(\hat{n}') + \mathcal{O}(\beta^2). \quad (30)
$$

This cancels the $-\beta \sin \partial_{\theta}$ acting on the spherical harmonic $Y_{lm}(\hat{n})$ in Eq. (29) to first order in $\beta$. Using

$$
\cos \theta = 2 \sqrt{\frac{\pi}{3}} Y_{10}(\hat{n}), \quad (31)
$$

we are left with

$$
a'_{l'm'} = a_{l'm'} \tag{32}
$$

$$
+ 2\sqrt{\frac{\pi}{3}} \mathcal{D}_{\text{kin}}(-1)^{m'} \\
\times \sum_{lm} \int_{S_2} d^2 \hat{n} Y_{10}(\hat{n}) Y_{lm}(\hat{n}) Y_{l'-m'}(\hat{n}). \quad (33)
$$

The integrals involving three spherical harmonics are
Gaunt coefficients, which satisfy certain selection rules. In
particular, the three $l$s that must satisfy the triangle
condition, i.e., $|l_1 - l_2| \le l_3 \le (l_1 + l_2)$. Then they are
given by (see Appendix 4 of [44] for more details)

$$
\begin{align}
& \int_{S_2} \sin \theta d\theta d\phi Y_{l_1 m_1}(\theta, \phi) Y_{l_2 m_2}(\theta, \phi) Y_{l_3 m_3}(\theta, \phi) \nonumber \\
& = \sqrt{\frac{(2l_1 + 1)(2l_2 + 1)(2l_3 + 1)}{4\pi}} \begin{pmatrix} l_1 & l_2 & l_3 \\ 0 & 0 & 0 \end{pmatrix} \begin{pmatrix} l_1 & l_2 & l_3 \\ m_1 & m_2 & m_3 \end{pmatrix}, \tag{34}
\end{align}
$$

otherwise, these integrals vanish. The (3 × 2) matrices are 3-j symbols that are related to the Clebsch-Gordan coefficients. They are nonvanishing only if the sum $m_1 + m_2 + m_3$ vanishes and the triangle inequality between the $l_i$s is satisfied, see Ref. [44] or some text on quantum mechanics. Since in (33) there is a sum over $l$ and $m$ and since the triangle condition must hold, we have schematically

$$
\begin{align}
& \sum_{lm} \int_{S_2} \sin \theta d\theta d\phi Y_{lm}(\theta, \phi) Y_{10}(\theta, \phi) Y_{l'(-m')}(\theta, \phi) \nonumber \\
& = \sqrt{\frac{(2l' + 3)(2l' + 1)}{4\pi}} \begin{pmatrix} l' + 1 & 1 & l' \\ 0 & 0 & 0 \end{pmatrix} \begin{pmatrix} l' + 1 & 1 & l' \\ m' & 0 & -m' \end{pmatrix} \nonumber \\
& \quad + \sqrt{\frac{(2l' - 1)(2l' + 1)}{4\pi}} \begin{pmatrix} l' - 1 & 1 & l' \\ 0 & 0 & 0 \end{pmatrix} \begin{pmatrix} l' - 1 & 1 & l' \\ m' & 0 & -m' \end{pmatrix}, \tag{35}
\end{align}
$$

where we have also used another important rule for the 3-j
symbols

$$
\begin{pmatrix} l_1 & l_2 & l_3 \\ m_1 & m_2 & m_3 \end{pmatrix} = 0 \quad \text{if} \quad m_1 + m_2 \neq -m_3. \qquad (36)
$$

$$
a'_{l'm} = \sum_{l=1}^{+\infty} K_{l'lm} a_{lm}, \quad (37)
$$

with the kernel

$$
K_{l'l'm} = \delta_{ll'} + \delta_{l(l'+1)} A_{lm} \mathcal{D}_{\text{kin}} + \delta_{l(l'-1)} A_{l+1m} \mathcal{D}_{\text{kin}}, \quad (38)
$$

We can then express Eq. (33) as

and define the functions
<figure><img src="image_2.png" /><figcaption>FIG. 1. We plot the coefficients $A_{lm}$ as a function of $l \in [1, 1000]$ and $m \in [-l, l]$. The $A_{lm}$s appear as proportionality factors in front of $D_{kin}$ in the correlation of the neighboring multipoles in Eqs. (42). For $m = |l|$, they vanish. For $m < |l|$, $A_{lm} > 0$, which allows for the use of the off diagonal multipole correlations to constrain the kinematic dipole $D_{kin}$.</figcaption></figure>

$$B_{lm}(D_{kin}/\beta) \equiv A_{lm} \frac{D_{kin}}{\beta} \equiv \sqrt{\frac{l^2 - m^2}{(2l+1)(2l-1)}} \cdot \frac{D_{kin}}{\beta}. \quad (39)$$

The functions $B_{lm}(D_{kin}/\beta)$ depend linearly on the ratio $D_{kin}/\beta$, while the correlation coefficients $A_{lm}$s are independent of $\beta$ and $D_{kin}$. For $|m| = l$ these coefficients

$$\langle a'_{lm} (a'_{l'm'})^* \rangle = C_l \delta_{ll'} \delta_{mm'} + D_{kin} \cdot \delta_{mm'} ([C_l + C_{(l-1)}] A_{lm} \delta_{l(l'+1)} + [C_l + C_{(l+1)}] A_{l+1m} \delta_{l(l'-1)}), \quad (42)$$

$$= C_l \delta_{ll'} \delta_{mm'} + \beta \cdot \delta_{mm'} ([C_l + C_{(l-1)}] B_{lm}(D_{kin}/\beta) \delta_{l(l'+1)} + [C_l + C_{(l+1)}] B_{l+1m}(D_{kin}/\beta) \delta_{l(l'-1)}). \quad (43)$$

This is no longer proportional to $\delta_{ll'}\delta_{mm'}$ but correlations between neighboring $a'_{lm}$s appear. This implies that an observer that is moving with respect to the statistically isotropic perturbed universe observes a statistically anisotropic distribution of sources. These deviations from statistical isotropy are encoded in the neighboring multipoles. Intuitively, aberration squashes the perturbation in the direction of motion according to Eq. (4), which naturally leads to a preferred direction, breaking statistical isotropy. These neighboring multipole correlations involving $C_l$ and $C_{l\pm 1}$ are proportional to $D_{kin}$. This suggests that correlations of $a_{lm}$ and $a_{l\pm 1,m}$ of a sky map of number counts may be used to constrain the amplitude of $D_{kin}$ or of $\beta$, if the ratio $D_{kin}/\beta$ is known. In the following section, we lay out the procedure to estimate $D_{kin}$ or $\beta$ from an observed map $dN'/d\Omega'[\hat{n}', S > S_*(\nu_o)]$.

## V. QUADRATIC ESTIMATORS

In this section, we derive a quadratic estimator for the kinematic dipole $D_{kin}$, we check that it is unbiased, and we

vanish, while for fixed $|m| \ll l$, they are typically of order $1/2$ and they have the large $l$ behavior

$$\lim_{l \to +\infty} A_{lm} = \frac{1}{2}. \quad (40)$$

We plot $A_{lm}$ in Fig. 1 for relevant values of $l$ and $m$. As one sees in the figure, the correlation coefficients are always positive and largest contributions to the cross-correlation coefficients $A_{lm}$ come from the smaller values of $m$.

In the next section, we compute the correlation between different $a'_{lm}$s, assuming that the Universe is statistically isotropic and homogeneous for the observer $O$.

## IV. CORRELATION OF NEIGHBORING MULTIPOLES

We assume that $O$ is a comoving observer for which the Universe appears statistically isotropic. This implies that the two-point correlation function measured by $O$ satisfies

$$\langle a_{lm} a_{l'm'}^* \rangle = C_l \delta_{ll'} \delta_{mm'}. \quad (41)$$

Here and in what follows, brackets denote ensemble averages. For the boosted observer, $O'$, to linear order in $\beta$ we find, using (37),

compute its variance and signal-to-noise ratio. It is straightforward to extend this estimator to the redshift dependent $D_{kin}(z)$ by replacing $a'_{lm} \to a'_{lm}(z)$. For simplicity, we assume a catalog with full sky coverage. In practice, number count catalogs cover only a fraction of the sky, which induces the important limitations that we discuss in Sec. VII. First one needs to determine the boosted coefficients using Eq. (28). We can estimate the variance $C_l$ of the $a'_{lm}$s, with the following quadratic estimator

$$\hat{C}_l = \frac{1}{2l+1} \sum_{m=-l}^{l} |a'_{lm}|^2, \quad (44)$$

which is boost independent to linear order in $\beta$. The predictions for the correlators of $l$ and $l+1$ is

$$\langle a'_{lm} (a'_{(l+1)m})^* \rangle = [C_{l+1} + C_l] A_{l+1m} D_{kin}. \quad (45)$$

We can estimate $D_{kin}$ with all products of neighboring $a'_{lm}$s available,
$$ \hat{D}_{\text{kin}} = \frac{1}{l_{\text{max}}(l_{\text{max}} + 2)} \sum_{l=1}^{l_{\text{max}}} \sum_{m=-l}^{l} \frac{a'_{lm}(a'_{(l+1)m})^*}{D_{lm}}, \quad (46) $$

where

$$ D_{lm} \equiv (C_{l+1} + C_l)A_{l+1m}. \quad (47) $$

Similarly, the redshift dependent kinematic dipole $D_{\text{kin}}(z)$ can be estimated using Eq. (46) by replacing $a'_{lm} \rightarrow a'_{lm}(z)$ and their variance $C_l \rightarrow C_l(z)$. If the ratio $D_{\text{kin}}/\beta$, which appears in $B_{lm}(D_{\text{kin}}/\beta)$ is known, for example, from direct

number count dipole measurements (together with $x$ and $\alpha$), one can directly estimate $\beta$, using the following estimator

$$ \hat{\beta} = \frac{1}{l_{\text{max}}(l_{\text{max}} + 2)} \sum_{l=1}^{l_{\text{max}}} \sum_{m=-l}^{l} \frac{a'_{lm}(a'_{(l+1)m})^*}{(C_{l+1} + C_l)B_{l+1m}(D_{\text{kin}}/\beta)}. \quad (48) $$

One can easily check that these estimators are unbiased. For example, we compute the bias of $\hat{D}_{\text{kin}}$

$$ \begin{aligned} b_{D_{\text{kin}}}(\hat{D}_{\text{kin}}) &= \langle \hat{D}_{\text{kin}} \rangle - D_{\text{kin}}, \\ &= \left\langle \frac{1}{l_{\text{max}}(l_{\text{max}} + 2)} \sum_{l=1}^{l_{\text{max}}} \sum_{m=-l}^{l} \frac{a'_{lm}(a'_{(l+1)m})^*}{(C_{l+1} + C_l)A_{l+1m}} \right\rangle - D_{\text{kin}}, \\ &= \frac{1}{l_{\text{max}}(l_{\text{max}} + 2)} \sum_{l=1}^{l_{\text{max}}} \sum_{m=-l}^{l} \frac{\langle a'_{lm}(a'_{(l+1)m})^* \rangle}{(C_{l+1} + C_l)A_{l+1m}} - D_{\text{kin}}, \\ &= \left( \frac{1}{l_{\text{max}}(l_{\text{max}} + 2)} \sum_{l=1}^{l_{\text{max}}} \sum_{m=-l}^{l} D_{\text{kin}} \right) - D_{\text{kin}} = 0, \end{aligned} \quad (49) $$

which shows that our estimator is unbiased. Next, we compute the variance of $\hat{D}_{\text{kin}}$

$$ \mathrm{Var}(\hat{D}_{\text{kin}}) = \langle \hat{D}_{\text{kin}}^2 \rangle - [\langle \hat{D}_{\text{kin}} \rangle]^2 = \langle \hat{D}_{\text{kin}}^2 \rangle - D_{\text{kin}}^2. \quad (50) $$

The squared expectation value of $\hat{D}_{\text{kin}}$ is given by

$$ \langle \hat{D}_{\text{kin}} \rangle^2 = \left( \frac{1}{l_{\text{max}}(l_{\text{max}} + 2)} \right)^2 \sum_{l=1}^{l_{\text{max}}} \sum_{l'=1}^{l_{\text{max}}} \sum_{m=-l}^{l} \sum_{m'=-l'}^{l'} \frac{\langle a'_{lm}(a'_{(l+1)m})^* \rangle \langle a'_{l'm'}(a'_{(l'+1)m'})^* \rangle}{D_{lm}D_{l'm'}}. \quad (51) $$

The expectation value of $\hat{D}_{\text{kin}}^2$ reads

$$ \langle \hat{D}_{\text{kin}}^2 \rangle = \left( \frac{1}{l_{\text{max}}(l_{\text{max}} + 2)} \right)^2 \sum_{l=1}^{l_{\text{max}}} \sum_{l'=1}^{l_{\text{max}}} \sum_{m=-l}^{l} \sum_{m'=-l'}^{l'} \frac{\langle a'_{lm}(a'_{(l+1)m})^* a'_{l'm'}(a'_{(l'+1)m'})^* \rangle}{D_{lm}D_{l'm'}}. \quad (52) $$

We can use Isserli's theorem [45] (better known as Wick's theorem [46]) to express the expectation value of four Gaussian random variables as a sum of products of expectation values of two random variables. We have

$$ \begin{aligned} \langle a'_{lm}(a'_{(l+1)m})^* a'_{l'm'}(a'_{(l'+1)m'})^* \rangle &= \langle a'_{lm}(a'_{(l+1)m})^* \rangle \langle a'_{l'm'}(a'_{(l'+1)m'})^* \rangle + \langle a'_{lm} a'_{l'm'} \rangle \langle (a'_{(l+1)m})^* (a'_{(l'+1)m'})^* \rangle \\ &\quad + \langle a'_{lm}(a'_{(l'+1)m'})^* \rangle \langle (a'_{(l+1)m})^* a'_{l'm'} \rangle, \\ &= \langle a'_{lm}(a'_{(l+1)m})^* \rangle \langle a'_{l'm'}(a'_{(l'+1)m'})^* \rangle + (-1)^{m'} (-1)^m \langle a'_{lm}(a'_{l'(-m')})^* \rangle \langle a'_{(l+1)(-m)}(a'_{(l'+1)m'})^* \rangle \\ &\quad + \langle a'_{lm}(a'_{(l'+1)m'})^* \rangle \langle a'_{l'm'}(a'_{(l+1)m})^* \rangle, \end{aligned} \quad (53) $$

where we have used that $a'_{lm} = (-1)^m a_{l(-m)}$. The first term on the right-hand side ends up canceling with the $\langle \hat{D}_{\text{kin}} \rangle^2$ from (51) in (50). We are left with
$$ 
\begin{aligned}
\mathcal{V}ar(\hat{D}_{\text{kin}}) &= \frac{1}{l_{\text{max}}^2 (l_{\text{max}} + 2)^2} \sum_{l=1}^{l_{\text{max}}} \sum_{l'=1}^{l_{\text{max}}} \sum_{m=-l_{\text{max}}'-l''}^{l} \sum_{l'''} \left( \frac{(-1)^{m+m'}}{D_{lm} D_{l'm'}} \frac{C_l \delta_{l'l} \delta_{m'(-m)}}{C_{l+1} \delta_{(l+1)(l'+1)} \delta_{(-m)m'}} + \frac{C_l \delta_{l(l'+1)} \delta_{mm'}}{D_{lm} D_{l'm'}} \frac{C_{l'} \delta_{l'(l+1)} \delta_{m'm}}{C_{l+1} \delta_{(l+1)(l'+1)} \delta_{m'm'}} \right), \\
&= \frac{1}{l_{\text{max}}^2 (l_{\text{max}} + 2)^2} \sum_{l=1}^{l_{\text{max}}} \sum_{l'=1}^{l_{\text{max}}} \sum_{m=-l}^{l} \frac{\delta_{l'l} C_l C_{l+1}}{D_{lm} D_{l'(-m)}} = \frac{1}{l_{\text{max}}^2 (l_{\text{max}} + 2)^2} \sum_{l=1}^{l_{\text{max}}} \sum_{m=-l}^{l} \frac{C_l C_{l+1}}{D_{lm} D_{l(-m)}}.
\end{aligned}
\quad (54) $$

The signal-to-noise ratio reads

$$ S/N = \frac{\langle \hat{D}_{\text{kin}} \rangle}{\mathcal{V}ar(\hat{D}_{\text{kin}})} = \frac{D_{\text{kin}} \cdot l_{\text{max}}(l_{\text{max}} + 2)}{\sqrt{\sum_{l=1}^{l_{\text{max}}} \sum_{m=-l}^{l} \frac{C_l C_{l+1}}{D_{lm} D_{l(-m)}}}}. \quad (55) $$

Future number count experiments are expected to measure $dN'/d\Omega'[\hat{n}', S > S_*(\nu_o)]$ with an angular resolution below the arcmin scale (with, e.g., 30 galaxies/arcmin² for Euclid [47]), translating to $l_{max} \ge 10^4$. It is not relevant here that the $C_l$s cannot be calculated within linear perturbation theory. They can be extracted from the observations themselves even in the nonlinear regime. The only relevant assumptions are statistical isotropy for observer O, as we show in the Appendix and $\beta \ll 1$ so that an expansion to linear order in $\beta$ makes sense. To estimate how the signal-to-noise ratio scales with $l_{max}$, we estimate the denominator of (55) for a smooth variance of the source number counts so that we can set $C_l \approx C_{l+1}$, such that

$$ \sum_{l=1}^{l_{\text{max}}} \sum_{m=-l}^{l} \frac{C_l C_{l+1}}{D_{lm} D_{l(-m)}} \approx \sum_{l=1}^{l_{\text{max}}} \sum_{m=-l}^{l} \frac{1}{4 A_{l+1m} A_{l+1(-m)}}. \quad (56) $$

To estimate the denominator, we replace the sum over m by an integral and use $A_{l(-m)} = A_{lm}$

$$ \begin{aligned} \sum_{m=-l}^{l} \frac{1}{4[A_{l+1m}]^2} &= \sum_{m=-l}^{l} \frac{(2l+3)(2l+1)}{4[(l+1)^2 - m^2]}, \\ &\approx \frac{(2l+3)(2l+1)}{4} \int_{-l}^{l} \frac{dm}{[(l+1)^2 - m^2]}, \\ &= \frac{(2l+3)(2l+1) \log(2l+1)}{4(l+1)}. \end{aligned} \quad (57) $$

We now integrate this result over l to find

$$ \sum_{l=1}^{l_{\text{max}}} \sum_{m=-l}^{l} \frac{1}{4[A_{l+1m}]^2} \approx \int_{1}^{l_{\text{max}}} dl \frac{(2l+3)(2l+1)\log(2l+1)}{4(l+1)}, \quad (58) $$

$$ = \frac{1}{8} \left[ \begin{aligned} & (4l_{\text{max}}(l_{\text{max}} + 2) - 2\log(2l_{\text{max}} + 2))\log(2l_{\text{max}} + 1) \\ & + 2\text{Li}_2(-3) - 2\text{Li}_2(-2l_{\text{max}} - 1) + 8 - 2l_{\text{max}}(l_{\text{max}} + 3) \\ & - 15\log(3) + \log(4)\log(9) \end{aligned} \right], \quad (59) $$

$$ \le 6l_{\text{max}}^2, \quad (60) $$

where the validity of the last inequality can be checked numerically for $l_{max} \le 10^5$. Here $Li_2$ denotes the dilogarithm given by

$$ \text{Li}_2(z) = \int_z^0 \frac{\log(1-t)}{t} dt = \sum_{k=1}^{\infty} \frac{z^k}{k^2}. \quad (61) $$

With (60), we can set

$$ S/N \ge \frac{D_{\text{kin}} \cdot l_{\text{max}}}{\sqrt{6}}, \quad (62) $$

which is larger than 1 for $l_{max} \ge 497$, assuming $\hat{\beta}$ converges to $\beta_{dip}$ and $D_{kin} = 4\beta_{dip}$. Note that this is rather a conservative estimate, since $D_{kin}$ is rather observed to be twice as large as $4\beta_{dip}$ [16–18]. For future experiments like Euclid, which can measure the clustering of photometric sources up to $l_{max} = 10^4$, we can hope to constrain $D_{kin}$ (or $\beta$ if the ratio $D_{kin}/\beta$ is known) with less than 5% error. As stated earlier, it is irrelevant that the $C_l$s cannot be calculated within linear perturbation theory, as they can be extracted from the observations themselves even in the nonlinear regime. This rather suggests a clever use of the $C_l$s beyond the linear regime.

## VI. ORIENTATION

In our treatment so far, we have assumed that the direction of the peculiar velocity $\beta$ is known and we have chosen the $\hat{z}$ axis in its direction. This is especially relevant, since most radio surveys agree relatively well with the direction of the CMB dipole, but they find a much too large amplitude for the velocity. However, in general we want to determine both the amplitude and the direction of $\beta$. This can be achieved easily, remembering the transformation of the coefficients $a_{lm}$ under rotation.

We start from Eq. (37) which relates the $a_{lm}$s of the comoving observer to the $a'_{lm}$s of the observer boosted along the $\hat{z}$ axis. Let us assume that the velocity $\beta$ is not along the $\hat{z}$ axis, but along a direction which is rotated with respect to $\hat{z}$ by a rotation $R \in SO(3)$. If we rotate the coordinate system by $R^{-1} = R^T$, then $\beta$ points in the $\hat{z}$ direction with respect to the new coordinate system. In this rotated system, Eq. (37) is valid. Under a rotation by $R^{-1}$, the $a_{lm}$s transform as
$$a_{lm}^{\text{rot}} = \sum_{m'=-l}^{l} D_{mm'}^{(l)}(R) a_{lm'}. \qquad (63)$$

Here $D_{mm'}^{(l)}(R)$ is the representation matrix of the representation $D^{(l)}$ of SO(3), see Ref. [44], Appendix 4, for details. In the rotated system, we have

$$\begin{aligned} a_{l'm}^{\text{rot}} &= \sum_{m''=-l'}^{l'} D_{mm''}^{(l'')} (R) a_{l'm''} = \sum_{l} K_{l'l'm} a_{lm}^{\text{rot}}, \\ &= \sum_{lm'} K_{l'l'm} D_{mm'}^{(l)} (R) a_{lm'}. \end{aligned} \qquad (64)$$

Hence, the general relation for $\beta$ rotated by $R$ with respect to the $\hat{z}$ axis is

$$a_{l'm} = \sum_{lm'} K_{l'lmm'}^{\text{rot}} a_{lm'}, \qquad (65)$$

$$\begin{aligned} K_{l'lmm'}^{\text{rot}} &= \sum_{m''} (D_{mm''}^{(l'')} (R))^* K_{l'lmm''} D_{m''m'}^{(l)} (R), \qquad (66) \\ &= \delta_{ll'} \delta_{mm'} + \mathcal{D}_{\text{kin}} \left( \delta_{l(l'+1)} \sum_{m''} (D_{mm''}^{(l'')} (R))^* A_{lm''} D_{m''m'}^{(l)} (R) \right. \\ & \qquad \left. + \delta_{l(l'-1)} \sum_{m''} (D_{mm''}^{(l'')} (R))^* A_{(l+1)m''} D_{m''m'}^{(l)} (R) \right). \end{aligned} \qquad (67)$$

Here, we have used that the $D^{(l)}$ are unitary representations, $D_{mm''}^{(l)}(R^{-1}) = (D_{mm''}^{(l)})^*$.

In general, a rotation is given by three Euler angles. Since $R$ is the rotation that turns $\hat{\beta}$ with its polar angles $(\theta_{\beta}, \varphi_{\beta})$ into the $\hat{z}$ axis, we can simply choose the Euler angles $(0, -\theta_{\beta}, -\varphi_{\beta})$ and insert

$$D_{m''m'}^{(l)}(R) = D_{m''m'}^{(l)}(0, -\theta_{\beta}, -\varphi_{\beta}) = \sqrt{\frac{4\pi}{2l+1-m'}} Y_{lm''}(\theta_{\beta}, \varphi_{\beta}). \qquad (68)$$

Here $sY_{lm}$ is the spin-weighted spherical harmonic of spin $s$ (see Ref. [44], Appendix 4 for details). In addition to the amplitude dependence given in Eqs. (38) to (39), the coefficients $K_{l'lmm'}^{\text{rot}}$ now also depend on the orientation $(\theta_{\beta}, \varphi_{\beta})$. Denoting

$$A^{\text{rot}}(l, m, m', \theta_{\beta}, \varphi_{\beta}) = \sum_{m''} D_{mm''}^{(l-1)*}(0, -\theta_{\beta}, -\varphi_{\beta}) \times A_{lm''} D_{m''m'}^{(l)}(0, -\theta_{\beta}, -\varphi_{\beta}), \qquad (69)$$

we have

$$\begin{aligned} K_{l'lmm'}^{\text{rot}}(\mathcal{D}_{\text{kin}}, \theta_{\beta}, \varphi_{\beta}) &= \delta_{ll'}\delta_{mm'} + \delta_{l(l'+1)}A^{\text{rot}}(l, m, m', \theta_{\beta}, \varphi_{\beta}) \\ & \quad \times \mathcal{D}_{\text{kin}} + \delta_{l(l'-1)}A^{\text{rot}}(l+1, m, m', \theta_{\beta}, \varphi_{\beta})\mathcal{D}_{\text{kin}}. \end{aligned} \qquad (70)$$

Note that under rotations $A^{\text{rot}}(l, \dots)$ transforms like a rank $l-1$ tensor from the left and a rank $l$ tensor from the right or vice versa. Note that since $A(l, m, m', 0, 0) \propto \delta_{mm'}$, an action from the left or from the right cannot be distinguished. For $\theta_{\beta} = \varphi_{\beta} = 0$, the tensor $A^{\text{rot}}$ is diagonal,

$$A^{\text{rot}}(l, m, m', 0, 0) = \delta_{mm'} A_{lm}. \qquad (71)$$

The product $a_{lm}a_{(l+1)m'}^*$ actually yields an unbiased estimator for $A^{\text{rot}}(l, m, m', \theta_{\beta}, \varphi_{\beta}) \cdot \mathcal{D}_{\text{kin}}$. To extract the three unknowns $(\mathcal{D}_{\text{kin}}, \theta_{\beta}, \varphi_{\beta})$, which allow us to reconstruct the vector $\mathcal{D}_{\text{kin}} = \mathcal{D}_{\text{kin}} \cdot \hat{\beta}$, one computes for each $l, m$, and $m'$, the products $a_{lm}a_{(l+1)m'}^*$ from observations in some basis. This three-dimensional array can be fitted by $A^{\text{rot}}(l, m, m', \theta_{\beta}, \varphi_{\beta}) \cdot \mathcal{D}_{\text{kin}}$ with the three parameters $\mathcal{D}_{\text{kin}}, \theta_{\beta}$, and $\varphi_{\beta}$, which determine the kinematic dipole due to the observer's peculiar velocity. While the amplitude $\mathcal{D}_{\text{kin}}$ enters only linearly, the dependence on the angles, $\theta_{\beta}$ and $\varphi_{\beta}$, is more complicated. The three parameters can be extracted, e.g., via a Markov chain Monte Carlo fitting procedure, assuming the $C_{ls}$ to be known and inserting the theoretical expressions for $A_{lm}$, which are independent of $\mathcal{D}_{\text{kin}}$.

## VII. DISCUSSION

In this paper, we have worked out the correlation of neighboring multipoles in the number count spherical harmonic coefficients due to a boost with velocity $||\hat{\beta}|| \ll 1$, in an otherwise isotropic Universe. We have found that these correlations are proportional to the kinematic dipole $\mathcal{D}_{\text{kin}}$. We derived an unbiased estimator $\hat{\mathcal{D}}_{\text{kin}}$ of $\mathcal{D}_{\text{kin}}$ in terms of the observed (boosted) coefficients $a'_{lm}$, $l \in \mathbb{N}$, $m \in [-l, \dots, l]$. The same estimator with redshift dependent coefficients can be used to measure $\mathcal{D}_{\text{kin}}(z)$, if the sources are arranged in redshift bins. Of course, the statistics for each redshift bin then decline. We computed the variance of the estimator $\hat{\mathcal{D}}_{\text{kin}}$ and have shown that for reasonably smooth variance of the source number counts, the signal-to-noise ratio scales (up to log $l_{\text{max}}$ corrections) as $\mathcal{O}(l_{\text{max}}\mathcal{D}_{\text{kin}})$, which becomes larger than 1 for $l_{\text{max}} \ge 497$ for the expected peculiar velocity $\beta_{\text{dip}}$ and kinematic dipole $\mathcal{D}_{\text{kin}} \sim 4\beta_{\text{dip}}$. This implies that order 5% precision may be achieved on $\beta$ if $l_{\text{max}} \ge 10^4$. Note that this is a conservative estimate since the kinematic dipole is rather observed to be twice larger than $4\beta_{\text{dip}}$ [16–18]. We have also derived an unbiased estimator $\hat{\beta}$ of $\beta$, which requires the additional knowledge of the ratio $\mathcal{D}_{\text{kin}}/\beta$. The latter can be estimated via Eq. (16) by assuming a cosmological model and by measuring $x(z)$. Alternatively, assuming that the redshift
evolutions of $\alpha$ and $x$ are uncorrelated (see Ref. [14]), one can measure them from the data and use Eq. (26) to obtain $D_{\text{kin}}$. A last option if the redshift evolutions of $\alpha$ and $x$ are of concern, is to build a parametric model for them and fit them together with $\beta$ on tomographic measurements of $D_{\text{kin}}$ with a sufficient number of redshift bins. We studied the effects of the orientation of $\hat{\beta}$ on the $l$ and $l \pm 1$ correlations and lay out the procedure that allows to determine also the orientation of $\hat{\beta}$ from observations.

While our results are optimistically set in full sky, they will be affected by an incomplete sky coverage, which also breaks statistical isotropy. In practice, number count catalogs cover only a fraction of the sky, the rest being masked for a number of reasons that include but are not limited to the footprint of the detector, the Milky Way, obscuration by dust, turbulence from the atmosphere, bad detectors, and instability of noise. One might then worry that by naively using the estimator given in (46), one may find deviations from statistical isotropy, which actually come from the mask rather than from the peculiar velocity of the observer. One way to handle this is to account for the mask by multiplying the underlying number counts by a mask function $W(\hat{n})$ that depends on the direction and on redshift and varies between zero and one, depending on the completeness of the survey. This multiplication by the mask function makes it a convolution in harmonic space that, in principle, allows us to disentangle neighboring multipole correlations generated by the mask from the ones, generated by the peculiar velocity of the observer. Masking in harmonic space has been successfully carried out in the analysis of microwave background fluctuations and of galaxy number counts [41,42].

Future work is needed to incorporate the effect of a survey mask on the estimator given in Eq. (46) and the impact on the estimator's bias, variance, and signal-to-noise ratio. The current work lays the foundation to use source number counts as an assessment of how we observe the Universe moving with our special velocity, if not from a special place.

## ACKNOWLEDGMENTS

We thank Camille Bonvin, Louis Legrand, Nastassia Grimm, Blake Sherwin, and specially Azadeh Moradinezhad Dizgah for useful conversations. This work is supported by the Swiss National Science Foundation. F. L. acknowledges support from the Swiss National Science Foundation through Grant No. IZSEZ0_207059, and from Progetto di Eccellenza 2022 of the Physics and Astronomy department at the University of Padova “Physics of the

Universe,” within the grant “CosmoGraLSS—Cosmology with Gravitational waves and Large Scale Structure.”

## APPENDIX: STATISTICAL ISOTROPY

In this appendix, we show for completeness that Eq. (41) relies only on statistical isotropy. In particular, it does not rely on the validity of perturbation theory or on Gaussianity of the fluctuations. Consider an observable, which is a real valued function, defined on the sky, i.e., $O: S_2 \to \mathbb{R}$. This can be the number count density or the temperature fluctuations of the CMB. Statistical isotropy implies that the correlation function depends only on the angle $\theta$ between $\hat{n}$ and $\hat{n}'$,

$$ \langle O(\hat{n})O(\hat{n}') \rangle = C(\mu), \quad (A1) $$

where $\mu = \cos(\theta) = \hat{n} \cdot \hat{n}'$. This function can be expanded on the interval $-1 \le \mu \le 1$ in Legendre polynomials

$$ C(\mu) = \sum_{l=0}^{+\infty} a_l P_l(\mu) = 4\pi \sum_{l=0}^{+\infty} \frac{a_l}{2l+1} \sum_{m=-l}^{l} Y_{lm}(\hat{n}) Y_{lm}^{*}(\hat{n}'), \quad (A2) $$

where we have used the addition theorem for spherical harmonics in the second equality. On the other hand, the function on the sphere $O(\hat{n})$ can be expanded in spherical harmonics

$$ O(\hat{n}) = \sum_{l=0}^{+\infty} \sum_{l'=0}^{l} a_{lm} Y_{lm}(\hat{n}) = \sum_{l=0}^{+\infty} \sum_{l'=0}^{l} a_{lm}^{*} Y_{lm}^{*}(\hat{n}), \quad (A3) $$

where the second equality holds because $O(\hat{n})$ is real. Plugging this expansion in Eq. (A1), one obtains

$$ C(\mu) = \sum_{l=0}^{+\infty} \sum_{l'=0}^{+\infty} \sum_{m=-l}^{l} \sum_{m'=-l'}^{l'} \langle a_{lm} a_{l'm'} \rangle Y_{lm}(\hat{n}) Y_{l'm'}(\hat{n}'). \quad (A4) $$

As the functions $Y_{lm}(\hat{n})Y_{l'm'}(\hat{n}')$ form an orthonormal basis on the Hilbert space of square integrable function on the cross product of two-spheres, $L^2(S_2 \times S_2)$, one can identify the coefficients in (A2) and (A4) to conclude that

$$ \langle a_{lm} a_{l'm'}^* \rangle = \delta_{ll'} \delta_{mm'} \frac{4\pi a_l}{2l+1} \equiv \delta_{l'l} \delta_{mm'} C_l. \quad (A5) $$
[1] N. Aghanim et al., Planck 2018 results. VI. Cosmological parameters, Astron. Astrophys. **641**, A6 (2020).

[2] D. J. Fixsen, E. S. Cheng, D. A. Cottingham, Jr. Eplee, R. E., R. B. Isaacman, J. C. Mather, S. S. Meyer, P. D. Noerdlinger, R. A. Shafer, R. Weiss, E. L. Wright, C. L. Bennett, N. W. Boggess, T. Kelsall, S. H. Moseley, R. F. Silverberg, G. F. Smoot, and D. T. Wilkinson, Cosmic microwave background dipole spectrum measured by the COBE FIRAS instrument, Astrophys. J. **420**, 445 (1994).

[3] D. J. Fixsen, E. S. Cheng, J. M. Gales, John C. Mather, R. A. Shafer, and E. L. Wright, The cosmic microwave background spectrum from the full COBE FIRAS data set, Astrophys. J. **473**, 576 (1996).

[4] N. Aghanim et al., Planck 2013 results. XXVII. Doppler boosting of the CMB: Eppur si muove, Astron. Astrophys. **571**, A27 (2014).

[5] D. W. Sciama, Peculiar Velocity of the Sun and the Cosmic Microwave Background, Phys. Rev. Lett. **18**, 1065 (1967).

[6] P. J. E. Peebles and David T. Wilkinson, Comment on the anisotropy of the primeval fireball, Phys. Rev. **174**, 2168 (1968).

[7] Guillem Domènech, Roya Mohayaee, Subodh P. Patil, and Subir Sarkar, Galaxy number-count dipole and superhorizon fluctuations, J. Cosmol. Astropart. Phys. **10** (2022) 019.

[8] Ivan Agullo, Dimitrios Kranas, and V. Sreenath, Anomalies in the CMB from a cosmic bounce, Gen. Relativ. Gravit. **53**, 17 (2021).

[9] Peter E. Freeman, C. R. Genovese, C. J. Miller, R. C. Nichol, and L. Wasserman, Examining the effect of the map-making algorithm on observed power asymmetry in wmap data, Astrophys. J. **638**, 1 (2006).

[10] P. Naselsky, W. Zhao, J. Kim, and S. Chen, Is the cosmic microwave background asymmetry due to the kinematic dipole?, Astrophys. J. **749**, 31 (2012).

[11] Jacques Colin, Roya Mohayaee, Mohamed Rameez, and Subir Sarkar, Evidence for anisotropy of cosmic acceleration, Astron. Astrophys. **631**, L13 (2019).

[12] Roya Mohayaee, Mohamed Rameez, and Subir Sarkar, Do supernovae indicate an accelerating universe?, Eur. Phys. J. Spec. Top. **230**, 2067 (2021).

[13] Prabhakar Tiwari and Adi Nusser, Revisiting the NVSS number count dipole, J. Cosmol. Astropart. Phys. **03** (2016) 062.

[14] Charles Dalang and Camille Bonvin, On the kinematic cosmic dipole tension, Mon. Not. R. Astron. Soc. **512**, 3895 (2022).

[15] G. F. R. Ellis and J. E. Baldwin, On the expected anisotropy of radio source counts, Mon. Not. R. Astron. Soc. **206**, 377 (1984).

[16] Nathan J. Secrest, Sebastian von Hausegger, Mohamed Rameez, Roya Mohayaee, Subir Sarkar, and Jacques Colin, A test of the cosmological principle with quasars, Astrophys. J. **908**, L51 (2021).

[17] Nathan Secrest, Sebastian von Hausegger, Mohamed Rameez, Roya Mohayaee, and Subir Sarkar, A Challenge to the Standard Cosmological Model, Astrophys. J. Lett. **937**, L31 (2022).

[18] C. A. P. Bengaly, R. Maartens, and M. G. Santos, Probing the cosmological principle in the counts of radio galaxies

at different frequencies, J. Cosmol. Astropart. Phys. **04** (2018) 031.

[19] T. M. Siewert, M. Schmidt-Rubart, and D. J. Schwarz, Cosmic radio dipole: Estimators and frequency dependence, Astron. Astrophys. **653**, A9 (2021).

[20] Caroline Guandalin, Jade Piat, Chris Clarkson, and Roy Maartens, Theoretical systematics in testing the Cosmological Principle with the kinematic quasar dipole, arXiv:2212.04925.

[21] Lawrence Dam, Geraint F. Lewis, and Brendon J. Brewer, Testing the Cosmological Principle with CatWISE Quasars: A Bayesian Analysis of the Number-Count Dipole, arXiv:2212.07733.

[22] Jeremy Darling, The universe is brighter in the direction of our motion: Galaxy counts and fluxes are consistent with the CMB dipole, Astrophys. J. Lett. **931**, L14 (2022).

[23] Nick Horstmann, Yannic Pietschke, and Dominik J. Schwarz, Inference of the cosmic rest-frame from supernovae Ia, Astron. Astrophys. **668**, A34 (2022).

[24] Yin-Zhe Ma, Christopher Gordon, and Hume A. Feldman, Peculiar velocity field: Constraining the tilt of the Universe, Phys. Rev. D **83**, 103002 (2011).

[25] A. Kashlinsky and F. Atrio-Barandela, Probing the rest-frame of the Universe with the near-IR cosmic infrared background, Mon. Not. R. Astron. Soc. **515**, L11 (2022).

[26] S. Mastrogiovanni, C. Bonvin, G. Cusin, and S. Foffa, Detection and estimation of the cosmic dipole with the Einstein Telescope and Cosmic Explorer, Mon. Not. R. Astron. Soc. **521**, 984 (2023).

[27] Adrian Ka-Wai Chung, Alexander C. Jenkins, Joseph D. Romano, and Mairi Sakellariadou, Targeted search for the kinematic dipole of the gravitational-wave background, Phys. Rev. D **106**, 082005 (2022).

[28] Tobias Nadolny, Ruth Durrer, Martin Kunz, and Hamsa Padmanabhan, A new way to test the cosmological principle: Measuring our peculiar velocity and the large-scale anisotropy independently, J. Cosmol. Astropart. Phys. **11** (2021) 009.

[29] Nidhi Pant, Aditya Rotti, Carlos A. P. Bengaly, and Roy Maartens, Measuring our velocity from fluctuations in number counts, J. Cosmol. Astropart. Phys. **03** (2019) 023.

[30] Anthony Challinor and Floor van Leeuwen, Peculiar velocity effects in high resolution microwave background experiments, Phys. Rev. D **65**, 103001 (2002).

[31] Liang Dai and Jens Chluba, New operator approach to the CMB aberration kernels in harmonic space, Phys. Rev. D **89**, 123504 (2014).

[32] J. Chluba, Aberrating the CMB sky: Fast and accurate computation of the aberration kernel, Mon. Not. R. Astron. Soc. **415**, 3227 (2011).

[33] Sayan Saha, Shabbir Shaikh, Suvodip Mukherjee, Tarun Souradeep, and Benjamin D. Wandelt, Bayesian estimation of our local motion from the Planck-2018 CMB temperature map, J. Cosmol. Astropart. Phys. **10** (2021) 072.

[34] Pedro da Silveira Ferreira and Miguel Quartin, First Constraints on the Intrinsic CMB Dipole and Our Velocity with Doppler and Aberration, Phys. Rev. Lett. **127**, 101301 (2021).
[35] Dominik J. Schwarz, Craig J. Copi, Dragan Huterer, and Glenn D. Starkman, CMB anomalies after planck, Classical Quantum Gravity **33**, 184001 (2016).

[36] A. Kashlinsky, F. Atrio-Barandela, D. Kocevski, and H. Ebeling, A measurement of large-scale peculiar velocities of clusters of galaxies: results and cosmological implications, Astrophys. J. Lett. **686**, L49 (2009).

[37] Siavash Yasini and Elena Pierpaoli, Beyond the Boost: Measuring the Intrinsic Dipole of the Cosmic Microwave Background Using the Spectral Distortions of the Monopole and Quadrupole, Phys. Rev. Lett. **119**, 221102 (2017).

[38] John David Jackson, *Classical Electrodynamics* (Wiley & Sons, New York, 1975).

[39] J. DeRose et al., Dark energy survey year 3 results: Cosmology from combined galaxy clustering and lensing validation on cosmological simulations, Phys. Rev. D **105**, 123520 (2022).

[40] Tilman Tröster et al., KiDS-1000 cosmology: Constraints beyond flat ΛCDM, Astron. Astrophys. **649**, A88 (2021).

[41] H. Camacho et al., Dark energy survey year 1 results: Measurement of the galaxy angular power spectrum, Mon. Not. R. Astron. Soc. **487**, 3870 (2019).

[42] F. Andrade-Oliveira et al., Galaxy clustering in harmonic space from the dark energy survey year 1 data: Compatibility with real-space results, Mon. Not. R. Astron. Soc. **505**, 5714 (2021).

[43] Roy Maartens, Chris Clarkson, and Song Chen, The kinematic dipole in galaxy redshift surveys, J. Cosmol. Astropart. Phys. **01** (2018) 013.

[44] Ruth Durrer, *The Cosmic Microwave Background* (Cambridge University Press, Cambridge, England, 2020).

[45] L. Isserlis, On a formula for the product.moment coefficient of any order of a normal frequency distribution in any number of variables, Biometrika **12**, 134 (1918).

[46] G. C. Wick, The evaluation of the collision matrix, Phys. Rev. **80**, 268 (1950).

[47] R. Laureijs, J. Amiaux, S. Arduini, J. L. Auguères, J. Brinchmann, R. Cole, M. Cropper, C. Dabin et al., Euclid definition study report, arXiv:1110.3193.