# From the Lab to the Theater: An Unconventional Field Robotics Journey

Ali Imran1, Vivek Shankar Varadharajan2, Rafael Gomes Braga1, Yann
Bouteiller2,
Abdalwhab Bakheet Mohamed Abdalwhab1, Matthis Di-Giacomo1, Alexandra
Mercader1,
Giovanni Beltrame2 and David St-Onge1 *This work was supported by FRQ PRISME-
ART program.1 are with the Lab of INIT Robots, Mechanical Engineering, Ecole
de technologie supérieure, 1100 Notre-Dame W., Canada<EMAIL_ADDRESS>are with MISTLab, Computer and Software Engineering, Polytechnique Montréal,
Canada<EMAIL_ADDRESS>
###### Abstract

Artistic performances involving robotic systems present unique technical
challenges akin to those encountered in other field deployments. In this
paper, we delve into the orchestration of robotic artistic performances,
focusing on the complexities inherent in communication protocols and
localization methods. Through our case studies and experimental insights, we
demonstrate the breadth of technical requirements for this type of deployment,
and, most importantly, the significant contributions of working closely with
non-experts.

##  I Introduction

The intersection of art and technology has given rise to a burgeoning field of
robotic artistic performances (such as [1, 2, 3]). From synchronized
choreographies to interactive installations, these performances push the
boundaries of traditional artistic expression, inviting audiences to engage
with the fusion of human creativity and robotic technology. However, behind
the curtain of artistic innovation lies a host of technical challenges that
rival those encountered in more conventional field deployments of robotic
systems.

Artistic deployment of robotics shares many common aspects with the more
conventional fields of robotics. First and foremost, accurate localization
methods are crucial in any deployment of robots outside a controlled lab
environment. The operating environment could be outdoors with access to GPS,
or it might be a GNSS-denied setting [4], such as in mining applications, cave
exploration, or industrial environments.

Over the last decade, major field deployment targeted robot swarms, such as
the DARPA OFFSET [5, 6]. In these settings, communication is one of the most
critical aspects for coordination and information sharing [7, 8].

In this paper, we explore the technical challenges of orchestrating robotic
artistic performances based on the use case of an artistic project called
DESSAIM111DESSAIM web page. DESSAIM stems from an interdisciplinary research
project combining engineering and performing arts, exploring the
expressiveness of three types of robotic swarms. Recycled materials transform
the robots’ appearance and the environment in which they operate. From these
imaginary territories, the poetry of material and form in motion emerges.

As this project calls for the deployment of self-organised multi-robot systems
in close proximity to the audience, our core technical challenges were posed
by communication protocols and localization methods in dynamic performance
spaces as well as safety. While the primary focus of robotic performances may
be on aesthetic expression and emotive storytelling, the underlying technical
infrastructure is critical to the seamless execution of representations.

This paper first describes the approach used to design swarm behaviors, then
the robotic systems deployed in this field experiment, and finally the
technical challenges and their solutions in terms of communication and
localization.

##  II Behavior composition

The team of artists delved into swarm behaviors and their inspiration from
nature to create the performances. We explored together the well-known
aggregation, diffusion, flocking [9], Lennard-Jones [10] and pursuit [11]
behaviors. In the end, our team created a library of twelve emergent
behaviors, to be applied to three swarms of robots at different scales: small
tabletop, medium flying robots, and wheeled bases with manipulators. They each
call upon swarm behavior primitives, but with variations in dynamic states and
combinations that extend their expressivity. All these behaviors were
programmed in Buzz [12], a domain-specific language for robot swarms, and
validated in a low-fidelity fast simulator, ARGoS3 [13]. The use of these two
tools was instrumental in sharing knowledge on swarm intelligence and
iterating on their behavior compositions.

Figure 1: Behavior composition: A “firework” composed of 4 behaviors in ARGoS

As an example, Fig. 1 presents the firework behavior where all robots are
attracted towards a point and then spread simultaneously in opposite
directions like an expanding circle. Velocities and timing in both directions
are different and finely tuned to create the desired effect.

(a) Tabletop swarm: Sushis (b) Flying swarm: CrazyCognies (c) Human-scale
swarm: Doodies

Figure 2: The three DESSAIM robotic swarms seen in their performing settings.

The three platforms introduced in the following section run the exact same
Buzz scripts, robot-agnostic, in a Python wrapper adapted to their hardware
called (Py)Buzz. (Py)Buzz is our swarm controller developed for artistic
performances; it implements the Buzz Virtual Machine (BVM) [12] to run
collective behaviors through Buzz scripts.

##  III DESSAIM swarms

As shown in Fig. 2, three swarms of various sizes and dynamics performed in
DESSAIM: tiny tabletop Sushis, flying caged CrazyCognies, and human-size
mobile manipulators, the Doodies.

###  III-A Tabletop swarm: Sushis

Our Sushis are small tabletop robots derived from Stanford’s Zooids [14].

Figure 3: Sushi with its top circuit including RF communication chip and
photodiode for localization.

A 48MHz ARM microcontroller (STM32F051C8) manages onboard control and
communication using a 2.4GHz radio chip (NRF24L01+). For localization, two
photodiodes are placed on top to detect and decode gray-coded patterns
projected from the top by a 3000Hz projector (DLP LightCrafter 4500). Two
small DC motors actuate their motion in differential drive mode. Their small
100 mAh LiPo battery can withstand at least one hour of continuous
performance.

Our revision of the original design includes a larger enclosure to facilitate
their assembly and maintenance as well as to accommodate a larger gearbox for
the drive motors. The expressivity of their motion is highly coupled to their
dynamics (velocity and acceleration), which required a larger band of possible
velocities, provided with the new gearbox. In our setup, their control is now
tied to (Py)Buzz.

###  III-B Flying swarm: CrazyCognies

Our Crazycognies are hybrid indoor aerial robots made from the combination of
Cogniflies [15], an open-sourced collision-resistant quadcopter, and the
Crazyflie quadcopter commercialized by BitCraze. We developed these platforms
for safe deployments in tight and cluttered environments - a requirement of
indoor artistic performances. The Crazyflies were not powerful enough for our
needs (flight time and payload), whereas the Cogniflies were not ready for
advanced indoor positioning. The resulting system shown in Fig. 4 weighs 310
grams and can fly reliably for about 7 minutes.

Figure 4: CrazyCognie with its Crazyflie on top, a Raspberry Pi just
underneath, and the iNav flight controller below, just above the battery pack.

Figure 5 shows a schematic of the software components deployed on the
Crazycognies. The driver software (cognifly-python222GitHub code repository
for Cognifly-Python) communicates with the Flight Controller (FC) through
UART, and that software is designed to enable the deployment of custom pose
estimation. The custom estimator communicates with the Crazyflie through USB
for horizontal pose estimates (X, Y) and sends it to the FC to achieve stable
flight in the theater room. As for the altitude, the highly unreliable nature
of the altitude estimates from UWB required us to use TFmini lidar
measurements instead.

Figure 5: Software schematic of a CrazyCogni

###  III-C Human-scale swarm: Doodies

Our swarm of Doodies consists of a fleet of six Clearpath Dingo [16]
differential drive robots, each equipped with a Kinova Gen3 Lite [17]
manipulator. Similar to the Crazycognies, the Doodies are also equipped with a
loco Crazyflie deck for position estimation. They run (Py)Buzz, our custom
Kalman filter, and all hardware drivers on an onboard Nvidia Xavier. For
safety, the platform uses small 12V Lead-acid batteries, one for the base and
another for the arm, granting autonomy for continuous motion of more than 3
hours. As shown in Fig. 6, the custom turret brings the whole robotic platform
to a human scale and includes an emergency stop in the back to increase
safety.

Figure 6: Doody

Figure 7 shows the software architecture on each robot built over modified
versions of the drivers for the Clearpath Dingo base, Kinova Kortex for the
robotic arm, Crazyflie [18] firmware and node to leverage UWB localization,
and drivers for the cameras to provide interaction capabilities. On the
Doodys, (Py)Buzz is bridged with the ROS ecosystem. Localization of the
Doodies is achieved with a Kalman filter fusing the wheel odometry, the IMU
measurements, and the 2D pose estimates from the UWB setup. The resulting pose
is sent to (Py)Buzz for the swarm-level collective behaviors.

Figure 7: Software Schematic of a Doody

##  IV Artistic performance challenges

The three swarms undergo five iterative performances in front of an audience
of 30 people throughout the creation process. Each performance concludes with
a mediation session with the public to describe the robot functionalities and
gather their feedback and perception. The creation process culminates with the
presentation of five identical performances in front of audiences of 45 people
as part of a local festival in Montréal, Canada.

###  IV-A High-frequency critical communication for coordination

Coordination among robots to synchronize their states and achieve collective
behaviors relies heavily on inter-robot communication.

Reliable communication in the application layer of the ISO/OSI model can be
achieved using TCP/IP protocols; some existing packages like nimbro_network
[19] enable it with ROS. However, the configuration task complexity increases
rapidly with the growing number of robots and is challenging for large swarms.
We thus selected Zenoh [20], a publish/subscribe/query protocol that operates
with a set of abstractions in the 4 and 5 layers of the OSI model.

Figure 8: The figure reports the bandwidth used to operate and maintain a
swarm of five Doodys, five CrazyCognies, a mobile marker node (used as an
attractor/repulsor point for behaviors), and 2 ground stations. The vertical
dotted lines (Doodys) and dashed lines (CrazyCognies) correspond to specific
events: a launch command (red), behavioral state switch command (green), and
behavioral stop command (purple).

In the application layer, we implement publishers and subscribers for each
robot’s gossip messages through a single topic, enabling us to add and remove
robots seamlessly. Zenoh requires only a small communication overhead of 5
bytes per packet. Figure 8 shows the recorded communication bandwidth during
one of the artistic performances. The bandwidth spikes to over 4 MB/s when
connecting to the Doodys and launching the software around 250s, similarly
around 500s when launching the Crazycognies. During the normal operation of
both swarms, the bandwidth stays below 2 MB/s since each node in the swarm
sends gossip below 250 bytes. When both the swarms are shut down, the
bandwidth drops to below 0.5 MB/s. The communication setup introduced in this
work was proven to be stable to swarms with up to 13 nodes and exhibits
graceful communication scalability.

###  IV-B Minimal infrastructure global positioning

Relative localization is vital for robots to achieve collective behaviors.
Range and bearing information between neighboring robots is often enough for
swarm behavior deployment, but most technologies (IR, vision, etc.) require
line-of-sight; a complex thing to achieve with tens of people walking among
the robots. Furthermore, in a public performance, the robots always need to
undergo safe motion with respect to the theater equipment (speakers,
lightning, curtains) and the people. However, deploying a vision-based motion
capture system is not an option: the artistic shows change venues often. Each
performance has a varying stage lighting design. Moreover, the rooms are
rather large to cover (for instance, 6m x 12m). We opt for a UWB-based
strategy, using preconfigured anchors in the environment, to derive a rough
global estimation. The anchors were configured in TDOA3 [21] mode to ping the
robots’ UWB tags. The resulting estimates obtained from the UWB undergo
another layer of filtering, which is particularly important for robots that
are more sensitive to noise - such as the Crazycognis.

Figure 9: Comparison of our UWB pose estimation (orange) and ground truth from
an Optitrack motion capture (blue) system during a pursuit behavior executed
by 5 CrazyCognies.

