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

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$\sigma$. 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, that 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 percent compared to number counting of binary black holes alone. We find that a few $10^6$ 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^5$ 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.


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. 1f 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 et al. 2018;Secrest et al. 2021;Siewert et al. 2021;Secrest et al. 2022) at redshifts 0 ≤  ≲ 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).Supernovae type Ia light curves (see e.g.Riess et al. (1995)) provide an alternative means of assessing the dipole anisotropy (Bonvin et al. 2006b).Indeed, measurements from supernovae catalogs performed in Singal (2022); Horstmann et al. (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 et al. (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. 2022), 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. (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 Stiskalek et al. (2021); Essick et al. (2023); Payne et al. (2020), 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, followed by the number counts of BBHs and BNSs.This is not surprising given the predominant impact of this parameter on the BNS inspiral waveform; the possibility of exploiting this fact has already been subject of various investigations (Finke et al. 2022;Chernoff & Finn 1993;Taylor et al. 2012;Taylor & Gair 2012).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 6 events, achieved with 5 years 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 5 events already, achieved in a year of observation.
The rest of the paper is structured as follow: in Sec. 2 we derive a theoretical modelling of the luminosity distance and chirp mass dipole.In Sec. 3 we define estimators of these dipoles and we compute their variance and covariance.In Sec. 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 Sec. 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 Sec. 6.

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  0 / in the derivation.

Luminosity distance
The luminosity distance of a source situated at conformal distance  from the observer is defined as the ratio between the intrinsic luminosity of the source, , and the observed flux, where 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  0 / the luminosity distance is given by (see Bonvin et al. (2006a) 2 ) where d () 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% 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/  .Combining measurements from different interferometers allows one to measure   for each binary system.We can then compute the mean luminosity distance from a population of binary systems in direction n. Here, denotes the number of events detected at distance  with a signalto-noise ratio (SNR)  above a threshold value  * , summed over all intrinsic binary masses  1 and  2 3 .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) where M is the redshifted chirp mass as measured in the detector frame (in the following we drop the "redshifted" for simplicity).In Eq. ( 5),  is the total detector frame mass of the system and is the GW frequency corresponding to the innermost stable circular orbit.As to Θ 2 , it is a geometrical factor that accounts for the binary inclination angle and average detector's antenna patterns: where  is the angle formed by the normal vector of the orbital plane and the line of sight.The numerical prefactor in the above equation represents an average value of all the detector's antenna patterns.The function F is calculated from the integral of the power spectral density   (  ) 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  * 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  0 / is given by 3 Note that here we have dropped the dependence on the SNR threshold,  * , in the left-hand side of Eq. ( 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.
A (,  1,2 ) ≡ 1 Inserting Eqs. ( 8) and (2) into Eq.( 3) and keeping only linear terms in  0 / we finally obtain for the mean luminosity distance where the monopole is given by The fractional dipole is given by where the prefactor    takes the form with (,  1,2 ,  * ) being the mean probability density function of detected sources, Note that in Eq. ( 13) we have introduced the parameter    factorising out a minus sign in such a way that    = 1 in the case where there are no threshold effects, i.e. when  = 0. We see from Eq. ( 14) that if there are threshold effects but  and A are constant in , then the integral vanishes and    = 1.We can indeed rewrite the monopole as 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 ( 0 /) 2 ).Inserting this into Eq.( 14) we see that the integral exactly vanishes if  and A are independent on .Since these functions are expected to evolve slowly with , we expect threshold effects to be partially suppressed, meaning that    will be close to 1, even in the case where  is non-zero.In Sec.4.4, we compute    from our population model of BNSs (where threshold effects are present), and compare this with measurements of    from our synthetic catalogues.We find that    can be well predicted and, as expected, that it is very close to 1.This is important, because    is fully degenerated with the observer velocity  0 / and if we do not know it, we cannot extract information on  0 / from a population of sources with threshold effects.

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 velocity, to first order, the redshifted chirp mass is given by where 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 Eq. ( 3), we can compute the mean detected chirp mass in direction n, averaged over all binary systems (i.e. over all masses  1,2 ) as This can be expanded in powers of n • v 0  as a monopole, M (0) , and dipole term: where Combining Eqs. ( 8) and ( 17) allows to write the fractional dipole term as where we introduced  M as the chirp mass analogue of    : As for the luminosity distance, we see that threshold effects are nonzero only if  and A vary with distance .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 Eq. ( 5).This has an impact on the modelling of threshold effects, leading to a different expression for the dipole signal.

Number counts
The dipole in the GWs number count has been derived in Mastrogiovanni et al. (2023) and follows directly from Eq. ( 8).The total number of GW events detected in direction n is given by where with

STATISTICAL ESTIMATORS OF THE DIPOLES
We now build statistical estimators for the three dipole signals defined in Sec. 2.

Luminosity distance
For each direction n ′ in the sky, we can build the following observable where  ′ is the angle between the (a priori unknown) dipole velocity direction and 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,  0 /, at the maximum.
Let us now build a statistical estimator for the observable    −n ′ .We divide the sky in  sky angular pixels of same solid angle, and associate a vector n  pointing to the center of each pixel .Within any pixel , there are   det detected events (labeled by ), whose corresponding luminosity distance is labeled by (  )   .The estimator is then defined as 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.
We then compute the variance of this estimator.It can be written as ratio of two variables / with and  = d defined in Eq. ( 28).The variance of such a ratio can be written as (if the variances of  and  are significantly smaller than their mean) var We can show that the third term in Eq. ( 30) is vanishing and the second one is smaller than the first one by a factor ( 0 /) 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  and can be written as To compute this, we divide each solid angle  into   bins in distance  of size Δ, and we rewrite the sum over  as a sum over the -bins.The sky is therefore divided into pixels of angular size ΔΩ and radial size Δ.We denote by    the number of events detected in the pixel centered in direction n  and distance   and by (  )   the corresponding luminosity distance.Due to shot noise, the number of events fluctuates around the mean from pixel to pixel, Note that here we neglect the fluctuations in    due to the uncertainty on sky localisation.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 where (Δ  )   is the error in the measurement of   from one binary system, in angular bin  and radial bin .Note that in Eqs. ( 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) where  tot is the total number of events detected and ⟨(Δ  )  ⟩ denotes the typical error on luminosity distance measurement at distance   .The term in the first line of Eq. ( 35) is the contribution from luminosity distance uncertainty to the variance of the estimator.If ⟨(Δ  )  ⟩ is independent on distance, we see that the variance from this term scales as 1/ tot , as expected for  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, that 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 where (,  * ) is the radial distribution of detected sources integrated over masses.In Fig. 1, we plot the variance (36) for BBHs as a function of  tot for different relative errors on the luminosity distance measurements, ranging from 10% to 100% (and constant in ).We use the radial distribution of sources obtained from our simulated catalogues (see Sec. 4 for more detail).We see that both the shot noise contribution and the distance-uncertainty contribution scale as 1/ tot .The relative importance of the two terms depends therefore only on the uncertainty on   .For an uncertainty of up to 20%, the shot noise contribution completely dominates in the variance.If the uncertainty reaches 50% 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, that can be observed at higher redshift than BNSs).In the following we will assume a 20% uncertainty on the measurement of   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.

