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

# Fast and spurious: a robust determination of our peculiar velocity with future galaxy surveys

Fabien Lacasa, Camille Bonvin, Charles Dalang, Ruth Durrer

► To cite this version:

Fabien Lacasa, Camille Bonvin, Charles Dalang, Ruth Durrer. Fast and spurious: a robust determination of our peculiar velocity with future galaxy surveys. Journal of Cosmology and Astroparticle Physics, 2024, 06, pp.045. 10.1088/1475-7516/2024/06/045. hal-04523232

HAL Id: hal-04523232

https://hal.science/hal-04523232v1

Submitted on 18 Jul 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>
PAPER • OPEN ACCESS

# Fast and spurious: a robust determination of our peculiar velocity with future galaxy surveys

To cite this article: Fabien Lacasa et al JCAP06(2024)045

View the [article online](#) for updates and enhancements.

You may also like

- Real-time tracking of source movement by Cherenkov emission imaging for high-dose-rate brachytherapy

K. Yogo, Y. Tatsuno, T. Shimo et al.

- Rapid State Transitions in the Galactic Black Hole Candidate Source GRS 1915+105

A. R. Rao, J. S. Yadav and B. Paul

- On stability of 1D exact inverse methods

P M Carrion
RECEIVED: March 13, 2024

ACCEPTED: May 27, 2024

PUBLISHED: June 19, 2024

# Fast and spurious: a robust determination of our peculiar velocity with future galaxy surveys

Fabien Lacasa ID, a,b Camille Bonvin ID, a Charles Dalang ID, c,d and Ruth Durrer ID a

E-mail: fabien.lacasa@unige.ch, camille.bonvin@unige.ch, c.dalang@qmul.ac.uk, ruth.durrer@unige.ch

**ABSTRACT:** To date, the most precise measurement of the observer's peculiar velocity comes from the dipole in the Cosmic Microwave Background (CMB). This velocity also generates a dipole in the source number counts, whose amplitude is governed not only by the observer velocity, but also by specific properties of the sources, that are difficult to determine precisely. Quantitative studies of the source number counts currently give dipoles which are reasonably well aligned with the CMB dipole, but with a significantly larger amplitude than that of the CMB dipole. In this work, we explore an alternative way of measuring the observer velocity from the source number counts, using correlations between neighboring spherical harmonic coefficients, induced by the velocity. We show that these correlations contain both a term sensitive to the source properties and another one directly given by the observer velocity. We explore the potential of a *Euclid*-like survey to directly measure this second contribution, independently of the characteristics of the population of sources. We find that the method can reach a precision of 4%, corresponding to a detection significance of 24σ, on the observer velocity. This will settle with precision the present “dipole tension”.

**KEYWORDS:** galaxy clustering, cosmological parameters from LSS, redshift surveys

ARXIV EPRINT: 2402.18438

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

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

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

<footnote><sup>c</sup>Queen Mary University of London, Mile End Road, E1 4NS, London, U.K.</footnote>

<footnote><sup>d</sup>Institute of Cosmology and Gravitation, University of Portsmouth, Burnaby Road, Portsmouth PO1 3FX, U.K.</footnote>
# Contents

<table><tr><td>1</td><td>Introduction</td><td>1</td></tr><tr><td>2</td><td>Effect of the observer velocity on the perturbed galaxy number counts</td><td>4</td></tr><tr><td></td><td>2.1 Boosted galaxy number counts</td><td>4</td></tr><tr><td></td><td>2.2 Spherical harmonic coefficients</td><td>8</td></tr><tr><td>3</td><td>Impact on observable statistics</td><td>10</td></tr><tr><td></td><td>3.1 Cosmic dipole</td><td>10</td></tr><tr><td></td><td>3.2 Correlation between multipoles</td><td>11</td></tr><tr><td>4</td><td>Estimators for β and F<sub>kin</sub></td><td>12</td></tr><tr><td></td><td>4.1 Single tracer case</td><td>12</td></tr><tr><td></td><td>4.2 Simple examples for the estimator and the signal</td><td>15</td></tr><tr><td></td><td>4.3 Generalisation to multi-tracer or multiple redshift bins</td><td>17</td></tr><tr><td>5</td><td>Forecast for future surveys</td><td>19</td></tr><tr><td></td><td>5.1 Individual redshift bins</td><td>20</td></tr><tr><td></td><td>5.2 Combining 10 redshift bins</td><td>21</td></tr><tr><td></td><td>5.3 Adding priors on F<sub>kin</sub></td><td>23</td></tr><tr><td></td><td>5.4 Adding priors on the magnification bias and evolution bias</td><td>23</td></tr><tr><td>6</td><td>Discussion and conclusions</td><td>24</td></tr><tr><td>A</td><td>Auxiliary calculations on spherical harmonics</td><td>27</td></tr><tr><td></td><td>A.1 Dipole times spherical harmonic</td><td>27</td></tr><tr><td></td><td>A.2 Derivative of spherical harmonic</td><td>29</td></tr><tr><td>B</td><td>Derivation of the quadratic estimators</td><td>30</td></tr><tr><td></td><td>B.1 Single parameter case</td><td>30</td></tr><tr><td></td><td>B.2 Two parameters case</td><td>31</td></tr><tr><td></td><td>B.3 Multi parameter case</td><td>32</td></tr><tr><td>C</td><td>Quadratic estimator and partial sky coverage</td><td>34</td></tr><tr><td>D</td><td>Revising assumptions: Gaussian covariance and known power spectrum</td><td>36</td></tr><tr><td></td><td>D.1 Gaussian covariance</td><td>36</td></tr><tr><td></td><td>D.2 Known power spectrum C<sub>ℓ</sub></td><td>36</td></tr><tr><td></td><td>D.3 Non-linear growth and bias</td><td>39</td></tr></table>

## 1 Introduction

The dipole of the Cosmic Microwave Background (CMB) has been measured with outstanding precision. Since its amplitude is about 100 times larger than all higher multipoles, it is
usually interpreted as due to our peculiar motion with respect to the cosmic rest frame. Such a motion does indeed generate a dipole $D_{\text{kin}}^{\text{CMB}} = v_o/c$ at linear order in the velocity $v_o$ [1]. Attributing the observed dipole entirely to our motion with respect to the cosmic rest frame, one infers a velocity amplitude and direction of [2, 3]

$$|v_0| = (369.82 \pm 0.11) \text{ km/s}, \quad (l, b) = (264.021^{\circ} \pm 0.011^{\circ}, 48.253^{\circ} \pm 0.005^{\circ}), \quad (1.1)$$

where $(l, b)$ denote the longitude and latitude of the velocity direction in galactic coordinates. Assuming that the kinematic dipole is contaminated by an intrinsic dipole of a similar amplitude as the higher multipoles induces an error of about 1% in the above determination of $v_o$.

In addition to the dipole in the CMB temperature, our motion is expected to generate a dipole in the distribution of cosmological sources around us. In recent years, several groups were able to detect this dipole in cosmic sources, observing radio galaxies [4–9] and, more recently, quasars [10–12]. Neglecting the redshift evolution of the sources, the kinematic dipole of radio galaxies and quasars, averaged over all redshifts (which is denoted here only, with a bar), can be written as [13]

$$\bar{D}_{\text{kin}} = [2 + x(1 + \alpha)] v_o/c, \quad (1.2)$$

where $x = -d \ln N(S > S_*)/d \ln S_*$ is the slope of the cumulative number of sources at the limiting flux of the survey $S_*$, and $\alpha = -d \ln S/d \ln f|_{f_o}$ is the spectral index of the sources evaluated at the observed frequency $f_o$. Eq. (1.2) is the so-called Ellis & Baldwin formula [13]. Surprisingly, measurements of the kinematic dipole of sources, while not very well agreeing with each other (see e.g. [14]), find a dipole with a direction roughly compatible with that of the CMB, but with an amplitude significantly larger than expected [4–8, 15, 16]. Recently, the authors of [10, 11], analysing a large sample of quasars, have obtained a result which is $5\sigma$ away from the value predicted from the CMB dipole. In particular, using eq. (1.2), they inferred an observer’s velocity which is twice as large as in eq. (1.1). These results hold when one accounts for the presence of the Galactic plane mask, which leads to a leakage of power from the dipole into higher multipoles, further increasing the tension with the CMB [12]. The authors consider that their finding is a true challenge for the cosmological standard model. This conclusion, as well as the significance of the tension with respect to the CMB, relies however entirely on the ability to predict the value of the prefactor $[2 + x(1 + \alpha)]$ in eq. (1.2), since $v_o$ cannot be measured directly, but only inferred from eq. (1.2). In particular, it has been shown that the redshift dependence of $\alpha$ and $x$ might significantly affect the expected kinematic dipole if these quantities are correlated [17, 18]. Of course, the importance of this effect depends on the strength of the correlation between the parameters [12], which is difficult to assess. However, it was recently shown that one can account for evolution effects by measuring $\alpha$ and $x$ on a subset of sources, close to the flux limit [19]. In the near future, measurements of the near-IR Cosmic Infrared Background (CIB) with Euclid promise to constrain the direction of the dipole to sub-degree precision [20]. This is the reason why we consider this direction fixed in our article. However, the amplitude of this CIB dipole is not only linked to $v_o$ but also to astrophysical factors which need to be predicted or measured accurately.
An independent avenue to shed light on this tension is to use the fact that our peculiar velocity introduces a preferred direction of observation in the sky and consequently breaks statistical isotropy. This effect leads to a mode coupling in the CMB harmonic coefficients leading to non-vanishing correlations between $a_{\ell,m}$ and $a_{\ell\pm1,m}$. These correlations have been detected in the CMB with high significance [21–23]. Using all available multipoles, [21] inferred a velocity $v_o = (384\pm78(\text{stat.})\pm115(\text{syst.}))$ km/s, which is compatible with eq. (1.1) although with much larger error bars.

The same effect can also be searched for in galaxy surveys, and ref. [24] showed that, compared to the CMB, the higher angular resolution of galaxy surveys yields better prospects for precise measurements of the observer velocity with this method. In [24], it was found that the correlations between $a_{\ell,m}$ and $a_{\ell\pm1,m}$ in galaxy surveys are governed by the same combination $D_{\text{kin}}$ as the dipole. In this article, we revisit this computation. We show that three terms were missing in the calculation of [24]. The first two terms have the same multipole dependence as in [24], they thus renormalise $D_{\text{kin}}$ into a new quantity that we call $F_{\text{kin}}$. The last term is very interesting, since it is directly proportional to the observer’s velocity $v_o$, with no dependence on the sources properties. Moreover it does not scale in the same way with multipole $\ell$ as the contribution proportional to $F_{\text{kin}}$. This allows us to efficiently break the degeneracy between $v_o$ and the source population parameters. As a consequence, the correlations between $a_{\ell,m}$ and $a_{\ell\pm1,m}$ (that we call hereafter $(\ell, \ell \pm 1)$ correlations) provide a robust way of measuring $v_o$ without needing a prior knowledge of the properties of the sources. This is particularly important, because precisely measuring the source parameters is not trivial. Our calculation accounts for the evolution of the $(\ell, \ell \pm 1)$ correlations with redshift and is therefore applicable to the case of a tomographic survey, with various redshift bins.

We forecast the detectability of $v_o$ using the $(\ell, \ell \pm 1)$ correlations, for a Euclid-like photometric galaxy survey. Combining the 10 redshift bins of Euclid, we find that the observer velocity can be detected at nearly $3\sigma$ if its amplitude is as predicted from the CMB dipole. The significance of the detection is limited by the degeneracy between the contribution proportional to $v_o$ and the contribution proportional to $F_{\text{kin}}$. We show that this degeneracy can be strongly reduced if we assume some prior knowledge on $F_{\text{kin}}$. Interestingly, with a prior of only 30% on the source parameters that enter into $F_{\text{kin}}$, we reach a detection significance of $9\sigma$, corresponding to a precision of 11% on the measurement of $v_o$. If we tighten the priors to 10%, we reach a significance of $24\sigma$, corresponding to a precision of 4%. Importantly, by building an optimal estimator for $v_o$, we show that our method does not rely on the modelling of the unboosted angular power spectrum of galaxies, $C_\ell$. The estimator is instead directly built from the measured $C_\ell$. For this reason, the determination of $v_o$ does not depend on our ability to model non-linear physics.

The present article is structured as follows: in section 2, we derive the impact of the observer velocity on the perturbed galaxy number counts, accounting for all terms of the form $v_o \times \delta_g$, where $\delta_g$ denotes the galaxy density perturbation. In section 3, we derive the correlations between the $\ell$ and $\ell \pm 1$ spherical harmonic coefficients. In particular, we show the appearance of an extra term directly proportional to $v_o$ in the $(\ell, \ell + 1)$ correlators which was missed in [24]. In section 4, we derive the maximum likelihood estimator for $F_{\text{kin}}(z)$ and $v_o$ using the $(\ell, \ell + 1)$ correlators, and we study the degeneracy between $F_{\text{kin}}(z)$ and
$v_o$. We then show how this degeneracy can be lifted by the study of correlations of different redshift bins or of a multi-tracer analysis. Finally, in section 5, we forecast the prospects for a Euclid-like galaxy survey and potential strategies to further improve the measurement precision. Several more technical aspects are relegated to the appendices.

We work in units in which the speed of light $c = 1$ such that $\beta \equiv v_o/c = v_o$. We consider a spatially flat perturbed Friedmann-Lemaître Universe with background metric $ds^2 = a^2(\tau)[-d\tau^2 + \delta_{ij}dx^i dx^j]$. An overdot denotes a derivative with respect to conformal time $\tau$, which relates to cosmic time $t$ via the scale factor $a(\tau)$: $a(\tau)d\tau = dt$. The conformal Hubble rate is defined as $\mathcal{H} = \dot{a}/a = H_a$, where $H$ is the physical Hubble rate. We denote the Kronecker delta by $\delta_{ij}^K$.

## 2 Effect of the observer velocity on the perturbed galaxy number counts

In this section, we derive the relation between the boosted and unboosted galaxy number counts. We adapt the calculations presented in [17, 25] and work out the details of the calculation sketched in [24], which misses an important term. We then compute the harmonic coefficients of the observed boosted number counts.

### 2.1 Boosted galaxy number counts

Let us denote by $\tilde{N}(\tilde{z}, \tilde{n})$ the unboosted number density of galaxies per redshift bin $d\tilde{z}$ and per solid angle $d\tilde{\Omega}$ in direction $\tilde{n}$ and at redshift $\tilde{z}$. Correspondingly, we denote the boosted, observed quantity $N(z, \tilde{n})$. We choose coordinates such that the boost $\beta \equiv v_o/c$ is aligned with the $+e_z$ axis. As shown in [17, 25, 26], the dipole in the distribution of galaxies is due to the impact of the boost on the monopole, the isotropic part of the distribution of galaxies $\tilde{N}(\tilde{z})$. The correlations between the harmonic coefficients, $a_{\ell,m}$ and $a_{\ell\pm1,m}$, on the other hand are due to the effect of the boost on the perturbed anisotropic distribution of galaxies $\tilde{N}(\tilde{z}, \tilde{n})$. In our derivation, we keep therefore terms that are linear in $v_o$ and that multiply the perturbed galaxy distribution. These terms are formally beyond linear order. As shown in [27], various terms contribute to the perturbed galaxy distribution at linear order. Calculating the impact of the boost on all those terms is beyond the scope of this article. Instead we focus on the dominant contribution to the perturbed galaxy distribution, namely the density fluctuation $\delta_g = b \cdot \delta_m$, where $b$ is the galaxy bias, and $\delta_m$ is the matter density contrast. We neglect redshift-space distortions, since we will use photometric galaxy surveys with relatively wide redshift bins, where redshift-space distortions are strongly suppressed. We also neglect the effect of lensing magnification. This term is completely negligible at low redshift, but it starts to become relevant at high redshift [28]. Neglecting it may therefore slightly affect the constraints in the highest redshift bins of Euclid, but we do not expect this to change the overall results. Finally, we also neglect the terms that are quadratic in the observer velocity, and which are suppressed with respect to terms of the form $v_o \cdot \delta_m$, since $v \sim \frac{H}{k}\delta_m \ll \delta_m$ for sub-horizon scales.

We now proceed to the calculation of the boosted number counts in terms of the unboosted number counts. The boosted redshift $z$, direction $\tilde{n} = (\theta, \phi)$ and solid angle $d\tilde{\Omega}$ are related
to their unboosted counterparts, denoted with tilde, by the following equations, see e.g. [17]

$$1 + \tilde{z} = (1 + z)(1 + \beta \cdot \mathbf{n}), \qquad (2.1)$$

$$\tilde{\mathbf{n}} = (1 + \mathbf{n} \cdot \beta)\mathbf{n} - \beta, \qquad (2.2)$$

$$\sin \tilde{\theta} = (1 + \mathbf{n} \cdot \beta) \sin \theta, \qquad (2.3)$$

$$\tilde{\theta} = \theta + \beta \sin \theta + \mathcal{O}(\beta^2), \qquad (2.4)$$

$$d\tilde{\Omega} = d\Omega \times (1 + 2\mathbf{n} \cdot \beta + \mathcal{O}(\beta^2)). \qquad (2.5)$$

The number of galaxies is conserved under the effect of the boost: all boosted galaxies observed in a solid angle $d\Omega$ and redshift bin $dz$ with flux above $S_*$, would be observed in solid angle $d\tilde{\Omega}$ and redshift bin $d\tilde{z}$ if the observer were not moving. Moreover, these galaxies all have an intrinsic luminosity larger than $L_*(z, \mathbf{n}) = 4\pi d_L^2(z, \mathbf{n})S_*$, where $d_L$ denotes the luminosity distance. As a consequence, we can write

$$N(z, \mathbf{n}, S_*) \, dz \, d\Omega = N(z, \mathbf{n}, L_*(z, \mathbf{n})) \, dz \, d\Omega = \tilde{N}(\tilde{z}, \tilde{\mathbf{n}}, L_*(z, \mathbf{n})) \, d\tilde{z} \, d\tilde{\Omega}, \qquad (2.6)$$

where we shortened notations for simplicity with $N(z, \mathbf{n}, S_*) \equiv N(z, \mathbf{n}, S > S_*)$, and similarly for $N(z, \mathbf{n}, L_*(z, \mathbf{n}))$. Dividing both sides of the equation by $d\Omega dz$ and Taylor expanding $\tilde{N}$ around the unboosted luminosity threshold $\tilde{L}_*$ (i.e. the luminosity threshold that would correspond to $S_*$ if the observer was not moving), we obtain

$$N(z, \mathbf{n}, S_*) = \left[ \tilde{N}(\tilde{z}, \tilde{\mathbf{n}}, \tilde{L}_*(\tilde{z}, \tilde{\mathbf{n}})) + (L_*(z, \mathbf{n}) - \tilde{L}_*(\tilde{z}, \tilde{\mathbf{n}})) \frac{\partial}{\partial \tilde{L}_*} \tilde{N}(\tilde{z}, \tilde{\mathbf{n}}, \tilde{L}_*) \right] \frac{d\tilde{\Omega}}{d\Omega} \frac{d\tilde{z}}{dz}. \qquad (2.7)$$

The luminosity derivative reads to zeroth order in $\beta$

$$\frac{\partial}{\partial \tilde{L}_*} \tilde{N}(\tilde{z}, \tilde{\mathbf{n}}, \tilde{L}_*) = - \frac{x(\tilde{z})}{\tilde{L}_*} \tilde{N}(\tilde{z}, \tilde{\mathbf{n}}, \tilde{L}_*) + \frac{\tilde{N}(\tilde{z}, \tilde{L}_*)}{\tilde{L}_*} \frac{\partial}{\partial \ln \tilde{L}_*} \left( \frac{\delta \tilde{N}(\tilde{z}, \tilde{\mathbf{n}}, \tilde{L}_*)}{\tilde{N}(\tilde{z}, \tilde{L}_*)} \right), \qquad (2.8)$$

where $\tilde{N}$ denotes the average number of galaxies at redshift $\tilde{z}$, $\delta \tilde{N}$ is the departure from the average, and the average magnification bias is defined as $x(\tilde{z}) \equiv -\partial \ln \tilde{N}/\partial \ln \tilde{L}_*$.

As the flux cut $S_*$ is a property of the survey, independent of the boost, it corresponds to different luminosity cuts for the boosted and unboosted observer,

$$S_* = \frac{\tilde{L}_*(\tilde{z}, \tilde{\mathbf{n}})}{4\pi d_L^2(\tilde{z}, \tilde{\mathbf{n}})} = \frac{L_*(z, \mathbf{n})}{4\pi d_L^2(z, \mathbf{n})}, \qquad (2.9)$$

from which we obtain to first order in perturbations

$$\frac{L_*(z, \mathbf{n}) - \tilde{L}_*(\tilde{z}, \tilde{\mathbf{n}})}{\tilde{L}_*(\tilde{z}, \tilde{\mathbf{n}})} = \frac{2(d_L(z, \mathbf{n}) - \tilde{d}_L(\tilde{z}, \tilde{\mathbf{n}}))}{\tilde{d}_L(\tilde{z}, \tilde{\mathbf{n}})}. \qquad (2.10)$$

The impact of the observer's velocity on the luminosity distance at fixed observed redshift is given by [29]

$$\frac{(d_L(z, \mathbf{n}) - \tilde{d}_L(\tilde{z}, \tilde{\mathbf{n}}))}{\tilde{d}_L(\tilde{z}, \tilde{\mathbf{n}})} = \frac{\mathbf{n} \cdot \beta}{r\mathcal{H}}, \qquad (2.11)$$
where $r$ is the comoving distance at redshift $z$ and $\mathcal{H}$ is the conformal Hubble parameter. Combining eqs. (2.1), (2.5), (2.7), (2.8) and (2.11), we obtain to first order in $\beta$

$$ N(z, \mathbf{n}, S_*) = (1 + 3\mathbf{n} \cdot \mathbf{\beta}) \tilde{N}(\tilde{z}, \tilde{\mathbf{n}}, \tilde{L}_*(\tilde{z}, \tilde{\mathbf{n}})) + \frac{2\mathbf{n} \cdot \mathbf{\beta}}{r\mathcal{H}} \left[ -x(\tilde{z})\tilde{N}(\tilde{z}, \tilde{\mathbf{n}}, \tilde{L}_*(\tilde{z}, \tilde{\mathbf{n}})) + \tilde{N}(\tilde{z}, \tilde{L}_*) \frac{\partial}{\partial \ln \tilde{L}_*} \left( \frac{\delta \tilde{N}(\tilde{z}, \tilde{\mathbf{n}}, \tilde{L}_*)}{\tilde{N}(\tilde{z}, \tilde{L}_*)} \right) \right]. \quad (2.12) $$

The quantity $\tilde{N}(\tilde{z}, \tilde{\mathbf{n}}, \tilde{L}_*(\tilde{z}, \tilde{\mathbf{n}}))$, which appears on the right-hand side of eq. (2.12), is the number of galaxies detected by an unboosted observer, at rest with respect to the CMB frame. From eq. (2.9), we see that all galaxies $\tilde{N}$ that have intrinsic luminosity larger than $\tilde{L}_*$ in the unboosted CMB frame, would be observed with a flux $S > S_*$ in that frame. Hence we can rewrite $\tilde{N}(\tilde{z}, \tilde{\mathbf{n}}, \tilde{L}_*(\tilde{z}, \tilde{\mathbf{n}})) = \tilde{N}(\tilde{z}, \tilde{\mathbf{n}}, S_*)$. This quantity depends on the unboosted redshift $\tilde{z}$ and the unboosted direction $\tilde{\mathbf{n}}$, that are not observable: one only measures the boosted redshift $z$ and angular direction $\mathbf{n}$. Hence to predict the value of $N(z, \mathbf{n}, S_*)$ in terms of observable quantities, it is necessary to further Taylor-expand $\tilde{N}(\tilde{z}, \tilde{\mathbf{n}}, S_*)$ around $z$ and $\mathbf{n}$. To first order in $\beta$, we have

$$ \tilde{N}(\tilde{z}, \tilde{\mathbf{n}}, S_*) = \tilde{N}(z, \mathbf{n}, S_*) + (\tilde{\mathbf{n}} - \mathbf{n}) \frac{\partial \tilde{N}(z, \mathbf{n}, S_*)}{\partial \mathbf{n}} + (\tilde{z} - z) \frac{\partial \tilde{N}(z, \mathbf{n}, S_*)}{\partial z}. \quad (2.13) $$

Using eq. (2.1), the third term can be written as

$$ (\tilde{z} - z) \frac{\partial \tilde{N}(z, \mathbf{n}, S_*)}{\partial z} = \mathbf{n} \cdot \mathbf{\beta} \frac{\partial \tilde{N}(z, \mathbf{n}, S_*)}{\partial \ln(1+z)}. \quad (2.14) $$

We first write

$$ \tilde{N}(z, \mathbf{n}, S_*) = \tilde{\tilde{N}}(z, S_*) \left( 1 + \frac{\delta \tilde{N}(z, \mathbf{n}, S_*)}{\tilde{\tilde{N}}(z, S_*)} \right), \quad (2.15) $$

and then we express the background number density (per unit solid angle and per redshift bin) in terms of the background comoving number density $n_c(z, S_*)$. Using that $d\Omega dz = (1+z) (\mathcal{H}/r^2) dV_c$, we obtain

$$ \tilde{\tilde{N}}(z, S_*) = \frac{r^2}{\mathcal{H}(1+z)} n_c(z, S_*) . \quad (2.16) $$

Evaluating the logarithmic redshift derivative of eq. (2.15), using eq. (2.16) we obtain

$$ \frac{\partial \tilde{N}(z, \mathbf{n}, S_*)}{\partial \ln(1+z)} = \left[ \frac{2}{r(z)\mathcal{H}(z)} - 1 + \frac{\dot{\mathcal{H}}}{\mathcal{H}^2} - f_{\text{evol}}(z) \right] \tilde{N}(z, \mathbf{n}, S_*) + \frac{\partial}{\partial \ln(1+z)} \left( \frac{\delta \tilde{N}(z, \mathbf{n}, S_*)}{\tilde{\tilde{N}}(z, S_*)} \right) \tilde{\tilde{N}}(z, S_*) . \quad (2.17) $$

Here, the overdot denotes the derivative with respect to conformal time and we have introduced the so-called evolution bias

$$ f_{\text{evol}}(z) = - \frac{\partial \ln n_c(z, S_*)}{\partial \ln(1+z)}. \quad (2.18) $$
Inserting this into eq. (2.13) and using eq. (2.4) to rewrite $(\tilde{n} - n)$, we obtain

$$ N(z, \mathbf{n}, S_*) = [1 + \cos \theta D_{\text{kin}}(z)] \tilde{N}(z, S_*) \\ + [1 + \cos \theta F_{\text{kin}}(z, \mathbf{n}) + \beta \sin \theta \partial_{\theta}] \tilde{N}(z, S_*) \delta_g(z, \mathbf{n}, S_*), \quad (2.19) $$

where

$$ \delta_g(z, \mathbf{n}, S_*) = \frac{\delta \tilde{N}(z, \mathbf{n}, S_*)}{\tilde{N}(z, S_*)}, \quad (2.20) $$

and

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

$$ F_{\text{kin}}(z, \mathbf{n}) = D_{\text{kin}}(z) + \left[ \frac{\partial \ln \delta_g}{\partial \ln(1+z)} + \frac{2}{r\mathcal{H}} \frac{\partial \ln \delta_g}{\partial \ln S_*} \right] \beta. \quad (2.22) $$

At linear order, the galaxy perturbations, $\delta_g$, can be written as the product of a scale-independent bias, which depends on redshift and flux threshold, and the matter density fluctuations: $\delta_g = b(z, S_*) \delta_m(z, \mathbf{n})$. In this regime $F_{\text{kin}}$ depends only on redshift and can be written as

$$ F_{\text{kin}}(z) = D_{\text{kin}}(z) + \left[ -f(z) + \frac{\partial \ln b}{\partial \ln(1+z)} + \frac{2}{r\mathcal{H}} \frac{\partial \ln b}{\partial \ln S_*} \right] \beta, \quad (2.23) $$

where $f \equiv d \ln \delta_m(a, \mathbf{n})/d \ln a$ is the growth rate of structure, which is scale-independent in the linear regime.

In the highly non-linear regime, the bias and growth rate would get scale-dependent corrections. We do not consider these in the main text, although we explain in appendix D.3 how one could consider them and show that they should not impact the conclusions of the article.

The first contribution in eq. (2.19), proportional to the background number of galaxies $\tilde{N}(z, S_*)$ generates only a dipolar modulation in the boosted number of galaxies. As already shown in [17, 25, 26], the amplitude of the dipole is governed by the combination $D_{\text{kin}}(z)$ given in (2.21), which depends on the magnification bias and evolution bias of the population of galaxies. The second contribution in eq. (2.19), proportional to the galaxy fluctuations, $\delta_g(z, \mathbf{n})$, generates correlations between $a_{\ell,m}$ and $a_{\ell \pm 1, m}$. We see that this contribution contains two terms: one proportional to $\cos \theta$, and a second term which depends on the variation of the galaxy fluctuations with $\theta$. This second term is a remapping contribution: it follows directly from the remapping of $\tilde{N}$ from the unboosted direction $\tilde{n}$ to the boosted observed direction $\mathbf{n}$.

Comparing the contribution proportional to $\tilde{N}(z, S_*)$ with the contribution proportional to $\delta_g(z, \mathbf{n})$ in eq. (2.19), we see that the second one is strongly suppressed with respect to the first one, since $\delta_g \ll 1$. Hence, the dipolar modulation is completely dominated by the first contribution, proportional to $D_{\text{kin}}$. The correlations between $a_{\ell,m}$ and $a_{\ell \pm 1, m}$, on the other hand, are only due to the second contribution. What is very interesting, is that this second contribution has two terms, $F_{\text{kin}}$ and the remapping, which do not have the same
angular dependence in $\theta$. As we will see in section 2.2, this will allow us to separate these two contributions, and measure separately $F_{\text{kin}}(z)$ and $\beta$. This provides a strong advantage with respect to measurements of the dipole, since $\beta$ can be measured without requiring a model for the characteristics of the sources. Finally, let us note that $D_{\text{kin}}$ and $F_{\text{kin}}$ are not exactly the same. $F_{\text{kin}}$ contains additional contributions, which depend on the evolution of the galaxy fluctuations with redshift and with flux cut.

## 2.2 Spherical harmonic coefficients

In order to find the effect of a boost on the statistical distribution of galaxies, we decompose the number density of galaxies in spherical harmonics

$$a_{\ell,m}(z) = \int_{S_2} d^2\mathbf{n} N(z, \mathbf{n}, S_*) Y_{\ell,m}^*(\mathbf{n}). \qquad (2.24)$$

In the following, we leave the dependence on the flux threshold and redshift implicit. We also decompose the unboosted number counts in spherical harmonics, with coefficients $\tilde{a}_{\ell,m}$ indicated with a tilde. The statistics of $\tilde{a}_{\ell,m}$ can be predicted from a theoretical model of galaxy clustering. Our goal is to relate the observed coefficients $a_{\ell,m}$ to the underlying $\tilde{a}_{\ell,m}$. We have

$$a_{\ell,m} = \int_{S_2} d^2\mathbf{n} N(z, \mathbf{n}) Y_{\ell,m}^*(\mathbf{n}) \qquad (2.25)$$

$$ \begin{aligned} & = D_{\text{kin}}(z) \tilde{N}(z) \int_{S_2} d^2\mathbf{n} \cos\theta Y_{\ell,m}^*(\mathbf{n}) \\ & \quad + \tilde{N}(z) \int_{S_2} d^2\mathbf{n} [1 + \cos\theta F_{\text{kin}}(z) + \beta \sin\theta \partial_\theta] \delta_g(z, \mathbf{n}) Y_{\ell,m}^*(\mathbf{n}) \\ & = \tilde{a}_{\ell,m} + \sqrt{\frac{4\pi}{3}} D_{\text{kin}}(z) \tilde{N}(z) \delta_{\ell,1}^K \delta_{m,0}^K + F_{\text{kin}}(z) \tilde{N}(z) \int_{S_2} d^2\mathbf{n} \delta_g(z, \mathbf{n}) \cos\theta Y_{\ell,m}^*(\mathbf{n}) \quad (2.26) \\ & \quad - \beta \tilde{N}(z) \int_{S_2} d^2\mathbf{n} \delta_g(z, \mathbf{n}) (\sin\theta \partial_\theta Y_{\ell,m}^*(\mathbf{n}) + 2 \cos\theta Y_{\ell,m}^*(\mathbf{n})) \end{aligned} $$

where we used that $\cos\theta = \sqrt{\frac{4\pi}{3}} Y_{1,0}(\mathbf{n})$ and we have performed an integration by part in the last line. Therefore, we write

$$a_{\ell,m} = \tilde{a}_{\ell,m} + \delta a_{\ell,m}^{\text{Dip}} + \delta a_{\ell,m}^{\text{F}} - \delta a_{\ell,m}^{\text{Remap}}, \qquad (2.27)$$

with

$$\delta a_{\ell,m}^{\text{Dip}} = \sqrt{\frac{4\pi}{3}} D_{\text{kin}}(z) \tilde{N}(z) \delta_{\ell,1}^{K} \delta_{m,0}^{K}, \qquad (2.28)$$

$$\delta a_{\ell,m}^{\text{F}} = F_{\text{kin}}(z) \tilde{N}(z) \int_{S_2} d^2\mathbf{n} \, \delta_g(z, \mathbf{n}) \cos\theta Y_{\ell,m}^*(\mathbf{n}), \qquad (2.29)$$

$$\delta a_{\ell,m}^{\text{Remap}} = \beta \tilde{N}(z) \int_{S_2} d^2\mathbf{n} \, \delta_g(z, \mathbf{n}) \left( \sin\theta \partial_\theta Y_{\ell,m}^*(\mathbf{n}) + 2 \cos\theta Y_{\ell,m}^*(\mathbf{n}) \right). \qquad (2.30)$$

The minus sign in eq. (2.27) has been introduced for convenience. The superscript “Dip” denotes the dipole contribution, which affects only the multipole $\ell = 1, m = 0$; the “F” superscript refers to the term proportional to $F_{\text{kin}}$, while the “Remap” refers to remapping
contribution, directly proportional to $\beta$. Using that the monopole $\int d^2n \tilde{N}(z, n) Y_{0,0}(n) = \sqrt{4\pi} \tilde{N}(z) = \tilde{a}_{0,0}^{\text{BG}}$, the dipole term in eq. (2.28) becomes

$$ \delta a_{\ell,m}^{\text{Dip}} = \frac{D_{\text{kin}}(z)}{\sqrt{3}} \tilde{a}_{0,0}^{\text{BG}} \delta_{\ell 1}^{\text{K}} \delta_{m 0}^{\text{K}}, \qquad (2.31) $$

where we have introduced the upper script “BG” to emphasise that this term is the only one that depends on the background number of galaxies $\tilde{N}(z)$, averaged over direction. All the other terms depends on the galaxy fluctuations $\delta_g$. The F term in eq. (2.29) can be rewritten, using calculations from appendix A.1

$$ \delta a_{\ell,m}^{\text{F}} = F_{\text{kin}}(z) \tilde{N}(z) \int_{S_2} d^2n \, \delta_g(z, n) [A_{\ell,m} Y_{\ell-1,m}^*(n) + A_{\ell+1,m} Y_{\ell+1,m}^*(n)] \quad (2.32) $$

$$ = F_{\text{kin}}(z) [A_{\ell,m} \tilde{a}_{\ell-1,m} + A_{\ell+1,m} \tilde{a}_{\ell+1,m}], \qquad (2.33) $$

where

$$ A_{\ell,m} \equiv \sqrt{\frac{\ell^2 - m^2}{(2\ell - 1)(2\ell + 1)}}. \qquad (2.34) $$

The Remap term requires identities derived in appendices A.1 & A.2, from which we obtain

$$ \delta a_{\ell,m}^{\text{Remap}} = \beta \tilde{N}(z) \int_{S_2} d^2n \, \delta_g(z, n) \left[ -(\ell + 1) A_{\ell,m} Y_{\ell-1,m}^*(n) + \ell A_{\ell+1,m} Y_{\ell+1,m}^*(n) \right] \quad (2.35) $$

$$ + 2\beta \tilde{N}(z) \int_{S_2} d^2n \, \delta_g(z, n) [A_{\ell,m} Y_{\ell-1,m}^*(n) + A_{\ell+1,m} Y_{\ell+1,m}^*(n)] $$

$$ = \beta [-(\ell - 1) A_{\ell,m} \tilde{a}_{\ell-1,m} + (\ell + 2) A_{\ell+1,m} \tilde{a}_{\ell+1,m}]. $$

In summary, the observed harmonic coefficient relate to the unboosted harmonic coefficients $\tilde{a}_{\ell,m}$ by

$$ a_{\ell,m} = \tilde{a}_{\ell,m} + \frac{D_{\text{kin}}}{\sqrt{3}} \tilde{a}_{0,0}^{\text{BG}} \delta_{\ell 1}^{\text{K}} \delta_{m 0}^{\text{K}} + F_{\text{kin}} [A_{\ell,m} \tilde{a}_{\ell-1,m} + A_{\ell+1,m} \tilde{a}_{\ell+1,m}] \\ - \beta [-(\ell-1) A_{\ell,m} \tilde{a}_{\ell-1,m} + (\ell+2) A_{\ell+1,m} \tilde{a}_{\ell+1,m}], \qquad (2.36) $$

which may be written alternatively as

$$ a_{\ell,m} = \tilde{a}_{\ell,m} + \frac{D_{\text{kin}}}{\sqrt{3}} \tilde{a}_{0,0}^{\text{BG}} \delta_{\ell 1}^{\text{K}} \delta_{m 0}^{\text{K}} + A_{\ell,m} (F_{\text{kin}} + \beta(\ell-1)) \tilde{a}_{\ell-1,m} \\ + A_{\ell+1,m} (F_{\text{kin}} - \beta(\ell+2)) \tilde{a}_{\ell+1,m}. \qquad (2.37) $$

The terms proportional to $\beta$ were previously missed in the literature [24]. These terms open up the possibility to measure both, $\beta$ and $F_{\text{kin}}$, independently. We investigate this possibility in detail in the following.

Eq. (2.37) has been computed for a fixed redshift $z$. In practice however, one always averages over redshift bins. Labelling by $i$ the $i$th redshift bin, and denoting by $N_i$ the total
number of galaxies in bin $i$ we can write the averaged boosted $a_{\ell,m}$ in that bin

$$ 
\begin{aligned}
a_{\ell,m}^{\text{av}}(z_i) &\equiv \frac{1}{N_i} \int dz \, a_{\ell,m}(z) \qquad (2.38) \\
&= \tilde{a}_{\ell,m}^{\text{av}}(z_i) - \beta \left[ -(\ell-1) A_{\ell,m} \tilde{a}_{\ell-1,m}^{\text{av}}(z_i) + (\ell+2) A_{\ell+1,m} \tilde{a}_{\ell+1,m}^{\text{av}}(z_i) \right] \\
&\quad + \frac{1}{N_i} \int dz \, F_{\text{kin}}(z) \left[ A_{\ell,m} \tilde{a}_{\ell-1,m}(z) + A_{\ell+1,m} \tilde{a}_{\ell+1,m}(z) \right] \\
&\quad + \frac{1}{N_i} \int dz \, D_{\text{kin}}(z) \tilde{a}_{0,0}^{\text{BG}}(z) \delta_{\ell,1}^{K} \delta_{m,0}^{K}.
\end{aligned}
$$

In the following, for simplicity we will assume that $F_{\text{kin}}(z)$ and $D_{\text{kin}}(z)$ are constant within the redshift bin, such that we can take them out of the integral over $z$ and evaluate them at the mean redshift $z_i$. We recover therefore eq. (2.37), but with the $a_{\ell,m}(z)$ replaced by $a_{\ell,m}^{\text{av}}(z_i)$. To simplify the notation, we drop the superscript “av” in the rest of the derivation.

# 3 Impact on observable statistics

In this section, we study the impact of the boost on the statistical properties of the spherical harmonic coefficients. We first recover the leakage of the monopole into the dipole. Then, we construct an observable which is proportional to the boost for which we express the proportionality factor in terms of the clustering properties of the galaxies.

## 3.1 Cosmic dipole

For most cosmological signals, the average on the sky is much larger than the fluctuations around the mean. This hierarchy between the monopole ($\propto \tilde{a}_{0,0}^{\text{BG}}$) and the fluctuations ($\tilde{a}_{\ell,m}$ for $\ell \ge 1$) implies that any leakage of the monopole to the higher order multipoles may be important. The effect of the observer's velocity, as described in the previous section, provides such a leakage from the monopole into the dipole, leading to

$$ 
\begin{aligned}
a_{1,0} &= \tilde{a}_{1,0} + \frac{D_{\text{kin}}}{\sqrt{3}} \tilde{a}_{0,0}^{\text{BG}} + F_{\text{kin}} A_{1,0} \tilde{a}_{0,0} + \beta (1-1) A_{1,0} \tilde{a}_{0,0} + \mathcal{O}(\tilde{a}_{2,0} \beta) \\
&\approx \tilde{a}_{1,0} + \frac{D_{\text{kin}}}{\sqrt{3}} \tilde{a}_{0,0}^{\text{BG}},
\end{aligned}
\qquad (3.1) $$

where we have used that the spherical average of $\delta_g(z, \mathbf{n})$ vanishes, since $\delta_g$ is customarily defined as the deviation from the measured mean number of galaxies, leading to $\tilde{a}_{0,0} = 0$ in the first line of eq. (3.1). Furthermore, we have neglected the contribution from the quadrupole $\mathcal{O}(a_{2,0} \beta)$, as it depends on $\delta_g$ and is therefore expected to be smaller than $\tilde{a}_{0,0}^{\text{BG}} D_{\text{kin}}$. Under this assumption, we see that the remapping term does not affect the cosmic dipole, which is therefore only sensitive to $D_{\text{kin}}$.

For a more intuitive interpretation, it helps to write this result in real space. Using eq. (2.28) and assuming $|\tilde{a}_{\ell,m}| \ll \tilde{a}_{0,0}^{\text{BG}} \quad \forall \ell \ge 1$, we obtain

$$ 
\begin{aligned}
N(z, \mathbf{n}) &\approx \sum_{\ell=0}^{1} \sum_{m=-\ell}^{\ell} a_{\ell,m} Y_{\ell,m}(\mathbf{n}) = \tilde{N}(z) \left[ 1 + D_{\text{kin}} \sqrt{\frac{4\pi}{3}} Y_{1,0}(\mathbf{n}) \right] + \sum_{m=-1}^{1} \tilde{a}_{1,m} Y_{1,m}(\mathbf{n}) \quad (3.2) \\
&\approx \tilde{N}(z) [1 + D_{\text{kin}} \cos(\theta)],
\end{aligned}
\qquad (3.3) $$
where in the second line, we further assumed that $\sqrt{3}|\tilde{a}_{1,m}| \ll D_{\text{kin}} \tilde{a}_{0,0}^{\text{BG}} \sim D_{\text{kin}} \tilde{N}(z)$. This amounts to neglecting the intrinsic dipole in comparison to the kinematically generated dipole. This assumption may not always be a good assumption, especially at low redshift, where the intrinsic dipole due to density fluctuations may be non-negligible. See for example [14, 26] for a study of the intrinsic dipole of radio galaxies and [30] for the dipole of type 1a Supernovae. Eq. (3.3) is the usual expression for the kinematic dipole, previously derived in [17, 25, 26]. Let us now study the impact of the boost on the correlation between nearest multipoles.

## 3.2 Correlation between multipoles

We want to compute the correlations between different multipoles due to the observer velocity. We assume that the galaxy density $\delta_g(z, \mathbf{n})$ is statistically isotropic. This translates to the 2-point correlation function of the spherical harmonic coefficients obeying

$$ \mathbb{E}[\tilde{a}_{\ell,m} \tilde{a}_{\ell',m'}^*] = C_{\ell} \delta_{\ell\ell'}^K \delta_{mm'}^K. \qquad (3.4) $$

Note that we use the notation $\mathbb{E}$ for expectation values (or ensemble averages), and we reserve brackets for the definition of an inner product later in section 4. The situation becomes different for the boosted spherical harmonic coefficients $a_{\ell,m}$. Indeed, the boost breaks statistical isotropy by introducing a preferred direction in the sky: the direction of $v_o$. Setting $v_o \propto e_z$, we have found that $a_{\ell,m}$ depend on $\tilde{a}_{\ell\pm1,m}$, which leads to cross-correlations between different $\ell$'s that no longer vanish. To measure this departure from isotropy, we introduce the following observable, which depends on $m$ and $\ell$

$$ O_{\ell,m} \equiv a_{\ell,m} a_{\ell+1,m}^* . \qquad (3.5) $$

At first order in $\beta$, the statistical average of $O_{\ell,m}$ is given by

$$ \begin{aligned} \mathbb{E}[O_{\ell,m}] &= \mathbb{E}[\tilde{a}_{\ell,m} \times (\delta a_{\ell+1,m}^{\text{F}} - \delta a_{\ell+1,m}^{\text{Remap}})^*] + \mathbb{E}[(\delta a_{\ell,m}^{\text{F}} - \delta a_{\ell,m}^{\text{Remap}}) \times \tilde{a}_{\ell+1,m}^*] \\ &= \mathbb{E}[\tilde{a}_{\ell,m} \times (F_{\text{kin}} A_{\ell+1,m} \tilde{a}_{\ell,m}^* + \beta \ell A_{\ell+1,m} \tilde{a}_{\ell,m}^*)] \\ &\quad + \mathbb{E}[(F_{\text{kin}} A_{\ell+1,m} \tilde{a}_{\ell+1,m} - \beta(\ell+2) A_{\ell+1,m} \tilde{a}_{\ell+1,m}) \times \tilde{a}_{\ell+1,m}^*] \\ &= (F_{\text{kin}} + \beta \ell) A_{\ell+1,m} C_{\ell} + (F_{\text{kin}} - \beta(\ell+2)) A_{\ell+1,m} C_{\ell+1} \\ &= A_{\ell+1,m} \cdot [F_{\text{kin}}(C_{\ell} + C_{\ell+1}) - \beta((\ell+2)C_{\ell+1} - \ell C_{\ell})]. \end{aligned} \qquad (3.6) $$

For the statistical analysis in section 4 and 5, we also need the covariance matrix of the $O_{\ell,m}$, which is defined as

$$ \text{Cov}(O_{\ell_1, m_1}, O_{\ell_2, m_2}) = \mathbb{E}[O_{\ell_1, m_1} O_{\ell_2, m_2}^*] - \mathbb{E}[O_{\ell_1, m_1}] \mathbb{E}[O_{\ell_2, m_2}^*]. \qquad (3.7) $$
We derive this covariance using the $f_{SKY}$ approximation<fnref n="1" /> and assuming gaussianity for the galaxy counts.<fnref n="2" /> To leading order in $\beta$, it gives

$$ \mathrm{Cov}(O_{\ell_1, m_1}, O_{\ell_2, m_2}) = \frac{C_{\ell_1} C_{\ell_1+1} \delta_{\ell_1 \ell_2}^{\mathrm{K}} \delta_{m_1 m_2}^{\mathrm{K}} + \mathcal{O}(\beta)}{f_{\mathrm{SKY}}} \quad (3.8) $$

As expected from statistical isotropy, we see that for a non-boosted observer ($\beta = 0$), the statistical average of $O_{\ell,m}$ in eq. (3.6) vanishes: there are no correlations between different multipoles. The boost, however, breaks statistical isotropy. To intuitively understand why this leads to correlations between different $\ell$'s we can look, for example, at the effect of aberration. Perturbations in the hemisphere in the direction of motion subtend a slightly smaller solid angle than in the opposite hemisphere as a result of aberration of the solid angle in which one counts $N(z, \mathbf{n})$. The consequence is that perturbations of different scales, parametrized by the $\ell$ index appear to a moving observer to be correlated. From eq. (3.6), we see that the observable $O_{\ell,m}$ contains two contribution: one proportional to $F_{\mathrm{kin}}$, and another one, the remapping term, directly proportional to $\beta$. Since these two contributions scale differently with $\ell$, as we see from eq. (3.6), we can in principle measure them separately by combining different $\ell$'s. This is particularly interesting, since it allows us to measure $\beta$ without needing prior knowledge of the galaxy specificities, parametrised by $x(z)$ and $f_{\mathrm{evol}}(z)$ which enter in $F_{\mathrm{kin}}$, see eq. (2.21).<fnref n="3" /> From eq. (3.6) we see that our ability to separate the two contributions is directly related to the shape of the $C_\ell$'s. A spectrum with structures will help breaking the degeneracies between the two terms. We shall discuss this in more detail in the next section.

# 4 Estimators for $\beta$ and $F_{\mathrm{kin}}$

In the previous section, we have introduced observables, the $O_{\ell,m}$'s which have a contribution proportional to $F_{\mathrm{kin}}$ and another one directly proportional to the observer velocity $\beta$. Since some part of the coefficients in front of each of these contributions are known, such as the $A_{\ell,m}$ via eq. (2.34), or can be estimated from the measured galaxy number counts, such as the $C_\ell$'s, we can build estimators for $\beta$ and $F_{\mathrm{kin}}$ from the observable $O_{\ell,m}$'s. In the first subsection, we detail how much can be achieved in an agnostic single tracer case. We then generalize the estimators to the case of multitracers and redshift tomography and finally, we study the case where we have some prior knowledge on $F_{\mathrm{kin}}$.

## 4.1 Single tracer case

In this section, we lay out the procedure to estimate $\beta$ and $F_{\mathrm{kin}}$ separately from the observed $O_{\ell,m}$'s.

<footnote n="1">¹This is a popular approximate treatment of partial sky coverage. In appendix C, we sketch how to treat exactly partial sky. We leave the implementation of this exact treatment for future work, as we expect partial sky coverage to have a minimal impact on the high multipoles where the bulk of our statistical power lies.</footnote>

<footnote n="2">²This sounds like a strong assumption when pushing the analysis to high multipoles. However, we explain in appendix D.1 why this can be assumed: because non-Gaussian contributions have a negligible impact for our specific target, $\beta$, to the level of precision we can claim with this method.</footnote>

<footnote n="3">³Note that the parameters that enter in $F_{\mathrm{kin}}(z)$ are not the same as those that contribute to the integrated dipole $\bar{D}_{\mathrm{kin}}$ defined in eq. (1.2). As demonstrated in [17], the integrated dipole depends on magnification bias and on the spectral index of the sources, whereas the dipole at fixed redshift $z$, as well as the correlations between different multipoles, both depend on magnification bias and evolution bias.</footnote>
We start by rewriting eq. (3.6) in vector form

$$ \mathbb{E}[\mathbf{O}] = F_{\text{kin}} \alpha + \beta \gamma, \quad (4.1) $$

where $\mathbf{O}$ is the vector collecting all $O_{\ell,m}$, and the same goes for $\alpha$ and $\gamma$, with

$$ \alpha_{\ell,m} = A_{\ell+1,m} (C_{\ell} + C_{\ell+1}) = \sqrt{\frac{(\ell+1)^2 - m^2}{(2\ell+1)(2\ell+3)}} (C_{\ell} + C_{\ell+1}), \quad (4.2) $$

$$ \gamma_{\ell,m} = -A_{\ell+1,m} ((\ell+2)C_{\ell+1} - \ell C_{\ell}) = -\sqrt{\frac{(\ell+1)^2 - m^2}{(2\ell+1)(2\ell+3)}} ((\ell+2)C_{\ell+1} - \ell C_{\ell}). \quad (4.3) $$

The coefficients $\alpha_{\ell,m}$ and $\gamma_{\ell,m}$ depend on the power spectra $C_{\ell}$. They can be computed either from measured power spectra or from theoretically modelled power spectra. The advantage of using measured power spectra is that we can include highly non-linear scales, where we do not have a reliable theoretical modelling. On the other hand, the theoretical power spectra have the advantage of not being affected by measurement uncertainties. In appendix D.2, we show that the optimal solution is to use a theoretical model at low $\ell$, where the measurement uncertainties are large (due to cosmic variance) and perturbation theory is reliable. At large $\ell$ on the other hand, it is better to use the measured power spectra that are very well measured and whose errors give negligible contributions to the variance of the estimators. In the following we thus neglect the spectra errors and consider them fixed.

There is a scenario where we can easily derive optimal (meaning unbiased and with minimal variance) Quadratic Estimators (QEs): when one of the parameters dominate the observable, i.e.

$$ \text{either} \quad F_{\text{kin}} \alpha_{\ell,m} \ll \beta \gamma_{\ell,m} \quad \text{or} \quad F_{\text{kin}} \alpha_{\ell,m} \gg \beta \gamma_{\ell,m}. \quad (4.4) $$

In that case we estimate only one parameter, the dominant one. This a problem that has been solved in the CMB lensing literature, see e.g. [31]. As shown in appendix B.1 the optimal QE is

$$ \text{either} \quad \hat{\beta} = \frac{1}{\|\gamma\|^2} \langle \mathbf{O}, \gamma \rangle \quad \text{or} \quad \hat{F}_{\text{kin}} = \frac{1}{\|\alpha\|^2} \langle \mathbf{O}, \alpha \rangle, \quad (4.5) $$

where we defined the inner product

$$ \langle \mathbf{x}, \mathbf{y} \rangle = \sum_{\ell_1, m_1; \ell_2, m_2} x_{\ell_1, m_1}^* (\Sigma^{-1})_{\ell_1, m_1; \ell_2, m_2} y_{\ell_2, m_2} = \sum_{\ell, m} \frac{x_{\ell, m}^* y_{\ell, m}}{C_{\ell} C_{\ell+1}} f_{\text{SKY}}, \quad (4.6) $$

with $\Sigma_{\ell_1, m_1; \ell_2, m_2} = \text{Cov}(O_{\ell_1, m_1}, O_{\ell_2, m_2})$ and where $\|\mathbf{x}\|^2 = \langle \mathbf{x}, \mathbf{x} \rangle$.

This scenario, where one parameter dominates, leads to simple estimators but the assumptions are not expected to be verified in practice since $\beta$ and $F_{\text{kin}}$ are generally of the same order of magnitude, and $\alpha$ and $\gamma$ are also of the same order magnitude. One can show that the estimators of eq. (4.4) can still be used (simultaneously) if the weights are orthogonal in the sense of the inner product, i.e.

$$ \text{if} \quad \langle \alpha, \gamma \rangle = 0. \quad (4.7) $$
In that case, both parameters are estimated independently. Unfortunately, condition (4.7) is also not expected to hold in practice. An interesting counter-example is that for $C_{\ell+1} = C_\ell$, we have $\gamma_{\ell,m} = \alpha_{\ell,m}$. So if all $C_\ell$'s are equal (e.g. in the shot-noise dominated regime), then $\beta$ and $F_{\text{kin}}$ are perfectly degenerate. A fortiori they cannot be estimated independently.

Instead, both parameters have to be estimated jointly. This requires generalising the derivation of Quadratic Estimators as explained in appendix B.2. Applying to the current case, we find that the optimal QEs are the estimators $(\hat{F}_{\text{kin}}, \hat{\beta})$ that solve the linear system

$$ \begin{pmatrix} \|\boldsymbol{\alpha}\|^2 & \langle \boldsymbol{\alpha}, \boldsymbol{\gamma} \rangle \\ \langle \boldsymbol{\alpha}, \boldsymbol{\gamma} \rangle & \|\boldsymbol{\gamma}\|^2 \end{pmatrix} \times \begin{pmatrix} \hat{F}_{\text{kin}} \\ \hat{\beta} \end{pmatrix} = \begin{pmatrix} \langle \boldsymbol{O}, \boldsymbol{\alpha} \rangle \\ \langle \boldsymbol{O}, \boldsymbol{\gamma} \rangle \end{pmatrix}. \qquad (4.8) $$

This system is invertible unless the Cauchy-Schwarz inequality is saturated $|\langle \boldsymbol{\alpha}, \boldsymbol{\gamma} \rangle| = \|\boldsymbol{\alpha}\| \|\boldsymbol{\gamma}\|$. In that case $\boldsymbol{\alpha}$ and $\boldsymbol{\gamma}$ are collinear, meaning that $\beta$ and $F_{\text{kin}}$ are degenerate and cannot be measured separately because they have the same impact on the data $\boldsymbol{O}$. We examine the degenerate case in more details in section 4.2.

Inverting eq. (4.8) we obtain

$$ \begin{pmatrix} \hat{F}_{\text{kin}} \\ \hat{\beta} \end{pmatrix} = \frac{1}{1 - r_{\alpha, \gamma}^2} \begin{pmatrix} \|\boldsymbol{\alpha}\|^{-2} & -\frac{r_{\alpha, \gamma}}{\|\boldsymbol{\alpha}\|\|\boldsymbol{\gamma}\|} \\ -\frac{r_{\alpha, \gamma}}{\|\boldsymbol{\alpha}\|\|\boldsymbol{\gamma}\|} & \|\boldsymbol{\gamma}\|^{-2} \end{pmatrix} \times \begin{pmatrix} \langle \boldsymbol{O}, \boldsymbol{\alpha} \rangle \\ \langle \boldsymbol{O}, \boldsymbol{\gamma} \rangle \end{pmatrix}, \qquad (4.9) $$

where we have introduced the correlation coefficient

$$ r_{\alpha, \gamma} \equiv \frac{\langle \boldsymbol{\alpha}, \boldsymbol{\gamma} \rangle}{\|\boldsymbol{\alpha}\| \|\boldsymbol{\gamma}\|}. \qquad (4.10) $$

From eq. (4.9), we see that $\hat{F}_{\text{kin}}$ and $\hat{\beta}$ can be directly inferred from the measurements of the observable $O_{\ell,m}$ and the modelling of the coefficients $\alpha_{\ell,m}$ and $\gamma_{\ell,m}$.

To determine how well $\hat{F}_{\text{kin}}$ and $\hat{\beta}$ can be measured, we need to compute the variance of these estimators, as well as their covariance. Those depend on the variance and covariance of $\langle \boldsymbol{O}, \boldsymbol{\alpha} \rangle$ and $\langle \boldsymbol{O}, \boldsymbol{\gamma} \rangle$. At zeroth order in $\beta$, we obtain

$$ \text{Var}(\langle \boldsymbol{O}, \boldsymbol{\alpha} \rangle) = \|\boldsymbol{\alpha}\|^2, \qquad (4.11) $$

$$ \text{Var}(\langle \boldsymbol{O}, \boldsymbol{\gamma} \rangle) = \|\boldsymbol{\gamma}\|^2, \qquad (4.12) $$

$$ \text{Cov}(\langle \boldsymbol{O}, \boldsymbol{\alpha} \rangle, \langle \boldsymbol{O}, \boldsymbol{\gamma} \rangle) = \langle \boldsymbol{\alpha}, \boldsymbol{\gamma} \rangle. \qquad (4.13) $$

With this, after some algebra, we obtain for the variance and covariance of $\hat{F}_{\text{kin}}$ and $\hat{\beta}$

$$ \text{Var}(\hat{F}_{\text{kin}}) = \frac{1}{\|\boldsymbol{\alpha}\|^2(1 - r_{\alpha, \gamma}^2)}, \qquad (4.14) $$

$$ \text{Var}(\hat{\beta}) = \frac{1}{\|\boldsymbol{\gamma}\|^2(1 - r_{\alpha, \gamma}^2)}, \qquad (4.15) $$

$$ \text{Cov}(\hat{F}_{\text{kin}}, \hat{\beta}) = -\frac{r_{\alpha, \gamma}}{\|\boldsymbol{\alpha}\| \|\boldsymbol{\gamma}\|(1 - r_{\alpha, \gamma}^2)}. \qquad (4.16) $$

The Pearson correlation coefficient for $\hat{F}_{\text{kin}}$ and $\hat{\beta}$ is directly related to $r_{\alpha, \gamma}$

$$ \frac{\text{Cov}(\hat{F}_{\text{kin}}, \hat{\beta})}{\sqrt{\text{Var}(\hat{F}_{\text{kin}}) \times \text{Var}(\hat{\beta})}} = -r_{\alpha, \gamma} = -\frac{\langle \boldsymbol{\alpha}, \boldsymbol{\gamma} \rangle}{\|\boldsymbol{\alpha}\| \|\boldsymbol{\gamma}\|}. \qquad (4.17) $$
From this, we see that the degeneracy between $\hat{F}_{\text{kin}}$ and $\hat{\beta}$ is directly linked to the degeneracy between the fitting templates $\alpha$ and $\gamma$, which increases the error bars on $\hat{F}_{\text{kin}}$ and $\hat{\beta}$ through the denominators of eq. (4.14)–(4.15).

In this section, we have derived analytical expressions for optimal QEs to estimate $\beta$ and $F_{\text{kin}}$. We now want to gain analytical intuition on the potential of these estimators and their meaning.

## 4.2 Simple examples for the estimator and the signal

In this section we consider some simplified cases in order to gain intuition on the potential to measure $\beta$ and $F_{\text{kin}}$ and break their degeneracy.

### 4.2.1 One parameter case

From previous considerations, if $\hat{F}_{\text{kin}} \cdot \alpha_{\ell,m} \gg \hat{\beta} \cdot \gamma_{\ell,m}$ (i.e. if the $F_{\text{kin}}$ effect dominates the equations) or if $|r_{\alpha,\gamma}| \ll 1$ (i.e. the two parameters are decorrelated), then the $\hat{F}_{\text{kin}}$ estimator reduces to

$$ \hat{F}_{\text{kin}} = \frac{1}{\|\boldsymbol{\alpha}\|^2} \langle \mathcal{O}, \boldsymbol{\alpha} \rangle, \quad (4.18) $$

and its variance is given by $\text{Var}(\hat{F}_{\text{kin}})^{-1} = \|\boldsymbol{\alpha}\|^2$. We can input the analytical expressions for $\alpha$ and $\Sigma$, to obtain an analytical estimate for this variance in terms of $\ell_{\text{max}}$

$$ \begin{aligned} \frac{\text{Var}(\hat{F}_{\text{kin}})^{-1}}{f_{\text{SKY}}} &= \sum_{\ell,m} A_{\ell+1,m}^2 \frac{(C_{\ell} + C_{\ell+1})^2}{C_{\ell} C_{\ell+1}} \approx \sum_{\ell,m} A_{\ell+1,m}^2 \times 4 = \frac{4}{3} \sum_{\ell=2}^{\ell_{\text{max}}} (\ell+1) \\ &= \frac{2}{3} (\ell_{\text{max}}^2 + 3\ell_{\text{max}} - 4), \end{aligned} \quad (4.19) $$

where we assumed $C_{\ell} \simeq C_{\ell+1}$ and used eq. (A.15) in appendix A to compute $\sum_m A_{\ell+1,m}^2$. This means that $F_{\text{kin}}$ can be measured with an absolute precision $\sim \sqrt{3/2 f_{\text{SKY}} \ell_{\text{max}}^{-1}}$ in this single parameter case. Since $F_{\text{kin}} \propto \beta \sim 10^{-3}$ we see that we need $\ell_{\text{max}} \sim 10,000$ to reach a relative precision of 10 percent with coming large surveys.

In the opposite limit, when $\hat{F}_{\text{kin}} \cdot \alpha_{\ell,m} \ll \hat{\beta} \cdot \gamma_{\ell,m}$ or if $|r_{\alpha,\gamma}| \ll 1$, the $\hat{\beta}$ estimator reduces to

$$ \hat{\beta} = \frac{1}{|\gamma|^2} \langle \mathcal{O}, \gamma \rangle, \quad (4.20) $$

and its variance follows

$$ \frac{\text{Var}(\hat{\beta})^{-1}}{f_{\text{SKY}}} = \|\gamma\|^2 = \sum_{\ell,m} A_{\ell+1,m}^2 \frac{((\ell+2)C_{\ell+1} - \ell C_{\ell})^2}{C_{\ell} C_{\ell+1}} = \sum_{\ell} \frac{\ell+1}{3} \frac{((\ell+2)C_{\ell+1} - \ell C_{\ell})^2}{C_{\ell} C_{\ell+1}}. \quad (4.21) $$

This can be simplified under two conditions: $C_{\ell} \simeq C_{\ell+1}$, but also $\ell|C_{\ell} - C_{\ell+1}| \ll C_{\ell+1}$, which is a stronger condition. In this case

$$ \frac{\text{Var}(\hat{\beta})^{-1}}{f_{\text{SKY}}} \approx \sum_{\ell,m} A_{\ell+1,m}^2 \times 4 = \frac{4}{3} \sum_{\ell=2}^{\ell_{\text{max}}} (\ell+1) = \frac{2}{3} (\ell_{\text{max}}^2 + 3\ell_{\text{max}} - 4). \quad (4.22) $$
If the conditions are not satisfied (which is usually the case), we have to use a numerical computation of eq. (4.21). In this case, the scaling of the variance with $\ell_{\max}$ is not directly straightforward as it depends on the scaling of $C_{\ell}$ with $\ell$. We will see in section 5 that for a Euclid-like survey, the variance of $\hat{\beta}$ can become significantly larger than that of $\hat{F}_{\text{kin}}$, due to the partial cancellation of $(\ell + 2) C_{\ell+1}$ with $\ell C_{\ell}$ in eq. (4.21).

Eqs. (4.19) and (4.21) hold in the case where one of the term dominates the signal (or when the two parameters are decorrelated). In the general case, the variance of each of the estimator is boosted by the factor $(1 - r_{\alpha, \gamma}^2)^{-1}$. In other words, if we note

$$M \equiv \begin{pmatrix} \|\alpha\|^2 & \langle\alpha, \gamma\rangle \\ \langle\alpha, \gamma\rangle & \|\gamma\|^2 \end{pmatrix} \qquad (4.23)$$

we have that $(M_{i,i})^{-1}$ is the unmarginalised (or single parameter) error bar, while $(M^{-1})_{i,i}$ is the marginalised error bar.

We now turn to simple cases where one can estimate the parameter degeneracy and its meaning.

### 4.2.2 Power law angular power spectrum

As a second example we consider a simple case for the behaviour of the signal, assuming that the angular power spectrum scales as a power law, i.e. $C_{\ell} \propto 1/\ell^n$. Recall the expectation value of the observable $O_{\ell,m}$,

$$E [O_{\ell,m}] = \alpha_{\ell,m} F_{\text{kin}} + \gamma_{\ell,m} \beta \qquad (4.24)$$

$$\alpha_{\ell,m} = A_{\ell+1,m} (C_{\ell} + C_{\ell+1}) \qquad (4.25)$$

$$\gamma_{\ell,m} = -A_{\ell+1,m} \underbrace{((\ell+2) C_{\ell+1} - \ell C_{\ell})}_{=2 \frac{d \ell(\ell+1)}{d \ell(\ell+1)} \sim \frac{1}{\ell} \frac{d \ell^2 C_{\ell}}{d \ell}} \qquad (4.26)$$

For a power law behaviour, we see that at high $\ell$ we have $\alpha_{\ell,m} \simeq -\gamma_{\ell,m} \frac{2}{2-n}$ and thus the $\beta$ and $F_{\text{kin}}$ parameters are essentially degenerate. Thus one can only measure the combination $G = F_{\text{kin}} - \beta \frac{2-n}{2}$. Using eq. (4.21) and a derivation similar to the one presented in eq. (4.19), one finds that $\text{Var}(\hat{G})^{-1} \sim \frac{2}{3} \ell_{\max}^2$.

Three particular values of $n$ are interesting for us. For a typical galaxy sample, at large $\ell$ but in the regime where the clustering signal still dominates over shot-noise, we have $n \simeq 1$, i.e. $\ell \times C_{\ell} \simeq \text{cst}$. In that case, one can measure the combination $F_{\text{kin}} - \beta/2$, with an absolute precision $\sim \sqrt{3/2}/\ell_{\max}$. This is in agreement with the signal to noise ratio presented in [24] (which only considered the $\tilde{D}_{\text{kin}}$ contribution). On the other hand, in the shot-noise dominated regime, we have $C_{\ell} = \text{cst}$, i.e. $n = 0$. In that case, the $\hat{F}_{\text{kin}}$ and $\hat{\beta}$ estimators are perfectly degenerate ($\alpha_{\ell,m} = -\gamma_{\ell,m}$) and one can measure the combination $\hat{F}_{\text{kin}} - \hat{\beta}$. This shows that the noise-dominated part of the spectrum is interesting and should not be thrown away since it allows us to measure this specific combination. In particular, combining the clustering part of the spectrum, where $n = 1$, with the shot-noise dominated part, where $n = 0$, helps us to break the degeneracy between $\hat{F}_{\text{kin}}$ and $\hat{\beta}$. Finally, the case $n = 2$ is also interesting: it corresponds asymptotically to a scale-invariant spectrum with $\ell(\ell+1)C_{\ell} = \text{cst}$. In that case,
the remapping term $\beta$ does not contribute to the signal, and one can only measure $\hat{F}_{\text{kin}}$. This is not surprising: since the remapping term redistributes in an anisotropic way fluctuations with different sizes, it has no impact if the fluctuations of all sizes have the same amplitude.

With the examples above, we see that in several simple cases, $\hat{F}_{\text{kin}}$ and $\hat{\beta}$ are degenerate. To lift this degeneracy, it is important that the $C_{\ell}$'s exhibit some non-trivial structure, like e.g. acoustic peaks, or that they change from one power law to another. In the next section, we generalise the approach to a multi-tracer analysis, or similarly an analysis with multiple redshift bins, and we show that in this case the degeneracy between $\hat{F}_{\text{kin}}$ and $\hat{\beta}$ can be significantly reduced.

## 4.3 Generalisation to multi-tracer or multiple redshift bins

Section 4.1 dealt with the case of a single tracer in a single redshift bin, for which there was only one $F_{\text{kin}}$ to estimate. Modern surveys however observe the galaxy distribution in sufficiently large redshift intervals that several bins are necessary. Since $F_{\text{kin}} = F_{\text{kin}}(z)$ is a function of redshift (due to the redshift evolution of $r(z)$, $\mathcal{H}(z)$ as well as the magnification bias $x(z)$ and the evolution bias $f_{\text{evol}}(z)$, see eq. (2.23)), one may want to estimate it in each redshift bin. Note that this helps to disentangle the intrinsic from the kinematic dipole in the number counts, as shown in [26]. There may also be different surveys observing different sources (spectroscopic galaxies, photometric galaxies, LRGs, Has, ...) for which the quantities $F_{\text{kin}}(z)$ defined in eq. (2.23) differ because the evolution or magnification biases are population-dependent. In contrast, the remapping contribution proportional to $\beta$ is the same for all the samples and all redshifts.

In the following, we consider $n_s$ different galaxy samples. Those can either be different redshift bins, or different tracers, or both. We call $F_{\text{kin}}^j$ with an index $j \in \{1, \dots, n_s\}$ the coefficient $F_{\text{kin}}$ corresponding to the sample $j$. We have therefore $n_s+1$ quantities to estimate: $\beta$ and $F_{\text{kin}}^j$ for $j \in \{1, \dots, n_s\}$. The observables $O_{\ell,m}$ introduced in eq. (3.5) now depend on a pair of probes $(j_1, j_2)$:

$$O_{\ell,m}^{(j_1, j_2)} = a_{\ell,m}(j_1) a_{\ell+1,m}(j_2)^*. \qquad (4.27)$$

Some algebra gives the mean and covariance of the observables

$$\mathbb{E}[O_{\ell,m}^{(j_1, j_2)}] = A_{\ell+1, m} \times [F_{\text{kin}}^{(j_2)} C_{\ell}^{(j_1, j_2)} + F_{\text{kin}}^{(j_1)} C_{\ell+1}^{(j_1, j_2)} - \beta ((\ell+2) C_{\ell+1}^{(j_1, j_2)} - \ell C_{\ell}^{(j_1, j_2)})], \qquad (4.28)$$

$$\text{Cov}(O_{\ell,m}^{(j_1, j_2)}, O_{\ell',m'}^{(j_3, j_4)}) = \frac{C_{\ell}^{(j_1, j_3)} C_{\ell+1}^{(j_2, j_4)}}{f_{\text{SKY}}} \delta_{\ell\ell'} \delta_{mm'} + \mathcal{O}(\beta). \qquad (4.29)$$

Note that while $\mathbb{E}[O_{\ell,m}^{(j_1, j_2)}]$ is real, the measured $O_{\ell,m}$ will in general be a complex number. The covariance is again diagonal in $(\ell, m)$, but in general, not in $(j_1, j_2)$, since different redshift bins can be correlated, as well as different populations of galaxies. It is convenient to write the mean in eq. (4.28) as the Euclidean scalar product between two vectors

$$\mathbb{E}[O_{\ell,m}^{(j_1, j_2)}] = \mathbf{f}_{\ell,m}^{(j_1, j_2)} \cdot \boldsymbol{\theta}, \qquad (4.30)$$
where $\theta \in \mathbb{R}^{n_s+1}$ is the vector of parameters to estimate and for each values of $\ell, m, j_1, j_2$, $f_{\ell,m}^{(j_1,j_2)}$ is an auxiliary vector of $n_s + 1$ elements. They are explicitly given by

$$ \theta_0 = \beta, \quad (f_{\ell,m}^{(j_1,j_2)})_0 = A_{\ell+1,m} \times \left[ -(\ell+2) C_{\ell+1}^{(j_1,j_2)} + \ell C_{\ell}^{(j_1,j_2)} \right], \quad (4.31) $$

$$ \theta_j = F_{\text{kin}}^{(j)}, \quad (f_{\ell,m}^{(j_1,j_2)})_j = A_{\ell+1,m} \times [C_{\ell}^{(j_1,j_2)} \delta_{j,j_2}^K + C_{\ell+1}^{(j_1,j_2)} \delta_{j,j_1}^K], \quad \forall j \in \{1, \dots, n_s\}. \quad (4.32) $$

In the following, we note the index of the vector of parameters with a Greek index, e.g. $\mu, \nu = 0, \dots, n_s$. We consider objects like $f$ as a vector whose elements are the $(f_{\ell,m}^{(j_1,j_2)})_\mu$. In analogy to tensorial calculations, we shall sometimes act with linear algebra operators on the $\mu$ index and sometimes on the $(j_1,j_2)_{\ell,m}$ indices, where the context and dimensionality of the operators should be sufficient to understand unambiguously what is meant.

As an application, we define a new inner product between objects having both $(\ell, m)$ and $(j_1, j_2)$ indices

$$ \langle \mathbf{x}, \mathbf{y} \rangle = \sum_{\ell,m} \left( \sum_{j_1,j_2,j_3,j_4} (x_{\ell,m}^{(j_1,j_2)})^* \operatorname{Cov}(O_{\ell,m}^{(j_1,j_2)}, O_{\ell,m}^{(j_3,j_4)})^{-1} y_{\ell,m}^{(j_3,j_4)} \right) \quad (4.33) $$

$$ = (x^A)^* (\Sigma^{-1})_{AB} y^B, \quad (4.34) $$

where in the second line we used Einstein summation convention and the definition

$$ \Sigma^{AB} = \operatorname{Cov}(O^A, O^B). \quad (4.35) $$

Using this inner product we can define a vector of size $n_s + 1$ that we note $\langle \mathbf{O}, \mathbf{f} \rangle$ and a matrix of shape $(n_s + 1) \times (n_s + 1)$ that we note $\langle \mathbf{f}, \mathbf{f} \rangle$:

$$ \langle \mathbf{O}, \mathbf{f} \rangle_\mu = \langle \mathbf{O}, \mathbf{f}_\mu \rangle = (O^A)^* (\Sigma^{-1})_{AB} f_\mu^B, \quad (4.36) $$

$$ \langle \mathbf{f}, \mathbf{f} \rangle_{\mu, \nu} = \langle \mathbf{f}_\mu, \mathbf{f}_\nu \rangle = (f_\mu^A)^* (\Sigma^{-1})_{AB} f_\nu^B. \quad (4.37) $$

Using this definitions and the derivation of multi-parameters QEs from appendix B.3, we find that optimal QEs for the $\theta$ components is given by the solution to the following linear system:

$$ \langle \mathbf{f}, \mathbf{f} \rangle \cdot \hat{\theta} = \langle \mathbf{O}, \mathbf{f} \rangle. \quad (4.38) $$

The system is invertible unless there is a perfect degeneracy between some parameters, which we found not to happen in practice. Inverting the system leads to the estimators $\hat{\theta} = \langle \mathbf{f}, \mathbf{f} \rangle^{-1} \cdot \langle \mathbf{O}, \mathbf{f} \rangle$. In the two parameters case of section 4.1, we could invert the matrix $\langle \mathbf{f}, \mathbf{f} \rangle$ analytically. This is no longer possible in the generalised case, although there is no obstacle to numerical inversion in practice.

We now compute the covariances of the $n_s + 1$ estimators. To do so, we start by writing the components of the estimator $\hat{\theta}^\mu$

$$ \hat{\theta}^\mu = (\langle \mathbf{f}, \mathbf{f} \rangle^{-1})^{\mu\nu} (O^A)^* (\Sigma^{-1})_{AB} f_\nu^B. \quad (4.39) $$
Using linearity, the covariance between the different estimators $\hat{\theta}^{\mu}$ becomes

$$ 
\begin{aligned}
\mathrm{Cov}(\hat{\theta}^{\mu_1}, \hat{\theta}^{\mu_2}) &= (\langle \mathbf{f}, \mathbf{f} \rangle^{-1})^{\mu_1 \nu_1} \Sigma^{AC} (\Sigma^{-1})_{AB} f_{\nu_1}^B (\langle \mathbf{f}, \mathbf{f} \rangle^{-1})^{\mu_2 \nu_2} (\Sigma^{-1})_{CD} (f_{\nu_2}^D)^* \\
&= \delta_B^C (\Sigma^{-1})_{CD} f_{\nu_1}^B (f_{\nu_2}^D)^* (\langle \mathbf{f}, \mathbf{f} \rangle^{-1})^{\mu_1 \nu_1} (\langle \mathbf{f}, \mathbf{f} \rangle^{-1})^{\mu_2 \nu_2} \\
&= \langle \mathbf{f}, \mathbf{f} \rangle_{\nu_1 \nu_2} (\langle \mathbf{f}, \mathbf{f} \rangle^{-1})^{\mu_1 \nu_1} (\langle \mathbf{f}, \mathbf{f} \rangle^{-1})^{\mu_2 \nu_2} \\
&= (\langle \mathbf{f}, \mathbf{f} \rangle^{-1})^{\mu_1 \mu_2},
\end{aligned}
\quad (4.40) $$

where we have also used the symmetry of the covariance matrix $\Sigma^{AB}$ and of the matrix $\langle \mathbf{f}, \mathbf{f} \rangle$.

In particular the variance of the estimator $\hat{\beta}$ is $(\langle \mathbf{f}, \mathbf{f} \rangle^{-1})^{00}$.

Thus, to forecast the uncertainties on the parameters of interest, we need to compute the inverse of the symmetric matrix $\langle \mathbf{f}, \mathbf{f} \rangle$. This calculation can be slightly simplified by performing analytically the sum over $m$. We use that $\Sigma^{AB}$ is diagonal in $\ell$ and $m$ and independent of $m$, and we denote by $\Sigma_{\ell}$ the matrix $\Sigma_{(\ell m)(\ell m)}^{(j_1, j_3), (j_2, j_4)} = C_{\ell}^{(j_1, j_3)} C_{\ell+1}^{(j_2, j_4)}/f_{\text{SKY}}$. With this, we find

$$ 
\begin{aligned}
\langle \mathbf{f}, \mathbf{f} \rangle_{00} &= \sum_{\ell} \left( -(\ell + 2) \bar{C}_{\ell+1} + \ell \bar{C}_{\ell} \right)^T \cdot \Sigma_{\ell}^{-1} \cdot \left( -(\ell + 2) \bar{C}_{\ell+1} + \ell \bar{C}_{\ell} \right) \underbrace{\sum_m A_{\ell+1, m}^2}_{= \frac{\ell+1}{3}} \\
&= \sum_{\ell} \frac{\ell+1}{3} \left( -(\ell + 2) \bar{C}_{\ell+1} + \ell \bar{C}_{\ell} \right)^T \cdot \Sigma_{\ell}^{-1} \cdot \left( -(\ell + 2) \bar{C}_{\ell+1} + \ell \bar{C}_{\ell} \right),
\end{aligned}
\quad (4.41) $$

$$ 
\begin{aligned}
\langle \mathbf{f}, \mathbf{f} \rangle_{0j} &= \sum_{\ell} \frac{\ell+1}{3} \left( -(\ell + 2) \bar{C}_{\ell+1} + \ell \bar{C}_{\ell} \right)^T \cdot \Sigma_{\ell}^{-1} \cdot \left( \bar{C}_{\ell} \delta_{j,j_4}^K + \bar{C}_{\ell+1} \delta_{j,j_3}^K \right), \\
\forall j \in \{1, \dots, n_s\},
\end{aligned}
\quad (4.42) $$

$$ 
\begin{aligned}
\langle \mathbf{f}, \mathbf{f} \rangle_{ij} &= \sum_{\ell} \frac{\ell+1}{3} \left( \bar{C}_{\ell} \delta_{i,j_2}^K + \bar{C}_{\ell+1} \delta_{i,j_1}^K \right)^T \cdot \Sigma_{\ell}^{-1} \cdot \left( \bar{C}_{\ell} \delta_{j,j_4}^K + \bar{C}_{\ell+1} \delta_{j,j_3}^K \right), \\
\forall i, j \in \{1, \dots, n_s\}.
\end{aligned}
\quad (4.43) $$

where $\bar{C}_{\ell}$ is understood as a vector with size $n_s(n_s + 1)/2$ containing $C_{\ell}^{(j_1, j_2)}$ for all pairs of samples and the indices $(j_1, j_2)$ and $(j_3, j_4)$ are also suppressed in the matrix $\Sigma_{\ell}$.

We now have all the tools to forecast the uncertainties on $\hat{\beta}$ and $\hat{F}_{\text{kin}}^{(j)}$.

## 5 Forecast for future surveys

In this section, we forecast how well a Euclid-like survey may measure $\beta$ and $F_{\text{kin}}$, assuming a fiducial cosmological model. We first consider the case of a single tracer for a single bin. We then show how the measurement of $\beta$ is improved when we combine 10 redshift bins. Finally, we study the impact of prior knowledge of $F_{\text{kin}}$ through a direct measurement of the dipole, and most importantly, the impact of prior knowledge on the source evolution and magnification biases.

We follow the survey specifications of [32] for a photometric Euclid-like survey. We choose our fiducial cosmology to be the flat $\Lambda$CDM model with parameters from Planck [33], namely $(\Omega_b, \Omega_m, h, n_s, \sigma_8, \sum m_{\nu}/[\text{eV}]) = (0.05, 0.32, 0.67, 0.96, 0.816, 0.06)$. The matter power
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Left: forecasted un marginalised absolute errors on ($\beta, F_{\text{kin}}$) from a single tracer analysis using the median redshift bin of a Euclid-like galaxy survey, as a function of $\ell_{\text{max}}$. Center: marginalised error. Right: degeneracy coefficient between the two parameters. All quantities are plotted as a function of the resolution of the analysis, parameterized by the maximal multipole $\ell_{\text{max}}$ included in the data vector. Note that the y-axis of the left and central plots are different, since the marginalisation significantly inflates the errors. In the central plot, the dashed line corresponds to $\beta = 1.23 \cdot 10^{-3}$, which is the expectation from CMB observations. The observer’s peculiar velocity falls short of being detectable with such an analysis.</figcaption></figure>

spectrum $P(k)$ on non-linear scales follows the prescriptions of [34, 35]. We consider a photometric mock galaxy survey of 15,000 deg$^2$ survey with a redshift distribution given by

$$n(z) \propto \left(\frac{z}{z_0}\right)^2 \exp\left[-\left(\frac{z}{z_0}\right)^{3/2}\right], \qquad (5.1)$$

where $z_0 = 0.9/\sqrt{2}$. The normalisation is set by requiring a total density of 30 galaxies per arcmin$^2$. We consider 10 redshift bins with boundaries defined such that they are equi-populated. Their boundaries are given by

$$z_i = (0.001, 0.42, 0.56, 0.68, 0.79, 0.90, 1.02, 1.15, 1.32, 1.58, 2.50).$$

The redshifts of all sources are then perturbed with photometric errors following the specifics from [32]. Galaxy clustering is modelled as a biased tracer of matter with a scale-independent bias $b(z) = \sqrt{1+z}$. The galaxy angular power spectrum is considered on a range of multipoles from $\ell_{\text{min}} = 10$ to $\ell_{\text{max}} = 10,000$. Since we do not have measurements for Euclid power spectra yet, we use theoretical spectra at all scales. Concerning the error budget, we do include cosmic variance and shot noise in our error budget but neglect additional instrumental errors.

## 5.1 Individual redshift bins

We first consider the case of individual redshift bins. In figure 1 we plot the errors on $\hat{\beta}$ and $\hat{F}_{\text{kin}}$ as a function of the maximal multipole $\ell_{\text{max}}$, for the fifth redshift bin of Euclid $z \in [0.79, 0.9]$. This is the median redshift bin, as a representative case. In figure 1, the left panel shows the un marginalised errors, which are the inverse square root of the diagonal elements of the $\langle \mathbf{f}, \mathbf{f} \rangle$ matrix, while the middle panel shows the marginalized errors, which are the diagonal elements of $\langle \mathbf{f}, \mathbf{f} \rangle^{-1}$. We also plot in the right panel, the correlation coefficient defined in
eq. (4.10), which governs the covariance between the estimators $\hat{\beta}$ and $\hat{F}_{\text{kin}}$. Looking at the unmarginalised errors, we see that $\hat{F}_{\text{kin}}$ is better measured than $\hat{\beta}$, especially for $\ell_{\text{max}} \lesssim 1000$. This is due to the partial cancellation of the contribution proportional to $\beta$ in eq. (3.6), when $(\ell+2)C_{\ell+1} \sim \ell C_{\ell}$. This in turns increases the variance of $\hat{\beta}$, as we see from eq. (4.21). In particular, when $n=1$ (clustering signal at large $\ell$), we have $(\ell+2)C_{\ell+1} - \ell C_{\ell} \sim 2/\ell$, leading to a large variance of $\hat{\beta}$. Adding more multipoles helps since it adds a part of the spectrum where shot-noise is non-negligible. In particular, in the regime where shot-noise dominates ($n=0$) we have $(\ell+2)C_{\ell+1} - \ell C_{\ell} \sim 2$, which decreases the variance of $\hat{\beta}$.

From the middle panel of figure 1, we see that the marginalised errors are significantly larger than the unmarginalised one, by almost a factor 10 for high $\ell_{\text{max}}$. This is due to the large degeneracy between $\hat{\beta}$ and $\hat{F}_{\text{kin}}$. The correlation coefficient plotted in the right panel of figure 1 is indeed close to 1 at large $\ell_{\text{max}}$. It is only around the BAO scales, where the $C_{\ell}$'s have a lot of structure, that the two contributions become less correlated. As seen from eqs. (4.14) and (4.15), this strong degeneracy between $\hat{\beta}$ and $\hat{F}_{\text{kin}}$ increases the variance of the estimators and consequently prevents us from detecting a $\beta$ of the order of $10^{-3}$, even if one pushes the analysis to $\ell_{\text{max}} = 10^4$.

In figure 2, we present the same results for the 10 photometric redshift bins of Euclid, when each of them is analysed separately. More precisely, we use eqs. (4.14) and (4.15) to determine the variance on $\hat{F}_{\text{kin}}^{(i)} = \hat{F}_{\text{kin}}(z_i)$ and $\hat{\beta}$ for each bin $i$. In the left panel of figure 2 we show the unmarginalised errors on $\hat{\beta}$ from each of the bins, whereas in the right panel we plot the marginalised errors. For large values of $\ell_{\text{max}}$, the error on $\hat{\beta}$ decreases when the redshift bin increases. This may seem surprising, since the angular power spectrum decreases at high redshift, where galaxy clustering is smaller than at low redshift. However, shot noise is larger at high redshift, where less galaxies are detected. As discussed before, since shot noise is affected by the remapping term, proportional to $\beta$, it helps us in measuring $\beta$. Even though the errors are smaller at high redshift, we see from the right panel of figure 2 that the degeneracy between $\hat{\beta}$ and $\hat{F}_{\text{kin}}^i$ is still too large to allow for a $1\sigma$ detection at $\ell_{\text{max}} = 10,000$.

To conclude, in this section we artificially considered each of the redshift bins individually and found that $\beta$ could reasonably not be detected. However, the remapping contribution proportional to $\beta$ is the same for all redshifts. Therefore, one can hope that combining redshift bins will improve the measurement of $\beta$ and reduce the degeneracies. This is what we turn to in the next section.

## 5.2 Combining 10 redshift bins

To perform a proper combined analysis of the 10 redshift bins, we need to use the estimators $\hat{\theta}$ obtained from solving eq. (4.38), where $\hat{\theta}^0 = \hat{\beta}$ and $\hat{\theta}^i = \hat{F}_{\text{kin}}^{(i)}$. The variance and covariance of these estimators are given in eqs. (4.41) to (4.43). To compute them, we need to account for the covariance of the $C_{\ell}$'s between different redshift bins. In figure 3, we show the unmarginalised (left panel) and marginalised (right panel) errors on $\hat{\beta}$ as well as on $\hat{F}_{\text{kin}}^{(i)}$ for four representative redshift bins. As expected, we see that the unmarginalised error on $\hat{\beta}$ combining the 10 bins are improved by a factor $\sim 3$ with respect to the single bin case. The improvement on the marginalised error on $\hat{\beta}$ is however more significant: the error decreases by a factor $\sim 5$. This is due to the fact that combining different redshift bins helps breaking
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Forecasted absolute error on $\beta$ as a function of $\ell_{\text{max}}$, for the ten redshift bins of an Euclid-like survey, when each bin is analysed separately. We show the un marginalised error in the left panel and the marginalised errors in the right panel. Again, the $\beta = 1.23 \cdot 10^{-3}$ value derived from the CMB dipole falls short of being detectable.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 3.** Forecasted absolute errors on $\beta$ and $F_{\text{kin}}^{(i)}$ as a function of $\ell_{\text{max}}$ for a few redshift bins $i \in \{1, 4, 7, 10\}$ for a multitracer analysis using the full power of the ten redshift bins of the Euclid-like survey. We show the un marginalised error in the left panel and the marginalised errors in the right panel. The CMB value $\beta = 1.23 \cdot 10^{-3}$ would be detected at $2.4\sigma$ at $\ell_{\text{max}} = 10,000$.</figcaption></figure>

the degeneracy between $\hat{\beta}$ and $\hat{F}_{\text{kin}}^{(i)}$ since the contribution from $\hat{\beta}$ to the observable $O_{\ell,m}$ varies with redshift. As a consequence, the errors on $\hat{F}_{\text{kin}}^{(i)}$ and $\hat{\beta}$ all become smaller than $10^{-3}$ around $\ell_{\text{max}} \simeq 10^3$. We therefore obtain a $2.4\sigma$ significance of $\hat{\beta}$ at $\ell_{\text{max}} = 10,000$ if $\beta = 1.23 \cdot 10^{-3}$.
## 5.3 Adding priors on $F_{\text{kin}}$

The multi-tracer analysis presented in the previous section works reasonably well. We now investigate how the constraints improve if there is prior knowledge of some of the quantities involved. At first, we considered modelling the redshift dependence of the coefficients $\hat{F}_{\text{kin}}(z)$ with 2-4 parameters. However, we found that doing so does not improve the marginalized error on $\hat{\beta}$ by more than a few percent.

Another possibility is to add priors on $\hat{F}_{\text{kin}}^{(i)}$ using measurements of the dipole, i.e. the spherical harmonic coefficient $a_{1,0}$ in eq. (3.1). Measuring the dipole bin per bin provides direct measurements of $D_{\text{kin}}^{(i)}$, which is related to $F_{\text{kin}}^{(i)}$ by eq. (2.23). The precision on the measurement of $D_{\text{kin}}^{(i)}$ is limited by shot noise, which scales as $D_{\text{sn}} \approx 3\sqrt{f_{\text{sky}}/N_i}$ where $N_i$ is the number of sources in the bin $i$, see [26]. Current measurements have reached a precision of 10% using more than a million quasars [10]. Since the *Euclid* photometric survey [36] is expected to detect more than $10^9$ sources, we will have, in each redshift bin, about 100 million galaxies. We expect therefore the shot-noise uncertainty to be 10 times smaller than in [10] in each of the redshift bins, leading to a precision on $D_{\text{kin}}^{(i)}$ of the order of 1%.

Note however that the dipole $a_{1,0}$ in eq. (3.1) depends not only on $D_{\text{kin}}^{(i)}$ but also on the intrinsic clustering dipole $\tilde{a}_{1,0}$ that may be important at low redshift $z \le 0.1$ [5, 26] and contaminate the measurements of $D_{\text{kin}}^{(i)}$ in our lowest redshift bin. Furthermore, $F_{\text{kin}}^{(i)}$ contains additional contributions with respect to $D_{\text{kin}}^{(i)}$ (see eq. (2.23)), which are sensitive to the growth rate and the evolution of the galaxy bias. To account for these two effects, we enlarge our priors on $F_{\text{kin}}^{(i)}$, choosing two cases in our analysis: a 10% prior and a 30% prior.

In figure 4, we show the unmarginalised (left panel) and marginalised (right panel) errors on $\hat{\beta}$ and on one of the coefficient $\hat{F}_{\text{kin}}^{(5)}$ combining the 10 redshift bins of *Euclid*. For small values of $\ell_{\text{max}}$, the unmarginalised error on $\hat{F}_{\text{kin}}^{(5)}$ are dominated by the priors, i.e. our method does not bring information for this parameter. However when increasing $\ell_{\text{max}}$, our method gradually decreases the errors down to $2 \times 10^{-4}$. From the right panel of figure 4, we see that the priors on $\hat{F}_{\text{kin}}^{(5)}$ have a strong impact on the marginalised errors on $\hat{\beta}$. Indeed, a prior knowledge of $\hat{F}_{\text{kin}}^{(5)}$ breaks the strong degeneracy between $\hat{F}_{\text{kin}}^{(5)}$ and $\hat{\beta}$, allowing for a significant detection of $\beta$. We see that with a 30% prior, $\beta = 1.23 \cdot 10^{-3}$ is detectable at $3.1\sigma$ at $\ell_{\text{max}} = 10,000$. With a 10% prior, which is still very reasonable, we reach $3\sigma$ significance at $\ell_{\text{max}} = 3,569$, $5\sigma$ significance at $\ell_{\text{max}} = 7,224$, and the maximum significance reached at $\ell_{\text{max}} = 10,000$ is of $6\sigma$.

## 5.4 Adding priors on the magnification bias and evolution bias

Another way to break the degeneracy between $\hat{F}_{\text{kin}}$ and $\hat{\beta}$ is to model directly the quantities that enter into $F_{\text{kin}}$, see eq. (2.23). In particular, the magnification bias, evolution bias, as well as the evolution of the galaxy bias with redshift and luminosity threshold $L_*$ can be measured from the population of galaxies. The comoving distance, Hubble parameter and growth function, that also affect $F_{\text{kin}}$, can be calculated assuming a cosmological model. All this information can then be translated into a prior on the coefficients $A_i \equiv F_{\text{kin}}^{(i)}/\beta$. The error on $A_i$ are expected to be dominated by the errors on the magnification bias, evolution bias and galaxy bias evolution, since the cosmological parameters should be determined
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Forecasted absolute errors on $\beta$ and $F_{\text{kin}}^{(i)}$, using a prior on $F_{\text{kin}}^{(i)}$. We consider a multitracer analysis with 10 redshift bins of a Euclid-like survey and plot the errors as a function of $\ell_{\text{max}}$. As illustration we choose a fiducial value of $F_{\text{kin}}^{(i)} = 4 \cdot \beta$, with the CMB value $\beta = 1.23 \cdot 10^{-3}$, and we consider two different priors: a 10% prior and a 30% prior. We show the un marginalised error in the left panel and the marginalised errors in the right panel.</figcaption></figure>

with percent precision by future surveys. As shown in [37], the evolution and magnification bias can be measured with a precision of the order of 10% using $\sim 10^4$ objects. Realistic expectations for Euclid are similar [28]. In our analysis, we take conservative priors of 30% and 10% on the coefficients $A_i$, $i \in \{1, \dots, 10\}$.

The errors are plotted in figure 5. We see that in this case, the marginalised error on $\beta$ reaches $10^{-4}$ at $\ell_{\text{max}} = 10,000$ for the 30% priors on $A_i$ and $4 \cdot 10^{-5}$ for the 10% priors. In particular, with a 30% prior, a value $\beta \simeq 1.23 \cdot 10^{-3}$ consistent with the CMB is detectable at $5\sigma$ already for $\ell_{\text{max}} = 962$. Using all multipoles up to $\ell_{\text{max}} = 10,000$ we reach a maximum significance of $8.9\sigma$. With a 10% prior we reach $5\sigma$ significance at $\ell_{\text{max}} = 832$, $10\sigma$ significance at $\ell_{\text{max}} = 1,769$, and a maximum significance of $24\sigma$ at $\ell_{\text{max}} = 10,000$. This corresponds to a relative precision of 4.1% in the measurement of $\beta$. The strong improvement brought by the priors on $A_i$ is due to the fact that in this case, not only $\hat{\beta}$ is used to measure $\beta$, but the 10 values of $\hat{F}_{\text{kin}}^{(i)}$ are also used, since we now use the fact that $\hat{F}_{\text{kin}}^{(i)}$ are proportional to $\beta$ with coefficients $A_i$ that we model with some imprecision. This analysis shows therefore that a survey like Euclid should allow us to determine the observer velocity using correlations between $a_{\ell,m}$, with a sufficient precision to unambiguously determine if it is consistent with CMB observations or with quasars and radio surveys.

# 6 Discussion and conclusions

In this work, we revisited the computation of the correlation of the $(\ell, \ell \pm 1)$ harmonic coefficients of the redshift-dependent source number counts. We found a term, which had been previously missed in the literature [24], and which is directly proportional to the boost factor $\beta$, without any dependence on the source population properties. This term comes
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Forecasted absolute errors on $\beta$ and $F_{\text{kin}}^{(i)}$ as a function of $\ell_{\text{max}}$ for a multitracer analysis using the 10 redshift bins of the Euclid-like survey. Here, we add a prior on the amplitude $A_i = F_{\text{kin}}^{(i)}/\beta$, the proportionality factor given in eq. (2.22). Precise knowledge of this factor is mostly limited by external measurements of the galaxy evolution bias and magnification bias.</figcaption></figure>

from the remapping of the angular position of the sources on the sky from aberration of a
moving observer. It is particularly interesting because it opens up the possibility to measure
$\beta$ and the coefficient $F_{\text{kin}}(z)$ independently. Using the $f_{\text{SKY}}$ approximation and assuming
Gaussianity of the galaxies distribution, we developed optimal quadratic estimators for $\beta$
and $F_{\text{kin}}(z)$. We showed how measurements of the $(\ell, \ell \pm 1)$ correlations in a single redshift
bin leave a significant degeneracy between $\hat{\beta}$ and $\hat{F}_{\text{kin}}$. We then developed a multi-tracer
analysis, that can be used with different populations of galaxies or different redshift bins
and showed how this can be used to break the degeneracy.

Applying our technique to specifications of a *Euclid*-like photometric survey, we forecasted
how well $\beta$ can be measured with only one redshift bin and found that it falls short from being
detectable at $1\sigma$, even when using all multipoles up to $\ell_{\text{max}} = 10,000$. We then combined
information from the 10 redshift bins of a *Euclid*-like photometric sample and showed that,
in this case, one can detect $\beta$ at nearly $3\sigma$ for a fiducial value $\beta = 1.23 \cdot 10^{-3}$, as expected
from the CMB dipole. We then showed that the detectability of $\beta$ can be strongly improved
by assuming some prior knowledge on the coefficients $F_{\text{kin}}^{(i)}$ in each redshift bins. In particular,
assuming a 10% prior on $F_{\text{kin}}^{(i)}$ leads to a $6\sigma$ detection of $\beta$ at $\ell_{\text{max}} = 10,000$. Even more
promising, prior knowledge of the evolution and magnification biases allows us to model the
coefficients $F_{\text{kin}}^{(i)}$ as: $F_{\text{kin}}^{(i)} = A_i \beta$, and impose priors on the coefficients $A_i$. In this case, we
showed that $\beta$ can be measured at $24\sigma$ for a 10% prior on $A_i$, which is actually a conservative
assumption [38]. The results are summarized in figure 6.

These potential measurements of $\beta$ with *Euclid*, using $(\ell, \ell \pm 1)$ correlations in the galaxy
number counts, will be paramount to better understand the tension between the CMB dipole
and the dipole in radio galaxies and quasars. In particular, since the CatWISE quasar dipole
measurement is roughly a factor of two larger than the expectation from the CMB dipole [10],
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Summary plot: marginalised errors on $\beta$ as a function of $\ell_{\text{max}}$, depending on the type of multi-tracer analysis. Without any prior of the kinematic dipole, the constraint on $\beta$ is the weakest but remains detectable at nearly $3\sigma$ at $\ell_{\text{max}} = 10,000$ with the CMB value $\beta = 1.23 \cdot 10^{-3}$. If there is 30% or 10% prior knowledge of the kinematic dipole, through a direct measurement of the dipole, the constraints improve to $3.1\sigma$ and to $6\sigma$ at $\ell_{\text{max}} = 10,000$, respectively. In the situation where one knows to 30% or 10% precision the proportionality factor of eq. (2.21), i.e. $A_i = F_{\text{kin}}^{(i)}/\beta$, the constraints are excellent and one can measure $\beta = 1.23 \cdot 10^{-3}$ with $8.9\sigma$ or $24\sigma$.</figcaption></figure>

it can be interpreted as an observer's peculiar velocity measurement which is twice as large as the one from the CMB. That means that an $X\sigma$ measurement of $\beta$ can rule out one or the other interpretation at $X\sigma$, if centered on one or the other value. This technique is thus very likely to provide the most robust and competitive constraint on the observer velocity from source number counts, confirming or not the validity of the cosmological principle.

To reach this conclusion, we made a few assumptions which are worth discussing. First, we used the popular $f_{\text{SKY}}$ approximation to account for partial sky coverage. Techniques for a more sophisticated treatment have already been developed for the analysis of CMB lensing and can be imported to the our case. We give guidelines on what needs to be done in appendix C and leave the implementation on data to future work. Since most information comes from relatively high multipoles, which are least affected by partial sky coverage, we expect only small changes. Second, the measurements of $\hat{\beta}$ and $\hat{F}_{\text{kin}}^{(i)}$ require an input for the power spectrum $C_\ell$ up to high multipoles. In our forecasts we assumed that this was known without error. This is not the case in practice since we do not (yet) possess precise cosmological modeling of this power spectrum down to the non-linear regime. However, in appendix D.2 we show that this is not an issue: to the level of precision we find in this article, one can use the measured power spectrum itself (or even better any effective model, like a spline interpolation in log-log) with negligible impact on the bias and error bar of the
QE. Our analysis therefore holds even when pushing to high multipoles. Third, we used a Gaussian covariance matrix (eq. (3.8)), which might seem problematic at high multipoles in the non-linear regime. However, in appendix D.1, we show that this is not an issue either: the impact of non-Gaussian covariance terms is not so much a question of multipole but instead of signal to noise. With the detection significance figures that we claim in this article we show that non-Gaussianity would have a negligible impact.

To conclude, our analysis shows that the $(\ell, \ell \pm 1)$ correlations of the harmonic coefficients of the number counts, that can be measured with a Euclid-like survey, are extremely powerful to infer the amplitude of the observer velocity $\beta$ in an unbiased and robust way.

## Acknowledgments

We thank Sylvain Gouyou Beauchamps for providing us with the predictions for $C_{\ell}^{\text{gal}}$ used in the numerical computation of this article. FL and CB acknowledge support from the Boninchi foundation through the project CosmoBoost. CB further acknowledges financial support from the Swiss National Science Foundation and from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (Grant agreement No. 863929; project title “Testing the law of gravity with novel large-scale structure observables”). C.D. is supported by ERC Starting Grant SHADE (grant no. StG 949572). RD acknowledges support from the Swiss National Science Foundation.

The calculations in this article made use of `astropy` [39, 40], `scipy` [41], ` numpy` [42], and ` matplotlib` [43].

# A Auxiliary calculations on spherical harmonics

## A.1 Dipole times spherical harmonic

Using $Y_{10} = \sqrt{3/4\pi} \cos\theta$, the following function appears in eq. (2.29) and (2.30)

$$f(\mathbf{n}) = \sqrt{\frac{4\pi}{3}} Y_{1,0}^{*}(\mathbf{n}) Y_{\ell, m}^{*}(\mathbf{n}). \qquad (A.1)$$

It is a function defined on the sphere, so we can decompose it in spherical harmonics. More precisely,

$$f(\mathbf{n}) = \sum_{\ell', m'} a_{\ell', m'}^{f} Y_{\ell', m'}(\mathbf{n}), \qquad (A.2)$$

with

$$a_{\ell', m'}^{f} = \int_{S_2} d^2\mathbf{n} \, f(\mathbf{n}) Y_{\ell', m'}^{*}(\mathbf{n}) = \sqrt{\frac{4\pi}{3}} \int_{S_2} d^2\mathbf{n} Y_{\ell,m}^{*}(\mathbf{n}) Y_{\ell',m'}^{*}(\mathbf{n}) Y_{1,0}^{*}(\mathbf{n}) = \sqrt{\frac{4\pi}{3}} G_{m,m',0}^{\ell,\ell',1} \qquad (A.3)$$

where $G_{m_1,m_2,m_3}^{\ell_1,\ell_2,\ell_3}$ is a Gaunt coefficient. The triangular conditions imply that $a_{\ell',m'}^{f}$ is nonzero only if

$$m' = -m, \quad \text{and} \quad \ell' = \ell \pm 1. \qquad (A.4)$$
In this case, we get

$$a_{\ell-1, -m}^{f} = \sqrt{\frac{4\pi}{3}} G_{m, -m, 0}^{\ell, \ell-1, 1} = \sqrt{\frac{4\pi}{3}} \sqrt{\frac{(2\ell-1)(2\ell+1)3}{4\pi}} \begin{pmatrix} \ell & \ell-1 & 1 \\ 0 & 0 & 0 \end{pmatrix} \begin{pmatrix} \ell & \ell-1 & 1 \\ m & -m & 0 \end{pmatrix}, \quad (A.5)$$

$$a_{\ell+1, -m}^{f} = \sqrt{\frac{4\pi}{3}} G_{m, -m, 0}^{\ell, \ell+1, 1} = \sqrt{\frac{4\pi}{3}} \sqrt{\frac{(2\ell+1)(2\ell+3)3}{4\pi}} \begin{pmatrix} \ell & \ell+1 & 1 \\ 0 & 0 & 0 \end{pmatrix} \begin{pmatrix} \ell & \ell+1 & 1 \\ m & -m & 0 \end{pmatrix}. \quad (A.6)$$

Here, we use Wigner symbol formulas from [44] and find that

$$\begin{pmatrix} \ell & \ell-1 & 1 \\ 0 & 0 & 0 \end{pmatrix} = (-1)^{\ell} \sqrt{\frac{\ell}{(2\ell-1)(2\ell+1)}}, \quad (A.7)$$

$$\begin{pmatrix} \ell & \ell-1 & 1 \\ m & -m & 0 \end{pmatrix} = (-1)^{\ell+m} \sqrt{\frac{\ell^2 - m^2}{\ell(2\ell-1)(2\ell+1)}}, \quad (A.8)$$

$$\begin{pmatrix} \ell & \ell+1 & 1 \\ 0 & 0 & 0 \end{pmatrix} = (-1)^{\ell+1} \sqrt{\frac{\ell+1}{(2\ell+1)(2\ell+3)}}, \quad (A.9)$$

$$\begin{pmatrix} \ell & \ell+1 & 1 \\ m & -m & 0 \end{pmatrix} = (-1)^{\ell+1+m} \sqrt{\frac{(\ell+1)^2 - m^2}{(\ell+1)(2\ell+1)(2\ell+3)}}, \quad (A.10)$$

which gives

$$\begin{aligned} a_{\ell-1, -m}^{f} &= (-1)^{m} \sqrt{(2\ell-1)(2\ell+1)} \sqrt{\frac{\ell}{(2\ell-1)(2\ell+1)}} \sqrt{\frac{\ell^2 - m^2}{\ell(2\ell-1)(2\ell+1)}} \\ &= (-1)^{m} \sqrt{\frac{\ell^2 - m^2}{(2\ell-1)(2\ell+1)}}, \end{aligned} \quad (A.11)$$

$$\begin{aligned} a_{\ell+1, -m}^{f} &= (-1)^{m} \sqrt{(2\ell+1)(2\ell+3)} \sqrt{\frac{\ell+1}{(2\ell+1)(2\ell+3)}} \sqrt{\frac{(\ell+1)^2 - m^2}{(\ell+1)(2\ell+1)(2\ell+3)}} \\ &= (-1)^{m} \sqrt{\frac{(\ell+1)^2 - m^2}{(2\ell+1)(2\ell+3)}}. \end{aligned} \quad (A.12)$$

Hence, we can write

$$\begin{aligned} \sqrt{\frac{4\pi}{3}} Y_{1,0}^{*}(\boldsymbol{n}) Y_{\ell,m}^{*}(\boldsymbol{n}) &= f(\boldsymbol{n}) = (-1)^{m} \sqrt{\frac{\ell^2 - m^2}{(2\ell-1)(2\ell+1)}} Y_{\ell-1, -m} \\ &\quad + (-1)^{m} \sqrt{\frac{(\ell+1)^2 - m^2}{(2\ell+1)(2\ell+3)}} Y_{\ell+1, -m} \\ &= A_{\ell,m} Y_{\ell-1,m}^{*}(\boldsymbol{n}) + A_{\ell+1,m} Y_{\ell+1,m}^{*}(\boldsymbol{n}), \end{aligned} \quad (A.13)$$

where we introduced

$$A_{\ell,m} \equiv \sqrt{\frac{\ell^2 - m^2}{(2\ell-1)(2\ell+1)}}. \quad (A.14)$$
This result is used to obtain eq. (2.33) from eq. (2.29). The coefficients $A_{\ell,m}$ obey the following sum rule, which will prove useful

$$ 
\begin{aligned}
\sum_{m=-\ell}^{\ell} A_{\ell+1,m}^2 &= \sum_{m=-\ell}^{\ell} \frac{(\ell+1)^2 - m^2}{(2\ell+1)(2\ell+3)} = \frac{(\ell+1)^2(2\ell+1) - \ell(\ell+1)(2\ell+1)/3}{(2\ell+1)(2\ell+3)} \\
&= \frac{1}{3} \frac{(\ell+1)(2\ell+1)(2\ell+3)}{(2\ell+1)(2\ell+3)} = \frac{\ell+1}{3}.
\end{aligned} 
\quad (A.15) $$

## A.2 Derivative of spherical harmonic

The spherical harmonics can be defined in terms of the associated Legendre polynomials $P_{\ell}^{m}$ as

$$ Y_{\ell,m}(\theta, \phi) = \sqrt{\frac{2\ell+1}{4\pi}} \frac{(\ell-m)!}{(\ell+m)!} e^{im\phi} P_{\ell}^{m}(\cos\theta) \quad (A.16) $$

where the normalisation is chosen for orthonormality ($\int_{S_2} d^2\hat{n} Y_{\ell,m} Y_{\ell'm'}^* = \delta_{\ell,\ell'} \delta_{m,m'}$). We want to compute the derivative $\sin\theta \partial_{\theta} Y_{\ell,m}(\theta, \phi)$, which appears in eq. (2.30). To do so, we use a recurrence formula for the associated Legendre polynomials:

$$ (1-x^2) \frac{d}{dx} P_{\ell}^{m}(x) = \frac{1}{2\ell+1} \left[ (\ell+1)(\ell+m) P_{\ell-1}^{m}(x) - \ell(\ell+1-m) P_{\ell+1}^{m}(x) \right] \quad (A.17) $$

Setting $x = \cos\theta$ this gives

$$ -\sin\theta \partial_{\theta} P_{\ell}^{m}(\cos\theta) = \frac{1}{2\ell+1} \left[ (\ell+1)(\ell+m) P_{\ell-1}^{m}(\cos\theta) - \ell(\ell+1-m) P_{\ell+1}^{m}(\cos\theta) \right] \quad (A.18) $$

Thus

$$ \sin\theta \partial_{\theta} Y_{\ell,m} = \frac{-1}{2\ell+1} \left[ \sqrt{\frac{2\ell+1}{4\pi}} \frac{(\ell-m)!}{(\ell+m)!} (\ell+1)(\ell+m) e^{im\phi} P_{\ell-1}^{m}(\cos\theta) - \sqrt{\frac{2\ell+1}{4\pi}} \frac{(\ell-m)!}{(\ell+m)!} \ell(\ell+1-m) e^{im\phi} P_{\ell+1}^{m}(\cos\theta) \right] \quad (A.19) $$

$$ = \frac{1}{2\ell+1} \left[ -\sqrt{\frac{2\ell+1}{2\ell-1}} \frac{\ell-m}{\ell+m} (\ell+1)(\ell+m) Y_{\ell-1,m} + \sqrt{\frac{2\ell+1}{2\ell+3}} \frac{\ell+1+m}{\ell+1-m} \ell(\ell+1-m) Y_{\ell+1,m} \right] \quad (A.20) $$

$$ = -\sqrt{\frac{(\ell+1)^2 (\ell^2 - m^2)}{(2\ell-1)(2\ell+1)}} Y_{\ell-1,m} + \sqrt{\frac{\ell^2 ((\ell+1)^2 - m^2)}{(2\ell+1)(2\ell+3)}} Y_{\ell+1,m}, \quad (A.21) $$

and taking the complex conjugate, we find

$$ \sin\theta \partial_{\theta} Y_{\ell,m}^* = -\sqrt{\frac{(\ell+1)^2 (\ell^2 - m^2)}{(2\ell-1)(2\ell+1)}} Y_{\ell-1,m}^* + \sqrt{\frac{\ell^2 ((\ell+1)^2 - m^2)}{(2\ell+1)(2\ell+3)}} Y_{\ell+1,m}^* \quad (A.22) $$

which can be simplified using $A_{\ell,m}$'s introduced in eq. (2.34):

$$ \sin\theta \partial_{\theta} Y_{\ell,m}^{*}(\mathbf{n}) = -(\ell+1) A_{\ell,m} Y_{\ell-1,m}^{*}(\mathbf{n}) + \ell A_{\ell+1,m} Y_{\ell+1,m}^{*}(\mathbf{n}). \quad (A.23) $$

We use this expression to calculate (2.30).
## B Derivation of the quadratic estimators

We recall that the off-diagonal correlators defined in eq. (3.5) have the expectation value given in eq. (4.1),

$$ \mathbb{E}[O_{\ell,m}(F_{\text{kin}}, \beta)] = F_{\text{kin}} \alpha_{\ell,m} + \beta \gamma_{\ell,m}. \quad (B.1) $$

The boost effect impacts linearly the (off-diagonal) correlator between two harmonic coefficients $a_{\ell,m}$. It is thus natural to attempt to build an estimator for $(F_{\text{kin}}, \beta)$ that is quadratic in the $a_{\ell,m}$'s. This motivation is reinforced by the striking analogy with CMB lensing, which impacts the same off-diagonal correlators (for CMB temperature and polarisation instead of galaxy counts) and for which Quadratic Estimators (QEs) have been developed with care in the literature, for instance in [31].

### B.1 Single parameter case

First, we assume that one parameter dominates in eq. (B.1), e.g. $F_{\text{kin}} \alpha_{\ell,m} \ll \beta \gamma_{\ell,m}$, so that we have a single parameter to estimate: $\beta$. This is the case where the derivation is the most similar to [31], whereas we will need to adapt calculations in later sections. We first write a general form for the quadratic estimator:

$$ \hat{\beta} = \frac{1}{N} \sum_{\ell,m} g_{\ell,m} (O_{\ell,m})^*. \quad (B.2) $$

The normalisation $N$ and the weights $g_{\ell,m}$ have to be optimised for the estimator to fulfill two conditions: (i) be unbiased, and (ii) have minimal variance.

Condition (i), i.e. $\mathbb{E}[\hat{\beta}] = \beta$, sets the normalisation as

$$ N = \sum_{\ell,m} g_{\ell,m} (\gamma_{\ell,m})^*. \quad (B.3) $$

Condition (ii) is fulfilled by requiring

$$ \frac{\partial \ln \text{Var}(\hat{\beta})}{\partial g_{\ell,m}} = 0 \quad \forall (\ell, m), \quad (B.4) $$

which gives after some algebra

$$ \sum_{\ell',m'} g_{\ell',m'} \Sigma_{\ell,m;\ell',m'} = \gamma_{\ell,m} \times \text{cst}, \quad (B.5) $$

where $\Sigma_{\ell,m;\ell',m'} \equiv \text{Cov}(O_{\ell,m}, O_{\ell',m'})$. After remarking that a constant multiplicative factor in $g_{\ell,m}$ cancels out with the normalisation in the estimator, we see that we can set

$$ g_{\ell,m} = \sum_{\ell',m'} \Sigma_{\ell,m;\ell',m'}^{-1} \gamma_{\ell',m'}, \quad (B.6) $$

where $\Sigma^{-1}$ is the inverse matrix. This also sets the normalisation as

$$ N = \sum_{\ell,m;\ell',m'} \gamma_{\ell,m} \Sigma_{\ell,m;\ell',m'}^{-1} (\gamma_{\ell',m'})^*. \quad (B.7) $$
All these notations are simplified by introducing the inner product between objects with indices $(\ell, m)$

$$ \langle x, y \rangle \equiv \sum_{\ell, m; \ell', m'} x_{\ell, m} \Sigma_{\ell, m; \ell', m'}^{-1} (y_{\ell', m'})^*, \quad (B.8) $$

which yields a compact equation for the estimator:

$$ \hat{\beta} = \frac{1}{\|\gamma\|^2} \langle \mathbf{O}, \gamma \rangle, \quad (B.9) $$

where $\|\gamma\|^2 = \langle \gamma, \gamma \rangle$ is the norm associated with the inner product.

A posteriori, we can make the whole derivation simpler by starting with the following general form for the quadratic estimator, instead of eq. (B.2)

$$ \hat{\beta} = \frac{1}{N} \langle \mathbf{O}, g \rangle. \quad (B.10) $$

In that case, the condition of unbiasedness reads $N = \langle \gamma, g \rangle$. Some algebra gives the variance as

$$ \mathrm{Var}(\hat{\beta}) = \frac{1}{N^2} \langle g, g \rangle, \quad (B.11) $$

and minimising it yields $g \propto \gamma$ which gives again eq. (B.9).

## B.2 Two parameters case

In this case, we do not assume that any of the two parameters has a dominant impact on the observable, so we need to estimate them both jointly.

An easy scenario is when the templates $\alpha$ and $\gamma$ are orthogonal: $\langle \alpha, \gamma \rangle = 0$. In that case, one sees that the condition of unbiasedness can be applied separately to each parameter. Thus by minimizing each variance we have two independent single-parameter QEs:

$$ \hat{F}_{\mathrm{kin}} = \frac{1}{\|\alpha\|^2} \langle \mathbf{O}, \alpha \rangle \qquad \hat{\beta} = \frac{1}{\|\gamma\|^2} \langle \mathbf{O}, \gamma \rangle. \quad (B.12) $$

Then these estimators trivially solve the following linear system

$$ M \times \begin{pmatrix} \hat{F}_{\mathrm{kin}} \\ \hat{\beta} \end{pmatrix} = \begin{pmatrix} \langle \mathbf{O}, \alpha \rangle \\ \langle \mathbf{O}, \gamma \rangle \end{pmatrix} \quad \text{with} \quad M = \begin{pmatrix} \|\alpha\|^2 & \langle \alpha, \gamma \rangle \\ \langle \alpha, \gamma \rangle & \|\gamma\|^2 \end{pmatrix}. \quad (B.13) $$

The more complex scenario is the general one where $\langle \alpha, \gamma \rangle \neq 0$. There, the condition of unbiasedness for each parameter yields equations which are not separate anymore. Furthermore, trying to minimise the variance of each estimator separately is a bad idea because it does not account for the covariance between the two estimators, so we can end up with two estimators which are heavily correlated and thus poor marginalised error bars on $\beta$. Instead, the trick is to use transformations to reduce the computation to the independent case.

Eq. (B.1) can be rewritten in vector form as $\mathbb{E}[\mathbf{O}] = F_{\mathrm{kin}}\alpha + \beta\gamma$. Thanks to the previous paragraph, we would know how to define optimal estimators if we instead had

$$ \mathbb{E}[\mathbf{O}] = A_1 e_1 + A_2 e_2 \quad \text{with} \quad \langle e_1, e_2 \rangle = 0. \quad (B.14) $$
Indeed the estimators would be $\hat{A}_i = \frac{1}{\|e_i\|^2} \langle \mathbf{O}, e_i \rangle$. It is possible to find such vectors ($e_1, e_2$) by orthogonalising the base ($\alpha, \gamma$). There are (infinitely) many possible choices for this orthogonalisation, but there is one that leads to simpler calculations and has the added benefit of normalisation ($\|e_i\| = 1$): take as transfer matrix the square root of $M$. Indeed, $M$ is a real symmetric matrix so it admits a square root: a real symmetric matrix $P$ such that $M = P^2 = PP^T$.

Now we define ($A_1, A_2$) and ($e_1, e_2$) such that

$$ \begin{pmatrix} A_1 \\ A_2 \end{pmatrix} = P \times \begin{pmatrix} F_{\text{kin}} \\ \beta \end{pmatrix} \quad \text{and} \quad \begin{pmatrix} e_1 \\ e_2 \end{pmatrix} = P^{-1} \times \begin{pmatrix} \alpha \\ \gamma \end{pmatrix}. \qquad (B.15) $$

With a little bit of algebra, one can check that eqs. (B.14) are verified by this choice.

Following the previous paragraphs, the optimal quadratic estimators for $A_i$ are $\hat{A}_i = \frac{1}{\|e_i\|^2} \langle \mathbf{O}, e_i \rangle$. We see that they solve (trivially) the linear system

$$ \mathbb{1}_{2 \times 2} \times \begin{pmatrix} \hat{A}_1 \\ \hat{A}_2 \end{pmatrix} = \begin{pmatrix} \langle \mathbf{O}, e_1 \rangle \\ \langle \mathbf{O}, e_2 \rangle \end{pmatrix} \quad \text{with} \quad \mathbb{1}_{2 \times 2} = \begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix}. \qquad (B.16) $$