Figure 9 compares the pose estimates of a robot when 5 Cogniflies display a
pursuit behavior. The pose estimates were obtained using UWB and a motion
capture system installed at the venue. The UWB-based pose estimates show
errors in the range of tens of centimeters, with occasional fluctuations
caused by the abrupt drone movements. The accuracy of the UWB-based systems is
significantly influenced by factors such as the initial calibration, the
number of anchors, and their configuration (gain, bandwidth, etc.). To enhance
this, our team has developed advanced calibration strategies, inspired by land
surveying techniques, enabling us to install our affordable UWB system
(<$1500) in less than 4 hours.

##  V CONCLUSIONS

We describe the three robotic swarms performing in the DESSAIM artistic
representations. Each of the swarms was running the exact same behavioral
script, as designed (composed) by the artists and the engineers. About 400
people experienced this intimate and immersive robotic performance, with the
audience freely moving among the performing robots. To achieve robustness
(capable of repeated performance) and safety, this paper highlighted two key
aspects: the Zenoh communication layer integration and the UWB localization.
We believe this report serves as a compelling demonstration of the value of
artistic manifestations as challenging field robotic deployments.

## ACKNOWLEDGMENTS

The authors would like to acknowledge the creative and empathetic contribution
of the talentful artists of DESSAIM: Audrey Rochette, Hélène Duval, Marcelle
Hudon, Magali Babin and Danielle Lecourtois. Several other engineering
students support this team effort: Jean Mazerolle, Corentin Boucher, Hao-Fang
Cheng, Clement Gassmann-Prince, Guillaume Ricard and Radhouen Khmiri.

## References

  * [1] B. Vorn, “I Want to Believe—Empathy and Catharsis in Robotic Art BT - Robots and Art: Exploring an Unlikely Symbiosis,” D. Herath, C. Kroos, and Stelarc, Eds. Singapore: Springer Singapore, 2016, pp. 365–377. [Online]. Available: https://doi.org/10.1007/978-981-10-0321-9_18
  * [2] C. P. Lansky, “How quantum physicist ballerina merritt moore learned to dance with a robot (plus, her newest film),” _Pointe Magazine_ , 10 2020. [Online]. Available: https://www.pointemagazine.com/merritt-moore-ballet-2648967672.html
  * [3] C. Cuan, “Output: Choreographed and recomposed human and industrial robot bodies across artistic modalities,” _Frontiers in AI, The Art of Human-Robot Interaction: Creative Perspectives from Design and the Arts_ , 2020.
  * [4] J. Braga, A. P. Aguiar, and J. B. de Sousa, “Coordinated multi-uav exploration strategy for large areas with communication constrains,” in _Iberian Robotics conference_. Springer, 2017, pp. 149–160.
  * [5] T. Rouček, M. Pecka, P. Čížek, T. Petříček, J. Bayer, V. Šalanskỳ, D. Heřt, M. Petrlík, T. Báča, V. Spurnỳ _et al._ , “Darpa subterranean challenge: Multi-robotic exploration of underground environments,” in _Modelling and Simulation for Autonomous Systems: 6th International Conference, MESAS 2019, Palermo, Italy, October 29–31, 2019, Revised Selected Papers 6_. Springer, 2020, pp. 274–290.
  * [6] C. Ju, J. Kim, J. Seol, and H. I. Son, “A review on multirobot systems in agriculture,” _Computers and Electronics in Agriculture_ , vol. 202, p. 107336, 2022.
  * [7] A. I. Mourikis and S. I. Roumeliotis, “Performance analysis of multirobot cooperative localization,” _IEEE Transactions on robotics_ , vol. 22, no. 4, pp. 666–681, 2006.
  * [8] J. Cortés and M. Egerstedt, “Coordinated control of multi-robot systems: A survey,” _SICE Journal of Control, Measurement, and System Integration_ , vol. 10, no. 6, pp. 495–503, 2017.
  * [9] M. Schranz, M. Umlauft, M. Sende, and W. Elmenreich, “Swarm robotic behaviors and current applications,” _Frontiers in Robotics and AI_ , vol. 7, p. 36, 2020.
  * [10] D. St-Onge, V. S. Varadharajan, I. Švogor, and G. Beltrame, “From design to deployment: Decentralized coordination of heterogeneous robotic teams,” _Frontiers in Robotics and AI_ , vol. 7, pp. 1–16, 2020.
  * [11] D. St-Onge, C. Pinciroli, and G. Beltrame, “Circle formation with computation-free robots shows emergent behavioural structure,” 2019, pp. 5344–5349.
  * [12] C. Pinciroli and G. Beltrame, “Buzz: An extensible programming language for heterogeneous swarm robotics,” in _2016 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_. IEEE, 2016, pp. 3794–3800.
  * [13] C. Pinciroli, V. Trianni, R. O’Grady, G. Pini, A. Brutschy, M. Brambilla, N. Mathews, E. Ferrante, G. Di Caro, F. Ducatelle, M. Birattari, L. M. Gambardella, and M. Dorigo, “ARGoS: a modular, parallel, multi-engine simulator for multi-robot systems,” _Swarm Intelligence_ , vol. 6, no. 4, pp. 271–295, 2012.
  * [14] M. Le Goc, L. H. Kim, A. Parsaei, J.-D. Fekete, P. Dragicevic, and S. Follmer, “Zooids: Building blocks for swarm user interfaces,” in _Proceedings of the 29th annual symposium on user interface software and technology_ , 2016, pp. 97–109.
  * [15] R. de Azambuja, H. Fouad, Y. Bouteiller, C. Sol, and G. Beltrame, “When being soft makes you tough: A collision-resilient quadcopter inspired by arthropods’ exoskeletons,” in _2022 International conference on robotics and automation (ICRA)_. IEEE, 2022, pp. 7854–7860.
  * [16] Clearpath Robotics, “Dingo indoor mobile robot,” https://clearpathrobotics.com/dingo-indoor-mobile-robot/, 2024, accessed: 2024-03-19.
  * [17] Kinova Robotics, “Discover our gen3 lite robot,” https://www.kinovarobotics.com/product/gen3-lite-robots, 2024, accessed: 2024-03-19.
  * [18] Bitcraze, “Crazyflie platform overview,” https://www.bitcraze.io/documentation/system/platform/, 2024, accessed: 2024-03-19.
  * [19] Max Schwarz, “nimbro network github page,” https://github.com/AIS-Bonn/nimbro_network, 2024, accessed: 2024-03-19.
  * [20] A. Corsaro, L. Cominardi, O. Hecart, G. Baldoni, J. Enoch, P. Avital, J. Loudet, C. Guimarães, M. Ilyin, and D. Bannov, “Zenoh: Unifying communication, storage and computation from the cloud to the microcontroller,” vol. DSD 2023, 09 2023.
  * [21] W. Zhao, J. Panerati, and A. P. Schoellig, “Learning-based bias correction for time difference of arrival ultra-wideband localization of resource-constrained mobile robots,” _IEEE Robotics and Automation Letters_ , vol. 6, no. 2, pp. 3639–3646, 2021.


# Detection and estimation of the cosmic dipole with the Einstein Telescope
and Cosmic Explorer

S. Mastrogiovanni1,2, C. Bonvin3, G. Cusin3,4, S. Foffa3
1 INFN, Sezione di Roma, I-00185 Roma, Italy
2 Artemis, Université Côte d’Azur, Observatoire de la Côte d’Azur, CNRS,
F-06304 Nice, France
3Université de Genève, Département de Physique Théorique and Gravitational
Wave Science Center, 24 quai Ernest-Ansermet,
CH-1211 Genève 4, Switzerland
4Sorbonne Université, CNRS, UMR 7095, Institut d’Astrophysique de Paris, 75014
Paris, France
<EMAIL_ADDRESS>0000-0003-1606-4183 0000-0002-5318-4064 0000-0002-4530-3051

(Accepted XXX. Received YYY; in original form ZZZ
ET-0195A-22)

###### Abstract

One of the open issues of the standard cosmological model is the value of the
cosmic dipole measured from the Cosmic Microwave Background (CMB), as well as
from the number count of quasars and radio sources. These measurements are
currently in tension, with the number count dipole being 2-5 times larger than
expected from CMB measurements. This discrepancy has been pointed out as a
possible indication that the cosmological principle is not valid. In this
paper, we explore the possibility of detecting and estimating the cosmic
dipole with gravitational waves (GWs) from compact binary mergers detected by
the future next-generation detectors Einstein Telescope and Cosmic Explorer.
We model the expected signal and show that for binary black holes, the dipole
amplitude in the number count of detections is independent of the
characteristics of the population and provides a systematic-free tool to
estimate the observer velocity. We introduce techniques to detect the cosmic
dipole from number counting of GW detections and estimate its significance. We
show that a GW dipole consistent with the amplitude of the dipole in radio
galaxies would be detectable with $>3\sigma$ significance with a few years of
observation ($10^{6}$ GW detections) and estimated with a $16\%$ precision,
while a GW dipole consistent with the CMB one would require at least $10^{7}$
GW events for a confident detection. We also demonstrate that a total number
$N_{\rm tot}$ of GW detections would be able to detect a dipole with amplitude
$v_{o}/c\simeq 1/\sqrt{N_{\rm tot}}$.

###### keywords:

gravitational waves – cosmology: cosmic background radiation – galaxies:
active

††pubyear: 2022††pagerange: Detection and estimation of the cosmic dipole with
the Einstein Telescope and Cosmic Explorer–A

##  1 Introduction

Our motion through the Universe generates a dipole in both the temperature
anisotropies of the Cosmic Microwave Background (CMB) (Planck Collaboration et
al., 2020) and in the angular distribution of electromagnetic sources (Colin
et al., 2017; Bengaly et al., 2018; Secrest et al., 2021; Siewert et al.,
2021; Secrest et al., 2022). If the cosmological principle is valid, these two
measures should have consistent values. However, it is a longstanding problem
that number counts of radio sources and of quasars at low and intermediate
redshifts exhibit a dipole that is well aligned with that of the CMB but with
an amplitude which is 2-5 times larger than expected, leading to a tension
reaching up to $\sim 5\sigma$. In Dalang & Bonvin (2022) it is argued that
this tension might be alleviated once one takes into account the redshift
evolution of the population of sources, and the value that the evolution rate
should have in order to remove the tension is found.

Gravitational-wave (GW) sources observed at cosmological distances can shed
light on the cosmic dipole problem. Since the first historical detection of
Abbott et al. (2016) during the first scientific observation run (O1), the
rate of observed GW events has drastically increased, to reach roughly one
detection per week in the O3 run. In total, about 90 binary black hole (BBH)
coalescences, as well as two binary neutron star (BNS) and two neutron star -
black hole mergers have been detected so far (The LIGO Scientific
Collaboration et al., 2021b). Thanks to this abundance of new data, which will
further accumulate in the forthcoming future (Abbott et al., 2018), it is
possible to probe cosmology (Abbott et al., 2021a; The LIGO Scientific
Collaboration et al., 2021a; Mancarella et al., 2022; Leyde et al., 2022) and
astrophysical rates of compact binary coalescences (Abbott et al., 2019,
2021b; The LIGO Scientific Collaboration et al., 2021c) using GW observations.

