# Combining chirp mass, luminosity distance, and sky localization from gravitational wave events to detect the cosmic dipole

N. Grimm<sup>1,★</sup> M. Pijnenburg,<sup>1★</sup> S. Mastrogiovanni<sup>1,2</sup> C. Bonvin<sup>1,3</sup> S. Foffa<sup>1</sup> and G. Cusin<sup>1,3</sup>

Accepted 2023 October 2. Received 2023 October 2; in original form 2023 August 7

## ABSTRACT

A key test of the isotropy of the Universe on large scales consists in comparing the dipole in the cosmic microwave background (CMB) temperature with the dipole in the distribution of sources at low redshift. Current analyses find a dipole in the number counts of quasars and radio sources that is 2–5 times larger than expected from the CMB, leading to a tension reaching 5σ. In this paper, we derive a consistent framework to measure the dipole independently from gravitational wave (GW) detections. We exploit the fact that the observer velocity does not only change the distribution of events in the sky, but also the luminosity distance and redshifted chirp mass, which can be extracted from the GW waveform. We show that the estimator with higher signal-to-noise ratio is the dipole in the chirp mass measured from a population of binary neutron stars. Combining all estimators (accounting for their covariance) improves the detectability of the dipole by 30–50 per cent compared to number counting of binary black holes alone. We find that a few 10<sup>6</sup> events are necessary to detect a dipole consistent with the CMB one, whereas if the dipole is as large as predicted by radio sources, it will already be detectable with 10<sup>5</sup> events, which would correspond to a single year of observation with next-generation GW detectors. GW sources provide therefore a robust and independent way of testing the isotropy of the Universe.

**Key words:** gravitational waves – galaxies: active – cosmic background radiation.

## 1 INTRODUCTION

One of the basic assumptions of the ΛCDM cosmological model is that our Universe is homogeneous and isotropic on large scales. The latter follows from the high large-scale isotropy in observations of the cosmic microwave background (CMB) temperature and of the large-scale structure of the Universe. Combining this with the cosmological principle further leads to the homogeneity of the Universe. In the ΛCDM model, the observed dipole anisotropy in the CMB temperature (Aghanim et al. 2020) is due to the fact that we, as observer, are moving with respect to the homogeneous and isotropic background.<fnref n="1" />

If this picture is correct, then the same motion should induce a dipole in the distribution of sources, due to aberration and magnification effects (Ellis & Baldwin 1984). This idea has been put to the test in the past years, through measurement of the dipole in the number counts of quasars and radio sources (Colin et al. 2017; Bengaly, Maartens & Santos 2018; Secrest et al. 2021, 2022; Siewert, Schmidt-Rubart & Schwarz 2021) at redshifts 0 ≤ z ≲ 3. The direction of these dipoles is well aligned with that of the CMB, however, the amplitude is 2–5 times larger than expected, leading

to a tension with the CMB dipole reaching up to 5.1σ (Secrest et al. 2022). Supernova type Ia light curves [see e.g. Riess, Press & Kirshner (1995)] provide an alternative means of assessing the dipole anisotropy (Bonvin, Durrer & Kunz 2006b). Indeed, measurements from supernova catalogues performed in Horstmann, Pietschke & Schwarz (2022) and Singal (2022) have found a dipole again aligned with the CMB. However, they lead to an amplitude either compatible with the measurement from radio and quasar sources (Singal 2022) or even lower than the CMB value (Horstmann et al. 2022). Another measurement in Sorrenti, Durrer & Kunz (2022) shows an amplitude consistent with the CMB, but a strong tension in the direction. Hence, the inconclusive measurements of the dipole from type Ia supernovae do not resolve any tensions.

The discrepancy in the dipole amplitude could be due to systematic effects in the quasar and radio source data sets, to an imperfect theoretical modelling of the expected dipole in the number counts, which currently neglects evolution effects (Dalang & Bonvin 2022; Guandalin et al. 2023), or to a violation of isotropy in our Universe. In this last case, the large dipole in the quasars and radio sources would not be solely due to the observer velocity, but it would have an intrinsic part generated by a large local anisotropy in the Universe (inconsistent with the ΛCDM predictions).

One way to test these scenarios is to use other data sets to measure the dipole at low redshift, and see if the result is consistent with the CMB dipole, or with the dipoles from quasars and radio sources. A promising avenue is to use gravitational waves (GWs) from binary systems of black holes (BBH) or neutron stars (BNS). In Cai et al.

<footnote n="1"><sup>1</sup>Département de Physique Théorique and Gravitational Wave Science Center, Université de Genève, 24 quai Ernest-Ansermet, CH-1211 Genève 4, Switzerland</footnote>

<footnote n="2"><sup>2</sup>INFN, Sezione di Roma, I-00185 Roma, Italy</footnote>

<footnote n="3"><sup>3</sup>Institut d'Astrophysique de Paris, Sorbonne Université, CNRS, UMR 7095, F-75014 Paris, France</footnote>

<footnote>* E-mail: nastassia.grimm@unige.ch (NG); artin.pijnenburg@unige.ch (MP)</footnote>

<footnote n="1">¹Note that in the ΛCDM model, fluctuations in the matter density with respect to the homogeneous and isotropic background also generate a dipole anisotropy, but this contribution is negligible with respect to the one due to the kinematic dipole.</footnote>
(2018), it was proposed to look for a dipole in the luminosity distance measured from GWs. This paper does not model the kinematic contribution to the luminosity distance dipole, but it forecasts how large a dipole should be to be detected (independently of its origin). In Mastrogiovanni et al. (2023), some of us proposed to use the distribution of BBH to measure the cosmic dipole. We derived a modelling of the signal, taking into account aberration and threshold effects, and showed that with the next generation of interferometers (XG), like the Einstein Telescope (ET) and the Cosmic Explorer (CE) it will be possible to confidently detect the dipole. In Kashyap et al. (2023), the authors propose to use both the number counts of GWs and the chirp mass to measure the dipole. They apply their method to current data and show that no dipole anisotropy is detected. This is in agreement with Payne et al. (2020); Stiskalek, Veitch & Messenger (2021); Essick et al. (2023), who found no evidence for an anisotropy in the distribution of current GW sources.

In this paper, we derive a consistent framework to use GW detections to optimally measure the cosmic dipole. In particular, we exploit three of the quantities that can be extracted from the GW waveform and amplitude: the angular position of the binary system, its luminosity distance, and its redshifted chirp mass. These three quantities are all affected by the observer velocity and can therefore be combined to measure the velocity in an optimal way. We build on the results of Mastrogiovanni et al. (2023) to derive a modelling of the mean luminosity distance and the mean chirp mass per angular pixel. This modelling accounts not only for aberration, but also for threshold effects. We show that the distance and chirp mass dipoles are correlated with each other, but both are independent of the number count dipole. Combining them therefore requires taking into account this correlation, in order not to overestimate the constraints. We also show that threshold effects can be robustly modelled and computed, which is essential in order not to bias the measurement of the observer velocity.

We apply our framework to synthetic catalogues of binary systems (BBH and BNS) that would be observed by the next generation of interferometers, like ET and CE, and we forecast how well the dipole can be detected in these catalogues. We explore two scenarios: one where the dipole would be purely kinematic and therefore consistent with the CMB value, and a second one where the dipole is consistent with the results from radio sources. For this second scenario, we take the extreme case of an observer velocity that would be 5 times larger than the CMB one (Bengaly et al. 2018), and we call this the ‘AGN case’. Note that if we take the AGN case at face value, the dipole cannot be only due to the observer velocity (that would be inconsistent with the one extracted from the CMB), but it would have a large intrinsic part due to a strong anisotropy in the large-scale structure, as discussed above. However, when assessing the detectability of such a dipole, it does not matter if it is purely of kinematic origin, or if there is also an intrinsic contribution. Therefore we can simulate the AGN dipole as if it were due to an observer velocity which is 5 times larger than expected.

We find that the chirp mass from BNSs provides the estimator with lower variance, i.e. larger signal-to-noise ratio (SNR), followed by the number counts of BBHs and BNSs. This is not surprising given the pre-dominant impact of this parameter on the BNS inspiral waveform; the possibility of exploiting this fact has already been subject of various investigations (Chernoff & Finn 1993; Taylor & Gair 2012; Taylor, Gair & Mandel 2012; Finke et al. 2022). Our analysis further shows that, combining all estimators, we could detect a dipole consistent with the CMB one at > 1 σ significance with more than 10⁶ events, achieved with 5 yr of observations of ET and CE. If the dipole is as large as predicted in the AGN case, then we can

detect it at > 3 σ significance with 10⁵ events already, achieved in a year of observation.

The rest of the paper is structured as follow: in Section 2, we derive a theoretical modelling of the luminosity distance and chirp mass dipole. In Section 3, we define estimators of these dipoles and we compute their variance and covariance. In Section 4, we measure the dipoles from our synthetic catalogues of events and assess the detectability of the six estimators (three for BBHs and three for BNSs) and their combination. In Section 5, we forecast how well the observer velocity can be measured with XG detectors and we show that an imperfect knowledge of threshold effects does not degrade the constraints significantly. We conclude in Section 6.

## 2 THEORETICAL MODELLING OF THE LUMINOSITY DISTANCE AND CHIRP MASS DIPOLE

To compute the impact of the observer velocity on the luminosity distance and the chirp mass measured from GWs, we follow the same steps as in Mastrogiovanni et al. (2023), which computed the dipole in the number of events. Since the observer velocity is significantly smaller than the speed of light, we keep only terms at linear order in $v_0/c$ in the derivation.

### 2.1 Luminosity distance

The luminosity distance of a source situated at conformal distance $r$ from the observer is defined as the ratio between the intrinsic luminosity of the source, $L$, and the observed flux, $F$

$$d_L(r, \mathbf{n}) = \sqrt{\frac{L}{4\pi F(r, \mathbf{n})}}, \qquad (1)$$

where $\mathbf{n}$ denotes the direction in which the observer sees the source. The observer velocity impacts the flux measured by the observer and consequently modifies the measured luminosity distance. At first order in $v_0/c$, the luminosity distance is given by [see Bonvin, Durrer & Gasparini (2006a)<fnref n="2" />]

$$d_L(r, \mathbf{n}) = \bar{d}_L(r) \left( 1 - \mathbf{n} \cdot \frac{\mathbf{v}_0}{c} \right), \qquad (2)$$

where $\bar{d}_L(r)$ denotes the background luminosity distance in a homogeneous and isotropic universe. Note that other perturbations, in particular the peculiar velocity of the source, affect the luminosity distance besides the observer velocity, see Bonvin et al. (2006a) for the full expression. However, here we are interested in the dipole which is strongly dominated by the observer velocity, and other contributions are negligible compared to typical uncertainties of about 20 per cent for luminosity distance measurements from GW events (Callister et al. 2020; Iacovelli et al. 2022).

The luminosity distance enters directly in the amplitude of the GWs, which decays as $1/d_L$. Combining measurements from different interferometers allows one to measure $d_L$ for each binary system. We

<footnote n="2">²Note that here we are interested in the luminosity distance at fixed conformal distance $r$, which is given by equation (55) of Bonvin et al. (2006a). Equation (59) gives instead the luminosity distance at fixed observed redshift. Note also a sign difference due to the fact that in Bonvin et al. (2006a) $\mathbf{n}$ denotes the direction of propagation of the signal, whereas here $\mathbf{n}$ is the direction of observation, which points in the opposite direction.</footnote>
can then compute the mean luminosity distance from a population of binary systems in direction **n**.

$$d_{L}^{\text{mean}}(\mathbf{n}) \equiv \frac{\int dr \, dL(r, \mathbf{n}) \frac{dN_{\text{det}}}{d\Omega dr}(r, \mathbf{n})}{\int dr \, \frac{dN_{\text{det}}}{d\Omega dr}(r, \mathbf{n})} . \qquad (3)$$

Here,

$$\frac{dN_{\text{det}}}{d\Omega dr}(r, \mathbf{n}) = \int dm_{1,2} \frac{dN_{\text{det}}}{d\Omega dr dm_{1,2}}(r, \mathbf{n}, m_{1,2}, \rho > \rho_{*}) \qquad (4)$$

denotes the number of events detected at distance *r* with an SNR *ρ* above a threshold value *ρ*<sub>*</sub>, summed over all intrinsic binary masses *m*<sub>1</sub> and *m*<sub>2</sub>.<sup>3</sup> Following Mastrogiovanni et al. (2023), we adopt a simplified, zero post-Newtonian (0PN) order for the SNR, which can be computed as (Finn & Chernoff 1993)

$$\rho^2(r, \mathbf{n}, m, \mathcal{M}) = \frac{5}{96\pi^{4/3}} \frac{\Theta^2}{d_L^2(r, \mathbf{n})} (GM)^{5/3} \mathcal{F}(f_{\text{ISCO}}^z(m)), \quad (5)$$

where *M* is the redshifted chirp mass as measured in the detector frame (in the following we drop the 'redshifted' for simplicity). In equation (5), *m* is the total detector frame mass of the system and

$$f_{\text{ISCO}} \equiv \frac{1}{6\sqrt{6(2\pi)}} \frac{c^3}{Gm} \approx 2.2\,\text{kHz} \left(\frac{M_{\odot}}{m}\right) \qquad (6)$$

is the GW frequency corresponding to the innermost stable circular orbit. As to $\Theta^2$, it is a geometrical factor that accounts for the binary inclination angle and average detector's antenna patterns:

$$\Theta^2 = 0.2 \left[ \left( \frac{1 + \cos^2 \iota}{2} \right)^2 + \cos^2 \iota \right], \qquad (7)$$

where $\iota$ is the angle formed by the normal vector of the orbital plane and the line of sight. The numerical pre-factor in the above equation represents an average value of all the detector's antenna patterns. The function $\mathcal{F}$ is calculated from the integral of the power spectral density $S_n(f)$ following Maggiore (2007).