Chirp mass
The same procedure can be applied to the chirp mass.We first build an observable equivalent to (26), We repeat the process of dividing the detected observations in sky angular pixels, as in Eq. ( 27), to associate the following statistical estimator to As for the luminosity distance, this estimator is found to be unbiased assuming that shot noise and uncertainty in the measurement of the chirp mass are subdominant compared to the mean.For the variance, we generalise the computation of Sec.3.1 and Appendix A done for the luminosity distance to have data not only distributed in   radial bins, but also in   bins of the two source masses,  1,2 , which constitute the binary.This is needed since, unlike the luminosity distance, the redshifted chirp mass function depends on all (,  1,2 ).As for the luminosity distance, we neglect subdominant Variance of the luminosity distance estimator for BBHs, plotted as a function of the total number of events  tot .The different panels are for different values of the relative error on   (assumed to be independent of redshift).We show separately the shot noise contribution and the distance-uncertainty contribution, as well as the total.
contributions from the observer velocity in the variance.We obtain where the -index is the radial bin index, while the ℓ-index is the index over the various  1,2 bins.Mℓ is the expected redshifted chirp mass in such a bin, while N ℓ  is the expected number of detected events in that bin, i.e. such that N   =   ℓ=1 N ℓ  .The quantity ΔM ℓ denotes the typical uncertainty in the measurement of the chirp mass in a radial bin centered at   and in the mass bin labelled by ℓ.
Taking the continuous limit, the variance (40) becomes 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  1,2 .This effect is encoded in the dependence of the distribution function  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  tot assuming a 10% uncertainty in the measurement of the chirp mass (expected for BBHs) and a 1% 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 2 ⊙ , while BBHs are observed between 5 ⊙ and 100 ⊙ .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.