It is thus not surprising that recently, it has been proposed to use GW
sources as a new and independent probe to measure the cosmic dipole. A first
approach proposed to measure the cosmic dipole is by studying anisotropies in
the GW stochastic background (Cusin & Tasinato, 2022; Valbusa Dall’Armi et
al., 2022; Bartolo et al., 2022). However, stochastic GW backgrounds have not
been detected yet. Another possibility, is to study the sky distribution of
the transient GW sources currently detected (Cai et al., 2018; Kashyap et al.,
2022; Essick et al., 2022). Although GW sources are the central paradigm of
all these works, the methods employed significantly differ. In Cai et al.
(2018), the authors perform a forecast on the accuracy with which the cosmic
dipole would be detected by fitting for the (modified) luminosity distance
distribution of GW events. In Kashyap et al. (2022), the authors try instead
to constrain the dipole anisotropy using the BBHs mass distribution from
current GW events (The LIGO Scientific Collaboration et al., 2021b) finding
that the mean mass is higher in the direction of the CMB dipole and invoking
the need of further study to understand the origin of their findings. In
contrast to Stiskalek et al. (2021), Kashyap et al. (2022), Essick et al.
(2022) find no-evidence of anisotropies, using 63 GW sources distributed over
the sky, and implementing a hierachical Bayesian analysis that takes into
account selection biases.

In our work, we focus on the detection and estimation of the cosmic dipole
with Einstein Telescope (ET) (Punturo et al., 2010) and Cosmic Explorer (CE)
(Reitze et al., 2019). These detectors, along with LISA (Amaro-Seoane et al.,
2017), belong to the so-called next generation (XG) of GW detectors and
represent the aim of the community to substantially scale-up the experience of
the LIGO-Virgo era. In this paper, we show that the number counts of BBH
sources will likely offer an optimal tool to detect the cosmic dipole. We
build an estimator, that has the advantage of being independent of unknown
characteristics of the GW sources (or their evolution) and that will therefore
allow us to infer our motion with respect to the cosmological frame in a
robust way.

This paper is structured as follows. In Sec. 2, we derive the dipole
modulation in the GW sources number counts, due to the observer velocity with
respect to the cosmological frame. We show in this section that GW number
counting of BBHs offers a clean tool to evaluate anisotropies. In Sec. 3, we
simulate GW detections using a detector network composed by ET and two CEs
(ET+2CE). We discuss detection capabilities and the measurement process using
a frequentist and Bayesian approach. In Sec. 4, we discuss possible limiting
factors of our approach in light of the possibility of constraining the cosmic
dipole with GW sources. Finally, in Sec. 5 we draw our conclusions.

##  2 Modeling the impact of the observer velocity on GW number counts

We write down explicitly the expression for the cosmic dipole in the GW source
distribution. In Sec. 2.1, we derive the theoretical modeling for the dipole
induced by the observer velocity in the number counts of GW signals emitted by
compact binary coalescences. In Sec. 2.2, we define a statistical estimator
for the cosmic dipole.

###  2.1 Theoretical framework

We define the number of GW sources $N_{\rm det}$ detectable per unit solid
angle $d\Omega$ and distance bin $dr$, in direction ${\bf{n}}$ at comoving
distance $r$ and with signal-to-noise ratio (SNR) $\rho$ larger than a given
threshold $\rho_{*}$ as

$\frac{dN_{\rm det}}{d\Omega
dr}(r,{\bf{n}},\rho>\rho_{*})\equiv\int_{\rho_{*}}^{\infty}d\rho\frac{dN}{d\Omega
dr}(r,{\bf{n}},\rho)\,.$ (1)

This number depends on the direction of observation ${\bf{n}}$ due to three
effects. First, the GW sources are not perfectly homogeneously distributed:
they live in galaxies which follow the large-scale structure of the Universe.
Second, the propagation of GWs is affected by inhomogeneities along the
trajectory, that change the apparent distribution of sources over the sky.
Line of sight effects include lensing and local matter effects, and also
effects due to the source peculiar motion see e.g. Takahashi & Nakamura
(2003); Barausse et al. (2014); Bonvin et al. (2017); Sberna et al. (2022);
Cusin & Tamanini (2021); Toubiana et al. (2021); Cusin et al. (2021); Bonvin
et al. (2022); Toscani et al. (2023) and Cusin et al. (2017); Pitrou et al.
(2020) in the context of a stochastic background. Finally, the motion of the
observer with respect to the source rest frame generates a further anisotropy
in the observed distribution. In this work, we are interested in the latter
effect, which gives a dipolar modulation in the distribution of sources. The
other two effects will also have a dipolar contribution, but, as has been
shown for quasars (Secrest et al., 2021), this clustering dipole is expected
to be negligible compared to the kinematic dipole. We checked that this
hypothesis is valid also for the BBH and BNS merger rate model used in this
work (see Sec. 3) for more details.

The cosmic dipole averaged over all distances $r$ is obtained from the
difference between Eq. (1) and its angular average, integrated over $r$.
Namely,

$\displaystyle\Delta({\bf{n}})\equiv\mathcal{D}[{\bf{n}}\cdot\hat{{\bf{v}}}_{o}]$
$\displaystyle=\dfrac{\int dr\left[\frac{dN_{\rm det}}{d\Omega
dr}(r,{\bf{n}},\rho>\rho_{*})-\frac{d\bar{N}_{\rm det}}{d\Omega
dr}(r,\rho>\rho_{*})\right]}{\int dr\frac{d\bar{N}_{\rm det}}{d\Omega
dr}(r,\rho>\rho_{*})}\,,$ (2)

where ${\bf{v}}_{o}$ denotes the observer velocity and
$\hat{{\bf{v}}}_{o}\equiv{\bf{v}}_{o}/|{\bf{v}}_{o}|$ is its direction. The
angular average over the sky, denoted with a bar, is given by

$\frac{d\bar{N}_{\rm det}}{d\Omega
dr}(r,\rho>\rho_{*})\equiv\frac{1}{4\pi}\int d\Omega\frac{dN_{\rm
det}}{d\Omega dr}(r,{\bf{n}},\rho>\rho_{*})\,.$ (3)

Eq. (1) depends on the observer velocity through two effects. First, the
observed solid angle is affected by aberration:

$d\Omega=d\bar{\Omega}\left(1-2{\bf{n}}\cdot\frac{{\bf{v}}_{o}}{c}\right)\,.$
(4)

And second, the SNR threshold at the observer, $\rho_{*}$, (which is a fixed
number) corresponds to different emitted GW power for sources situated in
different directions. Since we want to relate the $\rho_{*}$ dependence of Eq.
(1) to the astrophysical distribution of GW sources and the effect of the
kinematic dipole, we factorize the SNR as follows

$\displaystyle\rho^{2}(r,{\bf{n}})=I\times P(r,{\bf{n}})\,.$ (5)

The part $I$ is an ‘intrinsic” part that accounts for the fact that the SNR
depends on the astrophysical properties of the sources (such as masses). The
“propagation” part, $P(r,{\bf{n}})$ on the other hand describes how the SNR
depends on the relative position of the source and the observer. Since
distances and redshifts are affected by the observer velocity, this
propagation part is not isotropic: it depends on ${\bf{n}}$ (more precisely it
depends on the angle between ${\bf{n}}$ and the observer velocity). We will
see later in this section how to define these quantities for GWs emitted from
inspiralling binaries.

A source sitting at position $(r,{\bf{n}})$ with an SNR above threshold,
$\rho(r,{\bf{n}})>\rho_{*}$, must have intrinsic properties
$I>I_{*}(r,{\bf{n}})$ where

$\displaystyle I_{*}(r,{\bf{n}})\equiv\frac{\rho_{*}^{2}}{P(r,{\bf{n}})}\,.$
(6)

We see that because $P(r,{\bf{n}})$ depends on ${\bf{n}}$, a fixed SNR
threshold $\rho_{*}$ corresponds to different intrinsic properties
$I_{*}(r,{\bf{n}})$ in different directions. Using that the observer velocity
is small compared to the Hubble flow, we can then Taylor expand Eq. (1) around
the homogeneous background. We obtain

$\displaystyle\frac{dN_{\rm det}}{d\Omega dr}(r,{\bf{n}},\rho>\rho_{*})$
$\displaystyle=\frac{dN_{\rm det}}{d\Omega
dr}(r,{\bf{n}},I>I_{*}(r,{\bf{n}}))$ $\displaystyle\simeq\frac{dN_{\rm
det}}{d\bar{\Omega}dr}(r,I>\bar{I}_{*}(r))\left(1+2{\bf{n}}\cdot\frac{{\bf{v}}_{o}}{c}\right)$
$\displaystyle+\frac{\partial}{\partial I_{*}}\left(\frac{dN_{\rm
det}}{d\Omega dr}(r,I>I_{*})\right)_{I_{*}=\bar{I}_{*}}\delta
I_{*}(r,{\bf{n}})\,,$ (7)

where we have used Eq. (4) in the second line.

To compute $\delta I_{*}(r,{\bf{n}})$ we need a model for the SNR. For a
binary system of compact objects, the SNR at the zero Post-Newtonian (0PN)
order is given by (Finn & Chernoff, 1993)

$\rho^{2}(r,{\bf{n}},m,\mathcal{M})=\frac{5}{96\pi^{4/3}}\frac{\Theta^{2}}{D_{L}^{2}(r,{\bf{n}})}(G\mathcal{M}_{z})^{5/3}\mathcal{F}(f^{z}_{\rm
ISCO}(m))\,,$ (8)

where $\mathcal{M}_{z}=\mathcal{M}(1+z)$ is the redshifted chirp mass of the
system, $m=m_{1}+m_{2}$ is the total mass, $\Theta^{2}$ is a geometrical
factor that depends on the inclination of the binary and on the antenna
pattern of the detector, and $D_{L}$ denotes the luminosity distance. The
$\mathcal{F}$ quantifies the sensitivity of the GW detector, namely

$\mathcal{F}(f^{z}_{\rm
ISCO}(m))\equiv\int_{0}^{2f^{z}_{\text{ISCO}}}df\left[f^{7/3}S_{n}(f)\right]^{-1}\,,$
(9)

where $S_{n}(f)$ is the detector Power Spectral Density (PSD). In Eq. (9) the
upper integration bound is given by the redshifted frequency corresponding to
the innermost stable circular orbit (ISCO) of the system, i.e. the frequency
at which we consider the inspiraling phase of the system to end in our
approximation. It is defined as (Maggiore, 2007)

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

and $f^{z}_{\text{ISCO}}=f_{\text{ISCO}}/(1+z)$. Note that at this frequency
the 0PN approximation in Eq. (8) used for the SNR estimate becomes inaccurate.

Referring to Eq. (5), we now define the intrinsic part of the SNR as the one
that depends only on the intrinsic properties of the source

$\displaystyle
I(\mathcal{M},m)\equiv\frac{5}{96\pi^{4/3}}\Theta^{2}(G\mathcal{M})^{5/3}\mathcal{F}(f_{\rm
ISCO}(m))\,,$ (11)

and the propagation part as the part that depends on the relative position of
the source and the observer, i.e. on the redshift and luminosity distance of
the source

$\displaystyle
P(r,{\bf{n}})\equiv\frac{(1+z(r,{\bf{n}}))^{5/3}}{D_{L}^{2}(r,{\bf{n}})}\frac{\mathcal{F}(f^{z}_{\rm
ISCO}(m,r,{\bf{n}}))}{\mathcal{F}(f_{\rm ISCO}(m))}\,.$ (12)