Now we can define natural estimators for $\beta$ and $F_{\text{kin}}$ as

$$ \begin{pmatrix} \hat{F}_{\text{kin}} \\ \hat{\beta} \end{pmatrix} = P^{-1} \times \begin{pmatrix} \hat{A}_1 \\ \hat{A}_2 \end{pmatrix}. \qquad (B.17) $$

If we input this equation into eq. (B.16), some linear algebra shows that these estimators obey the same linear system as in the independent case, just with $\langle \alpha, \gamma \rangle \neq 0$:

$$ \begin{pmatrix} \|\alpha\|^2 & \langle \alpha, \gamma \rangle \\ \langle \alpha, \gamma \rangle & \|\gamma\|^2 \end{pmatrix} \times \begin{pmatrix} \hat{F}_{\text{kin}} \\ \hat{\beta} \end{pmatrix} = \begin{pmatrix} \langle \mathbf{O}, \alpha \rangle \\ \langle \mathbf{O}, \gamma \rangle \end{pmatrix}. \qquad (B.18) $$

This corresponds to eq. (4.8), the result we were looking for: a definition of the estimators that does not refer to $P$ anymore and can be entirely computed from the observable $\mathbf{O}$ and the galaxy power spectrum $C_\ell$ (the latter defining $\alpha$ and $\gamma$).