Since the SNR (5) depends on the luminosity distance and on the chirp mass, that are both affected by the observer velocity, the number of detected events above threshold $\rho_*$ will be modified by the observer velocity. These threshold effects add to the effect of aberration, which modifies the number of events per solid angle. As shown in Mastrogiovanni et al. (2023), the final result at linear order in $v_0/c$ is given by

$$\frac{dN_{\text{det}}}{d\Omega dr dm_{1,2}}(r, \mathbf{n}, m_{1,2}, \rho > \rho_{*}) = \frac{d\bar{N}_{\text{det}}}{d\Omega dr dm_{1,2}}(r, m_{1,2}, \rho_{*}) \left\{ 1 + \left[ 2 + s \left( \frac{1}{3} + \mathcal{A} \right) \right] \mathbf{n} \cdot \frac{\mathbf{v}_0}{c} \right\}, (8)$$

with<fnref n="4" />

$$s(r, \rho_{*}, m_{1,2}) \equiv -\frac{\partial}{2\partial \ln \rho_{*}} \ln \left( \frac{dN_{\text{det}}(r, \rho > \rho_{*}, m_{1,2})}{d\Omega dr dm_{1,2}} \right), \quad (9)$$

$$A(r, m_{1,2}) \equiv \frac{1}{\mathcal{F}\left(\frac{f_{\text{ISCO}}}{(1+z)}\right)} \left(\frac{2f_{\text{ISCO}}}{(1+z)}\right)^{-7/3} S_n\left(\frac{2f_{\text{ISCO}}}{(1+z)}\right)^{-1} \frac{2f_{\text{ISCO}}}{(1+z)} \qquad (10)$$

Inserting equations (8) and (2) into equation (3) and keeping only linear terms in $v_0/c$ we finally obtain for the mean luminosity distance

$$d_{L}^{\text{mean}}(\mathbf{n}) = d_{L}^{(0)} \left( 1 + \mathcal{D} d_{L}(\mathbf{n}) \right), \qquad (11)$$

where the monopole is given by

$$d_{L}^{(0)} \equiv \frac{1}{4\pi} \int d\Omega \, dL_{\text{mean}}(\mathbf{n}). \qquad (12)$$

The fractional dipole is given by

$$\mathcal{D}d_L(\mathbf{n}) = -\alpha_{d_L} \left( \mathbf{n} \cdot \frac{\mathbf{v}_0}{c} \right), \qquad (13)$$

where the pre-factor $\alpha_{d_L}$ takes the form

$$\alpha_{d_L} \equiv 1 - \int dr \, dm_{1,2} s \left( \frac{1}{3} + \mathcal{A} \right) \left[ \frac{\bar{d}_L(r)}{d_L^{(0)}} - 1 \right] p(r, m_{1,2}, \rho_*) , \qquad (14)$$

with $p(r, m_{1,2}, \rho_*)$ being the mean probability density function of detected sources,

$$p(r, m_{1,2}, \rho_*) = \frac{\frac{d\bar{N}_{\text{det}}}{d\Omega dr dm_{1,2}}(r, m_{1,2}, \rho > \rho_*)}{\int dr \, dm_{1,2} \frac{d\bar{N}_{\text{det}}}{d\Omega dr dm_{1,2}}(r, m_{1,2}, \rho > \rho_*)} . \qquad (15)$$

Note that in equation (13) we have introduced the parameter $\alpha_{d_L}$ factorizing out a minus sign in such a way that $\alpha_{d_L} = 1$ in the case where there are no threshold effects, i.e. when $s = 0$.

We see from equation (14) that if there are threshold effects but $s$ and $\mathcal{A}$ are constant in $r$, then the integral vanishes and $\alpha_{d_L} = 1$. We can indeed rewrite the monopole as

$$d_{L}^{(0)} = \frac{1}{4\pi} \int d\Omega \left[ \frac{\int dr \, dm_{1,2} dL(r, \mathbf{n}) \frac{dN_{\text{det}}}{d\Omega dr dm_{1,2}}(r, \mathbf{n}, m_{1,2}, \rho > \rho_{*})}{\int dr \, dm_{1,2} \frac{dN_{\text{det}}}{d\Omega dr dm_{1,2}}(r, \mathbf{n}, m_{1,2}, \rho > \rho_{*})} \right] \\ = \int dr \, dm_{1,2} \bar{d}_L(r) p(r, m_{1,2}, \rho_{*}), \qquad (16)$$

since the dipole contributions in the luminosity distance and in the number of detected events vanish when integrated over direction (note that here we neglect terms or order $(v_0/c)^2$). Inserting this into equation (14), we see that the integral exactly vanishes if $s$ and $\mathcal{A}$ are independent on $r$. Since these functions are expected to evolve slowly with $r$, we expect threshold effects to be partially suppressed, meaning that $\alpha_{d_L}$ will be close to 1, even in the case where $s$ is non-zero. In Section 4.4, we compute $\alpha_{d_L}$ from our population model of BNSs (where threshold effects are present), and compare this with measurements of $\alpha_{d_L}$ from our synthetic catalogues. We find that $\alpha_{d_L}$ can be well predicted and, as expected, that it is very close to 1. This is important, because $\alpha_{d_L}$ is fully degenerated with the observer velocity $v_0/c$ and if we do not know it, we cannot extract information on $v_0/c$ from a population of sources with threshold effects.

## 2.2 Chirp mass

The chirp mass that can be extracted from the GW is the redshifted chirp mass, i.e. the product of the intrinsic chirp mass of the binary system and the redshift. Since the latter is affected by the observer

<footnote n="3">³Note that here we have dropped the dependence on the SNR threshold, $\rho_*$, in the left-hand side of equation (4) for simplicity. When we refer to 'detected' events, by definition we refer to events above the SNR threshold. We state the dependence explicitly only when it is not integrated over masses.</footnote>

<footnote n="4">⁴Note that these functions depend on $m_{1,2}$, while in Mastrogiovanni et al. (2023), $s$ and $\mathcal{A}$ refer to mass integrated quantities</footnote>
velocity, to first order, the redshifted chirp mass is given by

$$ \mathcal{M}(r, \mathbf{n}, m_{1,2}) = \bar{\mathcal{M}}(r, m_{1,2}) \left(1 - \mathbf{n} \cdot \frac{\mathbf{v}_0}{c}\right), \qquad (17) $$

where $\bar{\mathcal{M}}(r, m_{1,2})$ is the background (redshifted) chirp mass in a homogeneous and isotropic universe. Note that, as for the luminosity distance, other perturbations contribute to the redshifted chirp mass, that can be neglected since they have a negligible impact on the dipole. Analogously to equation (3), we can compute the mean detected chirp mass in direction $\mathbf{n}$, averaged over all binary systems (i.e. over all masses $m_{1,2}$) as

$$ \mathcal{M}^{\text{mean}}(\mathbf{n}) = \frac{\int dr \, dm_{1,2} \mathcal{M}(r, \mathbf{n}, m_{1,2}) \frac{dN_{\text{det}}}{d\Omega dr \, dm_{1,2}}(r, \mathbf{n}, m_{1,2})}{\int dr \, \frac{dN_{\text{det}}}{d\Omega dr}(r, \mathbf{n})} \qquad (18) $$

This can be expanded in powers of $\mathbf{n} \cdot \frac{\mathbf{v}_0}{c}$ as a monopole, $\mathcal{M}^{(0)}$, and dipole term:

$$ \mathcal{M}^{\text{mean}}(\mathbf{n}) = \mathcal{M}^{(0)}(1 + \mathcal{D}\mathcal{M}(\mathbf{n})), \qquad (19) $$

where

$$ \mathcal{M}^{(0)} = \frac{1}{4\pi} \int d\Omega \, \mathcal{M}^{\text{mean}}(\mathbf{n}). \qquad (20) $$

Combining equations (8) and (17) allows to write the fractional dipole term as

$$ \mathcal{D}\mathcal{M}(\mathbf{n}) = -\alpha_{\mathcal{M}} \left(\mathbf{n} \cdot \frac{\mathbf{v}_0}{c}\right), \qquad (21) $$

where we introduced $\alpha_{\mathcal{M}}$ as the chirp mass analogue of $\alpha_{d_L}$:

$$ \alpha_{\mathcal{M}} = 1 - \int dr \, dm_{1,2} s \left(\frac{1}{3} + \mathcal{A}\right) \left[ \frac{\bar{\mathcal{M}}(r, m_{1,2})}{\mathcal{M}^{(0)}} - 1 \right] p(r, m_{1,2}, \rho_s). \qquad (22) $$

As for the luminosity distance, we see that threshold effects are non-zero only if $s$ and $\mathcal{A}$ vary with distance $r$.

Kashyap et al. (2023) have also computed the dipole in the mean chirp mass (that they call mass intensity). Their modelling differs from ours, since they assume that all events above a given mass threshold are detected. In practice, this is however not the case. What determines if an event is detected or not is the SNR, which depends not only on the mass of the binary system, but also on its distance from the observer, as can be seen from equation (5). This has an impact on the modelling of threshold effects, leading to a different expression for the dipole signal.

## 2.3 Number counts

The dipole in the GWs number count has been derived in Mastrogi-vanni et al. (2023) and follows directly from equation (8). The total number of GW events detected in direction $\mathbf{n}$ is given by

$$ \begin{aligned} N_{\text{det}}(\mathbf{n}) &= \int dr \, dm_{1,2} \frac{dN_{\text{det}}}{d\Omega dr \, dm_{1,2}}(r, \mathbf{n}, m_{1,2}) \\ &= N_{\text{det}}^{(0)}(1 + \mathcal{D}N_{\text{det}}(\mathbf{n})), \end{aligned} \qquad (23) $$

where

$$ \mathcal{D}N_{\text{det}}(\mathbf{n}) = 2\alpha_N \mathbf{n} \cdot \frac{\mathbf{v}_0}{c}, \qquad (24) $$

with

$$ \alpha_N = 1 + \frac{1}{2} \int dr \, dm_{1,2} s \left(\frac{1}{3} + \mathcal{A}\right) p(r, m_{1,2}, \rho_s). \qquad (25) $$

# 3 STATISTICAL ESTIMATORS OF THE DIPOLES

We now build statistical estimators for the three dipole signals defined in Section 2.

## 3.1 Luminosity distance

For each direction $\mathbf{n}'$ in the sky, we can build the following observable

$$ v_{d_L-n'} = -3 \int \frac{d\Omega}{4\pi} (\mathbf{n} \cdot \mathbf{n}') \mathcal{D}d_L(\mathbf{n}) \approx \alpha_{d_L} \frac{v_0}{c} \cos \theta', \qquad (26) $$

where $\theta'$ is the angle between the (a priori unknown) dipole velocity direction and $\mathbf{n}'$. This observable is maximized when evaluated along the dipole direction and, in the absence of threshold effects, it is exactly equal to the observer velocity, $v_0/c$, at the maximum.