Since the redshift and the luminosity distance are affected by the observer
velocity, the propagation factor, $P$ depends directly on ${\bf{v}}_{0}$. Note
that the redshift perturbation also enters via the upper bound of the integral
$f^{z}_{\rm ISCO}$, which reflects the fact that the observer velocity shifts
the observed frequency of the ISCO. Inserting (12) into (6) and expanding at
linear order in the velocity we obtain

$\displaystyle\delta I_{*}(r,{\bf{n}})=\bar{I}_{*}(r)$
$\displaystyle\Bigg{[}2\frac{\delta
D_{L}}{\bar{D}_{L}}-\frac{5}{3}\frac{\delta z}{1+\bar{z}}$ (13)
$\displaystyle+\frac{2f^{\bar{z}}_{\rm ISCO}}{\mathcal{F}(f^{\bar{z}}_{\rm
ISCO})}\frac{\left(2f^{\bar{z}}_{\rm
ISCO}\right)^{-7/3}}{S_{n}\left(2f^{\bar{z}}_{\rm ISCO}\right)}\frac{\delta
z}{1+\bar{z}}\Bigg{]}\,.$

Using that the redshift pertubation and luminosity distance perturbations are
given by

$\frac{\delta
z}{1+\bar{z}}=-{\bf{n}}\cdot\frac{{\bf{v}}_{o}}{c}\,,\quad\mbox{and}\quad\frac{\delta
D_{L}}{\bar{D}_{L}}=-{\bf{n}}\cdot\frac{{\bf{v}}_{o}}{c}\,,$ (14)

we find

$\displaystyle\delta I_{*}(r,{\bf{n}})=-\bar{I}_{*}(r)$
$\displaystyle\Bigg{[}\frac{1}{3}+\frac{2f^{\bar{z}}_{\rm
ISCO}}{\mathcal{F}(f^{\bar{z}}_{\rm ISCO})}\frac{\left(2f^{\bar{z}}_{\rm
ISCO}\right)^{-7/3}}{S_{n}\left(2f^{\bar{z}}_{\rm
ISCO}\right)}\Bigg{]}{\bf{n}}\cdot\frac{{\bf{v}}_{o}}{c}\,.$ (15)

Finally, we need to compute the variation of the cumulative number of events
above threshold:

$\frac{\partial}{\partial I_{*}}\left(\frac{dN_{\rm det}(r,I>I_{*})}{d\Omega
dr}\right)_{I_{*}=\bar{I}_{*}}\,.$ (16)

Without loss of generality, we assume that GW events are distributed in a
window of $I$, i.e. $I\in[I_{\text{min}},I_{\text{max}}]$. The intrinsic part
of the SNR, $I$, is indeed directly related to the chirp mass of the system,
and to the total mass (through $f_{\rm ISCO}$). Both these quantities have a
given distribution with finite width. Therefore only a range of values of $I$
are physical. We can write

$\displaystyle\frac{\partial}{\partial I_{*}}\left(\frac{dN_{\rm
det}(r,I>I_{*})}{d\Omega dr}\right)_{I_{*}=\bar{I}_{*}}=$
$\displaystyle-\frac{s(I_{*})}{I_{*}}\frac{dN_{\rm det}(r,I>I_{*})}{d\Omega
dr}\,,$ (17)

where the parameter $s$ is defined through

$\displaystyle
s(I_{*})=s(\rho_{*},r)\equiv-\frac{\partial\ln\left(\frac{dN_{\rm
det}(r,I>I_{*})}{d\Omega dr}\right)}{\partial\ln I_{*}}\,.$ (18)

This parameter directly depends on the population of sources. It is non-zero
only if $I_{\rm min}<I_{*}<I_{\rm max}$, or similarly if $\rho_{\rm
min}(r)<\rho_{*}<\rho_{\rm max}(r)$, i.e. when we are dealing with a
population of sources with a significant fraction of events across threshold:

$\displaystyle
s(\rho_{*},r)\propto\Theta(\rho_{*}-\rho_{\text{min}})\Theta(\rho_{\text{max}}-\rho_{*})\,.$
(19)

Inserting Eqs. (17), (15) and (4) into (7), and subtracting the angular
average, we find for the dipole

$\displaystyle\mathcal{D}\left[{\bf n}\cdot{\bf v}_{0}\right]$
$\displaystyle={\bf{n}}\cdot\frac{{\bf{v}}_{o}}{c}\int_{0}^{\infty}drf(r)\left[2+s(r,\rho_{*})\left(\frac{1}{3}+\mathcal{A}(r)\right)\right]$
$\displaystyle\equiv\alpha\,{\bf{n}}\cdot\frac{{\bf{v}}_{o}}{c}\,,$ (20)

where we implicitly defined the parameter $\alpha$ that will be used later on
in the analysis. The function $f(r)$ denotes the radial distribution of
sources

$f(r)\equiv\frac{\frac{dN}{drd\Omega}}{\int_{0}^{\infty}dr\frac{dN}{d\Omega
dr}}\,,$ (21)

and $\mathcal{A}$ is given by the second term in Eq. (15), averaged over all
sources at distance $r$:

$\displaystyle\mathcal{A}(r)\equiv$ $\displaystyle\int
dm_{1}dm_{2}\,p(m_{1},m_{2})\frac{1}{\mathcal{F}\left(\frac{f_{\text{ISCO}}}{(1+z)}\right)}$
$\displaystyle\times\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)}\,.$
(22)

Here $p(m_{1},m_{2})$ is the probability density function (PDF) of the source-
frame masses and $\mathcal{A}$ depends on $r$ through the redshift $z=z(r)$.
We have checked that $\mathcal{A}$ is a quantity of order $\mathcal{O}(1)$.

We observe that in Eq. (20), the term proportional to $s$ is relevant only if
a significant part of the population is across threshold. This term is due to
the third line of Eq. (7). For a population of sources where $s\simeq 0$, i.e.
such that nearly all sources are detected, then $\alpha(r)\simeq 2$, i.e. it
is a fixed constant, independent on distance.

In this case, the dipole can be used as a direct estimator of the observer
velocity. More precisely, for a perfectly isotropic distribution of sources,
i.e. with PDF $P(\Omega)=1/4\pi$, we can build the following observable
(choosing $v_{0}$ aligned along the azimutal axis)

$v_{{\bf n}^{\prime}}\equiv\frac{3}{2}\int{\rm d}\Omega\ P(\Omega)\ \left({\bf
n}\cdot{\bf n}^{\prime}\right)\ \mathcal{D}\left[{\bf n}\cdot{\bf
v}_{0}\right]=\frac{\alpha
v_{0}}{2c}\cos\theta^{\prime}\simeq\frac{v_{0}}{c}\cos\theta^{\prime}\,,$ (23)

where ${\bf n}^{\prime}\equiv\left(\sin\theta^{\prime}\
\cos\phi^{\prime}\,,\sin\theta^{\prime}\
\sin\phi^{\prime}\,,\cos\theta^{\prime}\right)$ is a vector pointing towards a
generic fixed direction. This observable is maximized when evaluated along the
a priori unknown dipole direction and is exactly equal to the observer
velocity, $v_{0}/c$, at the maximum.

If the term $s$ is not negligible, it is necessary to estimate it accurately
in order not to bias the measurement of the observer velocity. One possibility
is to measure $s$ directly from the catalog of events. Since we work at linear
order in the observer velocity, $s$ depends only on the isotropic distribution
of sources. Using Eq. (18) we can therefore write at zeroth order in $v_{0}/c$

$\displaystyle s(\rho_{*},r)\equiv-\frac{\partial\ln\left(\frac{dN_{\rm
det}(r,\rho>\rho_{*})}{d\Omega dr}\right)}{2\partial\ln\rho_{*}}\,.$ (24)

To measure $s$ one can bin the events in bins of comoving distance (using a
fiducial cosmology to translate the measured luminosity distance into $r$),
and in bins of SNR, $\rho$. The quantity $s$ is then given by the slope of the
cumulative number of events above $\rho$, evaluated at the chosen value
$\rho_{*}$. This method has been used for example in Lepori et al. (2022) for
galaxy number counts, providing a measurement of $s$ from the Euclid flagship
simulation. Similarly the quantity $\mathcal{A}(r)$ can be measured from the
catalog of events, using a fiducial cosmology to compute the probability
distribution in the source frame mass from the distribution of masses in the
observer frame. This may however be non-trivial and it requires also to
account for non-trivial detector sensitivities as a function of frequency.

A complementary possibility would be to study the effect of $s$ with Monte
Carlo simulations. For instance, as we show in Sec. 4, even if we are not able
to detect all the BBHs simulated, the effect of $s$ does not seem to introduce
a bias in the estimation of ${\bf v}_{0}$. Using Monte Carlo simulations, we
could make sure that the effect of $s$ on $\alpha$ is lower than the typical
statistical uncertainties on the estimation of the dipole which should scale
as $1/\sqrt{N_{\rm tot}}$. Let us also recall that $s$ would impact the
estimation of the kinematic dipole, but not its detectability with a given
number of detections.

###  2.2 A statistical estimator for the cosmic dipole

Let us now build an estimator for the observable $v_{{\bf n}^{\prime}}$. We
divide the sky in $N_{\rm sky}$ pixels of same solid angle, and associate a
vector ${\bf{n}}_{i}$ pointing to the center of each pixel. The estimator is
given by

$\hat{v}_{{\bf n}^{\prime}}=\frac{3}{2N_{\rm tot}}\sum_{i=1}^{N_{\rm
sky}}N_{\rm det}^{i}\cdot({\bf n}_{i}\cdot{\bf n}^{\prime})\,,$ (25)

where $N_{\rm tot}$ is the total number of events and $N_{\rm det}^{i}$ is the
number of sources that falls in the sky pixel $i$. This number can be written
as

$\displaystyle N_{\rm det}^{i}=\bar{N}_{\rm
det}\left(1+\alpha\,{\bf{n}}_{i}\cdot\frac{{\bf{v}}_{o}}{c}\right)+\Delta
N^{i}\,,$ (26)

where $\bar{N}_{\rm det}=N_{\rm tot}/N_{\rm sky}$ is the mean number of events
per pixel, and $\Delta N^{i}$ accounts for the fact that the actual
distribution of detected GW signals will not be exactly isotropic (even in the
absence of a dipole), because of the stochastic nature of the coalescence
distribution and of the detection process. Inserting (26) in (25), the
expected value of the estimator becomes

$\displaystyle\langle\hat{v}_{{\bf n}^{\prime}}\rangle$
$\displaystyle=\frac{3}{2N_{\rm sky}}\alpha\frac{v_{o}}{c}\sum_{i=1}^{N_{\rm
sky}}({\bf n}_{i}\cdot\hat{{\bf{v}}}_{o})({\bf n}_{i}\cdot{\bf
n}^{\prime})=\frac{\alpha v_{o}}{2c}\cos(\theta^{\prime})=v_{\bf
n^{\prime}}\,,$ (27)

because the isotropic part of the actual sky distribution averages to zero due
to the $({\bf n}_{i}\cdot{\bf n}^{\prime})$ factor, and the expectation value
of the stochastic noise is, by definition, zero. Hence what survives is
exactly the dipole component discussed in the previous section, which reduces
to the observer velocity when $\alpha=2$, i.e. when threshold effects are
negligible.