## B.3 Multi parameter case

This section is essentially a reskin of section B.2 with the same derivation but with new notations that allow for an arbitrary number of $F_{\text{kin}}$ parameters, as needed for the multi-tracer case.

We start with a generic equation for how the observable $\mathbf{O}$ relates to the parameters we want to estimate

$$ \mathrm{E}[\mathbf{O}^A] = \mathbf{f}^A \cdot \boldsymbol{\theta} = \mathbf{f}_{\mu}^A \boldsymbol{\theta}^{\mu}, \qquad (B.19) $$

where we used Einstein summation convention and the position of indices (up or down) is purely cosmetic. In the equation above, $\boldsymbol{\theta}$ is the vector of parameters, indexed with the Greek letter $\mu$. For example in the main text when we apply this to $n_s$ galaxy samples, we
will have $\mu = 0 \cdots n_s$, $\theta_0 = \beta$ and $\theta_j = F_{\text{kin}}^{(j)} \forall j = 1 \cdots n_s$. The observable $O$ is indexed with a capital letter $A$. For instance with a single sample ($n_s = 1$), as in section B.2, $A$ would be the multi-index of spherical harmonics: $A = (\ell, m)$. In the multi-sample case, the observable can be measured for all possible pairs of samples ($j_1, j_2$), so $A$ is a bigger multi-index: $A = (\ell, m, j_1, j_2)$. Finally, $f^A$ dictates how the parameters impact the observable; for instance, in the case of section B.2, we have $f^{\ell,m} = \begin{pmatrix} \alpha_{\ell,m} \\ \gamma_{\ell,m} \end{pmatrix}$.