Let us now build a statistical estimator for the observable $v_{d_L-n'}$. We divide the sky in $N_{\text{sky}}$ angular pixels of same solid angle, and associate a vector $\mathbf{n}_i$ pointing to the centre of each pixel $i$. Within any pixel $i$, there are $N_{\text{det}}^i$ detected events (labelled by $j$), whose corresponding luminosity distance is labelled by $(d_L)_j^i$. The estimator is then defined as

$$ \hat{v}_{d_L-n'} \equiv \frac{-3}{\hat{d}_L N_{\text{sky}}} \sum_{i=1}^{N_{\text{sky}}} \sum_{j=1}^{N_{\text{det}}^i} \frac{(d_L)_j^i}{N_{\text{det}}^i} (\mathbf{n}_i \cdot \mathbf{n}') \qquad (27) $$

with

$$ \hat{d}_L \equiv \frac{1}{N_{\text{sky}}} \sum_{i=1}^{N_{\text{sky}}} \sum_{j=1}^{N_{\text{det}}^i} \frac{(d_L)_j^i}{N_{\text{det}}^i}. \qquad (28) $$

As shown in Appendix A, if shot noise and the uncertainty in the measurement of the luminosity distance are subdominant (i.e. smaller than their mean), the estimator (27) is unbiased, i.e. $\langle \hat{v}_{d_L-n'} \rangle = v_{d_L-n'}$.

We then compute the variance of this estimator. It can be written as ratio of two variables $X/Y$ with

$$ X = \frac{-3}{N_{\text{sky}}} \sum_{i=1}^{N_{\text{sky}}} \sum_{j=1}^{N_{\text{det}}^i} \frac{(d_L)_j^i}{N_{\text{det}}^i} (\mathbf{n}_i \cdot \mathbf{n}') \qquad (29) $$

and $Y = \hat{d}_L$ defined in equation (28). The variance of such a ratio can be written as (if the variances of $X$ and $Y$ are significantly smaller than their mean)

$$ \text{var}\left(\frac{X}{Y}\right) \approx \frac{\langle X \rangle^2}{\langle Y \rangle^2} \left[ \text{var}(X) + \frac{\text{var}(Y)}{\langle X \rangle^2} - 2 \frac{\text{cov}(X, Y)}{\langle X \rangle \langle Y \rangle} \right]. \qquad (30) $$

We can show that the third term in equation (30) is vanishing and the second one is smaller than the first one by a factor $(v_0/c)^2$ and can therefore be neglected (see Appendix A for a detailed derivation). The variance of (27) is therefore directly proportional to the variance of $X$ and can be written as

$$ \begin{aligned} \text{var}(\hat{v}_{d_L-n'}) &= \langle \hat{v}_{d_L-n'}^2 \rangle - \langle \hat{v}_{d_L-n'} \rangle^2 \\ &= \frac{9}{(d_L^{(0)})^2} N_{\text{sky}}^2 \sum_{i=1}^{N_{\text{sky}}} \sum_{a=1}^{N_{\text{sky}}} \sum_{j=1}^{N_{\text{det}}^i} \sum_{b=1}^{N_{\text{det}}^a} (\mathbf{n}_i \cdot \mathbf{n}') (\mathbf{n}_a \cdot \mathbf{n}') \\ &\quad \times \left( \left\langle \frac{(d_L)_j^i (d_L)_b^a}{N_{\text{det}}^i N_{\text{det}}^a} \right\rangle - \left\langle \frac{(d_L)_j^i}{N_{\text{det}}^i} \right\rangle \left\langle \frac{(d_L)_b^a}{N_{\text{det}}^a} \right\rangle \right). \end{aligned} \qquad (31) $$

To compute this, we divide each solid angle $i$ into $N_r$ bins in distance $r$ of size $\Delta r$, and we rewrite the sum over $j$ as a sum over the $r$-bins. The sky is therefore divided into pixels of angular size $\Delta\Omega$
<figure><img src="image_1.png" /><figcaption>**Figure 1.** Variance of the luminosity distance estimator for BBHs, plotted as a function of the total number of events $N_{\text{tot}}$. The different panels are for different values of the relative error on $d_L$ (assumed to be independent of redshift). We show separately the shot noise contribution and the distance-uncertainty contribution, as well as the total.</figcaption></figure>

and radial size $\Delta r$. We denote by $N_p^{ik}$ the number of events detected in the pixel centred in direction $\mathbf{n}_i$ and distance $r_k$ and by $(d_L)_k^i$ the corresponding luminosity distance. Due to shot noise, the number of events fluctuates around the mean from pixel to pixel,

$$N_p^{ik} = \frac{dN_{\text{det}}(r_k, \mathbf{n}_i)}{d\Omega dr} \Delta r \Delta \Omega = \bar{N}_p^k + \Delta N_p^{ik}, \quad (32)$$

$$N_{\text{det}}^i = \bar{N}_{\text{det}} + \sum_j \Delta N_p^{ij}. \quad (33)$$

Note that here we neglect the fluctuations in $N_p^{ik}$ due to the uncertainty on sky localization. We assume that this uncertainty is of ~3 degrees, significantly smaller than the size of the angular pixels that we consider. In addition to shot noise, the variance of the estimator is affected by uncertainty in the measurement of the luminosity distance, which we write as

$$(d_L)_k^i = (\bar{d}_L)_k + (\Delta d_L)_k^i, \quad (34)$$

where $(\Delta d_L)_k^i$ is the error in the measurement of $d_L$ from one binary system, in angular bin $i$ and radial bin $k$. Note that in equations (32)-(34) we neglect the contributions from the observer velocity that would lead to subdominant contributions to the variance. The main contribution comes indeed from a 'fluctuation' of the monopole (due to shot noise and measurement uncertainty) that would mimic a dipole. Using the fact that errors on count and on luminosity distance are uncorrelated and that number counts follow a Poissonian statistics, we obtain for the variance (see Appendix A for more detail)

$$\text{var}(\hat{v}_{d_L-n'}^2) = \frac{3N_{\text{sky}}}{(d_L^{(0)})^2 N_{\text{tot}}^2} \left\{ \sum_{k=1}^{N_r} \langle (\Delta d_L)_k^2 \bar{N}_p^k \rangle + \sum_{k=1}^{N_r} (\bar{d}_L)_k^2 \bar{N}_p^k - \frac{1}{\bar{N}_{\text{det}}} \left( \sum_{k=1}^{N_r} (\bar{d}_L)_k \bar{N}_p^k \right)^2 \right\}, \quad (35)$$

where $N_{\text{tot}}$ is the total number of events detected and $\langle (\Delta d_L)_k \rangle$ denotes the typical error on luminosity distance measurement at distance $r_k$.

The term in the first line of equation (35) is the contribution from luminosity distance uncertainty to the variance of the estimator. If $\langle (\Delta d_L)_k \rangle$ is independent on distance, we see that the variance from this term scales as $1/N_{\text{tot}}$, as expected for $N_{\text{tot}}$ measurement of uncorrelated quantities. The contributions in the second line are due to shot noise: the mean distance in a given solid angle depends on the radial distribution of events. Since shot noise generates fluctuations in the number of events, a given solid angle can have more events closer to the observer, whereas another solid angle can have more events further away. This generates fluctuations in the mean distance, which can mimic the presence of a dipole. Note that since we are measuring the mean distance per angular pixel and not the sum of distances per pixel, we are not sensitive to the impact of shot noise on the total number of events per pixel, but rather to its impact on the radial distribution of sources.

Taking the continuous limit, the variance can be rewritten as

$$\text{var}(\hat{v}_{d_L-n'}) = \frac{3}{(d_L^{(0)})^2 N_{\text{tot}}} \times \int dr p(r, \rho_*) \left[ \langle \Delta d_L^2(r) \rangle + \bar{d}_L^2(r) - (d_L^{(0)})^2 \right], \quad (36)$$

where $p(r, \rho_*)$ is the radial distribution of detected sources integrated over masses. In Fig. 1, we plot the variance (36) for BBHs as a function of $N_{\text{tot}}$ for different relative errors on the luminosity distance measurements, ranging from 10 per cent to 100 per cent (and constant in $r$). We use the radial distribution of sources obtained from our simulated catalogues (see Section 4 for more detail). We see that both the shot noise contribution and the distance-uncertainty contribution scale as $1/N_{\text{tot}}$. The relative importance of the two terms depends therefore only on the uncertainty on $d_L$. For an uncertainty of up to 20 per cent, the shot noise contribution completely dominates in the variance. If the uncertainty reaches 50 per cent however, its contribution to the variance is not negligible anymore and the SNR is degraded. Similar results are obtained for BNSs, see Fig. A1 in Appendix A. The only difference comes from the radial distribution of sources, which differs for BBHs and BNSs and leads to a slightly larger variance for BBHs (due to a wider redshift range of observation for BBHs, which can be observed at higher redshift than BNSs). In the following we will assume a 20 per cent uncertainty on the measurement of $d_L$ as a proxy of typical distance uncertainties for GW events (Callister et al. 2020; Iacovelli et al. 2022), meaning that we are in the regime where shot noise completely dominates.

## 3.2 Chirp mass

The same procedure can be applied to the chirp mass. We first build an observable equivalent to (26),

$$v_{M-n'} \equiv -3 \int \frac{d\Omega}{4\pi} (\mathbf{n} \cdot \mathbf{n}') \mathcal{D}\mathcal{M}(\mathbf{n}) \approx \alpha_M \frac{v_0}{c} \cos \theta'. \quad (37)$$

We repeat the process of dividing the detected observations in sky angular pixels, as in equation (27), to associate the following statistical estimator to $v_{M-n'}$:

$$\hat{v}_{M-n'} \equiv \frac{-3}{\hat{M} N_{\text{sky}}} \sum_{i=1}^{N_{\text{sky}}} \sum_{j=1}^{N_{\text{det}}} \frac{(M_j^i)^2}{N_{\text{det}}^i} (\mathbf{n}_i \cdot \mathbf{n}') \quad (38)$$

with

$$\hat{M} \equiv \frac{1}{N_{\text{sky}}} \sum_{i=1}^{N_{\text{sky}}} \sum_{j=1}^{N_{\text{det}}} \frac{(M_j^i)^2}{N_{\text{det}}^i}. \quad (39)$$

As for the luminosity distance, this estimator is found to be unbiased assuming that shot noise and uncertainty in the measurement of
<figure><img src="image_2.png" /><figcaption>**Figure 2.** Different contributions to the variance of the chirp mass estimator, for a relative measurement error $(\Delta M^2(r, m_{1,2})/\bar{M}^2(r, m_{1,2}) = 1$ per cent (top), 10 per cent (bottom). In any case, shot noise is the fully dominant variance contribution, the measurement error giving a smaller contribution by orders of magnitude.</figcaption></figure>

the chirp mass are subdominant compared to the mean. For the variance, we generalize the computation of Section 3.1 and Appendix A done for the luminosity distance to have data not only distributed in $N_r$ radial bins, but also in $N_m$ bins of the two source masses, $m_{1,2}$, which constitute the binary. This is needed since, unlike the luminosity distance, the redshifted chirp mass function depends on all $(r, m_{1,2})$. As for the luminosity distance, we neglect subdominant contributions from the observer velocity in the variance. We obtain

$$ \text{var}(\hat{v}_{M-n'}) = \frac{3N_{\text{sky}}}{(M^{(0)})^2 N_{\text{tot}}^2} \left\{ \sum_{k=1}^{N_r} \sum_{\ell=1}^{N_m} ((\Delta M_{k\ell})^2) \bar{N}_p^{k\ell} + \sum_{k=1}^{N_r} \sum_{\ell=1}^{N_m} (\bar{M}_{k\ell})^2 \bar{N}_p^{k\ell} - \frac{1}{N_{\text{det}}} \left[ \sum_{k=1}^{N_r} \sum_{\ell=1}^{N_m} \bar{M}_{k\ell} \bar{N}_p^{k\ell} \right]^2 \right\}, \quad (40) $$

where the k-index is the radial bin index, while the l-index is the index over the various $m_{1,2}$ bins. $\bar{M}_{k\ell}$ is the expected redshifted chirp mass in such a bin, while $\bar{N}_p^{k\ell}$ is the expected number of detected events in that bin, i.e. such that $\bar{N}_p^{k\ell} = \sum_{\ell=1}^{N_m} \bar{N}_p^{k\ell}$. The quantity $\Delta M_{k\ell}$ denotes the typical uncertainty in the measurement of the chirp mass in a radial bin centred at $r_k$ and in the mass bin labelled by $l$.

Taking the continuous limit, the variance (40) becomes

$$ \text{var}(\hat{v}_{M-n'}) = \frac{3}{(M^{(0)})^2 N_{\text{tot}}} \int dr \, dm_{1,2} \, p(r, m_{1,2}, \rho_*) \times \left[ \langle (\Delta M)^2(r, m_{1,2}) \rangle + \bar{M}^2(r, m_{1,2}) - (M^{(0)})^2 \right]. \quad (41) $$

As for the luminosity distance, we see that the variance contains two contributions: one from measurement uncertainty of the chirp mass, and the second one from shot noise. In addition to the fact that shot noise affects the radial distribution of events in a given solid angle, it will also affect the distribution in masses $m_{1,2}$. This effect is encoded in the dependence of the distribution function $p$ on the masses and it will also generate fluctuations in the mean chirp mass that can mimic a dipole. In Fig. 2, we plot the different contributions as a function of $N_{\text{tot}}$ assuming a 10 per cent uncertainty

in the measurement of the chirp mass (expected for BBHs) and a 1 per cent uncertainty (expected for BNSs), see Iacovelli et al. (2022). We see that in both cases shot noise completely dominates over mass uncertainty. Comparing the shot noise contribution for BBHs and BNSs, we see that it is significantly smaller for BNSs, by a factor 3. This is due to the fact that BNSs have a narrower mass range than BBHs: the mass distribution of BNSs only span $2M_\odot$, while BBHs are observed between 5 and $100M_\odot$. Since shot noise changes the mass distribution of events, it has more impact in the second case. For example, a shot noise fluctuation generating one more event at the higher end of the mass range will more drastically affect the mean mass of BBHs than the mean mass of BNSs.

### 3.3 Number counts

The dipole estimator for the number counts and its variance have been derived in detail in Mastroiovanni et al. (2023). Here we just state the final result for completeness:

$$ v_{N-n'} \equiv \frac{3}{2} \int \frac{d\Omega}{4\pi} (\mathbf{n} \cdot \mathbf{n}') \mathcal{D}N(\mathbf{n}) \approx \alpha_N \frac{v_0}{c} \cos \theta', \quad (42) $$

$$ \hat{v}_{N-n'} = \frac{3}{2N_{\text{tot}}} \sum_{i=1}^{N_{\text{sky}}} N_{\text{det}}^i (\mathbf{n}_i \cdot \mathbf{n}') . \quad (43) $$

The variance of the number count estimator is only due to shot noise and is given by

$$ \text{var}(\hat{v}_{N-n'}) = \frac{3}{4N_{\text{tot}}}. \quad (44) $$

### 3.4 Covariance of the estimators

In order to use the distance, mass and number counts estimators together to measure the dipole from the same population of GW sources, it is necessary to compute their covariance. Below we show that the distance and mass estimators are correlated, but that these estimators are both uncorrelated with the number count estimator. Naturally, the estimators for two different populations of sources (for example BBHs and BNSs) are uncorrelated.

#### 3.4.1 Mass-distance covariance

The covariance is calculated following the same steps as for the variance in Sections 3.1 and 3.2. Here we assume that the uncertainty in the measurement of the chirp mass is uncorrelated with the uncertainty in the measurement of the distance. The chirp mass is mainly measured from the phase, whereas the luminosity distance is measured through the amplitude of the wave. Hence, this is certainly a good approximation for BNSs, which have a long inspiral in the detectors' band and thus a good phase determination. Moreover, inspection of typical BBHs signals shows that the $d_L-\mathcal{M}_c$ entries of the normalized Fisher matrix are generally smaller than the diagonal ones, which we have already shown to give subdominant contributions to the total noise (see Figs 1 and 2). This justifies the assumption of uncorrelated chirp mass and distance measurements uncertainties as well for the BBH type of sources.

We obtain for the covariance

$$ \text{cov}(\hat{v}_{d_L-n'}, \hat{v}_{M-n'}) = \frac{3N_{\text{sky}}}{M^{(0)} d_L^{(0)} N_{\text{tot}}^2} \left\{ \sum_{k=1}^{N_r} \sum_{\ell=1}^{N_m} \bar{M}_{k\ell}(\bar{d}_L)_k \bar{N}_p^{k\ell} - \frac{1}{N_{\text{det}}} \left[ \sum_{k=1}^{N_r} \sum_{\ell=1}^{N_m} \bar{M}_{k\ell} \bar{N}_p^{k\ell} \right] \left[ \sum_{j=1}^{N_r} (\bar{d}_L)_j \bar{N}_p^j \right] \right\}. \quad (45) $$
<figure><img src="image_3.png" /><figcaption>**Figure 3.** Value of the dipole estimators as a function of the angle between the true dipole direction and a chosen direction **n'**. The plots are generated using the AGN value of the observer velocity, and with $10^6$ BBH detections (left panel) and $10^6$ BNS detections (right panel). The grey shaded areas are the $1\sigma$, $2\sigma$, $3\sigma$ values associated with the shot noise and measurement uncertainties in the absence of dipole, obtained by shuffling the sources isotropically over the sky one hundred times. The horizontal dashed lines mark the theoretical expectations for the variance obtained in Section 3. Top plot: Number count estimator. Middle plot: Chirp mass estimator. Bottom plot: Luminosity distance estimator.</figcaption></figure>

Taking the continuous limit we obtain

$$ \operatorname{cov}(\hat{v}_{d_L-n'}, \hat{v}_{M-n'}) = \frac{3}{M^{(0)} d_L^{(0)} N_{\text{tot}}} \int dr \, dm_{1,2} \, p(r, m_{1,2}, \rho_*) \times [\bar{M}(r, m_{1,2}) \bar{d}_L(r) - M^{(0)} d_L^{(0)}]. \quad (46) $$

From equation (46) we see that the two estimators are correlated through shot noise. Indeed, if in a given solid angle there are more events closer to the observer, the mean distance in that solid angle will be smaller than on average, and the mean chirp mass as well. Hence a shot noise fluctuation that can mimic a dipole in one estimator will automatically mimic a dipole in the other estimator as well.

### 3.4.2 Count-distance and count-mass covariance

We start by computing the covariance between the number counts and the distance estimators. Since shot noise is uncorrelated with the measurement uncertainty in the distance, we obtain, following the same steps as previously

$$ \operatorname{cov}(\hat{v}_{d_L-n'}, \hat{v}_{N-n'}) = \frac{9}{2} \frac{1}{d_L^{(0)} N_{\text{tot}} N_{\text{sky}}} \sum_{i=1}^{N_{\text{sky}}} \sum_{a=1}^{N_{\text{sky}}} \sum_{\ell=1}^{N_r} (\mathbf{n}_i \cdot \mathbf{n}') (\mathbf{n}_a \cdot \mathbf{n}') \\ (\bar{d}_L)_\ell \bar{N}_p^\ell \left[ \sum_{q=1}^{N_r} \left\langle \frac{\Delta N_p^{aq}}{\bar{N}_{\text{det}}} \frac{\Delta N_p^{il}}{\bar{N}_p^\ell} \right\rangle - \sum_{q=1}^{N_r} \sum_{t=1}^{N_r} \left\langle \frac{\Delta N_p^{aq}}{\bar{N}_{\text{det}}} \frac{\Delta N_p^{it}}{\bar{N}_{\text{det}}} \right\rangle \right]. \quad (47) $$

Using that the number of detected events follows a Poisson distribution, we can easily show that the two terms in equation (47) exactly vanish. A similar calculation shows that the covariance between mass

and number counts also vanishes. This is not surprising: contrary to the distance and mass estimators, the number count estimator is not sensitive to the radial distribution of events. Only fluctuations in the total number of events in a solid angle can mimic a dipole in the number count. On the contrary, the mass and distance estimators are not sensitive to the total number of events in a solid angle (since we divide by the total number of events to obtain the mean) but rather to their radial distribution. As a consequence, a shot noise fluctuation that mimics a dipole in the number count estimator does not necessarily mimic a dipole in the mass and distance estimator, and vice versa.

## 3.5 Optimal combination

In the case where there are no threshold effects, the six estimators (three for BBHs and three for BNSs) are estimators of the same quantity: $v_0/\text{ccos } \theta'$. We can therefore look for an optimal estimator of this quantity, i.e. an estimator that would maximize the SNR. In the following we drop the dependence in **n'**, since the optimal estimator can be defined in exactly the same way for each value of **n'**.

As a first step, we build combinations of estimators that are independent. Since only the distance and the mass estimators are correlated for the same population of sources, we simply need to diagonalize that part of the covariance matrix. We obtain two new estimators (per population) given by

$$ \hat{v}_{\pm} = \frac{X^{\pm} \hat{v}_{d_L} + \hat{v}_{M}}{1 + X^{\pm}}, \quad (48) $$

with

$$ X^{\pm} = \frac{1}{2\operatorname{cov}(\hat{v}_{d_L}, \hat{v}_{M})} \times \left\{ \operatorname{var}(\hat{v}_{d_L}) - \operatorname{var}(\hat{v}_{M}) \pm \sqrt{(\operatorname{var}(\hat{v}_{M}) - \operatorname{var}(\hat{v}_{d_L}))^2 + 4 \operatorname{cov}^2(\hat{v}_{d_L}, \hat{v}_{M})} \right\}. \quad (49) $$

The denominator in equation (48) insures that the new estimators have mean $v_0/\text{ccos } \theta'$.

We can now build an optimal estimator of $v_0/\text{ccos } \theta'$ by linearly combining the independent estimators in the following way

$$ \hat{v}_{\text{ideal}} = \frac{\sum_{i=1}^{6} c_i \hat{v}_i}{\sum_{i=1}^{6} c_i}, \quad c_i = (\operatorname{var}(\hat{v}_i))^{-1}, \quad (50) $$

with $i$ labelling all the possible observables (number count, plus and minus combinations) for both BBHs and BNSs. This combination is the one that maximizes the SNR (in every direction **n'**), defined as

$$ (\text{SNR})^2 = \frac{\langle \hat{v}_{\text{ideal}} \rangle^2}{\operatorname{var}(\hat{v}_{\text{ideal}})}. \quad (51) $$

Note that the normalization in equation (50) is there to keep the condition $\langle \hat{v}_{\text{ideal}} \rangle = v_0/c \cos \theta'$, but it is otherwise irrelevant.

If threshold effects are relevant, the mass, distance and number count estimators are not estimators of the same quantity, since the respective $\alpha$'s are different. In this case, we need to divide each estimator by its respective $\alpha$ before combining them. Since the $\alpha$'s can be determined only with limited precision (using population models, see Section 4.4 for more detail), this would induce additional contributions to the variance of the optimal estimator and degrades the SNR. In Section 5.2, we quantify this effect using a Fisher analysis.
# 4 MEASUREMENT OF THE DIPOLE FROM SYNTHETIC CATALOGUES OF GW SOURCES

To test our estimators we build synthetic catalogues of BBHs and BNSs and use our estimators on these simulated events. We compare the measurement with our theoretical modelling for the signal and the variance.

## 4.1 Simulating BBHs and BNSs

To build catalogues of BBH and BNS sources, we draw source frame masses for BBHs and BNSs from the same probabilistic models used in Mastrogiovanni et al. (2023; see appendix A therein). These mass models are consistent with current BBH and BNS detections (Iacovelli et al. 2022). The redshift distribution of GW sources is determined by the merger rate model. In our simulation, the merger redshift is distributed according to

$$p(z) \propto [1 + (1 + z_p)^{-\gamma-k}] \frac{(1+z)^\gamma}{1 + \left(\frac{1+z}{1+z_p}\right)^{\gamma+k}} \frac{1}{1+z} \frac{dV_c}{dz}, \quad (52)$$

where $\gamma$, $k$, and $z_p$ are parameters that control the merger rate model, while $\frac{dV_c}{dz}$ is the differential of the comoving volume. As fiducial values, we take $\gamma = 2.7$, $k = 3$, and $z_p = 2$. The distribution of $\cos \iota$ is chosen to be uniform. Once a set of BBHs and BNSs is drawn, we add the effect of the observer velocity. Aberration is included by shifting the angular position by $\theta' = \theta - v_0/\sin\theta$, where $\theta$ is the angle between the source position and the observer velocity. The luminosity distance and chirp mass are modified according to equations (2) and (17) respectively. We produce two copies of both the BBHs and BNSs catalogue, one with the ‘CMB value’ of the observer velocity, $v_0/c = 1.2 \times 10^{-3}$ and the other one with the ‘AGN value’ of the observer velocity, which is 5 times larger: $v_0/c = 6 \times 10^{-3}$.

We can then calculate the SNR of each event using equation (5). We consider a network of XG detectors including ET (Punturo et al. 2010) and two CE (Dwyer et al. 2015; Reitze et al. 2019). The power spectral density of ET is set to the one used in Iacovelli et al. (2022); Mastrogiovanni et al. (2023), while for CE we take the power spectral density taken from the CE consortium<fnref n="5" /> If the SNR from this network exceeds a detection threshold of $\rho_* = 9$, we label the binary as detected.

Finally, we add an extra step in the simulation to mimic the fact that we will not be perfectly able to measure the sky location, luminosity distance and redshifted chirp mass of the source. Once a binary is detected, we do not save its true values for the sky position, distance and chirp mass but instead, we register a scattered value around the true one. We include Gaussian scatter of the sky location by 3 degrees around the true sky position and of the luminosity distance by 20 per cent of its true value. For the chirp mass, we use a scattering of 1 per cent for BNSs and of 10 per cent for BBHs. These are typical errors that we might expect to obtain with XG detectors (Iacovelli et al. 2022). Note that the 20 per cent uncertainty on the luminosity distance implies that shot noise is the dominant contribution, as can be seen from Figs 1 and A1.

The luminosity distances, redshifted chirp masses and sky positions are then used to compute the estimators using equations (27), (38), and (42). The results are shown in Fig. 3 for the case of the AGN observer velocity. We show the values of the estimator obtained for $10^6$ BBH detections (left panel) and $10^6$ BBN detections (right panel) as a function of the angle $\theta'$ between the true direction of the dipole and a chosen direction $\mathbf{n}'$. The spread in the signal comes from the fact that for a given angle $\theta'$ we have pixels in different azimuthal directions, which give slightly different values for the dipole estimator (due to shot noise and measurement uncertainties). At the equator there are clearly more azimuthal pixels than at the pole (exactly at the pole there is just one), leading to a larger spread at the equator.

As discussed in Mastrogiovanni et al. (2023), for BBHs threshold effects do not contribute to the amplitude of the dipole, since all events are above the threshold [see fig. 2 of Mastrogiovanni et al. (2023)]. As a consequence, at the maximum, i.e. when $\mathbf{n}'$ coincides with the direction of the observer velocity, the estimators are roughly equal to $v_0/c = 6 \times 10^{-3}$. The peak of the dipole estimators could be slightly shifted from the true position due to the variance. For BNSs, we see that the amplitude of the dipole at the peak is also very close to the observer velocity. This is due to the fact that, as we will show in Section 4.4, the amplitude of threshold effects are actually small for BNSs, below 10 per cent, which is smaller than the spread in the signal.

Fig. 3 also reports the $1\sigma$, $2\sigma$, $3\sigma$ fluctuations (grey areas) of the estimators due to shot noise and measurement uncertainties on the luminosity distance and chirp mass. The fluctuation levels are generated by shuffling GW detections isotropically over the sky a hundred times. We see that the chirp mass from BNSs is the estimator with smaller variance, consistent with the theoretical results of Section 3. For all estimators, the fluctuation levels obtained through sky shuffling agree with the theoretical calculation of the variance from equations (36), (41) and (44), that are indicated with dashed lines on the plot (see Section 4.2 for a more detailed comparison). When the estimator values exceed a certain noise threshold, the cosmic dipole can be detected (see Section 4.3 for a more in-depth discussion on detectability).

## 4.2 Comparison of the theoretical variance and covariance with simulations

The first sanity check that we perform is to see if our predictions of the variance and covariance agree with the numerical simulations. As explained above, the numerical variance and covariance are simply obtained by shuffling the sources isotropically over the sky. This removes the true dipole signal, meaning that the remaining fluctuations are due to shot noise and measurement uncertainties of the luminosity distance and chirp mass. Here we use 2'000 sky shufflings to obtain an accurate numerical estimate of the variance and covariance.

Theoretically, the variance of the number count estimator is given by equation (44), and simply depends on the total number of detected events. On the other hand, the variance of the mass and distance estimators as well as the covariance between them has to be computed from the integrals (36), (41) and (46), which requires a model for the $(r, m_{1,2})$ source distributions. To estimate this, we numerically generate GW events following reference mass and redshift distributions models, and approximate the integrals in (36), (41) and (46) by sampling those distributions. The sampled number of detections is increased until we reach numerical convergence.

The top plots of Fig. 4 show the simulated variance compared with the theoretical one, for BBHs and BNSs, as a function of the total number of events $N_{\text{tot}}$. The agreement between the simulated and theoretical variances is excellent. As expected from equations (36), (41) and (44) the variances scale as $N_{\text{tot}}^{-1}$. For the number count estimator, the variance is the same for BBHs and BNSs: it depends

<footnote n="5">⁵https://cosmicexplorer.org/sensitivity.html</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 4.** The plots report the variance (first row) and correlation (second row) of the number count, mass, and distance dipole estimators as a function of the number of detections $N_{\text{tot}}$. The variances and correlations are obtained by reshuffling 2000 times isotropically a population of GW sources. The blue circles indicate the values obtained for the BBH population while the orange diamonds the values for the BNS population. The black lines (dashed for BBHs and dotted for BNSs) indicate the theoretical variance calculated in Section 3. We find that the maximum deviation from the expected value of zero for the *number count-chirp mass* and *number count-luminosity distance* is 0.05, arising from the limited number of sky shufflings.</figcaption></figure>

only on the total number of events. For the chirp mass estimator, on the other hand, the variance for the BNSs is smaller by a factor of 3. As already discussed in Section 3.2, this is due to the narrower mass range of BNSs. For the luminosity distance estimator we see that the variance is also slightly smaller for BNSs. Again, this is due to the slightly narrower radial distribution of BNS events compared to BBH events.

Comparing the different estimators, we see that the one with smaller variance is the mass estimator for BNSs. This is due to the fact that shot noise only affects the mass estimator by changing the radial distribution and $m_{1,2}$-distribution of sources. If the chirp mass would be constant in $m_{1,2}$ and in $r$, the last two terms in the second line of equation (41) would cancel each other. Although the chirp mass does depend on $m_{1,2}$ and $r$, the distribution of BNSs in masses and redshift is narrow enough for a reduction to occur. This leads to a shot noise contribution smaller than the one for the number count estimators, where there is no such reduction. For BBHs however, the wider range leads instead to an increased variance.

Finally, we see that the variance of the distance estimator is slightly larger than that of the number count both for BBHs and BNSs, due to the large variation of distance with $r$. The luminosity distance varies indeed as $(1+z)r$, while the redshifted chirp mass varies as $(1+z)$, i.e. significantly slower.

Since the means of the three estimators are very similar (see Section 4.4), the mass estimator for BNSs is optimal in terms of SNR. On the other hand, BNSs are affected by threshold effects, meaning that the $\alpha$ coefficients need to be modelled if one wants to measure the observer velocity. Including the uncertainty in the modelling of the $\alpha$'s in the analysis generates an extra contribution to the variance, which needs to be accounted for. In Section 5, we quantify this effect and show that, despite it, the BNSs mass estimator strongly contributes to the constraints on the observer velocity.

The bottom panels of Fig. 4 show the values of the correlations (i.e. the covariance divided by the square root of the respective variances) as a function of the number of BBH and BNS detections. As we can see from the plot, also in the simulations, we obtain that the number count estimator does not correlate with the mass and distance estimators. On the other hand, as expected, we find that the distance and mass estimators are positively correlated, in excellent agreement

with the theoretical calculation. As explained before, this is due to the fact that the two observables are similarly sensitive to the radial distribution of sources.

## 4.3 Detection efficiency of the dipole

We now assess the detection efficiency of the dipole from the three estimators, for the BBH and BNS populations. For this we report in Figs 5 and 6, the detection probability versus false alarm probability (FAP) for several cases. The FAP identifies a threshold for the dipole detection and it is defined as the probability that a random fluctuation in the absence of a dipole (due to shot noise or measurement uncertainty), would result in a false positive. The detection probability is defined as the probability that, in the presence of a dipole, the estimator for the dipole detection would exceed the FAP threshold. For example, if we find a detection probability of 0.8 at an FAP of 0.32, it means that we have a probability of 80 per cent to detect the dipole at more than $1\sigma$ significance. In the other 20 per cent of times, the dipole would be consistent with zero within the $1\sigma$ error bars. However, noise alone gives one a 32 per cent per cent chance of a detection. Clearly, the lower the FAP, the higher our detection confidence. On the other hand, a lower FAP requires the dipole signal to be larger in order to be detected, which therefore decreases the detection probability. In Figs 5 and 6, we also draw with a solid black curve the noise-dominated limiting case, where the FAP equals the detection probability, i.e. where all the dipole detections are due to noise fluctuations alone.

As we can see from the plots, a dipole consistent with the CMB value becomes detectable after collecting more than $10^6$ GW events, while a dipole consistent with the AGN value becomes detectable already after $10^5$ detections. This result is valid for both the BBH and BNS populations.

Comparing the detection efficiencies of the individual estimators, we see that, in the BBH case (Fig. 5), the number count estimator (square) is clearly the most efficient one. This directly follows from the amplitude of the variance, which as shown in Fig. 4 is smaller for the number count estimator, due to the large mass and radial range of BBH events that increase the variance of the mass and distance estimators. In Fig. 5, we also show the result obtained by combining
<figure><img src="image_5.png" /><figcaption>**Figure 5.** Detection probability for the cosmic dipole (vertical axis) versus false alarm probability (horizontal axis) for a population of BBH detections. The different columns consider different numbers of GW detections. The first row corresponds to a dipole consistent with the CMB cosmic dipole and the second to a dipole consistent with the AGN dipole. The different marker types (see the legend) indicate the three estimators and their combination. The solid curve marks the detection probability/FAP relation in the case of fluctuations arising from an isotropic background. The detection probability is calculated by generating 200 population realizations while the FAP threshold is calculated using 2000 noise realizations with the sky shuffling method. The vertical solid, dashed, and dotted lines mark the $1\sigma$, $2\sigma$, and $3\sigma$ false alarm probabilities, respectively.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 6.** Same as Fig. 5, but for BNSs instead of BBHs.</figcaption></figure>

the three estimators, accounting for the covariance between the mass and distance estimators. We see that in the case of BBHs, the combined estimator (stars) does not perform better than the number count estimator, meaning that the mass and distance estimators are irrelevant. Note that the fact that sometimes the combined estimator is below the number count estimator is due to numerical uncertainties, related to the fact that the variance and covariance are generated using only 2'000 sky shufflings.

In Fig. 6, we show the results for BNSs. We clearly see that in this case the mass estimator has the highest detection efficiency, better than the number count estimator. This is directly related to the shot noise suppression discussed in Section 4.2. The distance estimator performs better for the BNSs than for the BBHs, due to the smaller redshift range spanned by BNSs, which also reduces the shot noise contribution with respect to the BBH case. Despite this, the efficiency of the distance estimator remains below that of the mass and number

count. Combining the three estimators leads to a non-negligible gain in terms of detection efficiency.<fnref n="6" />

Finally, we also simulated a more realistic scenario for the detection of the cosmic dipole, where we do not perform separate analyses for BNSs and BBHs, but instead combine all the BBHs and BNSs detected in a given observing time $T_{\text{obs}}$. Note that in this combination, we still apply the estimators separately on the BNS and BBH populations, since we want to preserve the fact that BNSs have a smaller mass range and radial range than BBHs. We consider that in 1 yr, the network of ET + 2CE detectors would be able to observe

<footnote n="6">⁶Note that if we want to combine the estimators to extract the observer velocity $v_0/c$, we need to divide each estimator by its appropriate $\alpha$, as discussed in Section 3.5. On the other hand, if we just want to assess the detectability of the combined dipole, we do not need this extra step. In this case, the mean of the combined estimator does however not provide a measurement of $v_0/c$.</footnote>
<figure><img src="image_7.png" /><figcaption>**Figure 7.** Detection probability for the cosmic dipole (vertical axis) versus false alarm probability (horizontal axis) using the combined estimator for a mixed population of GW sources. The different columns consider different observation times. Following Iacovelli et al. (2022), for 1 yr of observations we have taken 7.5 × 10⁴ BBH detections and 10⁵ BNS detections. The number of events expected in 6 months and in 3 and 5 yr of observations are found by linearly scaling the previous fiducial values. The first row corresponds to a dipole consistent with the CMB cosmic dipole and the second to a dipole consistent with the AGN dipole. The solid curve marks the detection probability/FAP relation in the case of fluctuations arising from an isotropic background. The detection probability is calculated by generating 200 population realizations while the FAP threshold is calculated using 2000 noise realizations obtained with the sky shuffling method. The vertical solid, dashed, and dotted lines mark the 1σ, 2σ, and 3σ false alarm probabilities, respectively.</figcaption></figure>

7.5 × 10⁴ BBHs and 10⁵ BNSs. In Fig. 7, we show the detection efficiency for the combined estimator for 6 months, 1 yr, 3 yr, and 5 yr of observations. On one hand, we find that, with 5 yr of observations and a bit less than 10⁶ sources detected, we can detect the dipole with 1σ significance, but it will be unlikely to reach a high significance of 3σ. On the other hand, we find that a dipole consistent with the AGN one could be detected with high significance already with 1 yr of observations, thanks to the BNS population. Therefore, a non-detection of the cosmic dipole in the first year of XG detectors would automatically rule out the AGN value of the dipole, thus providing a strong indication for the presence of unmodelled systematic in the AGN measurements.

## 4.4 Modelling of the α's

The results of the previous section show that the mass estimator for the BNSs is better than the count estimator, both for BNSs and BBHs, in terms of detectability of the dipole signal. However, BNSs (contrary to BBHs) are affected by threshold effects. In order to use the BNSs mass estimator to measure the observer velocity, it is necessary to have a modelling of the coefficients α.

The parameters α depend on the population of sources, through the parameters s and A defined in equations (9) and (10). We use the mass model defined in appendix A of Mastrogiovanni et al. (2023) to describe the population of BNSs. We then bin the events in (r, m₁,₂, ρ*) around ρ* = 9, and compute with finite differences the derivative with respect to ρ* to obtain s through equation (9). This needs to be done for each (r, m₁,₂) bin, at ρ* = 9. We then use interpolation in r and m₁,₂ to promote the binned s values to a function. The function A depends on F, which quantifies the

<table><caption>**Table 1.** Expected values for the BNS α parameters, for a fiducial astrophysical population model.</caption><thead><tr><th></th><th>α<sub>N</sub></th><th>α<sub>dL</sub></th><th>α<sub>M</sub></th></tr></thead><tbody><tr><td>BNS</td><td>1.08</td><td>0.94</td><td>0.98</td></tr></tbody></table>

sensitivity of the detector. We first evaluate numerically F on a discrete set of points and then interpolate between them, allowing to attribute one A value to each event. We can then estimate the integrals (14), (22) and (25) to obtain α<sub>dL</sub>, α<sub>M</sub> and α<sub>N</sub>.

The results are shown in Table 1. We see that the values are very close to one in all three cases. Threshold effects seem to be slightly suppressed in the mass and distance case with respect to the number count, probably due to the fact that the integrals (14) and (22) vanish if s and A are constant in r.

Since we know the observer velocity in our synthetic catalogues of BNS events, we can use the dipole to measure the parameters α<sub>dL</sub>, α<sub>M</sub>, and α<sub>N</sub> and compare this with our modelling. In reality, we would not be able to do that and the only quantity that can be measured is the product of v₀/c and the respective α. In Figs 8 and 9, we show the values obtained for the α's assuming an observer velocity consistent with the AGN dipole and CMB dipole, respectively. For BBHs, the histograms are centred around 1, as expected. For BNSs, the histograms are slightly displaced due to threshold effects. The peak is in excellent agreement with the theoretical predictions for the α's. This is important for two reasons: first it shows that threshold effects are indeed small and should not spoil too much our measurement of the observer velocity. In particular, one of the goals of using GWs to measure the dipole is to determine if GWs are consistent with the AGN dipole or with the CMB dipole. Having
<figure><img src="image_8.png" /><figcaption>**Figure 8.** The plots indicate the distribution of the $\alpha$ parameters obtained for the three estimators with 200 population realizations. The vertical black lines indicate the theoretical prediction, see Table 1. The histogram width is due to shot noise and distance and mass uncertainty, which vary over the 200 realizations. The first row of plots considers $10^6$ detections and the second row $10^7$ detections. The first column is for BBHs, while the second for BNSs. The plots are generated with the AGN value of the observer velocity.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 9.** Same as Fig. 8, but with the CMB value of the observer velocity.</figcaption></figure>

threshold effects of the order of 10 per cent means that this test can be done robustly also using BNSs. Indeed, we can conclude that if we find a dipole consistent with the AGN one, it would be unlikely that this was due to very large threshold effects increasing the $\alpha$'s by a factor 5. Second, the plots show excellent agreement between the theoretical modelling and the measured $\alpha$'s. Since in practice the $\alpha$'s can only be modelled (and not measured), it is important to know that this can be done in a robust way.

In Fig. 10, we check the dependence of the $\alpha$'s on the population model. We vary the value of one of the parameter in the population model by 10 per cent around its fiducial value and compute the histograms for the $\alpha$'s. Comparing the width in the $\alpha$'s with the one only due to shot noise and measurement uncertainty (Fig. 8), we see that varying the model does not generate an additional spread in the $\alpha$'s. This means that the uncertainty from the choice of model is smaller than the variance of the dipole.

## 5 DETERMINATION OF THE AMPLITUDE OF THE OBSERVER VELOCITY

Let us now use the Fisher formalism to estimate how well we can measure the amplitude of the velocity, $v_0/c$, by combining the different estimators. We consider two cases, the first one where the only unknown is the observer velocity, i.e. we assume that the parameters $\alpha_N$, $\alpha_{d_L}$ and $\alpha_M$ are perfectly known; and the second

<figure><img src="image_10.png" /><figcaption>**Figure 10.** The plots indicate the distribution of the $\alpha$ parameters obtained for the three estimators with 200 population realizations. Each realization also considers a random realization of the merger rate model parameters $\gamma$ with an uncertainty of 10 per cent around its fiducial value of $\gamma = 2.7$. The vertical black lines indicate the theoretical prediction from Table 1 for the $\gamma$ fiducial value. The histogram width is due to shot noise, distance and mass uncertainty, and also variation of the merger rate parameter. The top plot is for $10^6$ detections while the second one for $10^7$ detections. The plots are generated with the AGN value of the observer velocity.</figcaption></figure>

one where these parameters are treated as free parameters, that can be determined (through additional measurements or through a theoretical modelling) with some uncertainty. In both cases, we assume that the direction of the dipole has already been determined by maximizing the estimators $\hat{v}_N$ with respect to $\mathbf{n}$.

### 5.1 Known $\alpha$'s

We use the Fisher formalism to estimate the uncertainty on $v_0/c$ obtained from the measurement of the dipole amplitude from the number count, chirp mass and luminosity distance of the BBH and the BNS populations. The signal contains therefore six measurements:

$$S = (\hat{v}_N^{\text{BBH}}, \hat{v}_{d_L}^{\text{BBH}}, \hat{v}_M^{\text{BBH}}, \hat{v}_N^{\text{BNS}}, \hat{v}_{d_L}^{\text{BNS}}, \hat{v}_M^{\text{BNS}}). \quad (53)$$

The Fisher for the (unique) parameter $v_0/c$ is then given by

$$\mathcal{F}_{v_0/c} = \sum_{ij} \frac{\partial S_i}{\partial(v_0/c)} \text{cov}(S_i, S_j)^{-1} \frac{\partial S_j}{\partial(v_0/c)} \quad (54)$$

$$= \sum_{X=\text{BBH}, \text{BNS}} \begin{pmatrix} \alpha_N^X \\ \alpha_{d_L}^X \\ \alpha_M^X \end{pmatrix} \text{cov}_X^{-1} \begin{pmatrix} \alpha_N^X \\ \alpha_{d_L}^X \\ \alpha_M^X \end{pmatrix}, \quad (55)$$

where

$$\text{cov}_X = \begin{pmatrix} \text{var}(\hat{v}_N^X) & 0 & 0 \\ 0 & \text{var}(\hat{v}_{d_L}^X) & \text{cov}(\hat{v}_{d_L}^X, \hat{v}_M^X) \\ 0 & \text{cov}(\hat{v}_{d_L}^X, \hat{v}_M^X) & \text{var}(\hat{v}_M^X) \end{pmatrix}, \quad (56)$$

for $X = \text{BBH, BNS}$. Here, we have used that the BBH and BNS measurements are uncorrelated (meaning that the Fisher matrix can be written as a sum over the two populations) and that within one population, the number counts are uncorrelated with the distance and the mass, as shown in Section 3.4.2. The error on $v_0/c$ is then given by

$$\sigma_{v_0/c} = \sqrt{\mathcal{F}_{v_0/c}^{-1}}. \quad (57)$$

To compute $\sigma_{v_0/c}$, we need to know the values of the coefficients $\alpha_N^X$ for the different populations and the different estimators. For the BBH, since threshold effects are negligible we have $\alpha_N^{\text{BBH}} = \alpha_{d_L}^{\text{BBH}} =$
<table><caption>**Table 2.** Error, $\sigma_{v_0/c}$, obtained from: the individual estimators (first six columns), combining the three BBH estimators, the three BNS estimators, all estimators, or using the top three combinations of mass and number count estimator for BNSs and number count estimator for BBHs. In the top three rows, we consider a dipole consistent with the CMB one, while in the bottom three rows, we consider a dipole consistent with the AGN one.</caption><thead><tr><th></th><th></th><th>v̂<sup>BBH</sup><sub>N</sub><br>(per cent)</th><th>v̂<sup>BBH</sup><sub>d<sub>L</sub></sub><br>(per cent)</th><th>v̂<sup>BBH</sup><sub>M</sub><br>(per cent)</th><th>v̂<sup>BNS</sup><sub>N</sub><br>(per cent)</th><th>v̂<sup>BNS</sup><sub>d<sub>L</sub></sub><br>(per cent)</th><th>v̂<sup>BNS</sup><sub>M</sub><br>(per cent)</th><th>BBH<br>(per cent)</th><th>BNS<br>(per cent)</th><th>All<br>(per cent)</th><th>Top 3<br>(per cent)</th></tr></thead><tbody><tr><td rowspan="3">CMB</td><td>N<sub>tot</sub> = 10<sup>5</sup></td><td>223</td><td>327</td><td>356</td><td>206</td><td>305</td><td>194</td><td>178</td><td>139</td><td>110</td><td>119</td></tr><tr><td>N<sub>tot</sub> = 10<sup>6</sup></td><td>70</td><td>103</td><td>113</td><td>65</td><td>96</td><td>61</td><td>56</td><td>44</td><td>35</td><td>38</td></tr><tr><td>N<sub>tot</sub> = 10<sup>7</sup></td><td>22</td><td>33</td><td>36</td><td>21</td><td>30</td><td>19</td><td>18</td><td>14</td><td>11</td><td>12</td></tr><tr><td rowspan="3">AGN</td><td>N<sub>tot</sub> = 10<sup>5</sup></td><td>46</td><td>67</td><td>73</td><td>42</td><td>63</td><td>40</td><td>36</td><td>29</td><td>23</td><td>24</td></tr><tr><td>N<sub>tot</sub> = 10<sup>6</sup></td><td>14</td><td>21</td><td>23</td><td>13</td><td>20</td><td>13</td><td>12</td><td>9</td><td>7</td><td>8</td></tr><tr><td>N<sub>tot</sub> = 10<sup>7</sup></td><td>5</td><td>7</td><td>7</td><td>4</td><td>6</td><td>4</td><td>4</td><td>3</td><td>2</td><td>2</td></tr></tbody></table>

$\alpha_{\Lambda}^{\text{BBH}} = 1$. For the BNS, we use the values calculated theoretically in Section 4.4 and reported in Table 1.

The error, $\sigma_{v_0/c}$, is reported in Table 2 for different cases. First, we compute the error from each estimator taken individually. We see that, as expected, the chirp mass estimator for BNSs gives the smallest uncertainty. The number count estimators for BNSs and for BBHs are the other two best ones. The number count for BNSs is very slightly better than for BBHs, due to threshold effects, that increase $\alpha_N$ from 1 to 1.08. Combining the three BBH estimators improves the constraints by 20 per cent compared to using the number counts alone, while combining the three BNS estimators improves the constraints by 32 per cent. Combining all estimators improves the constraints by 50 per cent compared the number counts of BBHs alone (studied in Mastroiovanni et al. (2023)). We also show the results for the top 3 estimators, i.e. mass estimator for BNSs and number count estimators for BBHs and BNSs. We see that the constraints are very similar to those obtained with all estimators.

The absolute error on the observer velocity is the same for the CMB and AGN case, however the relative error is reduced by a factor 5 for the AGN case, as can be seen from Table 2. Hence we see that a robust measurement of the observer velocity requires $10^6$ events if the observer velocity is consistent with the CMB dipole, but only $10^5$ events if the observer velocity is consistent with the AGN one, which is consistent with the detection efficiency results in Figs 5 and 6. This means in particular that if we do not detect a dipole with $10^5$ events, the GW dipole is in tension with the AGN dipole.

## 5.2 Adding uncertainties on $\alpha$'s

Whereas for the BBHs the $\alpha$ coefficients are known and equal to 1, for BNSs it is not the case. These coefficients are affected by threshold effects. As shown in Section 4.4, these coefficients can be computed, assuming a given model for the population of BNSs. The uncertainty on the model will directly impact the determination of the $\alpha$'s. To account for these uncertainties in the Fisher computation, we add the coefficients in the signal, and assign a corresponding error in the covariance matrix. We then compute how this uncertainty degrades the constraints on $v_0/c$. We consider the combination of all estimators and compute the error on $v_0/c$ for different values of the uncertainties on the $\alpha$'s. The signal is written as

$$S = (v_{N}^{\text{BBH}}, v_{d_{L}}^{\text{BBH}}, v_{\mathcal{M}}^{\text{BBH}}, v_{N}^{\text{BNS}}, v_{d_{L}}^{\text{BNS}}, v_{\mathcal{M}}^{\text{BNS}}, \alpha_{N}^{\text{BNS}}, \alpha_{d_{L}}^{\text{BNS}}, \alpha_{\mathcal{M}}^{\text{BNS}}) . (58)$$

The Fisher matrix contains now four parameters: $\theta_1 = v_0/c$, $\theta_2 = \alpha_N^{\text{BNS}}$, $\theta_3 = \alpha_{d_L}^{\text{BNS}}$ and $\theta_4 = \alpha_{\mathcal{M}}^{\text{BNS}}$ and it is given by

$$F_{ab} = \sum_{ij} \frac{\partial S_i}{\partial \theta_a} \text{cov}(S_i, S_j)^{-1} \frac{\partial S_j}{\partial \theta_b} . \quad (59)$$

<table><caption>**Table 3.** Fisher bound on the error $\sigma_{v_0/c}$, obtained from the combination of the six estimators assuming different uncertainties on the $\alpha$'s for the BNSs. In the top three rows, we consider a dipole consistent with the CMB one, while in the bottom three rows we consider a dipole consistent with the AGN one.</caption><thead><tr><th></th><th></th><th>10 per cent</th><th>20 per cent</th><th>50 per cent</th></tr></thead><tbody><tr><td rowspan="3">CMB</td><td>N<sub>tot</sub> = 10<sup>5</sup></td><td>110 per cent</td><td>110 per cent</td><td>112 per cent</td></tr><tr><td>N<sub>tot</sub> = 10<sup>6</sup></td><td>35 per cent</td><td>36 per cent</td><td>40 per cent</td></tr><tr><td>N<sub>tot</sub> = 10<sup>7</sup></td><td>12 per cent</td><td>13 per cent</td><td>16 per cent</td></tr><tr><td rowspan="3">AGN</td><td>N<sub>tot</sub> = 10<sup>5</sup></td><td>23 per cent</td><td>24 per cent</td><td>28 per cent</td></tr><tr><td>N<sub>tot</sub> = 10<sup>6</sup></td><td>8 per cent</td><td>9 per cent</td><td>11 per cent</td></tr><tr><td>N<sub>tot</sub> = 10<sup>7</sup></td><td>3 per cent</td><td>3 per cent</td><td>4 per cent</td></tr></tbody></table>

The covariance matrix has the form

$$\text{cov} = \begin{pmatrix} \text{cov}_{\text{BBH}} & 0 & 0 & 0 & 0 \\ 0 & \text{cov}_{\text{BNS}} & 0 & 0 & 0 \\ 0 & 0 & (x \cdot \alpha_N^{\text{BNS}})^2 & 0 & 0 \\ 0 & 0 & 0 & (x \cdot \alpha_{d_L}^{\text{BNS}})^2 & 0 \\ 0 & 0 & 0 & 0 & (x \cdot \alpha_{\mathcal{M}}^{\text{BNS}})^2 \end{pmatrix} \quad (60)$$

where the 3 × 3 blocks $\text{cov}_{\text{BBH}}$ and $\text{cov}_{\text{BNS}}$ are given by equation (56). In equation (60), $x$ denotes the relative uncertainty on the determination of the $\alpha$ parameters for the BNS population, that we assume here to be the same. We consider three cases for $x$: 10 per cent, 20 per cent and 50 per cent. The error on $v_0/c$ is then given by

$$\sigma_{v_0/c} = \sqrt{(\mathcal{F}^{-1})_{v_0/c} v_0/c} . \quad (61)$$

The results are reported in Table 3. Comparing with the combination of all estimators (second last column) in Table 2, we see that having a 10 per cent uncertainty on all three $\alpha^{\text{BNS}}$ degrades the constraints on $v_0/c$ by at most 1 per cent compared to the case where these parameters are assumed to be known perfectly. Increasing this uncertainty to 50 per cent degrades the constraints by at most 5 per cent. Hence, even if our modelling of the threshold effects is not very accurate, the degradation remains small.

Comparing Table 3 with the BBH column of Table 2, we see that even in the case where the uncertainty on the $\alpha$ is as large as 50 per cent, we still gain information by including BNSs. For example, for $10^6$ events and a dipole consistent with the CMB one, we gain 29 per cent in the measurement of $v_0/c$ by adding BNSs. If we can model the $\alpha$'s with a precision of 20 per cent, this gain increases to 36 per cent (38 per cent for a 10 per cent precision). Seen as a function of the $\alpha^{\text{BNS}}$ uncertainty, the constraints from all estimators are upper-bounded by the constraints from the combination of the three BBH estimators. At worst, if the uncertainty on the $\alpha^{\text{BNS}}$ is too large, no information is gained by adding BNSs.
# 6 CONCLUSIONS

In this paper we have developed a robust framework to measure the cosmic dipole using GW detections. Contrary to radio sources or quasars, for which only the sky position can be used, GWs have the advantage of providing three quantities that are affected by the observer velocity: sky position, luminosity distance and redshifted chirp mass. We have developed estimators of these three dipoles, and we have calculated their variance and covariance. We have found that the mass and distance estimators are partially correlated, but that they are both uncorrelated with the number count estimator. Combining the three of them does therefore increase the detectability of the dipole.

BBHs have the advantage over BNSs to be unaffected by threshold effects, since all sources within the frequency range of ET and CE will have SNR above threshold. On the other hand, a significant fraction of BNSs will have an SNR below threshold, meaning that threshold effects are relevant in this case. The dipole from BNSs can of course be detected even without knowing the amplitude of threshold effects. However, to interpret the results and determine if the amplitude is consistent or not with the CMB dipole, it is necessary to have a modelling of these effects. We have developed such a modelling and computed the amplitude of threshold effects. For our population model, we have found that these effects are small, of the order of 10 per cent at most for all three estimators. The amplitudes of these effects depend of course on the population model that is used, however we expect that the order of magnitude we estimated will not change when changing the details of the population model. This shows that it is worth including BNSs to measure the observer velocity and test the isotropy of the Universe.

Comparing the three BNS and three BBH estimators, we have found that the BNS chirp mass estimator is the one with higher detectability, i.e. lower variance. This is due to the fact that the variance is fully dominated by shot noise, which generates fluctuations in the radial distribution of sources, consequently changing the mean mass per pixel. Since the intrinsic mass distribution of BNSs is very narrow, this shot noise contribution is mainly due to the redshift dependence in the chirp mass, which is significantly smaller than the spread in luminosity distance. After the BNS chirp mass, the other two best estimators are the number counts of BBHs and BNSs.

Overall, we have found that combining all events, we need a few $10^6$ events to detect a dipole consistent with the CMB one. On the other hand, if the dipole is consistent with the AGN one, we should detect it with $10^5$ events. This can be achieved already after 1 yr of observation. In this context, the fact that threshold effects are small is crucial, since it ensures that they cannot boost the dipole by a factor 5, thus mimicking the amplitude of the AGN dipole (which is 5 times larger than the CMB one). Hence, if we see results consistent with the AGN dipole, we can robustly conclude that it is not due to threshold effects, but rather to a large intrinsic anisotropy of the large scale structure.

# ACKNOWLEDGEMENTS

NG and CB acknowledge funding from the European Research Council (ERC) under the European Union's Horizon 2020 - Research and Innovation Framework Programme (Grant agreement No. 863929; project title 'Testing the law of gravity with novel large-scale structure observables'). CB is also supported by the Swiss National Science Foundation. SF is supported by the Fonds National Suisse, grant 200020.191957, and by the National Centres of Competence in Research SwissMAP. The work of GC is supported by CNRS and GC and MP acknowledge support from the Swiss National Science Foundation (Ambizione grant, 'Gravitational wave propagation in

MNRAS 526, 4673–4689 (2023)

the clustered universe'). SM acknowledges the support of the computing facilities at INFN Rome of the Amaldi Research Center funded by the MIUR program 'Dipartimento di Eccellenza' (CUP: B81118001170001).

# DATA AVAILABILITY

The simulations and numerical code underlying this paper are available on GitHub https://github.com/simone-mastrogiovanni/cosmic_dipole_GW_3G.

# REFERENCES

* Aghanim N. et al., 2020, A&A, 641, A1

* Bengaly C. A., Maartens R., Santos M. G., 2018, J. Cosmo. Astropart. Phys., 2018, 031

* Bonvin C., Durrer R., Gasparini M. A., 2006a, Phys. Rev. D, 73, 023523

* Bonvin C., Durrer R., Kunz M., 2006b, Phys. Rev. Lett., 96, 191302

* Cai R.-G., Liu T.-B., Liu X.-W., Wang S.-J., Yang T., 2018, Phys. Rev. D, 97, 103005

* Callister T., Fishbach M., Holz D. E., Farr W. M., 2020, ApJ, 896, L32

* Chernoff D. F., Finn L. S., 1993, ApJ, 411, L5

* Colin J., Mohayaee R., Rameez M., Sarkar S., 2017, MNRAS, 471, 1045

* Dalang C., Bonvin C., 2022, MNRAS, 512, 3895

* Dwyer S., Sigg D., Ballmer S. W., Barsotti L., Mavalvala N., Evans M., 2015, Phys. Rev. D, 91, 082001

* Ellis G. F. R., Baldwin J. E., 1984, MNRAS, 206, 377

* Essick R., Farr W. M., Fishbach M., Holz D. E., Katsavounidis E., 2023, Phys. Rev. D, 107, 043016

* Finke A., Foffa S., Iacovelli F., Maggiore M., Mancarella M., 2022, Phys. Dark Univ., 36, 100994

* Finn L. S., Chernoff D. F., 1993, Phys. Rev. D, 47, 2198

* Guandalin C., Piat J., Clarkson C., Maartens R., 2023, AJ, 953, 144

* Horstmann N., Pietschke Y., Schwarz D. J., 2022, A&A, 668, A34

* Iacovelli F., Mancarella M., Foffa S., Maggiore M., 2022, ApJ, 941, 208

* Kashyap G., Singh N. K., Phukon K. S., Caudill S., Jain P., 2023, J. Cosmol. Astropart. Phys., 2023, 042

* Maggiore M., 2007, Gravitational Waves. Vol. 1: Theory and Experiments. Oxford Master Series in Physics, Oxford University Press

* Mastrogiovanni S., Bonvin C., Cusin G., Foffa S., 2023, MNRAS, 521, 984

* Payne E., Banagiri S., Lasky P., Thrane E., 2020, Phys. Rev. D, 102, 102004

* Punturo M. et al., 2010, Class. Quantum Gravity, 27, 194002

* Reitze D. et al., 2019, BAAS, 51, 035

* Riess A. G., Press W. H., Kirshner R. P., 1995, ApJ, 438, L17

* Secrest N. J., Hausegger S. v., Rameez M., Mohayaee R., Sarkar S., Colin J., 2021, ApJ, 908, L51

* Secrest N. J., von Hausegger S., Rameez M., Mohayaee R., Sarkar S., 2022, ApJ, 937, L31

* Siewert T. M., Schmidt-Rubart M., Schwarz D. J., 2021, A&A, 653, A9

* Singal A. K., 2022, MNRAS, 515, 5969

* Sorrenti F., Durrer R., Kunz M., 2022 preprint (arXiv:2212.10328)

* Stiskalek R., Veitch J., Messenger C., 2021, MNRAS, 501, 970

* Taylor S. R., Gair J. R., 2012, Phys. Rev. D, 86, 023502

* Taylor S. R., Gair J. R., Mandel I., 2012, Phys. Rev. D, 85, 023535

# APPENDIX A: COMPUTATION OF THE MEAN AND VARIANCE

## A1 Mean of the luminosity distance estimator

The mean of the luminosity distance estimator (26) is given by

$$ \langle \hat{v}_{d_L - \mathbf{n}'} \rangle = \frac{-3}{N_{\text{sky}}} \sum_{i=1}^{N_{\text{sky}}} \sum_{j=1}^{N_{\text{det}}} \left\langle \frac{(d_L)_j^i}{\hat{d}_L N_i^{\text{det}}} \right\rangle (\mathbf{n}_i \cdot \mathbf{n}') . \quad (A1) $$
The number of detections in angular pixel *i* is given by

$$N_i^{\text{det}} = \int dr \frac{dN_{\text{det}}}{d\Omega dr} (\mathbf{r}, \mathbf{n}_i) \Delta\Omega + \Delta N^i, \qquad (\text{A2})$$

where $\Delta N^i$ is the shot noise contribution. The luminosity distance can be written as

$$(d_L)_j^i = d_L(r_j, \mathbf{n}_i) + (\Delta d_L)_j^i, \qquad (\text{A3})$$

where $(\Delta d_L)_j^i$ denotes the error in the measurement of the distance. Assuming that the shot noise is small compared to the number of events in the angular pixel, and that the error in the distance measurement is small compared to the distance, we can expand each term in equation (A1) keeping only linear contributions in $\Delta N^i$ and $(\Delta d_L)_j^i$. Since the mean of these quantities is zero: $\langle \Delta N^i \rangle = \langle (\Delta d_L)_j^i \rangle = 0$, equation (A1) becomes

$$\langle \hat{v}_{d_L - \mathbf{n}} \rangle = \frac{-3 \sum_{i=1}^{N_{\text{sky}}} (\Delta \Omega \int dr \frac{dN_{\text{det}}}{d\Omega dr}(\mathbf{r}, \mathbf{n}_i))^{-1} \sum_{j=1}^{N_{\text{det}}^i} d_L(r_j, \mathbf{n}_i)(\mathbf{n}_i \cdot \mathbf{n}')}{\sum_{m=1}^{N_{\text{sky}}} (\Delta \Omega \int dr \frac{dN_{\text{det}}}{d\Omega dr}(\mathbf{r}, \mathbf{n}_m))^{-1} \sum_{n=1}^{N_{\text{det}}^m} d_L(r_n, \mathbf{n}_m)}. \qquad (\text{A4})$$

To compute this, we divide each solid angle $i$ into $N_r$ bins in $r$ of size $\Delta r$, and we rewrite the sum over $j$ as

$$\sum_{j=1}^{N_{\text{det}}^i} d_L(r_j, \mathbf{n}_i) = \sum_{k=1}^{N_r} \frac{d N_{\text{det}}}{d \Omega dr} (r_k, \mathbf{n}_i) \Delta r \Delta \Omega d_L(r_k, \mathbf{n}_i), \qquad (\text{A5})$$

where $r_k$ denotes the centre of the bin number $k$ in $r$. In the continuous limit, i.e. taking $N_r \to \infty$ corresponding to $\Delta r \to 0$, this becomes

$$\Delta \Omega \int dr \frac{dN_{\text{det}}}{d\Omega dr}(\mathbf{r}, \mathbf{n}_i) d_L(r, \mathbf{n}_i). \qquad (\text{A6})$$

Similarly, in the continuous limit the sum over $N_{\text{sky}} \to \infty$ becomes

$$\sum_{i}^{N_{\text{sky}}} F(r, \mathbf{n}_i) = \frac{1}{\Delta \Omega} \int d\Omega F(r, \mathbf{n}). \qquad (\text{A7})$$

With this, the numerator in equation (A4) becomes

$$\begin{aligned} & \frac{-3}{\Delta \Omega} \int d\Omega \frac{\int dr \frac{dN_{\text{det}}}{d\Omega dr}(\mathbf{r}, \mathbf{n}) d_L(r, \mathbf{n})}{\int dr \frac{dN_{\text{det}}}{d\Omega dr}(\mathbf{r}, \mathbf{n})} (\mathbf{n} \cdot \mathbf{n}') \\ &= \frac{-3}{\Delta \Omega} \int d\Omega d_L^{\text{mean}}(\mathbf{n})(\mathbf{n} \cdot \mathbf{n}') \\ &= N_{\text{sky}} d_L^{(0)} \alpha_{d_L} \frac{v_0}{c} \cos \theta', \end{aligned} \qquad (\text{A8})$$

where we have used that $\Delta\Omega = 4\pi/N_{\text{sky}}$. Similarly, the denominator in equation (A4) becomes

$$\begin{aligned} & \frac{1}{\Delta \Omega} \int d\Omega \frac{\int dr \frac{dN_{\text{det}}}{d\Omega dr}(\mathbf{r}, \mathbf{n}) d_L(r, \mathbf{n})}{\int dr \frac{dN_{\text{det}}}{d\Omega dr}(\mathbf{r}, \mathbf{n})} \\ &= \frac{1}{\Delta \Omega} \int d\Omega d_L^{\text{mean}}(\mathbf{n}) \\ &= N_{\text{sky}} d_L^{(0)}. \end{aligned} \qquad (\text{A9})$$

Therefore, we obtain for the mean

$$\langle \hat{v}_{d_L - \mathbf{n}'} \rangle = \alpha_{d_L} \frac{v_0}{c} \cos \theta' = v_{d_L - \mathbf{n}'} . \qquad (\text{A10})$$

## A2 Variance of the luminosity distance estimator

To compute the variance of the luminosity distance, we use equation (30). We start by computing the variance of the variable $X$ defined

in (29)

$$\begin{aligned} \text{var}(X) &= \langle X^2 \rangle - \langle X \rangle^2 \\ &= \frac{9}{N_{\text{sky}}^2} \sum_{i=1}^{N_{\text{sky}}} \sum_{m=1}^{N_{\text{sky}}} \sum_{j=1}^{N_{\text{det}}^i} \sum_{n=1}^{N_{\text{det}}^m} (\mathbf{n}_i \cdot \mathbf{n}') (\mathbf{n}_m \cdot \mathbf{n}') \\ & \times \left( \left\langle \frac{\langle d_L \rangle_j^i \langle d_L \rangle_n^m}{N_{\text{det}}^i N_{\text{det}}^m} \right\rangle - \left\langle \frac{\langle d_L \rangle_j^i}{N_{\text{det}}^i} \right\rangle \left\langle \frac{\langle d_L \rangle_n^m}{N_{\text{det}}^m} \right\rangle \right). \end{aligned} \qquad (\text{A11})$$

As done for the mean, we divide each solid angle $i$ into $N_r$ bins in $r$ of size $\Delta r$, and we rewrite the sum over $j$ as a sum over the $r$-bins. We obtain

$$\begin{aligned} \text{var}(X) = & \frac{9}{N_{\text{sky}}^2} \sum_{i=1}^{N_{\text{sky}}} \sum_{m=1}^{N_{\text{sky}}} \sum_{k=1}^{N_r} \sum_{\ell=1}^{N_l} (\mathbf{n}_i \cdot \mathbf{n}') (\mathbf{n}_m \cdot \mathbf{n}') \\ & \times \left( \left\langle \frac{N_p^{ik}(d_L)_k^i N_p^{m\ell}(d_L)_\ell^m}{N_{\text{det}}^i N_{\text{det}}^m} \right\rangle - \left\langle \frac{N_p^{ik}(d_L)_k^i}{N_{\text{det}}^i} \right\rangle \left\langle \frac{N_p^{m\ell}(d_L)_\ell^m}{N_{\text{det}}^m} \right\rangle \right). \end{aligned} \qquad (\text{A12})$$

where $N_p^{ik}$ denotes the number of events in the pixel in direction $\mathbf{n}_i$ and centred at distance $r_k$ and $(d_L)_k^i = d_L(\mathbf{n}_i, r_k)$. Neglecting the dipole contribution, which does not contribute to the variance, we have

$$N_p^{ik} = \bar{N}_p^k + \Delta N_p^{ik} \qquad (\text{A13})$$

$$N_{\text{det}}^i = \bar{N}_{\text{det}} + \sum_j \Delta N_p^{ij} \qquad (\text{A14})$$

$$(d_L)_k^i = (\bar{d}_L)_k + (\Delta d_L)_k^i \qquad (\text{A15})$$

Here $(\Delta d_L)_k^i$ is the error for one measurement of $d_L$ in angular bin $i$ and radial bin $k$. Inserting this into equation (A12), keeping only the terms quadratic in $\Delta N_p^{ik}$ and $(\Delta d_L)_k^i$, and using that the two types of error are uncorrelated we obtain

$$\begin{aligned} \text{var}(X) = & \frac{9}{N_{\text{sky}}^2} \sum_{i=1}^{N_{\text{sky}}} \sum_{m=1}^{N_{\text{sky}}} \sum_{k=1}^{N_r} \sum_{\ell=1}^{N_l} (\mathbf{n}_i \cdot \mathbf{n}') (\mathbf{n}_m \cdot \mathbf{n}') \frac{\bar{N}_p^k \bar{N}_p^{\ell}(\bar{d}_L)_k(\bar{d}_L)_\ell}{\bar{N}_{\text{det}}^2} \\ & \times \left( \left\langle \frac{\Delta N_p^{ik}}{N_p^k} \frac{\Delta N_p^{m\ell}}{N_p^\ell} \right\rangle - \sum_n \left\langle \frac{\Delta N_p^{ik}}{N_p^k} \frac{\Delta N_p^{mn}}{N_{\text{det}}^m} \right\rangle - \sum_j \left\langle \frac{\Delta N_p^{m\ell}}{N_p^\ell} \frac{\Delta N_p^{ij}}{N_{\text{det}}^j} \right\rangle \right. \\ & \left. + \sum_{j,n=1}^{N_r} \left\langle \frac{\Delta N_p^{ij}}{N_{\text{det}}^k} \frac{\Delta N_p^{mn}}{N_{\text{det}}^m} \right\rangle + \left\langle \frac{(\Delta d_L)_k^i (\Delta d_L)_\ell^m}{(\bar{d}_L)_k(\bar{d}_L)_\ell} \right\rangle \right). \end{aligned} \qquad (\text{A16})$$

We then use that

$$\langle \Delta N_p^{ik} \Delta N_p^{m\ell} \rangle = \bar{N}_p^k \delta_{im} \delta_{k\ell}. \qquad (\text{A17})$$

With this, the first 4 terms in the parenthesis of equation (A16) become

$$\frac{1}{\bar{N}_p^k} \delta_{im} \delta_{k\ell} - \frac{1}{\bar{N}_{\text{det}}} \delta_{im}. \qquad (\text{A18})$$

To compute the error from the luminosity distance we need to rewrite the sum over $r$-bins as a sum over events:

$$\sum_{k=1}^{N_r} \sum_{\ell=1}^{N_r} \bar{N}_p^k \bar{N}_p^{\ell} \langle (\Delta d_L)_k^i (\Delta d_L)_\ell^m \rangle = \sum_{j=1}^{N_{\text{det}}^i} \sum_{n=1}^{N_{\text{det}}^m} \langle (\Delta d_L)_j^i (\Delta d_L)_n^m \rangle, \qquad (\text{A19})$$

and we use that the error for different events are uncorrelated

$$\langle (\Delta d_L)_j^i (\Delta d_L)_n^m \rangle = \langle (\Delta d_L(r_j))^2 \rangle \delta_{im} \delta_{jn}. \qquad (\text{A20})$$
<figure><img src="image_11.png" /><figcaption>**Figure A1.** Variance of the luminosity distance estimator for BNSs, plotted as a function of the total number of events $N_{tot}$. The different panels are for different values of the relative error on $d_L$ (assumed to be independent of redshift). We show separately the shot noise contribution and the distance-uncertainty contribution, as well as the total.</figcaption></figure>

Inserting this in equation (A19) we obtain

$$ 
\begin{aligned}
\sum_{k=1}^{N_r} \sum_{\ell=1}^{N_r} \bar{N}_p^k \bar{N}_p^{\ell} ((\Delta d_L)_k^j (\Delta d_L)_\ell^m)^n &= \sum_{j=1}^{N_{\text{det}}} ((\Delta d_L(r_j))^2) \delta_{im} \\
&= \sum_{k=1}^{N_r} \bar{N}_p^k ((\Delta d_L(r_k))^2) \delta_{im}.
\end{aligned}
\quad (A21) $$

Inserting this in the variance and using that $N_{\text{sky}} \bar{N}_{\text{det}} = N_{\text{tot}}$ we obtain

$$ \text{var}(X) = \frac{3N_{\text{sky}}}{N_{\text{tot}}^2} \left\{ \sum_{k=1}^{N_r} \bar{N}_p^k ((\Delta d_L)_k^2) \sum_{k=1}^{N_r} (\bar{d}_L)_k^2 \bar{N}_p^k - \frac{1}{\bar{N}_{\text{det}}} \left( \sum_{k=1}^{N_r} (\bar{d}_L)_k \bar{N}_p^k \right)^2 \right\}, \quad (A22) $$

where we have used that

$$ \sum_{i=1}^{N_{\text{sky}}} (\mathbf{n}_i \cdot \mathbf{n}')^2 = \frac{1}{\Delta\Omega} \int d\Omega (\mathbf{n} \cdot \mathbf{n}')^2 = \frac{N_{\text{sky}}}{3}. \quad (A23) $$

The calculation of the variance of $Y = \hat{d}_L$ is exactly the same, except for the factor -3 which is not in Y and the fact that $\hat{d}_L$ does not contain the product $\mathbf{n}_i \cdot \mathbf{n}'$ that is present in X, see equation (29). As a consequence, equation (A23) is replaced by

$$ \sum_{i=1}^{N_{\text{sky}}} = N_{\text{sky}}, \quad (A24) $$

leading to $\text{var}(Y) = 3\text{var}(X)$.

The calculation of the covariance between X and Y is also similar. In this case, equation (A23) is replaced by

$$ \sum_{i=1}^{N_{\text{sky}}} \mathbf{n}_i \cdot \mathbf{n}' = \frac{1}{\Delta\Omega} \int d\Omega \mathbf{n} \cdot \mathbf{n}' = 0, \quad (A25) $$

meaning that the covariance is exactly zero.

We can easily show that the mean X is proportional to $d_L^{(0)} \cdot v_0/c$, whereas the mean of Y is equal to $d_L^{(0)}$. As a consequence, the second term in equation (30) is suppressed by a factor $(v_0/c)^2$ with respect to the first one and we can neglect it. We therefore have

$$ \text{var}(\hat{v}_{d_L-n'}) = \frac{\text{var}(X)}{(d_L^{(0)})^2}. \quad (A26) $$

In Fig. A1, we show the variance of the luminosity distance for BNSs, as a function of $N_{\text{tot}}$ and for different values of the measurement uncertainty on $d_L$. Fig. 1 shows the same for BBHs.

## APPENDIX B: DERIVATION IN REDSHIFT SPACE

In equations (3), (18) and (23), we express luminosity distance, chirp mass, and number count dipoles in terms of integrals over r-dependent quantities. The comoving distance r is however not an observable quantity. On the other hand, the redshift z is observable. In this appendix, we first show how the dipoles at fixed z differ from the dipoles at fixed r. We then demonstrate that once we integrate over all sources, we obtain the same result, as expected.

We sketch the process as follows: consider any function $f(x, \epsilon)$ depending on one variable x which should be changed, and a small parameter $\epsilon$ (f may have other parameters that are neither small nor involved in the change of variable), such that expanding to first order we may write $f(x, \epsilon) = f^{(0)}(x) + \epsilon f^{(1)}(x)$. Assume further that f vanishes as $x \to 0$, $x \to \infty$. The goal is to perform the $\epsilon$-dependent change variable from x to $\tilde{x} = x + \epsilon g(x)$, i.e. $x = \tilde{x} - \epsilon g(\tilde{x})$. We obtain

$$ f(x, \epsilon) = f(\tilde{x} - \epsilon g(\tilde{x}), \epsilon) \approx f^{(0)}(\tilde{x}) + \epsilon \left( f^{(1)}(\tilde{x}) - \frac{d f^{(0)}(\tilde{x})}{d\tilde{x}} g(\tilde{x}) \right). \quad (B1) $$

Thus, in the change of variable, the same function f gains one extra term at first order.

In our case, x can be thought to be $\bar{z}$ (since a fixed value of r corresponds to a fixed value of the background redshift $\bar{z}$), and $\tilde{x} \equiv z \approx \bar{z} - (1 + \bar{z}) \mathbf{n} \cdot \frac{\mathbf{v}_0}{c}$. Hence we have: $\epsilon \equiv \mathbf{n} \cdot \frac{\mathbf{v}_0}{c}$ and $g(x) = -(1 + x)$.

With this knowledge, we can write the number count density of detected sources per redshift bin. Recalling that

$$ r(\bar{z}) = c \int_{0}^{\bar{z}} \frac{dz'}{H(z')} = r \left( z + (1+z)\mathbf{n} \cdot \frac{\mathbf{v}_0}{c} \right), \quad (B2) $$

we use equation (B1) with

$$ 
\begin{aligned}
f &\equiv \frac{dN_{\text{det}}}{d\Omega d\bar{z} dm_{1,2}} (\bar{z}, \mathbf{n}, m_{1,2}, \rho > \rho_*) \\
&= \frac{c}{H(\bar{z})} \frac{dN_{\text{det}}}{d\Omega dr dm_{1,2}} (r(\bar{z}), \mathbf{n}, m_{1,2}, \rho > \rho_*) ,
\end{aligned}
\quad (B3) $$

of which the functional form has been worked out in equation (8) until the dipolar order. Here H is the Hubble parameter. We end up with

$$ 
\begin{aligned}
& \frac{dN_{\text{det}}}{d\Omega dz dm_{1,2}}(z, \mathbf{n}, m_{1,2}, \rho > \rho_*) \\
&= \frac{d\bar{z}}{dz} \frac{dN_{\text{det}}}{d\Omega d\bar{z} dm_{1,2}} \left( z + (1+z)\mathbf{n} \cdot \frac{\mathbf{v}_0}{c}, \mathbf{n}, m_{1,2}, \rho > \rho_* \right),
\end{aligned}
\quad (B4) $$

from which equation (B1) allows us to write the monopole as

$$ 
\begin{aligned}
& \frac{d\bar{N}_{\text{det}}}{d\Omega dz dm_{1,2}}(z, m_{1,2}, \rho > \rho_*) = \\
& \frac{c}{H(z)} \frac{d\bar{N}_{\text{det}}}{d\Omega dr dm_{1,2}} \left( c \int_0^z \frac{dz'}{H(z')}, m_{1,2}, \rho > \rho_* \right),
\end{aligned}
\quad (B5) $$
while the dipole term is rewritten in terms of the monopole as

$$ \mathbf{n} \cdot \frac{\mathbf{v}_0}{c} \left\{ \frac{d\bar{N}_{\text{det}}}{d\Omega dz dm_{1,2}} \left[ 3 + s \left( \frac{1}{3} + \mathcal{A} \right) \right] + (1+z) \frac{d}{dz} \left( \frac{d\bar{N}_{\text{det}}}{d\Omega dz dm_{1,2}} \right) \right\}, \quad (B6) $$

where $\mathcal{A}$ and $\mathcal{A}$ are the functions defined in equations (9) and (10), evaluated at the arguments $(c \int_{0H(z')}^{z} \frac{dz'}{H(z')}, m_{1,2}, \rho_*)$ and $(c \int_{0H(z')}^{z} \frac{dz'}{H(z')}, m_{1,2})$ respectively.

Once equipped with the dipole expansion of detected events number count density per redshift, we can proceed analogously to Section 2 and introduce the dipoles in the luminosity distance and chirp mass distributions of detected events.

For the luminosity distance, the expansion in velocity at fixed $z$ has been done in Bonvin et al. (2006a) and it reads

$$ d_L(z, \mathbf{n}) = \bar{d}_L(z) \left( 1 + \frac{\mathbf{n} \cdot \mathbf{v}_0}{\mathcal{H}(z) \bar{r}(z)} \right), \quad (B7) $$

with $\mathcal{H}(z)$ the comoving Hubble parameter, and where $\bar{r}(z)$ is the monopole of the velocity expansion of $r$ at fixed $z$, which from the third expression in equation (B2) reads: $\bar{r}(z) = c \int_0^z \frac{dz'}{H(z')}$.

The redshifted chirp mass does not have a dipole with respect to fixed observed redshift $z$ slices, since it is simply the product of the source chirp mass with $(1+z)$, which is obviously constant on a slice of constant $z$.

As expected, the dipoles at fixed $r$ are different from the dipoles at fixed $z$. Let us now see what happens once we integrate over $r$ and $z$.

Integrating equation (B1) between 0 and $\infty$, we obtain

$$ \begin{aligned} & \int dx f^{(0)}(x) + \epsilon \int dx f^{(1)}(x) = \int dx f(x, \epsilon) = \int d\bar{x} f(x(\bar{x}), \epsilon), \epsilon \frac{dx}{d\bar{x}} \\ &= \int d\bar{x} f^{(0)}(\bar{x}) + \epsilon \int d\bar{x} \left( f^{(1)}(\bar{x}) - \frac{df^{(0)}(\bar{x})}{d\bar{x}} g(\bar{x}) - f^{(0)}(\bar{x}) \frac{dg(\bar{x})}{d\bar{x}} \right) \\ &= \int d\bar{x} f^{(0)}(\bar{x}) + \epsilon \int d\bar{x} f^{(1)}(\bar{x}), \end{aligned} \quad (B8) $$

where in the last step we used integration by parts and the assumed asymptotic behaviour for $f$ to make the boundary term vanish. The equivalence of the integrals holds both for the $\epsilon^0$ and $\epsilon^1$ terms.

Note that if the integration boundaries over $x$ are $[0, \infty[$, they should in principle become $[\epsilon g(0), \infty[$ for the $\bar{x}$ integration. However, at first order in $\epsilon$, we may take integrals over $[0, \infty[$ for $\bar{x}$ as well, the correction between the two being of order $\epsilon^2$.

With this, we can build the monopoles and dipoles in the distribution of detected luminosity distances and chirp masses, analogously to equations (3) and (18). The redshift integrals over the dipoles at fixed $z$ are equivalent to the $r$ integrals over dipoles at fixed $r$, by the above argument.

This paper has been typeset from a TeX/LA TeX file prepared by the author.