In order to assess the statistical detectability of the dipole, one must
determine the variance of the estimator, due to the stochastic fluctuations of
events, $\Delta N^{i}$, around the monopole. Since the pixels have associated
the same solid angle, the number of detections contained in every pixel is
drawn from the same distribution, that is a Poissonian with $\mu=\bar{N}_{\rm
det}$. One can then determine the variance of $\hat{v}_{{\bf n}^{\prime}}$
just by summing in quadrature the variances within each pixel:

$\langle\Delta\hat{v}_{{\bf n}^{\prime}}^{2}\rangle=\frac{9}{4N_{\rm
tot}^{2}}\sum_{i=1}^{N_{\rm sky}}\bar{N}_{\rm
det}\cos^{2}{\theta^{\prime}}_{i}\simeq\frac{3}{4N_{\rm tot}}\,,$ (28)

as the sum in $\cos^{2}\theta_{i}$ quickly converges to $N_{\rm sky}/3$,
already for $N_{\rm sky}<10$. We thus have the intuitive result that the
precision attainable in the dipole measurement scales down with the square
root of the total number of detections, as expected in an essentially
poissonian process.

From Eq. (28), one expects that in order to detect a velocity dipole $v_{0}/c$
of order $10^{-3}$, one needs a very large number of GW detections, $N_{\rm
tot}\simeq 10^{6}$, which can be reached only with XG detectors.

Another way to obtain the same result is to artificially generate a large
number of random sky distributions, for instance by randomly reshuffling the
positions of the detected GW events, and compute the variance of the estimator
on such a collection of sky distributions. As will be shown in the next
section, the result of this estimate is in very good agreement with the simple
poissonian computation outlined above.

##  3 Forecasts with XG GW detectors

In this section, we forecast the detectability of the cosmic dipole with GW
events. In Sec. 3.1 we start by discussing the detection prospects for compact
binary mergers with XG detectors. In Sec. 3.2 we study the significance of a
possible cosmic dipole detection using the estimator in Eq. (25). Finally, in
Sec. 3.3 we show an example of how to estimate the cosmic dipole using
Bayesian statistic.

###  3.1 Simulating compact binaries mergers with XG detectors

We consider a network of XG detectors composed by ET (Punturo et al., 2010)
and two CE (Dwyer et al., 2015; Reitze et al., 2019). For ET we assume the
same PSD used in Iacovelli et al. (2022), while for CE we use the PSD provided
by the CE consortium111https://cosmicexplorer.org/sensitivity.html. In this
exploratory approach, we detect a GW event using a SNR threshold calculated
with Eq. (8)222This is a reasonable approximation for BNS events, but not so
much for BBHs; however in the latter case threshold effects are less important
(see end of this section) and this crude approximation should not affect our
result.. As we detect the time-evolution of the antenna patterns of the GW
detectors, we set a lower frequency cut-off in Eq. (9) at $5$ Hz. Given the
sensitivity curves of the network and the two values of the source-frame
masses, Eq. (8) can be used to calculate the maximum redshift at which we will
be able to observe binaries with a given SNR as threshold.

Figure 1: Detection horizons for a ET+2CE detectors network. The redshift
horizons (vertical axis) are calculated as a function of total source-frame
mass (horizontal axis) for an equal mass binary “edge-on” with respect to the
observer (worst-case scenario). The different colors mark the horizons for 3
SNR thresholds of 6, 9 and 12. Detection horizons are calculated using a flat
$\Lambda$CDM cosmology with $H_{0}=67.7\mathrm{km\,s^{-1}\,Mpc^{-1}}$ and
$\Omega_{m}=0.3097$. Note that we only consider binaries with ISCO frequency
above $5$ Hz for our simulation.

Fig. 1 reports the maximum redshift, as a function of total source frame mass,
up to which we will be able to detect compact binary coalescences with a SNR
$>6,9$ and $12$. We note that the horizons reported in Fig. 1 are calculated
for edge-on binaries ($\cos\iota=0$) and therefore they represent the worst-
case scenario at which we will be able to detect compact binaries. Fig. 1
indicates that if we consider a SNR for detection of $9$, we will be able to
detect all the BNSs merging below redshift $\sim 0.8$ and all the BBHs merging
below redshift $\sim 2$.

BBHs are more massive than BNSs, which is why we expect them to have a louder
SNR and a higher detection range (see Eq. (8)). However, when a BBH is too
massive, it merges at low frequencies (see Eq. (10)) thus spending a small
number of cycles in the sensitivity band of GW detectors and collecting less
SNR. Moreover, low-frequency regions are not very sensitive for ground-based
GW detectors (Maggiore et al., 2020; Iacovelli et al., 2022). That is why in
Fig. 1, the detection range starts to decrease after $\sim 40M_{\odot}$. This
behavior mostly depends on the GW detector’s sensitivity as a function of
frequency. Ground-based GW detectors usually have the same design for
sensitivity as a function of frequency. While we might change BBHs population
models, detector sensitivities and other simulation prescriptions, recent
works indicate that we expect to detect almost all the population of BBHs
merging in the Universe with next-generation detectors (Maggiore et al., 2020;
Iacovelli et al., 2022).

We want to determine if these detection horizons can introduce a significant
selection bias when detecting populations of BNSs and BBHs. If a strong
selection bias is present, then the theoretical prediction for the dipole
depends on the parameter $s(\rho_{*},r)$ in Eq. (17), that we need to model.
If no strong selection bias is induced, i.e. if $s(\rho_{*},r)\simeq 0$, then
$\alpha\simeq 2$ and our estimator provides a model-independent measurement of
the observer velocity $v_{o}$.

Figure 2: Top panel: Total number (vertical axis) of GWs from BNS and BBH
mergers arriving on Earth within a given redshift shell (horizontal axis). The
shaded areas mark the 90% credible intervals associated to the population
model, while the solid lines their median values. Bottom panel: Total fraction
of GWs that will be detectable (with various SNR thresholds) within a given
redshift shell. The figures are generated using the same population model for
BNS and BBHs as in Iacovelli et al. (2022). The solid lines indicate the
fiducial model for the merger rates, while the shaded area the contours
identified by the 90% credible interval uncertainties on the rate models.

We simulate BNSs and BBHs following the same population models for BNSs and
BBHs reported in Iacovelli et al. (2022) and supported by current observations
(Abbott et al., 2021b). More details about the BNSs and BBHs mass
distributions are given in App. A. We use the same merger rate model as a
function of redshift for both BNSs and BBHs, given by, see Callister et al.
(2020)

$R(z)=R_{0}[1+(1+z_{p})^{-\gamma-k}]\frac{(1+z)^{\gamma}}{1+\left(\frac{1+z}{1+z_{p}}\right)^{\gamma+k}}\,,$
(29)

where $R_{0}$ is the merger rate today, and $\gamma$ and $k$ are two
parameters encoding the redshift evolution of the merger rate. Since $k$ and
$\gamma$ are not constrained by current observations, we assume fiducial
values of $k=3$ and $z_{p}=2$, which are consistent with the Star Formation
Rate (Madau & Dickinson, 2014). For $\gamma$ and $R_{0}$ we take values
consistent with the 90% credible intervals values found by The LIGO Scientific
Collaboration et al. (2021c). For BBHs we take $R_{0,\rm
BBH}=17^{+10}_{-6.7}{\rm Gpc^{-3}yr^{-1}}$ and $\gamma=2.7^{+1.7}_{-1.8}$,
while for BNSs we take $R_{0,\rm BNS}=106^{+190}_{-93}{\rm Gpc^{-3}yr^{-1}}$
and the same $\gamma$ as for BBHs (since it was not possible to constrain this
parameter from current observations). We take as fiducial rate models the ones
corresponding to the median values of the parameters. The number of GWs,
emitted within a redshift shell with $z^{\prime}<z$, crossing the Earth per
year is then given by

$\frac{dN}{dt_{d}}=\int_{0}^{z}R(z^{\prime})\frac{1}{1+z^{\prime}}\frac{dV_{c}}{dz^{\prime}}dz^{\prime}\,,$
(30)

where $dV_{c}/dz$ is the differential of the comoving volume and $t_{d}$ is
the detector time in years. In the top panel of Fig. 2 we report the total
number of GWs from BNS and BBH mergers arriving in 1 year of observing time.
With the assumed model for the merger rate of BNSs and BBHs we find that the
total number of BNS and BBH mergers saturates at redshift $\sim 2$. This is
expected from the assumed merger rate model that has a peak around redshift
$2$. If we consider all the observable Universe, we find that we might expect
to have between $(1-50)\cdot 10^{4}$ GWs from BBHs arriving on Earth per year
and $10^{4}-10^{7}$ GWs from BNSs per year. However, given the detection
horizons in Fig. 1 not all the GWs will be detectable.

The bottom panel of Fig. 2 shows the “detectable fraction” of BNS and BBH
mergers within a certain redshift. On one hand, the plot shows that we will be
basically able to detect all the BBH mergers in the Universe. This means that
all events are above threshold, leading to $s(\rho_{*},r)\simeq 0$ at all
distances. For BNSs, on the other hand, we see that above $z\simeq 1$, a
sizeable fraction of events will not be detectable. As a consequence
$s(\rho_{*},r)$ may be large, leading to a non-negligible contribution of
threshold effects to the parameter $\alpha$. Therefore, even though we expect
significantly less detections from BBHs than from BNSs, the former will offer
a cleaner measure of the observer velocity, which is independent of the
characteristics of the population.

Using the merger rate prescription in Eq. (29) we also estimate the clustering
dipole contribution following App. A of Secrest et al. (2021). To calculate
the clustering dipole contribution, we use nbodykit (Hand et al., 2018) to
obtain the matter density power spectrum at $z=0$ using the Halofit (Takahashi
et al., 2012) prescription to include non-linearities in the CLASS code
(Lesgourgues, 2011) for a Planck15 cosmology. By using the power spectrum at
$z=0$ we are overestimating the clustering dipole (since the power spectrum
decreases with redshift) and we therefore obtain a conservative estimation for
its importance with respect to the kinematic dipole. We obtain a value of the
clustering dipole of $\mathcal{D}_{\rm cl}\approx 1.2\cdot 10^{-4}$ which is
about one order of magnitude lower than the simulated cosmic kinematic dipole
for AGN and CMB.

Figure 3: Sky distribution of the estimator $\hat{v}_{\mathbf{n}}$ for the CMB
fiducial value (top map) and the AGN fiducial value (bottom map). The maps are
centered around the injected direction of the observer velocity. The figure is
generated by dividing the sky in equal size pixels of 53 deg2.

###  3.2 Detection and Estimation of the cosmic dipole

In this section we study the detectability of the GW dipole, considering two
fiducial values for its amplitude. The first one assumes that the GW dipole is
of purely kinematic nature, with a velocity consistent with the one inferred
from CMB observations (Planck Collaboration et al., 2020), i.e. $v^{\rm
CMB}_{\hat{{\bf{v}}}_{o}}=v_{o,{\rm CMB}}/c=1.2\cdot 10^{-3}$. Since
observations from radio sources and quasars find a dipole which is 2-5 times
larger than expected (Colin et al., 2017; Bengaly et al., 2018; Secrest et
al., 2021; Siewert et al., 2021; Secrest et al., 2022), we also consider a
second fiducial value for the GW dipole, that would be consistent with these
observations. Indeed, if the large dipole found in these studies is due to a
breaking of the Copernician principle (Secrest et al., 2022), i.e. to an
instrinsic large anisotropy in the distribution of structures, then this
anisotropy should be present also in the distribution of GW sources, that
follow the large-scale structure of the Universe. It is thus interesting to
assess the detectability of such a large dipole with GW sources. For this we
take the extreme case of a dipole which would be 5 times larger than the one
expected from CMB velocity, and with a direction aligned with it, i.e. $v^{\rm
AGN}_{\hat{{\bf{v}}}_{o}}=6.0\cdot 10^{-3}$.