We can now transcribe the derivation from the previous sections. First, we note $\Sigma_{A,B} \equiv \text{Cov}(O^A, O^B)$ the covariance matrix. This allows to generalise the inner product as

$$ \langle \mathbf{x}, \mathbf{y} \rangle = (x^A)^* (\Sigma^{-1})_{AB} y^B. \qquad (B.20) $$

We can use this inner product to define a matrix which generalises the $M$ of section B.2:

$$ M_{\mu,\nu} \equiv \langle \mathbf{f}_{\mu}, \mathbf{f}_{\nu} \rangle = (f_{\mu}^{A})^{*} (\Sigma^{-1})_{AB} f_{\nu}^{B}. \qquad (B.21) $$

In the following, we note simply $M \equiv \langle \mathbf{f}, \mathbf{f} \rangle$. This is a symmetric positive matrix so it admits a square root $P$. We use $P$ to define new vectors and tensors:

$$ \mathbf{A} = P \times \boldsymbol{\theta} \quad \text{and} \quad \mathbf{e} = P^{-1} \times \mathbf{f}, \qquad (B.22) $$

or more explicitly in components

$$ \mathbf{A}^{\mu} = P^{\mu}_{\nu} \boldsymbol{\theta}^{\nu} \quad \text{and} \quad \mathbf{e}_{\mu}^{A} = (P^{-1})^{\nu}_{\mu} \mathbf{f}_{\nu}^{A}. \qquad (B.23) $$