Number counts
The dipole estimator for the number counts and its variance have been derived in detail in state the final result for completeness: The variance of the number count estimator is only due to shot noise and is given by (44)

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.

Mass-distance covariance
The covariance is calculated following the same steps as for the variance in Secs.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   − M  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 Taking the continuous limit we obtain From Eq. ( 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.

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 Using that the number of detected events follows a Poisson distribution, we can easily show that the two terms in Eq. (47) exactly vanish.
A similar calculation shows that the covariance between mass and number counts also vanish.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.

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:  0 / cos  ′ .We can therefore look for an optimal estimator of this quantity, i.e. an estimator that would maximise the signal-tonoise ratio.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 diagonalise that part of the covariance matrix.We obtain two new estimators (per population) given by with The denominator in Eq. ( 48) insures that the new estimators have mean  0 / cos  ′ .
We can now build an optimal estimator of  0 / cos  ′ by linearly combining the independent estimators in the following way with  labeling 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 Note that the normalisation in Eq. ( 50) is there to keep the condition ⟨ videal ⟩ =  0 / cos  ′ , 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 's are different.In this case, we need to divide each estimator by its respective  before combining them.Since the 's can be determined only with limited precision (using population models, see Sec. 4.4 for more detail), this would induce additional contributions to the variance of the optimal estimator and degrades the SNR.In Sec.5.2, we quantify this effect using a Fisher analysis.

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.

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 where ,  and   are parameters that control the merger rate model, while d  d is the differential of the comoving volume.As fiducial values, we take  = 2.7,  = 3 and   = 2.The distribution of cos  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  ′ =  −  0 / sin , where  is the angle between the source position and the observer velocity.The luminosity distance and chirp mass are modified according to Eqs. ( 2) and ( 17) respectively.We produce two copies of both the BBHs and BNSs catalogue, one with the "CMB value" of the observer velocity,  0 / = 1.2 • 10 −3 and the other one with the "AGN value" of the observer velocity, which is 5 times larger:  0 / = 6 • 10 −3 .
We can then calculate the SNR of each event using Eq. ( 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 5 .If the SNR from this network exceeds a detection threshold of  * = 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% of its true value.For the chirp mass, we use a scattering of 1% for BNSs and of 10% for BBHs.These are typical errors that we might expect to obtain with XG detectors (Iacovelli et al. 2022).Note that the 20% 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 Eqs.( 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  ′ between the true direction of the dipole and a chosen direction n ′ .The spread in the signal comes from the fact that for a given angle  ′ 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   (2023)).As a consequence, at the maximum, i.e. when n ′ coincides with the direction of the observer velocity, the estimators are roughly equal to  0 / = 6 • 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 Sec.4.4, the amplitude of threshold effects are actually small for BNSs, below 10%, which is smaller than the spread in the signal.
Fig. 3 also reports the 1, 2, 3 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 Sec. 3.For all estimators, the fluctuation levels obtained through sky shuffling agree with the theoretical calculation of the variance from Eqs. ( 36), ( 41) and ( 44), that are indicated with dashed lines on the plot (see Sec. 4.2 for a more detailed comparison).When the estimator values exceed a certain noise threshold, the cosmic dipole can be detected (see Sec. 4.3 for a more in-depth discussion on detectability).

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 Eq. ( 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 (,  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  tot .The agreement between the simulated and theoretical variances is excellent.As expected from Eqs. ( 36), ( 41) and ( 44) the variances scale as  −1 tot .For the number count estimator, the variance is the same for BBHs and BNSs: it depends 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 Sec.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  1,2 -distribution of sources.If the chirp mass would be constant in  1,2 and in , the last two terms in the second line of Eq. ( 41) would cancel each other.Although the chirp mass does depend on  1,2 and , 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 .The luminosity distance varies indeed as (1 + ), while the redshifted chirp mass varies as (1 + ), i.e. significantly slower.
Since the means of the three estimators are very similar (see Sec. 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  coefficients need to be modelled if one wants to measure the observer velocity.Including the uncertainty in the modelling of the 's in the analysis generates an extra contribution to the variance, which needs to be accounted for.In Sec. 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.

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 esti- False Alarm Probability (FAP) 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 year of observations we have taken 7.5 • 10 4 BBH detections and 10 5 BNS detections.The number of events expected in 6 months and in 3 and 5 years 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 realisations while the FAP threshold is calculated using 2'000 noise realisations obtained with the sky shuffling method.The vertical solid, dashed and dotted lines mark the 1, 2, 3 false alarm probabilities.ties, 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 Sec.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 efficiency6 .
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  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 one year, the network of ET + 2CE detectors would be able to observe 7.5 • 10 4 BBHs and 10 5 BNSs.In Fig. 7 we show the detection efficiency for the combined estimator for 6 months, 1 year, 3 years and 5 years of observations.On one hand, we find that, with 5 years of observations and a bit less than 10 6 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 one year of observations, thanks to the BNS population.Therefore, a nondetection 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 un-modelled systematic in the AGN measurements.

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  and A defined in Eqs. ( 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 (,  1,2 ,  * ) around  * = 9, and compute with finite differences the derivative  1.The histogram width is due to shot noise and distance and mass uncertainty, which vary over the 200 realisations.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.with respect to  * to obtain  through Eq. ( 9).This needs to be done for each (,  1,2 ) bin, at  * = 9.We then use interpolation in  and  1,2 to promote the binned  values to a function.The function A depends on F , which quantifies the 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    ,  M and   .
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  and A are constant in .
Since we know the observer velocity in our synthetic catalogues of BNS events, we can use the dipole to measure the parameters    ,  M and   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  0 / 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 centered 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 threshold effects of the order of 10% 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 's by a factor 5. Second, the plots show excellent agreement between the theoretical modelling and the measured 's.Since in practice the '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 's on the population model.We vary the value of one of the parameter in the population model by 10% around its fiducial value and compute the histograms for the 's.Comparing the width in the '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 's.This means that the uncertainty from the choice of model is smaller than the variance of the dipole.

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,  0 /, 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   ,    and  M are perfectly known; and the second 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 maximising the estimators νn ′ with respect to n ′ .

Known 𝛼's
We use the Fisher formalism to estimate the uncertainty on  0 / 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: The Fisher for the (unique) parameter  0 / is then given by where for  = 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 Sec.3.4.2.The error on  0 / is then given by  1.
The error,   0 / , 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 estimator 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   from 1 to 1.08.Combining the three BBH estimators improves the constraints by 20% compared to using the number counts alone, while combining the three BNS estimators improves the constraints by 32%.Combining all estimators improves the constraints by 50% compared to the number counts of BBHs alone (studied in Mastrogiovanni 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.

Adding uncertainties on 𝛼's
Whereas for the BBHs the  coefficients are known and equal to 1, for BNSs it is not the case.These coefficients are affected by threshold effects.As shown in Sec.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 '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  0 /.We consider the combination of all estimators and compute the error on  0 / for different values of the uncertainties on the 's.The signal is written as The Fisher matrix contains now four parameters: and  4 =  BNS M and it is given by The covariance matrix has the form cov = (60) where the 3 × 3 blocks cov BBH and cov BNS are given by Eq. ( 56).In Eq. ( 60),  denotes the relative uncertainty on the determination of the  parameters for the BNS population, that we assume here to be the same.We consider three cases for : 10%, 20% and 50%.The error on  0 / is then given by 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% uncertainty on all three  BNS degrades the constraints on  0 / by at most 1% compared to the case where these parameters are assumed to be known perfectly.Increasing this uncertainty to 50% degrades the constraints by at most 5%.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  is as large as 50%, we still gain information by including BNSs.For example, for 10 6 events and a dipole consistent with the CMB one, we gain 29% in the measurement of  0 / by adding BNSs.If we can model the 's with a precision of 20%, this gain increases to 36% (38% for a 10% precision).Seen as a function of the  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  BNS is too large, no information is gained by adding BNSs.
Table 2. Error,   0 / , obtained from: the individual estimators (first 6 columns), combining the three BBH estimators, the three BNS estimators, all estimators, or using the top 3 combination 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.

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% 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 one year 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.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  has been done in Bonvin et al. (2006a)  where in the last step we used integration by parts and the assumed asymptotic behaviour for  to make the boundary term vanish.The equivalence of the integrals holds both for the  0 and  1 terms.Note that if the integration boundaries over  are [0, ∞[, they should in principle become [(0), ∞[ for the x integration.However, at first order in , we may take integrals over [0, ∞[ for x as well, the correction between the two being of order  2 .
With this, we can build the monopoles and dipoles in the distribution of detected luminosity distances and chirp masses, analogously to Eqs. (3) and ( 18).The redshift integrals over the dipoles at fixed  are equivalent to the  integrals over dipoles at fixed , by the above argument.

Figure 3 .
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, 2, 3 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 Sec. 3. Top plot: Number count estimator.Middle plot: Chirp mass estimator.Bottom plot: Luminosity distance estimator.

Figure 4 .
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  tot .The variances and correlations are obtained by reshuffling 2'000 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 Sec. 3. We find that the maximum deviation from the expected value of zero for the number countchirp mass and number count -luminosity distance is 0.05, arising from the limited number of sky shufflings.

Figure 8 .
Figure 8.The plots indicate the distribution of the  parameters obtained for the 3 estimators with 200 population realisations.The vertical black lines indicate the theoretical prediction, see Table1.The histogram width is due to shot noise and distance and mass uncertainty, which vary over the 200 realisations.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.

Figure 9 .
Figure 9. Same as Fig. 8, but with the CMB value of the observer velocity.

Figure 10 .
Figure 10.The plots indicate the distribution of the  parameters obtained for the 3 estimators with 200 population realisations.Each realisation also considers a random realisation of the merger rate model parameters  with an uncertainty of 10% around its fiducial value of  = 2.7.The vertical black lines indicate the theoretical prediction from Table1for the  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.
To compute   0 / , we need to know the values of the coefficients    for the different populations and the different estimators.For the BBH, since threshold effects are negligible we have  BBH  =  BBH   =  BBH M = 1.For the BNS, we use the values calculated theoretically in Sec.4.4 and reported in Table and it reads  (, n) = d () 1 + n • v 0 H () r () ,(B7)with H () the comoving Hubble parameter, and where r () is the monopole of the velocity expansion of  at fixed , which from the third expression in Eq. (B2) reads: r() =  ∫  0 d ′  ( ′ ).The redshifted chirp mass does not have a dipole with respect to fixed observed redshift  slices, since it is simply the product of the source chirp mass with (1 + ), which is obviously constant on a slice of constant .As expected, the dipoles at fixed  are different from the dipoles at fixed .Let us now see what happens once we integrate over  and .Integrating Eq. (B1) between 0 and ∞, we obtain∫ d  (0) () +  ∫ d  (1) () = ∫ d  (, ) = ∫ d x  ( x, ),  d d x = ∫ d x  (0) ( x) +  ∫ d x  (1) ( x) − d  (0) ( x) d x ( x) −  (0) ( x) d( x) d x = ∫ d x  (0) ( x) +  ∫ d x  (1) ( x) ,(B8) of Mastrogiovanni et al.

Table 1 .
Expected values for the BNS  parameters, for a fiducial astrophysical population model.
Table1for the  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.

Table 3 .
Fisher bound on the error   0 / , obtained from the combination of the six estimators assuming different uncertainties on the '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.