For each fiducial value of the dipole, we simulate $10^{4},10^{5},10^{6}$ and
$10^{7}$ BBHs detections $N_{\rm tot}$, using the ET+2CE network with a SNR
thresold of $9$ for detection and the population of BBHs described in the
previous section. The BBH detections are simulated as follow. Each simulated
BBH mass and redshift are drawn from the distributions described in the
previous section and in App. A. The original distribution of BBHs is isotropic
in sky. Then, for each BBHs, we introduce the effect of the observer velocity
by introducing an aberration following Eq. (4):
$\theta^{\prime}=\theta-(v_{o}/c)\sin\theta$, where $\theta$ is the angle
between the source position and the observer velocity 333Note that in the AGN
case, the dipole may not be due solely to the observer velocity, but rather to
a large anisotropy in the distribution of sources, as discussed above.
However, when assessing the detectability of such a dipole, it does not matter
if it is of kinematic origin or not. Hence, we can simply simulate it as if it
were due to a large velocity, 5 times larger than the one measured from the
CMB.. The detector-frame mass and luminosity distance change as described in
Eq. (14). For each BBH, we then calculate the SNR using Eq. (8). Then an
“observed” SNR is drawn from a $\chi^{2}$ distribution with 2 times number of
detectors d.o.f. If the “observed” SNR exceeds a SNR threshold of $9$, the
binary is labelled as detected. Finally, to mimic the sky localization
uncertainty given by the GW detection, we scatter the position of the GW
sources using a gaussian distribution with $3$ deg.

Once the total list of BBHs detections is obtained with their sky position, we
calculate the estimator of the GW dipole defined in Eq. (25). Fig. 3 shows the
skymap of $\hat{v}_{\mathbf{n}}$ calculated for the CMB and AGN fiducial
values using $10^{6}$ BBHs detections. As we can see from the figure, in the
AGN case, $\hat{v}_{\mathbf{n}}$ is maximized in the direction of the observer
velocity and it displays an amplitude of $\hat{v}_{\mathbf{n}}=5.73\cdot
10^{-3}$, similar to the one injected. For the CMB case on the other hand, the
estimator is maximized in a direction different from the one of the observer
velocity and it displays a maximum value of $\hat{v}_{\mathbf{n}}=1.3\cdot
10^{-3}$, higher than the one injected. The reason for this is that, in the
CMB case, the GW dipole will only be marginally detectable with $10^{6}$ BBHs.

Figure 4: Distribution of the estimator $\hat{v}_{\mathbf{n}}$ (red points)
for the CMB fiducial value (top plot) and the AGN fiducial value (bottom
plot), plotted as a function of the angle between ${\bf{n}}$ and the injected
velocity direction $\hat{{\bf{v}}}_{o}$. The shaded areas mark the background
$1,2,3\sigma$ confidence intervals generated by randomly shuffling the
detections over the sky. The horizontal dashed lines mark the standard
deviations generated with Eq. (28). The figure is generated by dividing the
sky equal size pixels of 53 deg2.

Fig. 4 shows another view of the values of the estimator reported in Fig. 3.
More precisely, we plot the value of the estimator as a function of the angle
between the direction of the observer velocity, $\hat{{\bf{v}}}_{o}$, and the
direction $\mathbf{n}$ at which the estimator is calculated. The figure also
displays the $1,2$ and $3$ $\sigma$ values of $\hat{v}_{\mathbf{n}}$ due to
the stochastic fluctuations of events around the monopole given by Eq. (28). A
$3\sigma$ detection of the GW dipole can be claimed if the
$\hat{v}_{\mathbf{n}}$ estimator exceeds the $3\sigma$ threshold from Poisson
noise. From the figure we can see that, with $10^{6}$ detections, a GW dipole
with amplitude compatible with that of the CMB dipole is marginally detectable
with $1\sigma$ confidence (the probability that it is generated from a
fluctuation of the monopole is not negligible), while a GW dipole with
amplitude compatible with that of the AGN dipole is clearly detectable.

Figure 5: Distribution of the estimator $\hat{v}_{\mathbf{n}}$ evaluated in
the injected velocity direction ${\bf{n}}=\hat{{\bf{v}}}_{o}$ with
$10^{4},10^{5},10^{6}$ and $10^{7}$ BBH detections. The histograms are
obtained by simulating 1000 populations of BBHs. The vertical dashed lines
indicate the fiducial values of the GW dipole in the CMB case (blue) and the
AGN case (orange). The dashed and dotted gray lines indicate the $1,2,3\sigma$
contribution from Poisson noise, generated using Eq. (28).

In order to better understand the detectability of the GW dipole, we repeat
1000 times the previous simulations for $10^{4},10^{5},10^{6}$ and $10^{7}$
detected BBHs and we calculate the fraction of cases for which the dipole will
be detected. Fig. 5 shows the distribution of the estimator
$\hat{v}_{\mathbf{n}}$, evaluated in the velocity direction
${\bf{n}}=\hat{{\bf{v}}}_{o}$, that we obtain for the CMB and AGN cases. From
the figure we can see that the distributions are centered around the injected
value of the dipole. This is a confirmation that, for BBHs, threshold effects
are negligible and $\alpha=2$ (such that $\alpha/2=1$ in Eq. (27)). From the
figure, we can see that when we only have $10^{4}$ detections, the
distribution of the dipole estimator for the AGN and the CMB cases almost
coincide. Moreover, the $1,2,3\sigma$ intervals of these distributions almost
coincide with the $1,2,3\sigma$ variance from Poisson noise. In other words,
with only $10^{4}$ detections the dipole is not detectable. As we collect more
and more detections, the dipole estimator distributions tend to become
separate from the Poisson noise.

Figure 6: FAP versus detection probability plots for the injected CMB and AGN
fiducial values (blue and orange lines). The vertical dashed lines indicate
the FAP at the standard $1,2,3\sigma$ credible intervals. The black dashed
line is the reference to indicate when the detection probability coincides
with the FAP (they are the same distribution).

In order to better quantify the detectability of the dipole, we plot False
Alarm probability (FAP) vs detection probability in Fig. 6. The FAP identifies
a threshold for the dipole detection and it is defined as the probability that
a random fluctuation of the number of GW detections in absence of dipole,
would result in a false positive. The detection probability is defined as the
probability that, in presence of a dipole, the estimator for the dipole
detection would exceed the FAP threshold. Decreasing the FAP allows us to be
more sure on the dipole nature of our detection but it decreases our
sensitivity for the dipole detection. FAP versus detection probabilities plots
can be used to: (i) estimate what are the detection prospects for a given
threshold (significance of the detection) and (ii) clearly check in what
regime the dipole is detectable. In fact, in the case where we are not able to
detect the dipole, we expect the detection probability to be equal to the FAP,
as the detection estimator would follow the same distribution.

As we can see from the top left panel of Fig. 5, with $10^{4}$ GWs detections,
the distribution of the estimator for the CMB case agrees with the confidence
intervals traced by the detection thresholds. This is reflected in the top
left panel of Fig. 6, where we see that the FAP and detection probability
follow the same distribution, thus indicating that the only dipole that can be
detected in this case, is the one from random fluctuations around the
isotropic background (false positive). For the AGN case, the detection
probability is slightly larger than the FAP, but still no robust detection
could be claimed in this case. With $10^{5}$ GWs detections ($\sim$ 2-3 year
of observation), we can see from Fig. 6 that there is a 75% probability of
detecting a dipole with AGN amplitude using a FAP of $2\sigma$. With $10^{6}$
detections (achievable in 10 years for an optimistic scenario), a dipole as
high as the AGN one would be detectable at 100%. On the other hand, if the GW
dipole is compatible with the one observed in the CMB, then there is a 50%
probability that we would detect it with a FAP of $2\sigma$.

###  3.3 A Bayesian study

Bayesian statistic can also be used to provide evidence for the GW dipole and
estimate its parameters. As discussed in Section 2.1, the number of observed
BBHs in a pixel situated in direction ${\bf{n}}_{i}$ is given by

$\displaystyle N_{\rm det}^{i}=\bar{N}_{\rm
det}\left(1+\alpha\,{\bf{n}}_{i}\cdot\frac{{\bf{v}}_{o}}{c}\right)+\Delta
N^{i}\,,$ (31)

where $\bar{N}_{\rm det}$ is the number of detections per pixel due to the
monopole, $\Delta N^{i}$ a statistical fluctuation , and $\alpha\approx 2$
since for BBHs threshold effects are negligible. The likelihood of obtaining
$k_{i}$ detections in one pixel ${\bf{n}}_{i}$ is then given by

$\mathcal{L}\left(k_{i}|\bar{N}_{\rm
det},\frac{{\bf{v}}_{o}\alpha}{c},{\bf{n}}_{i}\right)\propto e^{-N^{i}_{\rm
det}}\cdot\left(N^{i}_{\rm det}\right)^{k_{i}}\,.$ (32)

The overall likelihood of obtaining $\\{k\\}=\\{k_{1},\ldots,k_{N_{\rm
sky}}\\}$ detections when dividing the sky in $N_{\rm sky}$ equal area pixels
is

$\mathcal{L}\left(\\{k\\}|\bar{N}_{\rm
det},\frac{{\bf{v}}_{o}\alpha}{c}\right)=\prod_{i}^{N_{\rm
sky}}\mathcal{L}\left(k_{i}|\bar{N}_{\rm
det},\frac{{\bf{v}}_{o}\alpha}{c},{\bf{n}}_{i}\right)\,.$ (33)

Finally, by applying the Bayes theorem, we can obtain posterior distributions
on $\bar{N}_{\rm det},\alpha v_{o}/c$ and $\hat{{\bf{v}}}_{o}$ by calculating

$\displaystyle p\left(\bar{N}_{\rm
det},\frac{{\bf{v}}_{o}\alpha}{c}|\\{k\\}\right)$
$\displaystyle\propto\mathcal{L}\left(\\{k\\}|\bar{N}_{\rm
det},\frac{{\bf{v}}_{o}\alpha}{c},\right)\times\pi\left(\bar{N}_{\rm
det},\frac{{\bf{v}}_{o}\alpha}{c},\right)\,,$ (34)

where $\pi(\cdot)$ is a prior term.

We implement the likelihood in Eq. (33) in a nested sampling code to obtain
posterior distributions for $N_{\rm tot}=\bar{N}_{\rm det}N_{\rm sky}$ and
$\alpha v_{o}/c$, in the case studies of the CMB and AGN dipole amplitudes,
estimated with $10^{6}$ GWs. We use the python code bilby (Ashton et al.,
2019; Romero-Shaw et al., 2020) and its implementations of the nested sampling
algorithm dynesty (Higson et al., 2019). We use an isotropic prior for the
dipole direction, a flat in log prior for $\alpha v_{o}/c$ between $10^{-6}$
and $10^{-1}$ and a flat in log prior for the total number of events $N_{\rm
tot}$ between $10^{5}$ and $10^{8}$.