These tensors obey

$$ \mathbb{E}[\mathbf{O}] = \mathbf{A} \times \mathbf{e}, \quad \text{i.e.} \quad \mathbb{E}[\mathbf{O}^A] = \mathbf{A}^\mu \mathbf{e}_\mu^A, \quad \text{with} \quad \langle \mathbf{e}_\mu, \mathbf{e}_\nu \rangle = \delta_{\mu,\nu}^K. \qquad (B.24) $$

Since the $e$ basis is orthonormal, we have optimal independent estimators for the components of the vector of parameters $\mathbf{A}: \hat{\mathbf{A}}^{\mu} = \langle \mathbf{O}, \mathbf{e}_{\mu} \rangle$. These estimators trivially solve the linear system

$$ 1 \times \hat{\mathbf{A}} = \langle \mathbf{O}, \mathbf{e} \rangle. \qquad (B.25) $$

When we define natural estimators for $\boldsymbol{\theta}$ as

$$ \hat{\boldsymbol{\theta}} = P^{-1} \times \hat{\mathbf{A}}, \qquad (B.26) $$

and we input this definition into eq. (B.25), some linear algebra shows that these estimators obey

$$ \langle \mathbf{f}, \mathbf{f} \rangle \cdot \hat{\boldsymbol{\theta}} = \langle \mathbf{O}, \mathbf{f} \rangle, \qquad (B.27) $$