Figure 7: Sky area at 90% credible intervals for the dipole sky direction for
the AGN test case with $10^{6}$ detections. The red dot marks the direction of
the dipole. The sky localization has a radius of $\sim$ 13 deg.

For the case of a GW dipole with amplitude compatible with the AGN one, we
find that with $10^{6}$ GWs detections, we are able to estimate the direction
of the observer velocity with an uncertainty of $\sim 13$ deg ($\sim 2200$
deg2 of area) at 90% credible intervals. Fig. 7 shows the 90% credible
intervals in the sky identified by the posterior.

Figure 8: Plots of the posterior distribution (and its marginals) for the
total number of detections $N_{\rm tot}$ and the dipole amplitude $\alpha
v_{o}/c$ for the AGN test case with $10^{6}$ detections. The solid black lines
mark the injected values (assuming $\alpha=2$).

Fig. 8 shows instead the estimate of the total number of detections and the
value of the GW dipole amplitude, $\alpha v_{o}/c$. The value of the recovered
dipole amplitude is $\alpha v_{o}/c=1.15^{+0.18}_{-0.18}\cdot 10^{-2}$ at
68.3% symmetric credible intervals. The amplitude of a dipole consistent with
the AGN one can therefore be measured with a precision of 16%.

We can also calculate the Bayes factor to asses the detection of the dipole,
i.e.

$\mathcal{B}^{\rm dip}_{\rm mono}=\frac{p(\rm dip|\\{k\\})}{p(\rm
mono|\\{k\\})},$ (35)

where $p(\rm dip|\\{k\\})$ and $p(\rm mono|\\{k\\})$ are the evidences for the
dipole model (with $\alpha v_{o}/c\neq 0$) and the monopole model (with
$\alpha v_{o}/c=0$). For the AGN case, with $10^{6}$ GW detections we obtain a
$\log_{10}(\mathcal{B}^{\rm dip}_{\rm mono})=6.0$, so strong preference for
the presence of a dipole.

On the other hand, if the GW dipole has an amplitude consistent with that of
the CMB dipole, the situation is different. First, we find that with $10^{6}$
detections it will not be possible to constrain the sky location. This is
expected, since in Sec. 3.2 we showed that the GW dipole would be hardly
detectable in this case with $10^{6}$ events.

Figure 9: Plots of the posterior distribution (and its marginals) for the
total number of detections $N_{\rm tot}$ and the dipole amplitude $\alpha
v_{o}/c$ for the CMB test case with $10^{6}$ detections. The solid black lines
mark the injected values (assuming $\alpha=2$).

Moreover, if we look at the posteriors on the total number of events, $N_{\rm
tot}$, and on the dipole amplitude, $\alpha v_{o}/c$, plotted in Fig. 9, we
can see that the total number of detections is clearly constrained, while the
dipole amplitude cannot be constrained. However, it is possible to define an
upper-limit in this case that results in $\alpha v_{o}/c<3.0\cdot 10^{-3}$ at
95% credible intervals (the injected value was $2.4\cdot 10^{-3}$). We also
compute the Bayes factor for this case finding that
$\log_{10}(\mathcal{B}^{\rm dip}_{\rm mono})=-0.17$, indicating that there is
no clear preference for the presence of a dipole or not.

##  4 Discussion

In the previous section we have discussed several aspects related to the
detection and estimation of the cosmic dipole from GW counting using both
frequentists and bayesian techniques. We have shown that BBHs detected with
the ET+2CE might be a “clean” probe for estimation of the dipole, since
threshold effects are negligible in this case and the dipole is consequently
only due to aberration ($\alpha\approx 2$).

The first crucial task to consider when searching for a dipole with GW number
counts is the assessment of the detection significance. For the frequentists
techniques, we have shown that the significance can be evaluated using
standard p-value techniques for Poissonian statistic and by reshuffling GW
events in the sky to build a distribution of the noise due to the stochastic
distribution of sources. For Bayesian statistic instead, we have shown that
the detection can be evaluated with Bayes factors between the dipole and the
monopole models. In both cases, we have obtained that in order to detect a
dipole with amplitude $\alpha v_{o}/c$, one would need at least $\sim(c/\alpha
v_{o})^{2}$ detections, in order to be significantly confident that the
detected dipole is not a fluctuation of the monopole distribution. We stress
that, even with a few detections, it would not be surprising to find that the
distribution of events is not perfectly isotropic in the sky. However, a low
number of detections does not have the statistical significance to fit a
dipole distribution. For instance, in Fig. 3, we have shown that for a
simulated GW dipole consistent with the one of the CMB, our estimator shows a
peak over the sky which is not directly aligned with the simulated dipole.
However in this case, its value is not statistically significant and a
detection cannot be claimed.

This could explain the puzzling results of Kashyap et al. (2022), where using
LIGO/Virgo data the authors reconstruct a dipole direction orthogonal to the
CMB one. In that work the authors indeed do find an asymmetry in the
distribution of GW events over the sky, however they do not assign any
statistical significance to their results. If a significance was worked out,
it would very probably indicate that the dipole reconstruction is
statistically not significant due to the very low statistics of events used in
the analysis.

Figure 10: 95% confidence level upper-limits (blue dots) on the dipole
amplitude generated from populations of BBHs isotropically distributed over
the sky (no cosmic dipole present). The black dashed line indicates a fit for
the upper-limit scaling. The error bars are generated by repeating each
simulation 50 times.

To better demonstrate that the dipole detection capabilities scales as
$1/\sqrt{N_{\rm tot}}$, in Fig. 10, we provide upper limits on the dipole
amplitude obtained from BBHs distributions that are generated isotropically
over the sky. By using $10^{4},10^{5},10^{6}$ and $10^{7}$ BBH detections
simulated in the case where no cosmic dipole is present, we estimate the 95%
confidence level upper-limit on the amplitude of the cosmic dipole. The best
fit for the scaling of the 95% confidence level upper limit
$\mathrm{UL}_{95\%}$ on ${\bf\frac{v_{0}\alpha}{c}}$ is given them by

$\log_{10}[\mathrm{UL}_{95\%}]=a\log_{10}[N_{\rm tot}]+b,$ (36)

with $a=-0.47\pm 0.06$ and $b=0.4\pm 0.4$ therefore including the typical
scaling $\propto N_{\rm tot}^{-1/2}.$

From the figure, we can see that the 95% upper limit on the dipole amplitude
scales as expected with the number of detections.

In this study, we have focused on the dipole in GW number counting, showing in
particular that BBHs provide a very clean way of measuring the observer
velocity, since the theoretical prediction for the signal is very simple in
this case. This is in contrast with the dipole from radio galaxies and
quasars, for which the amplitude of the dipole depends directly on the
properties of the sources, namely their spectral index and flux distribution.
As shown in Dalang & Bonvin (2022), if the populations evolve with redshift,
then the theoretical prediction may change significantly and the tension with
the CMB dipole may disappear. In this context, BBHs will provide a robust way
to determine if the tension between the CMB dipole and the AGN dipole is due
to our imperfect knowledge of quasars and radio galaxies properties, to
systematic (that will necessarily be different with GW observations) or to a
breaking of the Copernician principle.

This being said, it is still interesting to investigate if other GW estimators
could help explaining the tension between the CMB and AGN dipole. In Kashyap
et al. (2022), the authors propose to measure the dipole in the distribution
over the sky of the detector-frame mass. Similarly to the case of number
counting, detecting a cosmic dipole with this method, would require the
detection of a relative discrepancy in detector-frame mass distribution of the
order of $v_{o}/c$. To determine if this method can work, it is therefore
crucial to assess the noise expected on such a measurement. Indeed, in
addition to Poisson noise that affects number counting, the width of the mass
distribution would induce additional contributions to the variance of the
estimator. For BBHs, for which the width is expected to be large, the variance
may be important 444Note that the mass distribution of BBHs is directly
estimated from observed data, with complex hierarchical analysis tools such as
the ones of Abbott et al. (2019, 2021b); The LIGO Scientific Collaboration et
al. (2021c). For BNSs, on the other hand, the variance may be significantly
smaller since the mass distribution is much more peaked. It may also be
interesting to combine mass and number counting, to determine if threshold
effects can be mitigated with a specific estimator. We will explore all these
aspects in a future work.

Finally, let us comment on some limitations of our prospect studies, that can
be tackled in future works. In our study, we used the 0PN approximation to
calculate the SNR of GW events. This is an approximation that underestimates
the SNR of the detections as it neglects the merging part, which is important
for BBHs. Therefore, we might expect threshold effects to be even less
important than what we discussed. Another crucial assumption that we made was
to neglect the sky direction dependence of the Antenna patterns of GW
detectors. Due to their geometrical configuration, GW detectors are not
equally sensitive to all the directions, and even for a detector network, the
antenna patterns averaged during an year of observation would not be
isotropic. For instance, for the simulated ET+2CE network, the antenna
patterns averaged over one year of observation is

$\langle F^{2}\rangle\simeq-0.097\cos^{4}\delta-0.040\cos^{2}\delta+0.878\,,$
(37)

where $\delta$ is the source declination. As ${\rm SNR}^{2}\propto\langle
F^{2}\rangle$, and $\langle F^{2}\rangle$ varies by about $10\%$ over the sky,
we expect the SNR to vary by about $3\%$ over the sky. This means that there
might be more or less sources detected over the sky due to this sensitivity
variation, that may contaminate the measurement of a dipole of the order of
$10^{-3}$. However, if threshold effects are negligible, as for the case of
ET+2CE and BBHs, this sky-dependent selection bias will be negligible. In
other cases, this effect can be calculated form first principles and used to
adjust the estimator. Finally, one should consider also the sky-localization
uncertainties associated to the GW detection, which could be higher than $3$
deg (although the results reported in Iacovelli et al. (2022) suggest this is
a fair approximation for most BBH events); this would require an analysis done
as in Essick et al. (2022).

##  5 Conclusions

In this paper we have discussed aspects and prospects for detecting and
estimating the cosmic dipole due to the observer velocity using GWs detected
with XG detectors.

In Sec. 2 we have introduced the theoretical framework to evaluate the effect
of the observer velocity on GW detections. We have shown that this velocity
introduces an aberration on the GW localization, and that it modifies the
number of detections above a given SNR threshold through the redshifted chirp
mass and the luminosity distance. For BBHs, we have demonstrated that these
threshold effects are negligible and that the amplitude of the GW dipole is
directly given by $2v_{o}/c$.

In Sec. 3 we have discussed several frequentists and bayesian techniques to
detect and estimate the presence of the cosmic dipole from GW counting. We
have shown that, with $10^{6}$ BBH detections, which would be observable in a
few years of observations with a ET+2CE network, it will be possible to detect
a cosmic dipole with amplitude similar to the one estimated from AGN, with a
precision of $\sim 16\%$. On the other hand, with $10^{6}$ detections, a GW
dipole with amplitude compatible with that of the CMB would only be marginally
detectable. With $10^{7}$ detections however, we would be able to
significantly detect a GW dipole with amplitude compatible with both the AGN
and the CMB one. If we include BNSs and model the threshold contributions to
the dipole, $10^{7}$ detections could be reachable in $\sim 10$ years of
observations of ET+2CE.

Finally, in Sec. 4 we have discussed critical aspects related to the detection
of the cosmic dipole using number counting and the detector-frame mass
distribution. Moreover, we have discussed the impact of some of the
assumptions that we made in this exploratory study.