which corresponds to eq. (4.38), the result that we were looking for: a definition of optimal estimators solely in terms of the starting quantities $\mathbf{O}$ and $\mathbf{f}$.
## C Quadratic estimator and partial sky coverage

In this appendix we show how to generalise the quadratic estimators derived in the main text in order to account for the survey partial sky coverage. To do so, we follow the formalism developed for Quadratic Estimators of CMB lensing, which translates straightforwardly to our case with minor modifications. Specifically we use methods and notations similar to [45]. The process is composed of two steps: (i) Wiener filtering, and (ii) construction of the QE and mean field subtraction. The steps are explained below, with a summary at the end.

### (i) Wiener filtering

We call $d$ the data, i.e. the masked map. It is considered as a vector of length set by the number of pixels. We call $s$ the signal, i.e. the full sky galaxy map, and $S$ its covariance (which is diagonal in harmonic space). Following [45], we then write

$$d = s + n, \qquad (C.1)$$

where $n$ is a fictious independent noise which has infinite variance ($N^{-1} = 0$) in masked pixels. The goal of this fictious noise is to downweight to zero all masked pixels in the Wiener filtering. $N$ is diagonal in pixel space: $N_{p,p'} = \delta_{p,p'}^K N_p$ where $p$ is the pixel index. What value should we give to $N_p$ in observed (non-masked) pixels? In principle it should be set to zero, but in practice it is useful to set it to a small value $\epsilon > 0$ so that we can invert and get $N^{-1}$ in all pixels. Then our Wiener-filtered data $d^{WF}$ (ref. [45] calls it $\bar{s}$) is the inverse covariance weighting of $d$:

$$d^{WF} = Wd, \quad \text{where} \quad W = [S+N]^{-1} = S^{-1}[S^{-1} + N^{-1}]^{-1}N^{-1}, \qquad (C.2)$$

where the rewriting on the right hand side allows to only need $N^{-1}$ and not $N$, as the latter formally has infinite values. The first multiplication $N^{-1}d$ downweights masked pixels to zero, which is what we want.

### (ii) Construction of the QE and mean field subtraction

We should now build our estimator only in terms of $d^{WF}$. So our observable is

$$O_{\ell,m} = d_{\ell,m}^{WF} (d_{\ell+1,m}^{WF})^*. \qquad (C.3)$$

Its expectation value is

$$\mathbb{E}[O_{\ell,m}] = \langle (Wd)_{\ell,m} (Wd)_{\ell+1,m}^* \rangle \qquad (C.4)$$

$$= \mathbb{E}[(Ws)_{\ell,m}(Ws)_{\ell+1,m}^*] + \mathbb{E}[(Wn)_{\ell,m}(Wn)_{\ell+1,m}^*] \qquad (C.5)$$

$$= \mathbb{E}[(Ws^{\text{boost}})_{\ell,m}(Ws^{\text{boost}})_{\ell+1,m}^*] + \mathbb{E}[(Ws^{\text{unboost}})_{\ell,m}(Ws^{\text{unboost}})_{\ell+1,m}^*] \\ + \mathbb{E}[(Wn)_{\ell,m}(Wn)_{\ell+1,m}^*] \qquad (C.6)$$

$$\approx F_{\ell} F_{\ell+1} \mathbb{E}[O_{\ell,m}^{\text{fullsky}}] + \mathbb{E}[O_{\ell,m}^{\text{unboost}}]_{\text{MC}}, \qquad (C.7)$$

where in the second line, we assumed that the signal and noise are uncorrelated. In the third line, we separated the galaxy clustering contributions coming from the boost (of order
$\mathcal{O}(\beta^1)$) and the unboosted classical part (of order $\mathcal{O}(\beta^0)$) that now leaks in the estimator due to isotropy breaking by the mask. Finally, in the fourth line, we (i) used the “diagonal” approximation from [45], $(Ws)_{\ell,m} \simeq F_\ell s_{\ell,m}$, i.e. we neglect the small amount of mask-induced mode-mixing, and (ii) took the unboosted contributions (from both the signal and noise) through Monte-Carlo (MC) average over (unboosted) simulations. This latter MC average is called the “mean-field” in the lensing literature. Following that tradition, we call

$$O_{\ell,m}^{MF} \equiv \mathbb{E}[O_{\ell,m}^{\text{unboost}}]_{\text{MC}}, \qquad (C.8)$$

and our new observable is the mean-field subtracted $O'_{\ell,m} \equiv O_{\ell,m} - O_{\ell,m}^{MF}$. Then we have

$$\mathbb{E}[O'_{\ell,m}] = F_{\ell} F_{\ell+1} A_{\ell+1,m} \times [F_{\text{kin}}(C_{\ell} + C_{\ell+1}) - \beta ((\ell+2)C_{\ell+1} - \ell C_{\ell})] \quad (C.9)$$

$$= {}^T \mathbf{f}_{\ell,m} \cdot \mathbf{\theta}, \qquad (C.10)$$

with the vector of parameters $\mathbf{\theta} = (\beta, F_{\text{kin}})$. Hence, we are in a position similar to the fullsky case and we need only a final element to build the estimator: the covariance of $O'$:

$$\text{Cov}(O'_{\ell,m}, O'_{\ell',m'}) = \mathbb{E}[(Wd)_{\ell,m} (Wd)_{\ell',m'}^* (Wd)_{\ell'+1,m's}]_c \quad (C.11)$$

$$\approx F_{\ell}^{2} F_{\ell+1}^{2} C_{\ell} C_{\ell+1} \delta_{\ell,\ell'}^{K} \delta_{m,m'}^{K} + \mathcal{O}(\beta), \qquad (C.12)$$

where we used the diagonal approximation again.

Then we can use the same developments as in section 4.1 and 4.3: we can define an inner product

$$\langle \mathbf{x}, \mathbf{y} \rangle = \sum_{\ell,m} \frac{(x_{\ell,m})^* y_{\ell,m}}{F_{\ell}^2 F_{\ell+1}^2 C_{\ell} C_{\ell+1}}. \qquad (C.13)$$

and the optimal quadratic estimators are given by

$$\hat{\mathbf{\theta}} = \langle \mathbf{f}, \mathbf{f} \rangle^{-1} \times \langle \mathbf{O}', \mathbf{f} \rangle. \qquad (C.14)$$

### (iii) Summary

In summary, compared to the full sky case the differences are

* the raw observable is measured on Wiener-filtered maps

* the mean-field has to be subtracted from the observable

* the average and covariance of $\mathbf{O}'$ are modified with Wiener filter factors leading to different inner products. For instance $\langle \mathbf{O}', \mathbf{f} \rangle$ writes explicitly as

$$\langle \mathbf{O}', \mathbf{f} \rangle = \sum_{\ell,m} \frac{(O'_{\ell,m})^* F_{\ell} F_{\ell+1} A_{\ell+1,m}}{F_{\ell}^2 F_{\ell+1}^2 C_{\ell} C_{\ell+1}} \left( -(\ell+2) C_{\ell+1} + \ell C_{\ell} \right) \quad (C.15)$$
# D Revising assumptions: Gaussian covariance and known power spectrum

In this section, we examine the validity of two assumptions that we made in the main text and that might seem questionable due to non-linearities on small scales: (i) using a Gaussian covariance, and (ii) that we know perfectly the angular power spectrum $C_\ell$ needed as an ingredient for the covariance and the templates $\alpha$ and $\gamma$. We do so by quantifying the impact of going beyond these assumptions, and we will show that this impact can be neglected to the degree of precision achievable by our method.

## D.1 Gaussian covariance

When we push the analyses to high $\ell_{max}$, the non-Gaussian contribution (coming from the trispectrum) to the covariance of a quantity quadratic in the galaxy density becomes inevitably more important. In fact this non-Gaussianity has a non-negligible impact on cosmological constraints using the galaxy power spectrum already for $\ell_{max} = 300$, as shown in [46]. This non-Gaussian contribution is complex, being composed of a large number of terms [47], although Super-Sample Covariance (SSC) is generally the term with the dominant impact [48].

However, in our case we are not interested in cosmological constraints using measurements of the $C_\ell$'s. Instead we are interested in obtaining a detection of $\beta$ with the $(\ell, \ell+1)$ correlator. That makes a large difference, because the impact of non-Gaussianity is not so much dictated by $\ell_{max}$ but much more by the statistical power of the measurement. Indeed, [49] showed that the impact of SSC on parameter constraints scale with the (square of the) signal to noise ratio (S/N) of the measurement. More precisely, in the low S/N regime the relative difference of the Fisher matrix element is [49]

$$ \left| \frac{\delta F}{F} \right| = \cos^2 \theta \times \frac{(S/N)^2}{(S/N)_{\text{max}}^2} \qquad (\text{D.1}) $$

where $\cos \theta$ is a quantity defined in [49] that quantifies how much the inferred parameter is degenerate with the impact of SSC. As hinted by the notation this number lies in the range $-1 \le \cos \theta \le 1$. Here, we are going to be conservative and take the upper limit $\cos^2 \theta = 1$. Then $(S/N)_{max}$ is the SSC-induced plateau of signal to noise. For the median Euclid redshift bin, ref. [49] found a value of $(S/N)_{max} \simeq 250$. We are going to be extremely conservative and use this value for our best possible scenario: the multi-tracer analysis of all 10 Euclid bins. This is extremely conservative because the impact of SSC scales down with the survey volume. With these assumptions, for the best possible scenario where we claim a detectability of $24\sigma$, the dominant non-Gaussian impact would increase our error bars by at most

$$ \left| \frac{\delta F}{F} \right| \le \left( \frac{24}{250} \right)^2 \sim 1\%. \qquad (\text{D.2}) $$

We thus conclude that the impact of non-Gaussianity can safely be neglected at the level of precision we claim in the main text.

## D.2 Known power spectrum $C_\ell$

The coefficients $\alpha_{\ell,m}$ and $\gamma_{\ell,m}$ that enter in our estimators depend on the power spectra $C_\ell$ (see eqs. (4.2)–(4.3)). As discussed in the main text, one can either use measured power
spectra or theoretically modelled power spectra. Here, we show that the optimal solution is a “hybrid model”, where one uses the theoretically modelled spectra at low $\ell$, in the perturbative regime, and the measured spectra at high $\ell$, when the measurement uncertainties are negligible, but no reliable model exists. More precisely, we use the fact that

* at the highest multipoles, $\ell_{max} \sim 10,000$ we enter the shot-noise dominated regime. In this regime, the power spectrum is known: $C_{\ell} = 1/N_{gal}$ where $N_{gal}$ is the actual number of galaxies (per steradian) present in the catalog. So this number is known exactly (see e.g. the shot-noise subtraction discussion in [47]). Hence, we have perfect knowledge of $C_{\ell}$ there,

* at low multipoles a perturbative model can be used for a faithful reproduction of the Baryonic Acoustic Oscillations (BAO). Such models are already available in the literature due to the cosmological importance of the BAO. They have sub-percent accuracy, which we will show below to be largely sufficient for our purpose,

* at intermediate multipoles which lie in the non-linear regime but are not yet dominated by shot-noise, for our purpose we do not need a cosmological model of $C_{\ell}$. Instead, the raw (binned) measurement of $C_{\ell}$ can be used. Indeed it can easily reach percent level accuracy, which we will show below to be largely sufficient for our purpose.

To justify our second and third points, we now quantify the issues raised by an incorrect $C_{\ell}$ and show each time that percent level precision is enough for our purpose. To do so, we first justify the second point by considering the impact of a fixed incorrect $C_{\ell}$, due to the (small) imprecision of the galaxy clustering model. Next we justify the third point by considering a stochastic $C_{\ell}$ measured from the galaxy clustering data itself. In both cases, we will show the calculations in the simpler case of a single tracer and single parameter estimate of $\beta$.

First is the case of a fixed incorrect $C_{\ell}$. Here the QE has two potential issues: bias and increase of error bar. Concerning the bias, one sees that an incorrect covariance matrix does not bias the QE; instead the bias comes exclusively from using an incorrect template $\gamma_{mod}$ different from the true template $\gamma_{true}$. Some algebra shows that the expectation value of the QE is

$$ \mathbb{E}[\hat{\beta}] = \frac{\langle \gamma_{\text{true}}, \gamma_{\text{mod}} \rangle}{\langle \gamma_{\text{mod}}, \gamma_{\text{mod}} \rangle} \times \beta \qquad (\text{D.3}) $$

Hence an x% precision on $C_{\ell}$ (and thus on $\gamma$) yields a bias on $\beta$ of at most x%. In the best scenario in the main text we had 4-5% precision on $\beta$, which means that we need a precision on $C_{\ell}$ of a few %. The latter is well within reach of current models, see e.g. [50].

Concerning the increase of error bar, one sees that it comes mostly from using an incorrect covariance matrix $\Sigma_{mod}$ different from the true covariance $\Sigma_{true}$. Some algebra shows that the variance of the QE is

$$ \text{Var}(\hat{\beta}) = \frac{\gamma^T \Sigma_{\text{mod}}^{-1} \Sigma_{\text{true}} \Sigma_{\text{mod}}^{-1} \gamma}{\|\gamma\|^4}. \qquad (\text{D.4}) $$
This does not simplify analytically in the general case, except if we take the relative variation of the power spectrum to be small and constant with multipole, in which case we get

$$ \frac{\delta \mathrm{Var}(\hat{\beta})}{\mathrm{Var}(\hat{\beta})} = 2 \frac{\delta C_{\ell}}{C_{\ell}}. \qquad (D.5) $$

This means that a precision on $C_\ell$ of a few % also gives a precision on the error bar on $\beta$ of a few percent, which is sufficient for our needs.

In conclusion to this first case, we see that it is a viable option for our purpose to use current non-linear galaxy modeling in order to produce a template for $C_\ell$, in the range of multipoles where the precision of the model is below a few percent.

Second, we examine the case where we use the data itself to estimate $C_\ell$. We consider this a superior alternative, in terms of both precision and practicality, in the high multipole regime outside of the BAO. This case is in principle more complex because we have to account for the fact that stochastic variations in the template $\gamma$ and the estimated covariance $\Sigma$ are correlated together as well as with fluctuations of the observable $O_{\ell,m}$ that we are fitting against. This affects the two parts of the estimator: the normalisation $\|\gamma\|^2$ and the inner product $\langle O, \gamma \rangle$. But we can argue that the relative variation of both these quantities is of the same order as the relative variation of the power spectrum and then determine conditions for the latter to be small enough.

We begin with the normalisation; with an estimated power spectrum $\hat{C}_\ell$, it becomes

$$ \|\gamma\|_{\mathrm{est}}^2 = \sum_{\ell} \frac{\ell+1}{3} \frac{((\ell+2)\hat{C}_{\ell+1} - \ell\hat{C}_{\ell})^2}{\hat{C}_{\ell}\hat{C}_{\ell+1}} \qquad (D.6) $$

$$ \approx \sum_{\ell} \frac{\ell+1}{3} \frac{((\ell+2)C_{\ell+1} - \ell C_{\ell})^2}{C_{\ell}C_{\ell+1}} \times \left[ 1 - \frac{\delta C_{\ell}}{C_{\ell}} - \frac{\delta C_{\ell+1}}{C_{\ell+1}} + 2 \frac{\delta ((\ell+2)C_{\ell+1} - \ell C_{\ell})}{(\ell+2)C_{\ell+1} - \ell C_{\ell}} \right], \quad (D.7) $$

where the variations $\delta X$ are now stochastic variables. Noting that for a local power law $C_\ell \propto \ell^{-n}$ we have

$$ (\ell + 2) C_{\ell+1} - \ell C_{\ell} \approx \frac{1}{\ell} \frac{d\ell^2 C_{\ell}}{d\ell} = (2-n)C_{\ell}, \qquad (D.8) $$

and thus that its relative variation is the same as $\frac{\delta C_\ell}{C_\ell}$, we can put the upper limit on the standard deviation

$$ \frac{\sigma(\|\gamma\|^2)}{\|\gamma\|^2} \le 4 \left| \frac{\sigma(C_\ell)}{C_\ell} \right|. \qquad (D.9) $$

Next we have the inner product; with an estimated power spectrum $\hat{C}_\ell$ it becomes

$$ \langle O, \gamma \rangle_{\mathrm{est}} = \sum_{\ell,m} -O_{\ell+1,m} A_{\ell+1,m} \frac{(\ell+2) \hat{C}_{\ell+1} - \ell \hat{C}_{\ell}}{\hat{C}_{\ell} \hat{C}_{\ell+1}} f_{\mathrm{SKY}} \qquad (D.10) $$

$$ \approx \sum_{\ell,m} -O_{\ell+1,m} A_{\ell+1,m} \frac{(\ell+2)C_{\ell+1} - \ell C_{\ell}}{C_{\ell}C_{\ell+1}} f_{\mathrm{SKY}} \times \left[ 1 - \frac{\delta C_{\ell}}{C_{\ell}} - \frac{\delta C_{\ell+1}}{C_{\ell+1}} + \frac{\delta ((\ell+2)C_{\ell+1} - \ell C_{\ell})}{(\ell+2)C_{\ell+1} - \ell C_{\ell}} \right], \qquad (D.11) $$
with the same reasoning as for $||\gamma||^2$, we get the upper limit

$$ \frac{\sigma(\langle \mathbf{O}, \gamma \rangle)}{\langle \mathbf{O}, \gamma \rangle} \le 3 \left| \frac{\sigma(C_\ell)}{C_\ell} \right| \qquad (D.12) $$

Thus, for both parts of the QE, the increase of error bar is under control for our precision requirements if the standard deviation of the estimated power spectrum is less than a percent: $\sigma(C_\ell)/C_\ell \le 1\%$.

The calculation above stands for the increase of error bar. Remains the potential issue of the bias. Calculating the bias requires to develop the normalisation and inner product to second order, since $E(\delta C_\ell) = 0$ but $E((\delta C_\ell)^2) \ne 0$. This is cumbersome and not particularly enlightening; instead it suffices to remark that the resulting relative biases will necessarily be proportional to $(\sigma(C_\ell)/C_\ell)^2$. Thus, if we apply the condition $\sigma(C_\ell)/C_\ell \le 1\%$ to keep control of the error bars, we will a fortiori keep control of the bias.

To close the argument we thus need to find the requirements for $\sigma(C_\ell)/C_\ell \le 1\%$ and show that they are realistic. For the input $C_\ell$ we simply use the actual binned power spectrum measurement.<fnref n="4" /> There, the standard $f_{\text{SKY}}$ formula for the cosmic variance gives

$$ \left( \frac{\sigma(C_\ell)}{C_\ell} \right)^2 = \frac{2}{(2\ell + 1) \Delta\ell f_{\text{SKY}}} \qquad (D.13) $$

For coming galaxy surveys like Euclid we have $f_{\text{SKY}} \sim 1/3$. Thus for a bin width $\Delta\ell = 100$, our requirement is fulfilled at a multipole as low as $\ell = 300$. Since the power spectrum is expected to be modeled accurately down to at least $\ell = 2000$ [32], we see that we can easily patch the two regimes: cosmological model at low multipoles and actual measurement at high multipoles.

We conclude that the main text assumption of a known power spectrum is largely valid for our precision requirements.

## D.3 Non-linear growth and bias

We start by recalling the equation for the boosted counts, eq. (2.19), isolating the $F_{\text{kin}}$ term:

$$ N(z, \mathbf{n}, S_*) = [1 + \cos \theta D_{\text{kin}}(z)] \tilde{N}(z, S_*) + [1 + \beta \sin \theta \partial_{\theta}] \tilde{N}(z, S_*) \delta_{\text{g}}(z, \mathbf{n}, S_*) \\ + \cos \theta F_{\text{kin}}(\mathbf{n}, z) \tilde{N}(z, S_*) \delta_{\text{g}}(z, \mathbf{n}, S_*), \qquad (D.14) $$

where

$$ F_{\text{kin}}(\mathbf{n}, z) = D_{\text{kin}}(z) + \left[ \frac{\partial \ln \delta_{\text{g}}}{\partial \ln(1+z)}(\mathbf{n}, z) + \frac{2}{rH} \frac{\partial \ln \delta_{\text{g}}}{\partial \ln S_{*}}(\mathbf{n}, z) \right] \beta \qquad (D.15) $$

In the main text, the log derivative reduced to the growth rate or bias derivatives which were assumed scale independent. This assumption cannot be sustained in practice in the non-linear regime. To study this, we rewrite the count equation as

$$ N(z, \mathbf{n}, S_*) = [1 + \cos \theta D_{\text{kin}}(z)] \tilde{N}(z, S_*) + [1 + \beta \sin \theta \partial_{\theta}] \tilde{N}(z, S_*) \delta_{\text{g}}(z, \mathbf{n}, S_*) \\ + \cos \theta \tilde{N}(z, S_*) \mathcal{D}_k \delta_{\text{g}}(z, \mathbf{n}, S_*), \qquad (D.16) $$

<footnote n="4">⁴We note that this could easily be refined. Indeed, in this regime the spectra are expected to be very smooth, close to a power law with a slowly varying spectral index. This could be used to refine a more precise estimate, e.g. with spline interpolation in the space $(\ln \ell, \ln C_\ell)$.</footnote>
where we defined the differential operator

$$ \mathcal{D}_k = D_{\text{kin}} + \beta \left[ \frac{\partial}{\partial \ln(1+z)} + \frac{2}{rH} \frac{\partial}{\partial \ln S_*} \right] = D_{\text{kin}} + \beta \mathcal{D}. \quad (D.17) $$

This can be propagated into the derivation of section 2, where it only changes eq. (2.33) into

$$ \delta a_{\ell,m}^{\text{F}} = \tilde{N}(z) \int_{\mathbb{S}_2} d^2 \boldsymbol{n} \mathcal{D}_k \delta_g(z, \boldsymbol{n}) [A_{\ell,m} Y_{\ell-1,m}^*(\boldsymbol{n}) + A_{\ell+1,m} Y_{\ell+1,m}^*(\boldsymbol{n})] \quad (D.18) $$

$$ = A_{\ell,m} a_{\ell-1,m} (\tilde{N} \mathcal{D}_k \delta_g) + A_{\ell+1,m} a_{\ell+1,m} (\tilde{N} \mathcal{D}_k \delta_g). \quad (D.19) $$

In section 3, the effect of non-linear growth and bias impacts only eq. (3.6), that it modifies into

$$ \mathbb{E}[O_{\ell,m}] = A_{\ell+1,m} \cdot \left[ \left( C_{\ell}^{\delta_g, \mathcal{D}_k \delta_g} + C_{\ell+1}^{\delta_g, \mathcal{D}_k \delta_g} \right) - \beta \left( (\ell+2) C_{\ell+1}^{\delta_g, \delta_g} - \ell C_{\ell}^{\delta_g, \delta_g} \right) \right], \quad (D.20) $$

while previously only $C_\ell = C_\ell^{\delta_g, \delta_g}$ appeared.

If we define $F_{\text{kin}}(\ell)$ as the ratio of the two different power spectra

$$ F_{\text{kin}}(\ell) = \frac{C_{\ell}^{\delta_g, \mathcal{D}_k \delta_g}}{C_{\ell}^{\delta_g, \delta_g}} = D_{\text{kin}} + \beta \frac{C_{\ell}^{\delta_g, \mathcal{D} \delta_g}}{C_{\ell}^{\delta_g, \delta_g}}, \quad (D.21) $$

we see that we get and equation similar to eq. (3.6)

$$ \mathbb{E}[O_{\ell,m}] = A_{\ell+1,m} \cdot \left[ (F_{\text{kin}}(\ell) C_{\ell} + F_{\text{kin}}(\ell+1) C_{\ell+1}) - \beta \left( (\ell+2) C_{\ell+1} - \ell C_{\ell} \right) \right]. \quad (D.22) $$

In practice, the ratio $C_\ell^{\delta_g, \mathcal{D} \delta_g} / C_\ell^{\delta_g, \delta_g}$ should be modeled using a non-linear description of the large scale structure, if a sufficiently precise one is available, or using simulations.<fnref n="5" /> Part of the ratio can also be determined from the survey itself, namely the derivative $\frac{\partial}{\partial \ln S_*}$ which may be the hardest to model or simulate accurately, by taking several flux cuts and estimating the derivative with finite difference. From this, we would have a parametric model for $F_{\text{kin}}(\ell)$, whose parameters are constrained by priors from the model/simulations/survey.

In order to study the impact of introducing such a parametric model for $F_{\text{kin}}(\ell)$, we choose here to test as toy model a simple agnostic form:

$$ F_{\text{kin}}(\ell) = F_{\text{kin}}^{[0]} + F_{\text{kin}}^{[1]} \times \frac{\ell}{\ell_0}, \quad (D.23) $$

with $\ell_0$ an arbitrary pivot scale. In our implementation, we chose $\ell_0 = \ell_{\text{max}} = 10,000$ so that both parameters have the same order of magnitude: $|F_{\text{kin}}^{[1]}| = |F_{\text{kin}}^{[0]}|$ corresponds to a 100% variation of $F_{\text{kin}}$ over the multipole range. For our test to be conservative, we decide to assume no prior on any of the two parameters, i.e. both $F_{\text{kin}}^{[0]}$ and $F_{\text{kin}}^{[1]}$ are free.

Some straightforward generalisation of the derivation in appendix B shows that we can build joint Quadratic Estimators for $\beta$ and the $F_{\text{kin}}$ parameters. For instance, in the

<footnote n="5">⁵A small number of simulations should be sufficient because the ratio will be prone to cosmic variance cancellation. Indeed upwards (downwards) fluctuations of matter density will increase (decrease) both the numerator and the denominator of the ratio.</footnote>
<table><caption>**Table 1.** First column: all cases studied in the main text, following the nomenclature in the legend of figure 6. Second column: detection significance for $\beta$ at $\ell_{max} = 10,000$ when $F_{kin}$ is constant with multipole, i.e. the results quoted in the main text. Third column: detection significance for $\beta$ at $\ell_{max} = 10,000$ when $F_{kin}$ is allowed to vary with $\ell$ through eq. (D.23), and we marginalise over that scale dependence. Fourth column: relative difference between the second and third columns.</caption><thead><tr><th>case</th><th>constant F<sub>kin</sub></th><th>scale-dependent F<sub>kin</sub></th><th>relative difference</th></tr></thead><tbody><tr><td>Original / free F<sub>kin</sub></td><td>2.5 σ</td><td>2.0 σ</td><td>19%</td></tr><tr><td>30% F<sub>kin</sub> prior</td><td>3.1 σ</td><td>2.6 σ</td><td>18%</td></tr><tr><td>10% F<sub>kin</sub> prior</td><td>6.0 σ</td><td>5.0 σ</td><td>16%</td></tr><tr><td>30% A<sub>i</sub> prior</td><td>8.9 σ</td><td>8.7 σ</td><td>2%</td></tr><tr><td>10% A<sub>i</sub> prior</td><td>24 σ</td><td>21 σ</td><td>11%</td></tr></tbody></table>

single tracer case, we have

$$ \begin{pmatrix} \|\alpha\|^2 & \langle\alpha, \gamma\rangle \langle\alpha, \delta\rangle \\ \langle\gamma, \alpha\rangle & \|\gamma\|^2 \langle\gamma, \delta\rangle \\ \langle\delta, \alpha\rangle & \langle\delta, \gamma\rangle \|\delta\|^2 \end{pmatrix} \times \begin{pmatrix} \hat{F}_{\text{kin}}^{[0]} \\ \hat{\beta} \\ \hat{F}_{\text{kin}}^{[1]} \end{pmatrix} = \begin{pmatrix} \langle\mathbf{O}, \alpha\rangle \\ \langle\mathbf{O}, \gamma\rangle \\ \langle\mathbf{O}, \delta\rangle \end{pmatrix}, \qquad (\text{D.24}) $$

where the weights are

$$ \alpha_{\ell,m} = A_{\ell+1,m} (C_{\ell} + C_{\ell+1}) \qquad \gamma_{\ell,m} = -A_{\ell+1,m} ((\ell+2)C_{\ell+1} - \ell C_{\ell}) \\ \delta_{\ell,m} = A_{\ell+1,m} (\ell C_{\ell} + (\ell+1)C_{\ell+1}). \qquad (\text{D.25}) $$

The multi-tracer case is obtained through a straightforward generalization.

We implemented these estimators and reran the analysis of section 5 with the new free parameter $F_{\text{kin}}^{[1]}$. We found that it yields a modest increase of the error bars on $\beta$, with the impact on the detection significance at $\ell_{max} = 10,000$ summarised in table 1.

We see that the impact of this conservative non-linear toy model is small, ranging between 2% and 20%. We found the impact to generally be even smaller when the multipole range is cut to a smaller $\ell_{max}$. We conclude that the impact of non-linear growth and bias will not alter significantly the conclusions in the main text, and in particular the prospect of a significant detection of $\beta$ with the proposed method.

## References

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

[2] A. Kogut et al., *Dipole anisotropy in the COBE DMR first year sky maps*, Astrophys. J. 419 (1993) 1 [astro-ph/9312056] [INSPIRE].

[3] PLANCK collaboration, Planck 2018 results. I. Overview and the cosmological legacy of Planck, Astron. Astrophys. 641 (2020) A1 [arXiv:1807.06205] [INSPIRE].

[4] C.A.P. Bengaly, R. Maartens and M.G. Santos, Probing the cosmological principle in the counts of radio galaxies at different frequencies, JCAP 04 (2018) 031 [arXiv:1710.08804] [INSPIRE].
[5] P. Tiwari and A. Nusser, *Revisiting the NVSS number count dipole*, JCAP 03 (2016) 062
[arXiv:1509.02532] [INSPIRE].

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

[7] J.D. Wagenveld, H.-R. Klöckner and D.J. Schwarz, *The cosmic radio dipole: bayesian estimators on new and old radio surveys*, Astron. Astrophys. 675 (2023) A72 [arXiv:2305.15335] [INSPIRE].

[8] P. Tiwari and P.K. Aluri, *Large angular-scale multipoles at redshift ~ 0.8*, Astrophys. J. 878 (2019) 32 [arXiv:1812.04739] [INSPIRE].

[9] J. Colin, R. Mohayaee, M. Rameez and S. Sarkar, *High redshift radio galaxies and divergence from the CMB dipole*, Mon. Not. Roy. Astron. Soc. 471 (2017) 1045 [arXiv:1703.09376] [INSPIRE].

[10] N.J. Secrest et al., *A test of the cosmological principle with quasars*, Astrophys. J. Lett. 908 (2021) L51 [arXiv:2009.14826] [INSPIRE].

[11] N.J. Secrest et al., *A challenge to the standard cosmological model*, Astrophys. J. Lett. 937 (2022) L31 [arXiv:2206.05624] [INSPIRE].

[12] L. Dam, G.F. Lewis and B.J. Brewer, *Testing the cosmological principle with CatWISE quasars: a bayesian analysis of the number-count dipole*, Mon. Not. Roy. Astron. Soc. 525 (2023) 231 [arXiv:2212.07733] [INSPIRE].

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

[14] Y.-T. Cheng, T.-C. Chang and A. Lidz, *Is the radio source dipole from NVSS consistent with the cosmic microwave background and ΛCDM?*, Astrophys. J. 965 (2024) 32 [arXiv:2309.02490] [INSPIRE].

[15] A.K. Singal, *Large peculiar motion of the solar system from the dipole anisotropy in sky brightness due to distant radio sources*, Astrophys. J. Lett. 742 (2011) L23 [arXiv:1110.6260] [INSPIRE].

[16] J. Colin, R. Mohayaee, M. Rameez and S. Sarkar, *Evidence for anisotropy of cosmic acceleration*, Astron. Astrophys. 631 (2019) L13 [arXiv:1808.04597] [INSPIRE].

[17] C. Dalang and C. Bonvin, *On the kinematic cosmic dipole tension*, Mon. Not. Roy. Astron. Soc. 512 (2022) 3895 [arXiv:2111.03616] [INSPIRE].

[18] C. Guandalin, J. Piat, C. Clarkson and R. Maartens, *Theoretical systematics in testing the cosmological principle with the kinematic quasar dipole*, Astrophys. J. 953 (2023) 144 [arXiv:2212.04925] [INSPIRE].

[19] S. von Hausegger, *The expected kinematic matter dipole is robust against source evolution*, arXiv:2404.07929 [INSPIRE].

[20] EUCLID collaboration, *Euclid preparation. The near-IR background dipole experiment with Euclid*, arXiv:2401.17945 [INSPIRE].

[21] PLANCK collaboration, *Planck 2013 results. XXVII. Doppler boosting of the CMB: eppur si muove*, Astron. Astrophys. 571 (2014) A27 [arXiv:1303.5087] [INSPIRE].

[22] S. Saha et al., *Bayesian estimation of our local motion from the Planck-2018 CMB temperature map*, JCAP 10 (2021) 072 [arXiv:2106.07666] [INSPIRE].

[23] P.S. Ferreira and M. Quartin, *First constraints on the intrinsic CMB dipole and our velocity with Doppler and aberration*, Phys. Rev. Lett. 127 (2021) 101301 [arXiv:2011.08385] [INSPIRE].
[24] C. Dalang, R. Durrer and F. Lacasa, Statistical effects of the observer's peculiar velocity on source number counts, Phys. Rev. D 107 (2023) 103514 [arXiv:2209.12812] [INSPIRE].

[25] R. Maartens, C. Clarkson and S. Chen, The kinematic dipole in galaxy redshift surveys, JCAP 01 (2018) 013 [arXiv:1709.04165] [INSPIRE].

[26] T. Nadolny, R. Durrer, M. Kunz and H. Padmanabhan, A new way to test the cosmological principle: measuring our peculiar velocity and the large-scale anisotropy independently, JCAP 11 (2021) 009 [arXiv:2106.05284] [INSPIRE].

[27] C. Bonvin and R. Durrer, What galaxy surveys really measure, Phys. Rev. D 84 (2011) 063505 [arXiv:1105.5280] [INSPIRE].

[28] EUCLID collaboration, Euclid preparation. XIX. Impact of magnification on photometric galaxy clustering, Astron. Astrophys. 662 (2022) A93 [arXiv:2110.05435] [INSPIRE].

[29] C. Bonvin, R. Durrer and M.A. Gasparini, Fluctuations of the luminosity distance, Phys. Rev. D 73 (2006) 023523 [Erratum ibid. 85 (2012) 029901] [astro-ph/0511183] [INSPIRE].

[30] F. Sorrenti, R. Durrer and M. Kunz, The dipole of the Pantheon+SH0ES data, JCAP 11 (2023) 054 [arXiv:2212.10328] [INSPIRE].

[31] W. Hu and T. Okamoto, Mass reconstruction with CMB polarization, Astrophys. J. 574 (2002) 566 [astro-ph/0111606] [INSPIRE].

[32] EUCLID collaboration, Euclid preparation. VII. Forecast validation for Euclid cosmological probes, Astron. Astrophys. 642 (2020) A191 [arXiv:1910.09273] [INSPIRE].

[33] PLANCK collaboration, Planck 2015 results. XIV. Dark energy and modified gravity, Astron. Astrophys. 594 (2016) A14 [arXiv:1502.01590] [INSPIRE].

[34] R. Takahashi et al., Revising the halofit model for the nonlinear matter power spectrum, Astrophys. J. 761 (2012) 152 [arXiv:1208.2701] [INSPIRE].

[35] S. Bird, M. Viel and M.G. Haehnelt, Massive neutrinos and the non-linear matter power spectrum, Mon. Not. Roy. Astron. Soc. 420 (2012) 2551 [arXiv:1109.4416] [INSPIRE].

[36] EUCLID collaboration, Euclid definition study report, arXiv:1110.3193 [INSPIRE].

[37] M.S. Wang, F. Beutler and D. Bacon, Impact of relativistic effects on the primordial non-Gaussianity signature in the large-scale clustering of quasars, Mon. Not. Roy. Astron. Soc. 499 (2020) 2598 [arXiv:2007.01802] [INSPIRE].

[38] R. Maartens et al., Magnification and evolution biases in large-scale structure surveys, JCAP 12 (2021) 009 [arXiv:2107.13401] [INSPIRE].

[39] ASTROPY collaboration, Astropy: a community python package for astronomy, Astron. Astrophys. 558 (2013) A33 [arXiv:1307.6212] [INSPIRE].

[40] ASTROPY collaboration, The Astropy project: building an open-science project and status of the v2.0 core package, Astron. J. 156 (2018) 123 [arXiv:1801.02634] [INSPIRE].

[41] E. Jones et al., SciPy: open source scientific tools for python, http://www.scipy.org/ (2001).

[42] T. Oliphant, NumPy: a guide to NumPy, Trelgol publishing, U.S.A. (2006).

[43] J.D. Hunter, Matplotlib: a 2D graphics environment, Comput. Sci. Eng. 9 (2007) 90 [INSPIRE].

[44] M. Abramowitz and I.A. Stegun, Handbook of mathematical functions, Dover (1972).

[45] PLANCK collaboration, Planck 2015 results. XV. Gravitational lensing, Astron. Astrophys. 594 (2016) A15 [arXiv:1502.01591] [INSPIRE].
[46] F. Lacasa and R. Rosenfeld, *Combining cluster number counts and galaxy clustering*, JCAP 08
(2016) 005 [arXiv:1603.00918] [INSPIRE].

[47] F. Lacasa, *Covariance of the galaxy angular power spectrum with the halo model*, Astron.
Astrophys. 615 (2018) A1 [arXiv:1711.07372] [INSPIRE].

[48] F. Lacasa, *The impact of braiding covariance and in-survey covariance on next-generation galaxy
surveys*, Astron. Astrophys. 634 (2020) A74 [arXiv:1909.00791] [INSPIRE].

[49] F. Lacasa and J. Grain, *Fast and easy super-sample covariance of large scale structure
observables*, Astron. Astrophys. 624 (2019) A61 [arXiv:1809.05437] [INSPIRE].

[50] J. Donald-McCann, F. Beutler, K. Koyama and M. Karamanis, *matryoshka: halo model
emulator for the galaxy power spectrum*, Mon. Not. Roy. Astron. Soc. 511 (2022) 3768
[arXiv:2109.15236] [INSPIRE].