Next generation GW detectors, that will give us access to at least hundreds of
thousands of GW detections per year, have therefore the potential to solve a
growing tension associated to the standard cosmological model.

## Acknowledgements

We thank Nicola Tamanini and Archisman Ghosh for discussions and exchanges.
S.M. is supported by the ANR COSMERGE project, grant ANR-20-CE31-001 of the
French Agence Nationale de la Recherche. S.F. is supported by the Fonds
National Suisse, grant $200020\\_191957$, and by the SwissMap National Center
for Competence in Research. C.B. acknowledges funding from the Swiss National
Science Foundation and from the European Research Council (ERC) under the
European Union’s Horizon 2020 research and innovation program (Grant agreement
No. 863929; project title “Testing the law of gravity with novel large-scale
structure observables"). The work of G.C. is supported by the CNRS and by
Swiss National Science Foundation (Ambizione grant, “Gravitational wave
propagation in the clustered universe").

## Data Availability

The simulations and numerical code underlying this paper are available on
GitHub .

## References

  * Abbott et al. (2016) Abbott B. P., et al., 2016, Phys. Rev. Lett., 116, 061102
  * Abbott et al. (2018) Abbott B. P., et al., 2018, LRR, 21, 3
  * Abbott et al. (2019) Abbott B. P., et al., 2019, ApJL, 882, L24
  * Abbott et al. (2021a) Abbott B. P., et al., 2021a, ApJ, 909, 218
  * Abbott et al. (2021b) Abbott R., et al., 2021b, ApJL, 913, L7
  * Amaro-Seoane et al. (2017) Amaro-Seoane P., et al., 2017, arXiv, p. arXiv:1702.00786
  * Ashton et al. (2019) Ashton G., et al., 2019, ApJS, 241, 27
  * Barausse et al. (2014) Barausse E., Cardoso V., Pani P., 2014, Phys. Rev. D, 89, 104059
  * Bartolo et al. (2022) Bartolo N., et al., 2022, arXiv, p. arXiv:2201.08782
  * Bengaly et al. (2018) Bengaly C. A. P., Maartens R., Santos M. G., 2018, JCAP, 2018, 031
  * Bonvin et al. (2017) Bonvin C., Caprini C., Sturani R., Tamanini N., 2017, Phys. Rev. D, 95, 044029
  * Bonvin et al. (2022) Bonvin C., Cusin G., Pitrou C., Mastrogiovanni S., Congedo G., Gair J., 2022
  * Cai et al. (2018) Cai R.-G., Liu T.-B., Liu X.-W., Wang S.-J., Yang T., 2018, PhRvD, 97, 103005
  * Callister et al. (2020) Callister T., Fishbach M., Holz D. E., Farr W. M., 2020, ApJL, 896, L32
  * Colin et al. (2017) Colin J., Mohayaee R., Rameez M., Sarkar S., 2017, MNRAS, 471, 1045
  * Cusin & Tamanini (2021) Cusin G., Tamanini N., 2021, Mon. Not. Roy. Astron. Soc., 504, 3610
  * Cusin & Tasinato (2022) Cusin G., Tasinato G., 2022, JCAP, 2022, 036
  * Cusin et al. (2017) Cusin G., Pitrou C., Uzan J.-P., 2017, Phys. Rev. D, 96, 103019
  * Cusin et al. (2021) Cusin G., Durrer R., Dvorkin I., 2021, Universe, 8, 19
  * Dalang & Bonvin (2022) Dalang C., Bonvin C., 2022, Mon. Not. Roy. Astron. Soc., 512, 3895
  * Dwyer et al. (2015) Dwyer S., Sigg D., Ballmer S. W., Barsotti L., Mavalvala N., Evans M., 2015, Phys. Rev. D, 91, 082001
  * Essick et al. (2022) Essick R., Farr W. M., Fishbach M., Holz D. E., Katsavounidis E., 2022, arXiv, p. arXiv:2207.05792
  * Finn & Chernoff (1993) Finn L. S., Chernoff D. F., 1993, Phys. Rev., D47, 2198
  * Hand et al. (2018) Hand N., Feng Y., Beutler F., Li Y., Modi C., Seljak U., Slepian Z., 2018, Astron. J., 156, 160
  * Higson et al. (2019) Higson E., Handley W., Hobson M., Lasenby A., 2019, S&C, 29, 891
  * Iacovelli et al. (2022) Iacovelli F., Mancarella M., Foffa S., Maggiore M., 2022, arXiv, p. arXiv:2207.02771
  * Kashyap et al. (2022) Kashyap G., Singh N. K., Phukon K. S., Caudill S., Jain P., 2022, arXiv, p. arXiv:2204.07472
  * Lepori et al. (2022) Lepori F., et al., 2022, Astron. Astrophys., 662, A93
  * Lesgourgues (2011) Lesgourgues J., 2011, arXiv, p. arXiv:1104.2932
  * Leyde et al. (2022) Leyde K., Mastrogiovanni S., Steer D. A., Chassande-Mottin E., Karathanasis C., 2022, arXiv, p. arXiv:2202.00025
  * Madau & Dickinson (2014) Madau P., Dickinson M., 2014, ARA&A, 52, 415
  * Maggiore (2007) Maggiore M., 2007, Gravitational Waves. Vol. 1: Theory and Experiments. Oxford Master Series in Physics, Oxford University Press, http://www.oup.com/uk/catalogue/?ci=9780198570745
  * Maggiore et al. (2020) Maggiore M., et al., 2020, JCAP, 2020, 050
  * Mancarella et al. (2022) Mancarella M., Genoud-Prachex E., Maggiore M., 2022, PhRvD, 105, 064030
  * Pitrou et al. (2020) Pitrou C., Cusin G., Uzan J.-P., 2020, Phys. Rev. D, 101, 081301
  * Planck Collaboration et al. (2020) Planck Collaboration et al., 2020, A&A, 641, A6
  * Punturo et al. (2010) Punturo M., et al., 2010, CQGra, 27, 194002
  * Reitze et al. (2019) Reitze D., et al., 2019, Bull. Am. Astron. Soc., 51, 035
  * Romero-Shaw et al. (2020) Romero-Shaw I. M., et al., 2020, MNRAS, 499, 3295
  * Sberna et al. (2022) Sberna L., et al., 2022, Phys. Rev. D, 106, 064056
  * Secrest et al. (2021) Secrest N. J., von Hausegger S., Rameez M., Mohayaee R., Sarkar S., Colin J., 2021, ApJL, 908, L51
  * Secrest et al. (2022) Secrest N., von Hausegger S., Rameez M., Mohayaee R., Sarkar S., 2022, arXiv, p. arXiv:2206.05624
  * Siewert et al. (2021) Siewert T. M., Schmidt-Rubart M., Schwarz D. J., 2021, A&A, 653, A9
  * Stiskalek et al. (2021) Stiskalek R., Veitch J., Messenger C., 2021, MNRAS, 501, 970
  * Takahashi & Nakamura (2003) Takahashi R., Nakamura T., 2003, Astrophys. J., 595, 1039
  * Takahashi et al. (2012) Takahashi R., Sato M., Nishimichi T., Taruya A., Oguri M., 2012, Astrophys. J., 761, 152
  * The LIGO Scientific Collaboration et al. (2021a) The LIGO Scientific Collaboration et al., 2021a, Constraints on the cosmic expansion history from GWTC-3, doi:10.48550/ARXIV.2111.03604, https://arxiv.org/abs/2111.03604
  * The LIGO Scientific Collaboration et al. (2021b) The LIGO Scientific Collaboration et al., 2021b, arXiv, p. arXiv:2111.03606
  * The LIGO Scientific Collaboration et al. (2021c) The LIGO Scientific Collaboration et al., 2021c, arXiv, p. arXiv:2111.03634
  * Toscani et al. (2023) Toscani M., Rossi E. M., Tamanini N., Cusin G., 2023
  * Toubiana et al. (2021) Toubiana A., et al., 2021, Phys. Rev. Lett., 126, 101105
  * Valbusa Dall’Armi et al. (2022) Valbusa Dall’Armi L., Ricciardone A., Bertacca D., 2022, arXiv, p. arXiv:2206.02747

##  Appendix A Mass population model

We use the Power Law+peak model from Abbott et al. (2021b) to describe the
source frame distribution of BBHs masses. This model is composed by two
statistical distributions: a truncated power law distribution

$\mathcal{P}(x|x_{\rm min},x_{\rm
max},\alpha)\propto\begin{cases}x^{\alpha}&\left(x_{\rm min}\leqslant
x\leqslant x_{\rm max}\right)\\\ 0&\mathrm{Otherwise}.\end{cases}$ (38)

and Gaussian distribution with mean $\mu$ and standard deviation $\sigma$.

$\mathcal{G}(x|\mu,\sigma)=\frac{1}{\sigma\sqrt{2\pi}}\exp{\left[-\frac{(x-\mu)^{2}}{2\sigma^{2}}\right]}\,.$
(39)

The distribution of the source frame masses $m_{1},m_{2}$ is factorized as

$\pi(m_{1},m_{2}|\Phi_{m})=\pi(m_{1}|\Phi_{m})\pi(m_{2}|m_{1},\Phi_{m}),$ (40)

where $\pi(m_{1}|\Phi_{m})$ is the Power Law+Peak model and
$\pi(m_{2}|m_{1},\Phi_{m})$ is a truncated power law distribution. The
secondary mass is conditioned to the constraint $m_{2}<m_{1}$, i.e.

$\pi(m_{2}|m_{1},m_{\rm min},\alpha)=\mathcal{P}(m_{2}|m_{\rm
min},m_{1},\beta)\,.$ (41)

The Power Law+Peak distribution is

$\displaystyle\pi(m_{1}|m_{\rm min},m_{\rm max},\alpha,\lambda_{\rm
g},\mu_{\rm g},\sigma_{\rm g})=$ $\displaystyle(1-\lambda_{\rm
g})\mathcal{P}(m_{1}|m_{\rm min},m_{\rm max},-\alpha)$
$\displaystyle+\lambda_{\rm g}\mathcal{G}(m_{1}|\mu_{\rm g},\sigma_{\rm g})\,$
(42)

where the power law part has slope $-\alpha$ between $m_{\rm min}$ and $m_{\rm
max}$, while the gaussian component has mean $\mu_{\rm g}$ and standard
deviation $\sigma_{\rm g}$ and accounts $\lambda_{g}$ total fraction of the
distribution. We also apply an additional smoothing at the lower edge of the
distribution

$\displaystyle\pi(m_{1},m_{2}|\Phi_{m})$ $\displaystyle=$
$\displaystyle[\pi(m_{1}|\Phi_{m})S(m_{1}|\delta_{m},m_{\rm min})]\times$
$\displaystyle[\pi(m_{2}|m_{1},\Phi_{m})S(m_{2}|\delta_{m},m_{\rm min})]$

where $S$ is a sigmoid-like window function as described in Abbott et al.
(2021b).

The parameters used for the simulation are $\alpha=3.4,\beta=1.1,m_{\rm
max}=87M_{\odot},m_{\rm
min}=5.1M_{\rm\odot},\delta_{m}=4.8M_{\odot},\sigma_{g}=3.6M_{\odot},\mu_{g}=34M_{\odot},\lambda_{g}=0.03$.