# Anisotropy of the astrophysical gravitational wave background: Analytic expression of the angular power spectrum and correlation with cosmological observations

Giulia Cusin, Cyril Pitrou, Jean-Philippe Uzan

► To cite this version:

Giulia Cusin, Cyril Pitrou, Jean-Philippe Uzan. Anisotropy of the astrophysical gravitational wave background: Analytic expression of the angular power spectrum and correlation with cosmological observations. Physical Review D, 2017, 96 (10), pp.103019. 10.1103/PhysRevD.96.103019. hal-01669773

HAL Id: hal-01669773

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

Submitted on 5 Oct 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Anisotropy of the astrophysical gravitational wave background: Analytic expression of the angular power spectrum and correlation with cosmological observations

Giulia Cusin*

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

Cyril Pitrou† and Jean-Philippe Uzan‡

Institut d'Astrophysique de Paris, CNRS UMR 7095, Université Pierre & Marie Curie—Paris VI, 98 bis Bd Arago, 75014 Paris, France and Sorbonne Universités, Institut Lagrange de Paris, 98 bis, Bd Arago, 75014 Paris, France

(Received 21 June 2017; published 28 November 2017)

Unresolved sources of gravitational waves are at the origin of a stochastic gravitational wave background. While the computation of its mean density as a function of frequency in a homogeneous and isotropic universe is standard lore, the computation of its anisotropies requires understanding the coarse graining from local systems, to galactic scales and then to cosmology. An expression of the gravitational wave energy density valid in any general spacetime is derived. It is then specialized to a perturbed Friedmann-Lemaître spacetime in order to determine the angular power spectrum of this stochastic background as well as its correlation with other cosmological probes, such as the galaxy number counts and weak lensing. Our result for the angular power spectrum also provides an expression for the variance of the gravitational wave background.

DOI: 10.1103/PhysRevD.96.103019

## I. INTRODUCTION

Astronomical observations include the detection of radiation from resolved objects as well as various stochastic backgrounds of radiations (electromagnetic, gravitational, ...), due to the superposition of the signals from all unresolved sources. The electromagnetic backgrounds include the cosmic microwave background with a black body spectrum at 2.725 K and the extragalactic background, made up of all the light emitted by stars, galaxies, quasars etc. since their formation. The analysis of these backgrounds provides information on the dynamics of our universe and on the distribution and evolution of its large scale structure. Similarly, there exists a neutrino background, the observation of which seems to be out of reach, and a background of gravity waves (GW).

The recent detection by the Advanced Laser Interferometric Gravitational-wave Observatory (LIGO) of the gravitational wave sources GW150914 [1] provided the first observation of the merging of a binary black hole (BBH) system. This first detection has been followed by other two observations of similar systems, i.e. GW151226 [2] and GW170104 [3], and by the very recent observation of a black hole merging from both the LIGO and VIRGO detectors [4]. Following these observations, the rate and mass of coalescing

binary black holes appear to be greater than many previous expectations. As a result, the stochastic background from unresolved compact binary coalescences is expected to be particularly loud. In [5] a search for the isotropic stochastic gravitational wave background is performed using data from Advanced LIGO's first observing run. The dimensionless energy density of gravitational waves is constrained to be $\Omega_{GW} < 1.7 \times 10^{-7}$ with 95% confidence, assuming a flat energy density spectrum in the most sensitive part of the LIGO band (20–86 Hz). This is a factor of ~33 times more sensitive than previous measurements [5,6], and this improves bounds on the stochastic background obtained from the analysis of big-bang nucleosynthesis [7,8] and of the cosmic microwave background [9,10] at 100 Hz. At low frequencies, Pulsar Timing Array (see below) gives a bound $\Omega_{GW} < 1.3 \times 10^{-9}$ for $f = 2.8 \times 10^{-9}$ Hz [11].

The possibility of measuring and mapping the gravitational wave background is discussed in Refs. [12–17] and a discussion of the different methods which can be used by LIGO and Laser Interferometer Space Antenna (LISA) to reconstruct an angular resolved map of the sky are presented in Ref. [18]. An analogous discussion for pulsar timing arrays can be found in Refs. [19–21].

Stochastic GW backgrounds are made up of the superposition of astrophysical signals from unresolved sources. The different backgrounds contribute at different frequencies and have distinct statistical properties, making it potentially possible to distinguish them [22]. In the standard cosmological model [23], the existence of a primordial GW

<footnote>* giulia.cusin@unige.ch</footnote>

<footnote>†pitrou@iap.fr</footnote>

<footnote>‡uzan@iap.fr</footnote>
background from the amplification of vacuum quantum fluctuations is a generic prediction of any inflationary phase. Gravitational waves may also be produced at the end of inflation during the reheating phase, see e.g. Ref. [24] for an analytic and numerical study of the process. More speculative sources of an early GW background include pre-big bang modes, cosmic strings [25–29], first-order phase transitions in the early universe [30,31], and magnetic fields [32]; see Ref. [33] for a review on those topics and Refs. [7,34] for general reviews. In addition, an astrophysical contribution (AGWB) results from the superposition of a large number of unresolved sources since the beginning of stellar activity. The nature of the AGWB may differ from its cosmological counterpart, expected to be stationary, unpolarized, statistically Gaussian and isotropic, by analogy with the cosmic microwave background. Many different sources may contribute to the AGWB, including black holes and neutron star mergers [35–41], supermassive black holes [42], exploding supernovae (SNe) [43], neutron stars [44–46], stellar core collapse [47], and population III binaries [48].

The observational landscape is also growing and covers large bands of frequencies; see e.g. Ref. [49] for a review.<fnref n="1" /> At extremely low frequencies ~$10^{-16}$ Hz bounds come mainly from the analysis of CMB B-modes while at low frequency of order $10^{-10} – 10^{-6}$ Hz, there are pulsar timing arrays such as the radio telescope Parks Pulsar Timing Array<fnref n="2" /> (PPTA), the Large European Array for Pulsar Timing<fnref n="3" /> (LEPTA), and the future International Pulsar Timing Array<fnref n="4" /> (IPTA). At low frequencies (typically $10^{-6} – 10^0$ Hz) detection relies on space-borne detectors, such as LISA<fnref n="5" /> and the evolved Laser Interferometer Space Antenna<fnref n="6" /> (eLISA) launched in 2016. High frequency (typically $10^0 – 10^5$ Hz) observations rely on ground-based detectors, such as LIGO and its advanced configuration (aLIGO), VIRGO,<fnref n="7" /> the Einstein Telescope<fnref n="8" /> (ET) or its american counterpart Cosmic Explorer (CE) [50]. This spectrum covers most of the theoretical predictions.

So far, constraints have been set on the total energy density of the GW background integrated over the sky, as a function of frequency, $\Omega_{GW}(\nu)$, even though some constraints on the anisotropy have been obtained by PTA [51,52], (see below for definition) and by the advanced LIGO first observing run [5]. On the theoretical side, the energy density of GW has been modeled and parametrized under the assumption that both our universe and the distribution of sources are

homogeneous and isotropic, see e.g. Refs. [22,38]. These assumptions can be relaxed in order to take into account that astrophysical sources are located in cosmic structures that indeed have a distribution that can be computed in a given cosmological model. Therefore, the flux of energy from all astrophysical sources (resolved and unresolved) is not constant across the sky and depends on the direction of observation. The goal of our work is to present an analytic framework to describe and compute the anisotropy in the observed energy density of the AGWB, taking into account the presence of inhomogeneities in the matter distribution and geometry in our observed universe.

To that goal, we need to distinguish different scales in our problem, represented in Fig. 1. The first and largest scale is the cosmological scale of the observer for whom galaxies are pointlike GW sources. The total flux of GW he or she receives is the sum of the fluxes of all individual galaxies. The flux from a given galaxy is related to the galaxy effective luminosity of gravitational waves. This is a physical quantity defined in the galaxy rest frame that takes into account contributions from the various sources inside the galaxy (galactic scale). Each of these sources has a peculiar motion with respect to the galaxy rest frame, inducing an extra dispersion in the signal. The luminosity of an astrophysical source in the galaxy rest frame can in turn be expressed in terms of local quantities characterizing the source in its rest frame (local astrophysical scale). In such a way, one can express cosmological observables characterizing GW at the observer scale in terms of quantities defined on local scales of single GW sources. In other words, we shall provide a final parameterization for the observed GW flux, which inherits two main contributions: the first has a cosmological origin, i.e. it is related to the specific cosmology we are looking at (i.e. geometry of the universe, distribution of the large scale structure). The second comes from local physics and encodes information on the specific processes of GW emission. This latter contribution depends on the details of the emission processes and on the intrinsic properties of GW sources inside a galaxy (e.g. rates and densities of the different kinds of processes).

It follows that the anisotropies of the AGWB encode information about (1) the different mechanisms of GW production, (2) the astrophysical distribution of GW sources in galaxies, (3) galaxy formation and distribution of the large scale structures of the universe, and (4) the spacetime geometry along the line of sight. The question is which of these effects dominates on different frequency bands. In particular, if local effects are dominant, the AGWB will teach us mostly about the distribution of GW sources, while if it is the cosmological contribution that dominates, we shall have a new cosmological observable robust to local physics unknown. The answer to this question will be discussed in our companion article [53], where we also analyze the shape and amplitude of the AGWB angular power spectra for different frequencies together with the correlation with other

<footnote n="1">¹The associated code http://rhcole.com/apps/GWplotter/ allows one to generate plots of noise curves for many detectors and associated target sources.</footnote>

<footnote n="2">²http://www.atnf.csiro.au/research/pulsar/ppta/</footnote>

<footnote n="3">³http://www.leap.eu.org</footnote>

<footnote n="4">⁴http://www.ipta4gw.org</footnote>

<footnote n="5">⁵www.lisamission.org</footnote>

<footnote n="6">⁶https://www.elisascience.org</footnote>

<footnote n="7">⁷https://www.ego-gw.it/public/about/whatIs.aspx</footnote>

<footnote n="8">⁸http://www.et-gw.eu</footnote>
<figure><img src="image_2.png" /><figcaption>FIG. 1. Schematic representation of our approach. Left: Spacetime representation. The observer detects a signal in the direction $e_o^\mu$, with solid angle $d^2\Omega_o$ around this direction. This defines a bundle of null geodesics along its past light cone. The physical volume element $dV$ is a 3-dimensional volume element, defined as the intersection of a 4-dimensional volume element with the observer past light cone. Right: The three scales in our problem: cosmological, galactic, and astrophysical.</figcaption></figure>

cosmological probes. Our theoretical framework can be straightforwardly generalized to discuss different cosmological scenarios: for example, if black hole binary systems are of stellar origin, then the AGWB is expected to be correlated with the galaxy distribution, while if dark matter is related to a 20–100 solar mass primordial BH component [54–56], the correlation between AGWB and the galaxy distribution will be smaller; see e.g. Ref. [57], where the possibility of distinguishing these models is discussed. The correlation of resolved BH binaries with galaxy clustering and weak lensing is also discussed in Ref. [58], which in particular discusses the possibility for second generation GW detectors to measure such correlations.

Our approach relies on integrating all sources. Alternatively, one can also derive an effective source term for a Boltzmann approach to the AGWB, as tried in Ref. [59]. For resolved sources, our analysis provides a modelization of the GW luminosity distance [60–62], see also Ref. [63] for the cross-correlation of these GW standard sirens to galaxy clustering. On the other side, the total energy density of the AGWB depends only on the ratio between the luminosity distance $D_L$ and the angular diameter distance $D_A$. This ratio is given by $(1+z)^2$ because of the reciprocity relation, which holds in any spacetime geometry [64,65].

This study is organized as follows. In Sec. II we introduce the line of sight approach, we define the various scales of the

problem and derive a parametrization for the energy density of the AGWB valid in a generic cosmology. We also describe the coarse graining necessary to link observable quantities to local quantities characterizing the processes of GW emission. In Sec. III we specialize this result to the case of a perturbed Friedmann-Lemaître-Robertson-Walker (FLRW) cosmology. This provides all the tools to compute the angular power spectrum of the AGWB energy density in Sec. IV. To conclude, in Secs. V and VI we compute the correlators with other cosmological probes, namely weak lensing by the large scale structure and galaxy number counts. The applicability of this formalism is discussed in Sec. VII, where we pave the way to the further numerical investigations presented in our companion article [53].

**Notations:** The cosmic scale factor is normalized to unity at the time of observation, $a_o = 1$, and we use units defined such that $c = \hbar = k_B = 1$. A dot corresponds to a derivative with respect to conformal time, i.e. $\dot{X} = \partial_\eta X$.

## II. GENERAL FRAMEWORK

The spectrum of the AGWB is characterized by the dimensionless density parameter

$$ \Omega_{\text{GW}}(\nu_o) = \frac{\nu_o \, d\rho_{\text{GW}}(\nu_o)}{\rho_c \, d\nu_o}, \qquad (1) $$
where $\rho_c \equiv 3H_0^2/(8\pi G)$ is the critical density, $\nu_o$ is the frequency as measured by the observer, and $\rho_{GW}$ is the AGWB energy density [66]. The corresponding spectrum $d\rho_{GW}/d\nu_o$ is related to the directional dependent spectrum, whose general derivation is the goal of this section, through

$$ \frac{d\rho_{GW}(\nu_o)}{d\nu_o} = \int \frac{d^3\rho_{GW}(\nu_o, e_o)}{d\nu_o d^2\Omega_o} d^2\Omega_o, \qquad (2) $$

where $e_o$ is a spatial 3-vector defining the direction of observation [see Eq. (6) below for definition].

## A. Line of sight approach to GW propagation

We work in the eikonal limit in which the propagation of GW is described in a similar way as light in the geometrical optic approximation [66–69]. The spacetime metric can be decomposed as

$$ (g_{\mu\nu})_{\text{tot}} = g_{\mu\nu} + h_{\mu\nu}, \qquad (3) $$

where $h_{\mu\nu}$ is the GW metric perturbation and we refer to $g_{\mu\nu}$ as the *general background metric* (e.g. in cosmology it will contain the background FLRW metric as well as cosmological perturbations).

Einstein equations imply that GWs follow null geodesics [66,67] $x^\mu(\lambda)$, where $\lambda$ is an affine parameter along the geodesic. Its tangent vector

$$ k^\mu = \frac{dx^\mu}{d\lambda} \qquad (4) $$

is a null vector ($k^\mu k_\mu = 0$) that satisfies the geodesic equation

$$ \frac{Dk^\mu}{D\lambda} = k^\nu \nabla_\nu k^\mu = 0, \qquad (5) $$

where $\nabla_\nu$ denotes the covariant derivative associated to the metric $g_{\mu\nu}$.

Consider an observer with 4-velocity $u^\mu$ ($u_\mu u^\mu = -1$). At any time, his worldline is the origin of the observer past light cone containing all observed GW rays, see Fig. 1. The 4-velocity $u^\mu$ defines a preferred spatial section and the spatial direction of GW propagation, defined as the opposite of the direction of propagation of the signal converging to the observer. It is spanned by the purely spatial unit vector $e^\mu$,

$$ e^\mu u_\mu = 0, \quad e^\mu e_\mu = 1, \qquad (6) $$

which provides the 3 + 1 decomposition of the wave 4-vector<fnref n="9" />

$$ k^\mu = E(u^\mu - e^\mu), \qquad (7) $$

where $E = 2\pi\nu = -u_\mu k^\mu$ is the cyclic frequency of the GW signal in the observer’s rest frame. We define the spatial projection of the wave 4-vector as

$$ p^\mu \equiv (g^{\mu\nu} + u^\mu u^\nu)k_\nu = -Ee^\mu. \qquad (8) $$

The redshift $z_G$ is defined from the ratio between the emitted frequency $\nu_G$ in the source’s rest frame and the observed frequency in the observer’s rest frame $\nu_o$ as

$$ 1 + z_G \equiv \frac{\nu_G}{\nu_o} = \frac{u_G^\mu k_\mu(\lambda_G)}{u_o^\mu k_\mu(\lambda_o)}, \qquad (9) $$

where $u_G^\mu$ is the 4-velocity of the source and $u_o^\mu$ is the 4-velocity of the observer. The source<fnref n="10" /> G located at a redshift $z_G$ is emitting GW with a given frequency spectrum. From the definition (9), it follows that the frequency measured in $O$, $\nu_o$ is related to the frequency at the emission $\nu_G$ by

$$ \nu_G = (1 + z_G)\nu_o. \qquad (10) $$

Let us conclude by stressing that $x^\mu(\lambda)$ is the worldline of a graviton which intersects the worldline of the observer at the time of observation. Then, by definition,

$$ x^\mu(\lambda_o) = x_o^\mu, \quad \left. \frac{dx^\mu(\lambda)}{d\lambda} \right|_{\lambda=\lambda_o} = E_o(u_o^\mu - e_o^\mu). \qquad (11) $$

Therefore, $x^\mu(\lambda)$ is a function of the direction of observation and of the 4-position of the observer, i.e. $x^\mu(\lambda) = x^\mu(\lambda, e_o^\mu, x_o^\mu)$. In the following, to make the notation compact, the dependence on $e_o^\mu$ and $x_o^\mu$ will be understood.

## B. Parametrization of the observed GW radiation

As discussed in the introduction, providing a general parametrization for the GW energy flux from unresolved astrophysical sources requires to consider 3 scales.

(1) **Cosmological scale.** The observer measures a GW flux in a direction $e_o$ and solid angle $\Omega_o$. The angular resolution of the observer is such that we assume galaxies to be pointlike sources emitting GWs and comoving with the cosmic flow.

(2) **Galactic scale.** A galaxy is described by a set of parameters such as its mass, mean metallicity, .... We associate with each galaxy an effective luminosity encoding the global effects of all the GW sources contained into it. These sources have a velocity $\Gamma$ in the galaxy rest frame. Therefore, to get the effective luminosity of each type of astrophysical source, we need to sum the luminosities of

<footnote n="9">⁹In other terms, introducing a space projector and a time projector, $S_{\mu\nu} = g_{\mu\nu} + u_\mu u_\nu$ and $T_{\mu\nu} = -u_\mu u_\nu$, respectively, we have $e^\mu = -S^\mu_\nu k^\nu / E$.</footnote>

<footnote n="10">¹⁰As discussed in the introduction, we identify sources with galaxies containing astrophysical sources of GW.</footnote>
all single sources of that type, averaged over their velocity distribution function in the galaxy rest frame.

(3) Astrophysical scale. This is the local scale of single GW sources. For any source (binary inspiraling, binary merging, supernova, ...) one can study the astrophysical process of GW production and calculate its energy spectrum and function of local parameters characterizing the source.

We denote as an i-type source a specific source of GWs. Each type of source is fully characterized by specific values for its local parameters $\theta^{(i)}$ (e.g. orbital parameters, masses, ...) and by the parameters characterizing its host galaxy (mass $M_G$, metallicity $z_G$, ...), which we will indicate collectively as $\theta_G$. These parameters are in full generality time-dependent functions. We consider $\theta_G$ and $\theta^{(i)}$ evaluated at the time at which the signal received by the observer O at the time of observation $t_o$ has been emitted.

Let us now turn to the definitions of the physical quantities on each scale and on their connections.

## 1. Cosmological scale

Since the intrinsic angular resolution of present observations does not allow us to resolve galaxies, we describe them as pointlike sources emitting GWs and comoving with the cosmic flow. Each galaxy is characterized by its effective luminosity, which depends on its parameters $\theta_G$ and that needs to be expressed in terms of the properties of the sources contained in the galaxy.

Consider a galaxy observed at a redshift $z_G$ in direction $e_o$. The total flux of energy received by the observer O in this direction is related to the absolute effective luminosity of the galaxy (i.e. the power radiated by the galaxy, in the galaxy rest frame) by the following relation:

$$ \Phi(z_G, e_o, \theta_G) \equiv \frac{L_G(\theta_G)}{4\pi D_L^2(z_G, e_o)}, \quad (12) $$

where $D_L(z_G, e_o)$ is the luminosity distance of the galaxy. The angular diameter distance is defined as

$$ D_A^2 \equiv \frac{A_G}{\Omega_o}, \quad (13) $$

where $A_G$ is the physical area of the galaxy and $\Omega_o \ll 1$ is its apparent observed angular size. Luminosity and angular diameter distances are related by the reciprocity relation (see e.g. [64,65,70,71] for a test of this relation)

$$ D_L = (1 + z_G)^2 D_A. \quad (14) $$

Let us introduce the effective frequency spectrum of GWs of the galaxy, $\mathcal{L}_G(\nu_G, \theta_G)$, normalized in such a way that

$$ \int_0^\infty \mathcal{L}_G(\nu_G, \theta_G) d\nu_G = \mathcal{L}_G(\theta_G). \quad (15) $$

Making use of Eq. (10), we find

$$ \mathcal{L}_G(\theta_G) = (1 + z_G) \int \mathcal{L}_G[\nu_o(1 + z_G), \theta_G] d\nu_o. \quad (16) $$

We can then define the flux measured by the observer in the frequency range ($\nu_o, \nu_o + d\nu_o$). Using the definition (12), it follows that

$$ \Phi(z_G, e_o, \nu_o, \theta_G) d\nu_o \equiv \frac{(1 + z_G)}{4\pi D_L^2(z_G, e_o)} \mathcal{L}_G(\nu_G, \theta_G) d\nu_o. \quad (17) $$

This quantity is called the specific flux of GW radiation.

## 2. Galactic scales

We now need to relate the former quantities to those associated to the astrophysical GW sources inside each given galaxy. In the galaxy rest frame, one can introduce a local system of coordinates (different from the cosmological coordinate system). In this system, a GW source characterized by $\theta^{(i)}$ has a velocity $\Gamma(\theta^{(i)}, \theta_G)$ that also depends on the parameters of the host galaxy (in particular, the velocity dispersion depends on the total mass).

We assume the velocities to have a distribution function $f(\Gamma, \theta_G)$ where the explicit dependence on $\theta_G$ reminds us that different galaxies have, in principle, associated different distribution functions. The function $f$ is normalized as

$$ \int d^3\Gamma f(\Gamma, \theta_G) = 1, \quad (18) $$

and we assume that each source type is globally at rest with respect to the galaxy, i.e.

$$ \int d^3\Gamma \Gamma f(\Gamma, \theta_G) = 0. \quad (19) $$

To describe the collective emission of all the sources of a galaxy, one should first perform a Lorentz transformation to shift from the reference frame of the source to the one of the galaxy (described e.g. by the normal vector to the galactic plane) to describe the emission of a single source in the galactic rest frame. Then, one shall perform a rotation to describe the emission in galactic coordinates. These transformations are described in the Appendix. However, these effects are linear in $\Gamma$ at lowest order. When averaging over all the sources of a given galaxy, this lowest order contribution will cancel thanks to Eq. (19) and higher contributions will be negligible since we expect $|\Gamma|/c \ll 1$.

We denote as $A_G^{(i)}$ the axis characterizing the orientation of the source in this galactic frame. In the source rest frame,
we expect the GW emission to be rotationally invariant around the axis of the source $A_S^{(i)}$, i.e. the spectrum depends on the direction $e_S^{(i)}$ only through $\mu_S^{(i)} = \cos \alpha_S^{(i)} = e_S \cdot A_S^{(i)}$. As explained in detail in the Appendix, the luminosity emitted by a given GW source in the galaxy rest frame, $L_G^{(i)}$, is a function of the velocity of the source and of its orientation in the galaxy rest frame, $A_G^{(i)}$, and of the line of sight in the galaxy frame $e_G$. Explicitly, the effective luminosity per unit solid angle $d^2\Omega_G$ (defined in the galaxy rest frame) is

$$ \frac{d^2 L_G^{(i)}}{d^2 \Omega_G} = \frac{d^2 L_G^{(i)}}{d^2 \Omega_G} (\Gamma, \theta^{(i)}, e_G, A_G^{(i)}, \theta_G), \quad (20) $$

where the subscript G indicates again that we are defining quantities in the galaxy rest frame.

In full generality, the computation of the total flux in the galaxy rest frame requires both the information on the orientation $A_G^{(i)}$ of each system and its velocity $\Gamma(\theta^{(i)}, \theta_G)$ in the galaxy rest frame. In the following, we assume that $A_G^{(i)}$ and $\Gamma(\theta^{(i)}, \theta_G)$ are uncorrelated random variables on which we average. Physically, the average on $A_G^{(i)}$ corresponds to an average over the possible orientations of a given source. This integration is justified by the fact that we do not have access to the information on the distribution of sources inside a galaxy. All this leads to the fact that the effective luminosity of the galaxy is isotropic in its rest frame, justifying our notation $\mathcal{L}_G(\theta_G)^{11}$. Under this hypothesis, each galaxy is a point source radiating isotropically in its rest frame.

The GW sources can be split in two categories, depending on the duration of the emission process with respect to the characteristic time of evolution of the system. We consider the following:

(a) *Inspiraling binary sources*, for which the evolution of the orbital parameters is slow compared to galactic time scales so that the emission of GWs can be averaged on many periods. To compute the luminosity associated with inspiraling systems, we need to compute the number of systems and the radiated power.

(b) *Merging binary sources and exploding SNe*, for which the duration of the emission is short compared to the typical galactic time scale. To compute the luminosity associated with these systems, we need to determine the rate of events and the total radiated energy.

We shall thus decompose the total effective luminosity of a galaxy as

$$ \mathcal{L}_G(\theta_G) = L_G^I(\theta_G) + L_G^M(\theta_G) + L_G^{SN}(\theta_G), \quad (21) $$

respectively indicating the contributions from inspiraling binaries (I), mergers (M), and supernovae (SN). As we shall see, the last two contributions have formally the same structure. However, to keep the notation clear, we write them separately.

The first contribution can be decomposed as

$$ L_G^I(\theta_G) = \int_0^\infty d\nu_G \mathcal{L}_G^I(\nu_G, \theta_G), \quad (22) $$

$$ \begin{aligned} & = \sum_{(i)} \int d\theta^{(i)} \mathcal{N}^{(i)}(\theta^{(i)}, \theta_G) \int d^3\Gamma f(\Gamma, \theta_G) \\ & \times \int \frac{d^2\Omega_G}{4\pi} d^2A_G^{(i)} \frac{d^2L_G^{(i)}}{d^2\Omega_G}(\Gamma, \theta^{(i)}, e_G, A_G^{(i)}, \theta_G), \end{aligned} \quad (23) $$

where the sum runs over different types of inspiraling binary systems (e.g. back hole–black hole, neutron star–neutron star, ...). In Eq. (23), $\mathcal{N}^{(i)}(\theta^{(i)}, \theta_G)$ is the number of sources in the inspiraling phase with parameters $\theta^{(i)}$ in a galaxy with parameters $\theta_G$, and $L_G^{(i)}$ is the effective luminosity of a single source inside the galaxy. Note that in full generality we should have considered $\mathcal{N}^{(i)}(\theta^{(i)}, \theta_G, A_G^{(i)})$, but we assumed the distribution of the axis to be isotropic. It follows that the dependence on the direction of emission drops out in the integration over $d^2\Omega_G$ so that the final result, $L_G^I(\theta_G)$, reduces to a monopole for the whole galaxy.

To grasp the physical meaning of the single source effective luminosity (per unit of solid angle $d^2\Omega_G$) and to understand how this can be related to the power radiated by the source, we express it as

$$ \frac{d^2 L_G^{(i)}}{d^2 \Omega_G} (\Gamma, \theta^{(i)}, e_G, A_G^{(i)}, \theta_G) = \int d\nu_G \frac{d^2 \mathcal{L}_G^{(i)}}{d^2 \Omega_G} (\nu_G, \Gamma, \theta^{(i)}, e_G, A_G^{(i)}, \theta_G) \quad (24) $$

$$ = \int d\nu_G \frac{d^4 E_G^{(i)}}{d^2 \Omega_G dt_G d\nu_G} (\nu_G, \Gamma, \theta^{(i)}, e_G, A_G^{(i)}, \theta_G), \quad (25) $$

where the power per unit of solid angle $dE_G^{(i)}/(dt_Gd\nu_Gd^2\Omega_G)$ in the galaxy rest frame is related to the emitted power per unit solid angle (in the source local frame) through a rotation and Lorentz transformation. Inserting Eq. (24) into Eq. (23), we can read the following parametrization for the effective luminosity per unit of effective frequency

<footnote n="11">¹¹Note that if this is not the case, one would also need to know the orientation of each galaxy with respect to the line of sight.</footnote>
$$ \mathcal{L}_G^I(\nu_G, \theta_G) = \sum_{(i)}^I \int d\theta^{(i)} \mathcal{N}^{(i)}(\theta^{(i)}, \theta_G) \int d^3\Gamma f(\Gamma, \theta_G) \\ \times \int \frac{d^2\Omega_G}{4\pi} d^2\mathbf{A}_G^{(i)} \frac{d^4 E_G^{(i)}}{d^2\Omega_G dt_G d\nu_G} \\ \times (\nu_G, \Gamma, \theta^{(i)}, \mathbf{e}_G, \mathbf{A}_G^{(i)}, \theta_G). \qquad (26) $$

The second contribution to Eq. (21) from mergers can be decomposed as

$$ L_G^M(\theta_G) = \int_0^\infty d\nu_G \mathcal{L}_G^M(\nu_G, \theta_G), \qquad (27) \\ = \sum_{(i)}^M \int d\theta^{(i)} \frac{d\mathcal{N}^{(i)}}{dt_G}(\theta^{(i)}, \theta_G) \int d^3\Gamma f(\Gamma, \theta_G) \\ \times \int d\nu_G \frac{d^2\Omega_G}{4\pi} d^2\mathbf{A}_G^{(i)} \frac{d^3 E_G^{(i)}}{d^2\Omega_G d\nu_G} \\ \times (\nu_G, \Gamma, \theta^{(i)}, \mathbf{e}_G, \mathbf{A}_G^{(i)}, \theta_G), \qquad (28) $$

where the sum runs over all types of mergers (i.e. back hole–black hole, neutron star–neutron star, ...) and $d\mathcal{N}^{(i)}/dt_G$ is the merging rate of systems of type i in the galaxy rest frame. $dE_G^{(i)}/d\nu_G$ is the energy spectrum per unit solid angle (in the galaxy rest frame) from a source of type i with parameters $\theta^{(i)}$ and velocity $\Gamma$ in the galaxy rest frame. The latter is related to the emitted spectrum in the source local frame through a Lorentz transformation.<fnref n="12" /> Comparing Eqs. (27) and (28), we can read the following expression for the effective luminosity associated to merging objects, per unit of effective frequency:

$$ \mathcal{L}_G^M(\nu_G, \theta_G) = \sum_{(i)}^M \int d\theta^{(i)} \frac{d\mathcal{N}^{(i)}}{dt_G}(\theta^{(i)}, \theta_G) \int d^3\Gamma f(\Gamma, \theta_G) \\ \times \int \frac{d^2\Omega_G}{4\pi} d^2\mathbf{A}_G^{(i)} \frac{d^3E_G^{(i)}}{d^2\Omega_G d\nu_G}(\nu_G, \Gamma, \theta^{(i)}, \mathbf{e}_G, \mathbf{A}_G^{(i)}, \theta_G). \qquad (29) $$

The decomposition of the third contribution to Eq. (21) from SNe follows the same steps described above for the case of merging binary systems, and the final result is

$$ \frac{d\mathcal{N}^{(i)}}{dt_G}(\theta^{(i)}, \theta_G) = \text{SFR} \times \mathcal{N}^{(i)}(\theta^{(i)}, \theta_G), $$

where the star formation rate (SFR) function depends on the details of local physics.

analogous to Eq. (29), with the sum running over the different types of SNe and with the quantity $d\mathcal{N}^{(i)}/dt_G$ denoting the explosion rate of SNe of type i in the galaxy rest frame.

Equations (26) for inspirals and (29) for merger and supernovae are the main results of this subsection.

### 3. Astrophysical scale

For each galaxy and each type of astrophysical source, one can calculate the GW spectrum or power emitted by the source in its local frame as a function of the source local parameters. Local quantities in the source frame, that we denote with a subscript S, are related to quantities in the galaxy rest frame through a Lorentz transformation, once the velocity of the source with respect to the galaxy rest frame is known. This is detailed in the Appendix. The effect of the velocity of the center of mass of a binary on its wavefront is discussed in Ref. [72].

At linear order in the velocity $\Gamma$, the effects of the peculiar motion of a source in its host galaxy can be neglected on average. Explicitly, the effective luminosity of a galaxy per unit of effective frequency, Eqs. (26) and (29), can be rewritten in terms of local quantities at the source since at linear order in $\Gamma$, $\nu_G = \nu_S$, $E_G = E_S$, and $t_G = t_S$, and the directional dependence reduces to a dependence in $\mu_G^{(i)} = \mathbf{e}_G \cdot \mathbf{A}_G^{(i)}$ if in the source rest frame the dependence on directions is only through $\mu_S^{(i)}$ (see Appendix).

### 4. Summary: From astrophysical to cosmological scale

We can now gather all these pieces to derive the expression of the flux of energy received by the observer O in the direction $e_o$, per unit of observed frequency in terms of quantities defined at the small scales (galactic and astrophysical). The definition (17) gives the specific flux

$$ \Phi(z_G, e_o, \nu_o, \theta_G) d\nu_o = \frac{(1+z_G)\mathcal{L}_G(\nu_G, \theta_G)}{4\pi D_L^2(z_G, e_o)} d\nu_o, \qquad (30) $$

where the total effective luminosity per unit of effective frequency $\mathcal{L}_G$ can be decomposed as the sum of contributions from inspiraling binaries, mergers, and SNe as

$$ \mathcal{L}_G \equiv \mathcal{L}_G^I + \mathcal{L}_G^M + \mathcal{L}_G^{SN}, \qquad (31) $$

where $\mathcal{L}_G^{I,M,SN}$ are respectively defined in Eqs. (26) and (29), and can be expressed in terms of local quantities defined in the local frame of the GW sources inside the galaxy.

### C. Final results

Equation (30) provides a general expression for the specific flux (i.e. flux per unit of observed frequency) of the gravitational waves received by an observer O in a direction

<footnote n="12">¹²The rate of mergers can be further parametrized as</footnote>
$e_o$ from a galaxy at redshift $z_G$, assumed to be a pointlike source and isotropically radiating in its rest frame.

In practice, because of limited spatial resolution, an instrument responds to the flux per solid angle, i.e. the intensity of GWs from the source.<fnref n="13" /> Our goal is to derive the total energy density of GWs received in a solid angle $d^2\Omega_o$ around the direction of observation $e_o$ per unit of observed frequency $\nu_o$. This requires adding the specific fluxes received from the various galaxies contained in this solid angle.

Formally, we need to consider the null geodesic $x^\mu(\lambda)$ passing through $O$ in the direction $e_o$, i.e. satisfying the conditions (11). At any event corresponding to a given value of $\lambda$ of the affine parameter, the local physical 3-volume element is

$$d^3V \equiv \sqrt{-g} \epsilon_{\mu\nu\alpha\beta} u^\mu dx^\nu dx^\alpha dx^\beta. \quad (32)$$

The total intensity of GW observed in the direction $e_o$ is then easily obtained by summing the contributions from all the galaxies along the line of sight. Explicitly, if we want the contribution between $\lambda$ and $\lambda + d\lambda$ we need to multiply the number of galaxies in $d^3V$ by the specific flux received from a single galaxy with parameters $\theta_G$, integrating over the range of possible values for $\theta_G$. One gets the final expression

$$\begin{aligned} \frac{d^3 \rho_{\text{GW}}}{d\nu_o d^2 \Omega_o} (\nu_o, e_o) \\ = \int d\lambda \int d\theta_G \Phi[x^\mu(\lambda), \nu_o, \theta_G] \frac{d^3 \mathcal{N}_G}{d\lambda d^2 \Omega_o} [x^\mu(\lambda), \theta_G], \end{aligned} \quad (33)$$

where $\Phi[x^\mu(\lambda), \nu_o, \theta_G]$ is the specific flux received at frequency $\nu_o$ in direction $e_o$ from a galaxy located at $x^\mu(\lambda)$, defined in Eq. (30), and the direction-dependent spectrum on the left-hand side is related to the GW energy density in Eq. (2). $d^3\mathcal{N}_G(x^\mu(\lambda), \theta_G)$ represents the number of galaxies with parameters $\theta_G$ contained in the physical volume $d^3V$ defined in Eq. (32), seen by the observer $O$ under the solid angle $d^2\Omega_o$. We introduce $n_G[x^\mu(\lambda), \theta_G]$, the physical density of galaxies with parameters $\theta_G$, defined by

$$d^3\mathcal{N}_G[x^\mu(\lambda), \theta_G] \equiv n_G[x^\mu(\lambda), \theta_G]d^3V[x^\mu(\lambda)]. \quad (34)$$

To simplify our final result, it is useful to rewrite Eq. (32) expressing the physical volume element as

$$d^3V[x^\mu(\lambda)] = d^2\Omega_o D_A^2(\lambda) \sqrt{p_\mu(\lambda) p^\mu(\lambda)} d\lambda, \quad (35)$$

where we used the fact that $d^3V$ is the volume with cross section $D_A^2$ and depth $\sqrt{p_\mu(\lambda)p^\mu(\lambda)}d\lambda = -(u_\mu k^\mu)d\lambda$ along the line of sight, $p^\mu$ being defined in Eq. (8).

Substituting Eqs. (34)–(35) into Eq. (33) and using the expression (30) for the specific flux, we get

$$\begin{aligned} \frac{d^3 \rho_{\text{GW}}}{d\nu_o d^2 \Omega_o} (\nu_o, e_o) \\ = \frac{1}{4\pi} \int d\lambda \int d\theta_G \frac{\sqrt{p_\mu(\lambda) p^\mu(\lambda)}}{[1 + z_G(\lambda)]^3} \\ \times n_G[x^\mu(\lambda), \theta_G] \mathcal{L}_G(\nu_G, \theta_G), \end{aligned} \quad (36)$$

where we have used the reciprocity relation (14). In this equation, the total effective luminosity per unit of effective frequency, $\mathcal{L}_G$ is defined in Eq. (31) as a sum of contributions from inspiraling binaries, mergers, and SNe. It depends on the local physics inside the galaxies while the other factors in the integral depend on the cosmology.

In order to compare our expression with standard results in the literature, we split the GW energy density as

$$\frac{d^3 \rho_{\text{GW}}}{d\nu_o d^2 \Omega_o} (\nu_o, e_o) \equiv \frac{d^3 \rho_{\text{GW}}^I}{d\nu_o d^2 \Omega_o} + \frac{d^3 \rho_{\text{GW}}^M}{d\nu_o d^2 \Omega_o} + \frac{d^3 \rho_{\text{GW}}^{SN}}{d\nu_o d^2 \Omega_o}, \quad (37)$$

respectively, for inspiraling binaries, mergers, and SNe. The explicit expression of each one of these terms can be straightforwardly derived comparing (37) to Eq. (36) with $\mathcal{L}_G \equiv \mathcal{L}_G^I + \mathcal{L}_G^M + \mathcal{L}_G^{SN}$.

We conclude observing that the expression (36) as well as the definition (37) are completely general in the sense that they do not assume any specific form for the metric. They can be specialized to any cosmology.

### III. COSMOLOGICAL FRAMEWORK

We now consider the standard cosmological framework [23] in which the Universe is modeled by a FLRW universe with Euclidean spatial sections and with scalar perturbations. In Newtonian gauge, the metric $g_{\mu\nu}$ is given by

$$ds^2 = a^2 [-(1 + 2\psi)d\eta^2 + (1 - 2\phi)\delta_{ij}dx^i dx^j], \quad (38)$$

where the metric of the constant time hypersurfaces is

$$\delta_{ij} dx^i dx^j = d\chi^2 + \chi^2 (d\theta^2 + \sin^2\theta d\phi^2) \quad (39)$$

in terms of the comoving radial distance $\chi$. The two Bardeen potentials are decomposed as

$$\psi = \Psi + \Pi, \quad \phi = \Psi - \Pi. \quad (40)$$

In the standard $\Lambda$CDM model, the matter content at late time is dominated by cold dark matter (CDM), described by

<footnote n="13">¹³What is actually measured in an image is the specific intensity $I$, i.e. the intensity in a specific frequency range and solid angle around the direction of observation.</footnote>
a pressureless fluid, and by the cosmological constant. It follows that the Bardeen potentials, $\phi$ and $\psi$, are equal, so that $\psi = \phi = \Psi$ and $\Pi = 0$.

We assume that the galaxies are all comoving with the cosmic flow.<fnref n="14" /> To first order in perturbations, the four-velocity of the cosmic fluid is given by

$$u^{\mu} \equiv \frac{1}{a}(1 - \psi, v^{i}) \equiv \bar{u}^{\mu} + \delta u^{\mu}, \qquad (41)$$

where $v^i$ is the peculiar velocity field. From the matter conservation equation, the galaxy peculiar velocity can be related to the gravitational potential through the Euler equation.

### A. Strictly spatially homogeneous and isotropic universe

We start by considering an unperturbed FLRW universe [i.e. $\phi = \psi = 0$ in Eq. (38)]. Homogeneity and isotropy imply that $\partial_i$ is the Killing vector associated with the invariance under spatial translations. Therefore, $g(\partial_i, k) = k_i$ is a constant of geodesic motion. Moreover, since $k$ is a null vector, $E^2 = g^{ij}k_i k_j$. From the geodesic equation, we get

$$\bar{k}^{\mu}(\lambda) = \frac{1}{a^2(\lambda)}(1, -e_o), \qquad (42)$$

where $e_o$ denotes the spatial direction of observation. From now on it will be denoted simply by $e$, and we denote with an overbar background quantities. It follows

$$\bar{E}(\lambda) = \frac{1}{a(\lambda)}. \qquad (43)$$

Further integrating Eq. (42), we find for the background geodesic

$$\bar{x}^0(\lambda) - \bar{x}^0(\lambda_o) = \eta - \eta_o \equiv \Delta\eta, \qquad (44)$$

$$\bar{x}^i(\lambda) - \bar{x}^i(\lambda_o) = e^i \Delta\eta. \qquad (45)$$

The spatial direction of propagation of GW is defined in Eq. (8) and now reads

$$\bar{p}^{\mu}(\lambda) = \frac{1}{a^2(\lambda)}(0, -e), \qquad (46)$$

and the physical volume defined in Eq. (35) reduces to

$$d^3V = \frac{1}{(1 + \bar{z})^3} \frac{\chi^2(\bar{z})}{H(\bar{z})} d^2\Omega_o d\bar{z} \qquad (47)$$

so that Eq. (36) reduces to

$$\frac{d^3 \rho_{GW}}{d\nu_o d^2 \Omega_o}(\nu_o) = \frac{1}{4\pi H_0} \int d\bar{z} \frac{1}{E(\bar{z})} \frac{1}{(1 + \bar{z})^4} \times \int d\theta_G \bar{n}_G(\bar{z}, \theta_G) \mathcal{L}_G(\nu_G, \theta_G). \qquad (48)$$

In Eq. (48), $\bar{n}_G(\bar{z}, \theta_G)$ is the background (homogeneous and isotropic) density of galaxies. We have used that in a FLRW universe $D_L(\bar{z}) = (1 + \bar{z})\chi(\bar{z})$ and we have defined $E(\bar{z}) \equiv H(\bar{z})/H_0$. As expected, this expression does not depend on the direction of observation. When specialized to merging binaries (i.e. setting $\mathcal{L}_G^{\text{L,SN}} = 0$ inside $\mathcal{L}_G$), Eq. (48) once replaced in Eqs. (2) and (1) is in agreement with the standard literature, e.g. with Eq. (3.18) of Ref. [38].<fnref n="15" />

### B. Perturbed FLRW universe

We now consider the effect of the scalar perturbations in the metric (38). Our goal is to obtain an expression for the GW background energy density at first order in the metric perturbations. Let us consider a graviton whose null worldline is intersecting the worldline of the observer at time $\eta_o$. We decompose it as

$$x^{\mu}(\lambda) = \bar{x}^{\mu}(\lambda) + \xi^{\mu}(\lambda), \qquad (49)$$

where $\bar{x}^\mu$ is the null geodesic of the FLRW background spacetime and $\xi^\mu$ represents the deviation from the background geodesic. To solve the geodesic equation we use that geodesics are the same in conformally related metrics. We introduce a metric $\tilde{g}_{\mu\nu}$ such that $ds^2 = a^2 d\bar{s}^2$ and

$$k^{\mu} = \frac{1}{a^2} \tilde{k}^{\mu}, \qquad \lambda = a^2 \tilde{\lambda}. \qquad (50)$$

The null geodesic of the metric $\tilde{g}_{\mu\nu}$ reads

$$\frac{d\delta\tilde{k}^{\mu}}{d\tilde{\lambda}} \equiv \frac{d^2\tilde{\xi}^{\mu}}{d\tilde{\lambda}^2} = -\delta\tilde{\Gamma}_{\alpha\beta}^{\mu} \tilde{k}^{\alpha} \tilde{\bar{z}}^{\beta}, \qquad (51)$$

where an overbar indicates background quantities and $\delta\tilde{\Gamma}_{\alpha\beta}^{\mu}$ denotes the Levi-Civita connection of the $\tilde{g}_{\mu\nu}$ metric at first

<footnote n="14">¹⁴The velocity of galaxies is not biased $v(z, e) = v_{\text{CDM}}(z, e)$.</footnote>

<footnote n="15">¹⁵In Ref. [38] the peculiar motion of sources inside galaxies is not taken into account. Therefore, to compare our result with the one of [38], we need to set $f(\Gamma, \theta_G) = \delta(\Gamma)$ in Eq. (48). As a consequence, in our result we need to identify the frequency $\nu_G$ with the frequency $\nu_S$. This explains the difference between our result (48) and the one of Ref. [38] in the expression of the local astrophysical contribution (i.e. inside $\mathcal{L}_G$ in our expression). Moreover, in Ref. [38] $n_G$ is the comoving density of the galaxy and not the physical one, and the merging rate is defined with respect to the observer time and not the time at the source. This explains the additional factor of $(1 + \bar{z})^{-4}$ in our expression.</footnote>
order in scalar perturbations. Solving Eq. (51) and going back to the original metric $g_{\mu\nu}$ gives

$$ \delta k^0(\lambda) = -\frac{2}{a^2}[\Psi + \Pi]_{\lambda_o}^{\lambda} + \frac{2}{a^2}\int_{\eta_o}^{\eta} d\eta' \dot{\Psi}, \quad (52) $$

$$ \delta k^i(\lambda) = \frac{2}{a^2}(\Pi - \Pi_o - \Psi)e^i - \frac{2}{a^2}\int_{\eta_o}^{\eta} d\eta' \partial_i \Psi, \quad (53) $$

where we have imposed initial conditions at $\lambda = \lambda_o$ in such a way that the null-geodesic condition $k_\mu k^\mu = 0$ is satisfied up to first order in perturbations. From Eq. (4), $dx^i/d\eta = k^i/k^0$; hence, we obtain

$$ \xi^i = -2 \int_{\eta_o}^{\eta} d\eta' [e^i \Psi + (\eta - \eta')(\delta_i^j - e_i e^j) \partial_j \Psi], \quad (54) $$

where we have set $\xi^i(\eta_o) = 0$.

The redshift of a graviton emitted from a galaxy at spacetime position $G$ and observed at $O$ is defined in Eq. (9). Since gravitons follow null geodesics, we find

$$ 1 + z_G = \frac{1}{a_G} \left[ 1 + (\Psi + \Pi - v \cdot e)|_G^0 + 2 \int_{\eta_o}^{\eta_G} d\eta \dot{\Psi} \right], \quad (55) $$

which we decompose as

$$ z_G \equiv \bar{z} + \delta z, \quad (56) $$

where the first-order correction $\delta z$ is due to the metric perturbations and peculiar velocities.

We now need to specialize Eq. (36) to the perturbed metric (38) and expand the various functions in perturbations in order to derive the expression of the energy density received up to first order in perturbations. The density of galaxies $n(x^\mu, \theta_G)$ can be expanded as

$$ n_G(x^\mu) = \bar{n}_G(\bar{x}^\mu) + \xi^\mu \nabla_\mu \bar{n}_G(\bar{x}^\mu) + \delta n_G(\bar{x}^\mu), \quad (57) $$

where we have omitted the dependence on $\theta_G$. The first term is the background contribution while the other two are first-order corrections. Note that since $\bar{n}_G$ is a background quantity, it does not depend on spatial coordinates, so $\xi^\mu \nabla_\mu \bar{n}_G = 0$. We define

$$ \frac{\delta n_G(\bar{x}^\mu, \theta_G)}{\bar{n}(\bar{x}^\mu)} \equiv b(\eta)\delta_{\text{CDM}}(\bar{x}^\mu), \quad (58) $$

where in full generality the bias $b$ depends on time coordinates; $\delta_{\text{CDM}}$ denotes the CDM density contrast. Similarly, the physical volume in $x^\mu(\lambda)$ can be expanded as

$$ d^3 V(x^\mu) = d^3 \bar{V}(\bar{x}^\mu) + \xi^\mu \nabla_\mu d^3 \bar{V}(\bar{x}^\mu) + \delta d^3 V(\bar{x}^\mu). \quad (59) $$

However, as explained in Sec. II C, by using the expression of the physical volume in terms of the observed solid angle and the angular diameter distance, Eq. (35)—the final

expression for the total flux received, Eq. (36)—depends on the angular diameter distance only through the ratio $D_A/D_L$, given by $(1+z)^{-2}$ by the reciprocity relation (14). We thus need to expand only the direction of propagation $p^\mu$ defined in Eq. (8) up to first order in perturbations. We get $p^\mu \equiv \bar{p}^\mu + \delta p^\mu$ with $\bar{p}^\mu$ given in Eq. (46) and

$$ \delta p^0 = -\frac{1}{a^2} v \cdot e, \quad (60) $$

$$ \delta p^i = -\frac{2}{a^2} \left( e^i (\Psi - \Pi + \Pi_o) + \int_{\eta_o}^{\eta} d\eta' \partial_i \Psi + \frac{v^i}{2} \right). \quad (61) $$

After simplifications, it can be shown that

$$ \sqrt{p^\mu p^\nu g_{\mu\nu}} = \frac{1}{a} + \frac{1}{a} \left[ \Psi - \Pi + 2\Pi_o + e \cdot v + 2 \int_{\eta_o}^{\eta} d\eta' e \cdot \nabla \Psi \right]. \quad (62) $$

This quantity can be related to the ratio between proper time $\tau$ of the observer and conformal time $\eta$ as

$$ \frac{d\tau}{d\eta} \equiv \sqrt{S_{\alpha}^{\mu} S_{\beta}^{\nu} \frac{dx^{\alpha}}{d\eta} \frac{dx^{\beta}}{d\eta} g_{\mu\nu}} = \sqrt{p^{\mu} p^{\nu} g_{\mu\nu}} \left(\frac{d\lambda}{d\eta}\right), \\ = a[1 + \Psi + \Pi + v \cdot e]. \quad (63) $$

Now, we need to substitute Eqs. (55), (57), and (62) for the redshift, galaxy number density, and the norm of the direction of propagation of GW in Eq. (36) and keep terms up to first order in scalar perturbations. For redshift and galaxy density, we get, respectively,

$$ 1 + z_G = \frac{1}{a} \left[ 1 + \Psi_o - \Psi + \Pi_o - \Pi - v_o \cdot e + v \cdot e + 2 \int_{\eta_o}^{\eta} d\eta' \dot{\Psi} \right], \quad (64) $$

$$ n_G(x^\mu) = \bar{n}_G(\eta) + \delta n_G. \quad (65) $$

As a cross-check, one can compute $\sqrt{p^\mu p^\nu g_{\mu\nu}} \equiv -(u^\mu k_\mu) = -(1+z_G)(u^\mu k_\mu)_o$ with $(1+z_G)$ given by Eq. (64) and verify that it reduces to the expression (62) for the norm of $p^\mu$. Gathering Eqs. (62), (65), and (64), we get
$$ 
\begin{aligned}
& \frac{\sqrt{p_{\mu}(\lambda)p^{\mu}(\lambda)}}{(1+z_G(\lambda))^3} n_G[x^{\mu}(\lambda), \theta_G] \\
& = a^2 \bar{n}_G \left[ 1 + b\delta_{\text{CDM}} + 2\Psi - \Psi_o + 2\Pi - \Pi_o \right. \\
& \qquad \left. - 2\mathbf{e} \cdot \nabla v + 3\mathbf{e} \cdot \nabla v_o - 4 \int_{\eta_o}^{\eta} d\eta' \dot{\Psi} \right].
\end{aligned}
\quad (66) $$

Then, inserting this result into Eq. (36) for the total energy density of GW and using the conformal time as integration variable, we get

$$ 
\begin{aligned}
& \frac{d^3 \rho_{\text{GW}}}{d\nu_o d^2\Omega_o}(\nu_o, \mathbf{e}) \\
& = \frac{1}{4\pi} \int d\eta a^4 \int d\theta_G \bar{n}_G \mathcal{L}_G(\nu_G, \theta_G) \\
& \quad \times \left[ 1 + b\delta_{\text{CDM}} + 4\Psi + 4\Pi - 3\Psi_o - 3\Pi_o - 2\mathbf{e} \cdot \nabla v \right. \\
& \quad \left. + 3\mathbf{e} \cdot \nabla v_o - 6 \int_{\eta_o}^{\eta} d\eta' \dot{\Psi} \right],
\end{aligned}
\quad (67) $$

where we have substituted $v = \nabla v$ since we are considering only the contribution of scalar perturbations. The integration runs along the unperturbed geodesic, parametrized as $\mathbf{x} = \mathbf{x}_o + \mathbf{e}(\eta_o - \eta)$.

We observe that in Eq. (67) the effective luminosity $\mathcal{L}_G(\nu_G, \theta_G)$ is not a background quantity since $\nu_G = (1+z_G)\nu_o$ where the observed redshift $z_G$ is given by Eq. (55) up to first order in perturbations. We write

$$ \nu_G = \bar{\nu}_G + \delta\nu_G, \quad (68) $$

with

$$ \bar{\nu}_G = (1 + \bar{z})\nu_o, \quad \delta\nu_G = \nu_o\delta z, \quad (69) $$

where $\bar{z}$ and $\delta z$ are defined in Eq. (55). We can therefore expand the effective luminosity in Eq. (67) around the background value of the effective frequency $\nu_G$ as

$$ \mathcal{L}_G(\nu_G, \theta_G) = \mathcal{L}_G(\bar{\nu}_G, \theta_G) + \frac{\partial \mathcal{L}_G}{\partial \nu_G}(\nu_G, \theta_G)|_{\bar{\nu}_G} \delta \nu_G. \quad (70) $$

We plug this expansion into Eq. (67), and we split the result into a homogeneous and isotropic background contribution and a first order contribution as

$$ 
\begin{aligned}
& \frac{d^3 \rho_{\text{GW}}}{d\nu_o d^2\Omega_o}(\eta_o, \mathbf{x}_o, \mathbf{e}, \nu_o) \\
& = \frac{d^3 \bar{\rho}_{\text{GW}}}{d\nu_o d^2\Omega_o}(\eta_o, \nu_o) + \mathcal{E}(\eta_o, \mathbf{x}_o, \mathbf{e}, \nu_o),
\end{aligned}
\quad (71) $$

where we have explicitly kept the dependence on the observer spacetime coordinates $(\eta_o, \mathbf{x}_o)$. The perturbation $\mathcal{E}(\eta_o, \mathbf{x}_o, \mathbf{e}, \nu_o)$ is explicitly given by

$$ 
\begin{aligned}
& \mathcal{E}(\eta_o, \mathbf{x}_o, \mathbf{e}, \nu_o) \\
& = \frac{1}{4\pi} \int d\eta a^4 \int d\theta_G \bar{n}_G \mathcal{L}_G(\bar{\nu}_G, \theta_G) \\
& \quad \times \left\{ b\delta_{\text{CDM}} + 4\Psi + 4\Pi - 2\mathbf{e} \cdot \nabla v \right. \\
& \quad \left. - 6 \int_{\eta_o}^{\eta} d\eta' \dot{\Psi} + \frac{1}{\mathcal{L}_G} \frac{\partial \mathcal{L}_G}{\partial \nu_G} \Big|_{\bar{\nu}_G} \frac{\nu_o}{a} \right. \\
& \quad \left. \times \left[ \mathbf{e} \cdot \nabla v - \Psi - \Pi + 2 \int_{\eta_o}^{\eta} d\eta' \dot{\Psi} \right] \right\},
\end{aligned}
\quad (72) $$

where we have set to zero all quantities at the observer, since they do not depend on $\mathbf{e}$ and thus do not contribute to the anisotropy and are a correction to the monopole of the GW background. This expression can be compared to Eq. (9) of Ref. [59] and gives the expression of the emissivity in terms of the astrophysical details of the galaxies and thus its frequency dependence.

In Eqs. (71) and (72) we distinguish two different types of contributions: (1) a contribution from local physics, proportional to the galaxy effective luminosity $\mathcal{L}_G$ and (2) a cosmological contribution which depends on metric perturbations, matter density contrast, and velocities. In the cosmological part, we can distinguish different types of terms: a contribution propositional to the CDM overdensity, coming from the perturbative expansion of the galaxy density $n_G$ in the covariant expression (36); a Doppler contribution due to the peculiar motion of galaxies with respect to the observer rest frame; local contributions and a contribution integrated along the line of sight. We expect the local contribution proportional to the CDM overdensity to be the dominant term, together with the Doppler one. A quantitative analysis of this expression and an estimate of its various contributions will be presented in [53].

We conclude this section observing that the cosmological factor in Eq. (72) has a structure very similar to the one of the Sachs-Wolfe formula for CMB temperature anisotropies. Nevertheless, this result is absolutely not trivial and could not be guessed from its electromagnetic counterpart. Indeed, it is true that in both cases (CMB and AGWB) a line of sight approach is exploited together with the eikonal approximation. However, for the electromagnetic case, photons are coming to the observer from a constant time hypersurface, while in the case under study, sources emitting GWs have a nontrivial redshift distribution. As a consequence, in this gravitational situation, when calculating (71), we had to take into account an effect of volume distortion, absent in the CMB case.<fnref n="16" />

<footnote n="16">¹⁶In deriving Eqs. (36) and (71) the reciprocity relation (14) has been exploited and it allowed us to get considerable computational simplifications.</footnote>
IV. ANGULAR POWER SPECTRUM

The function $\mathcal{E}(\eta_0, \mathbf{x}_0, \nu_0, \mathbf{e})$ defined in Eq. (72) is a stochastic variable which depends on the observer's position and direction of observation. The cosmological variables are correlated stochastic variables whose spectra are related to a scenario of structure formation [23]. It follows that $\mathcal{E}$ can be characterized by its angular correlation function

$$C(\theta) \equiv \langle \mathcal{E}(\eta_0, \mathbf{x}_0, \nu_0, \mathbf{e}_1) \mathcal{E}(\eta_0, \mathbf{x}_0, \nu_0, \mathbf{e}_2) \rangle. \quad (73)$$

Statistical isotropy implies that the correlation function depends only on the relative angle $\theta$; $\mathbf{e}_1 \cdot \mathbf{e}_2 = \cos\theta$. This correlation function can be expanded in Legendre polynomials

$$C(\theta) \equiv \sum_{\ell} \frac{2\ell + 1}{2\pi} C_{\ell} P_{\ell}(\mathbf{e}_1 \cdot \mathbf{e}_2). \quad (74)$$

A given multipole $\ell$ corresponds to the typical angular scale $\pi/\theta$ and $C_\ell$ is an estimate of the variance of GW energy density fluctuations on that scale. To get its expression, we first decompose $\mathcal{E}$ in spherical harmonics as

$$\mathcal{E}(\eta_0, \mathbf{x}_0, \nu_0, \mathbf{e}) = \sum_{\ell m} a_{\ell m}(\eta_0, \mathbf{x}_0, \nu_0) Y_{\ell m}(\mathbf{e}), \quad (75)$$

$$\hat{\mathcal{E}}_k(\eta_0, \mathbf{e}, \nu_0) = \frac{1}{4\pi} \int d\eta a^4 \int d\theta_G \mathcal{L}_G(\bar{\nu}_G, \theta_G) \bar{n}_G(\eta, \theta_G) \\ \times \left\{ e^{ik\mu\Delta\eta} \left[ 4\hat{\Psi}_k(\eta) + 4\hat{\Pi}_k(\eta) + b\hat{\delta}_k(\eta) - 2ik\mu\hat{v}_k(\eta) + \frac{1}{\mathcal{L}_G} \frac{\partial \mathcal{L}_G}{\partial \nu_G} \Big|_{\bar{\nu}_G} \frac{\nu_0}{a} (-\hat{\Psi}_k(\eta) - \hat{\Pi}_k(\eta) + ik\mu\hat{v}_k(\eta)) \right] \right. \\ \left. + 2 \left( \frac{1}{\mathcal{L}_G} \frac{\partial \mathcal{L}_G}{\partial \nu_G} \Big|_{\bar{\nu}_G} \frac{\nu_0}{a} - 3 \right) \int_{\eta_0}^{\eta} d\eta' e^{ik\mu\Delta\eta'} \hat{\Psi}_k(\eta') \right\}, \quad (79)$$

where $\mathbf{e} \cdot \mathbf{k} = k\mu$. Each stochastic variable $\hat{X}_k(\eta)$ in Eq. (79) can be decomposed as the product of a transfer function $X_k(\eta)$ and a unique random variable $\hat{a}_k$ as

$$\hat{X}_k(\eta) = \hat{X}_k(\eta)\hat{a}_k. \quad (80)$$

$\hat{X}_k(\eta)$ depends only on the modulus of $\mathbf{k}$ and $\hat{a}_k$ satisfies

$$\langle \hat{a}_k \hat{a}_k^* \rangle = (2\pi)^3 \delta^{(3)}(\mathbf{k} - \mathbf{k}'). \quad (81)$$

It follows that the mode function of $\hat{\mathcal{E}}$ can be decomposed as

$$\hat{\mathcal{E}}_k(\eta_0, \nu_0, \mathbf{e}) = \hat{\mathcal{E}}_k(\eta_0, \nu_0, \mu)\hat{a}_k. \quad (82)$$

where the coefficients of the development are given by

$$a_{\ell m}(\eta_0, \mathbf{x}_0, \nu_0) \equiv \int d^2 \mathbf{e} \mathcal{E}(\eta_0, \mathbf{x}_0, \mathbf{e}, \nu_0) Y_{\ell m}^*(\mathbf{e}). \quad (76)$$

Making use of the properties of Legendre polynomials, it is easy to show that

$$(2\ell + 1)C_{\ell} \equiv \sum_{m} \langle a_{\ell m} a_{\ell m}^{*} \rangle, \quad (77)$$

where the brackets stand for an ensemble average over the stochastic variables. Statistical homogeneity implies that this shall not depend on $\mathbf{x}_0$. From now on we will thus omit the dependence of $a_{\ell m}$ on quantities at the observer.

In order to calculate $C_\ell$, it is useful to first expand $\mathcal{E}$ in Fourier modes as

$$\mathcal{E}(\eta_0, \mathbf{x}_0, \nu_0, \mathbf{e}) = \int \frac{d^3 \mathbf{k}}{(2\pi)^3} \hat{\mathcal{E}}_k(\eta_0, \nu_0, \mathbf{e}), \quad (78)$$

where we have reabsorbed a phase $\exp(ik \cdot \mathbf{x}_0)$ into the definition of the field $\hat{\mathcal{E}}$. All variables are decomposed in a similar way and with the same convention for the phase $\exp(ik \cdot \mathbf{x}_0)$. Since $\mathcal{E}$ is an integral over the background geodesic $\mathbf{x}(\eta) = \mathbf{x}_0 + \mathbf{e}(\eta_0 - \eta)$, we get

We now expand the exponentials in spherical harmonics and spherical Bessel functions using the standard relation

$$e^{ip \cdot x} = 4\pi \sum_{\ell m} i^{\ell} j_{\ell}(px) Y_{\ell m}^{*}(\hat{p}) Y_{\ell m}(\hat{x}). \quad (83)$$

With the definition

$$\hat{\mathcal{E}}_k(\eta_0, \nu_0, \mu) = 4\pi \sum_{\ell m} i^{\ell} \hat{\mathcal{E}}_{\ell}(k, \nu_0) Y_{\ell m}^{*}(\hat{k}) Y_{\ell m}(\mathbf{e}), \quad (84)$$

we finally get
$$
\begin{equation}
\begin{split}
\hat{\mathcal{E}}_\ell(k, \nu_o) ={}& \frac{1}{4\pi} \int d\eta a^4 \int d\theta_G \mathcal{L}_G(\bar{\nu}_G, \theta_G) \bar{n}_G(\eta, \theta_G) \\
& \times \left\{ [4\hat{\Psi}_k(\eta) + 4\hat{\Pi}_k(\eta) + b\hat{\delta}_k(\eta)] j_\ell(k\Delta\eta) - 2k\hat{v}_k(\eta) j'_\ell(k\Delta\eta) - 6 \int_{\eta_o}^\eta d\eta' \hat{\Psi}_k(\eta') j_\ell(k\Delta\eta') \right. \\
& \left. + \frac{1}{\mathcal{L}_G} \frac{\partial \mathcal{L}_G}{\partial \nu_G} \Big|_{\bar{\nu}_G} \frac{\nu_o}{a} \left[ -\hat{\Psi}_k(\eta) j_\ell(k\Delta\eta) - \hat{\Pi}_k(\eta) j_\ell(k\Delta\eta) + k\hat{v}_k(\eta) j'_\ell(k\Delta\eta) + 2 \int_{\eta_o}^\eta d\eta' \hat{\Psi}_k(\eta') j_\ell(k\Delta\eta') \right] \right\},
\end{split}
\tag{85}
\end{equation}
$$

where the prime acting on the spherical Bessel denotes a
derivative with respect to its argument.

We can now go back to the expression (77) defining the
angular power spectrum $C_\ell$. Plugging in Eq. (76) defining
the multipoles $a_{\ell m}$, the Fourier transform (78) of $\mathcal{E}$, and
using Eqs. (82)–(85), we get

$$
C_{\ell}(\nu_{o}) = \frac{2}{\pi} \int dk k^{2} |\mathcal{E}_{\ell}(k, \nu_{o})|^{2}. \quad (86)
$$

This expression is similar to those of CMB temperature and
polarization angular power spectrum. Note, however, that
we have one angular spectrum per frequency band.

From the result (86) we can compute the variance of
$\Omega_{\text{GW}}$ due to the distribution of the large scale structures. For
each frequency band, the variance is defined by

$$
\sigma_{\text{GW}}^{2}(\nu_{o}) \equiv \langle \delta\rho_{\text{GW}}(\nu_{o}, e_{1}) \delta\rho_{\text{GW}}(\nu_{o}, e_{2}) \rangle |_{e_{1}, e_{2}=1}, \quad (87)
$$

so using Eq. (74), we get

$$
\sigma_{\text{GW}}^{2}(\nu_{o}) \equiv \sum_{\ell} \frac{(2\ell + 1)}{2\pi} C_{\ell}(\nu_{o}). \quad (88)
$$

The computation of $\sigma_{\text{GW}}^2(\nu_o)$ requires solving the evolution of the perturbations and of the galaxy number counts (via a merger-tree analysis). The evolution of the cosmological perturbations can be split in terms of a transfer function and an initial power spectrum. When the dominant contribution to (86) is proportional to $\delta_{\text{CDM}}$ and the initial power spectrum is scale invariant, it takes the form

$$
\sigma_{\text{GW}}^2(\nu_o) = \sum_{\ell} \frac{2\ell + 1}{16\pi^4} \int \frac{dk}{k} \times \left[ \int d\eta a^4(\eta) \int d\theta_G \mathcal{L}(\bar{\nu}_G, \theta_G) \bar{n}_G(\eta, \theta_G) \times b(\eta) T_{\delta}(k, \eta) j_{\ell}(k\Delta\eta) \right]^2, \quad (89)
$$

where $T_{\delta}(k, \eta)$ is the transfer function of matter density contrast. Approximate analytic expressions for $\sigma_{\text{GW}}^2(\nu_o)$ and its numerical computation will be presented in our companion article [53].

V. CORRELATION WITH LENSING

The GW energy density (33) depends on the local
properties of galaxies as well as on cosmological pertur-
bations. This means that it will be correlated with other
probes and, in particular, weak lensing of background
galaxies and galaxy number counts.

Weak lensing describes the deformation of the shape of
background galaxies by the gravitational potential of the
large scale structure. Its description is now standard; see
e.g. Refs. [23,73–75]. The observed distortion in a given
direction is characterized by the Jacobi matrix $\mathcal{D}$ and the
associated amplification matrix $\mathcal{A}$, which relate the
observed and intrinsic shapes.

In a perturbed FLRW spacetime, the amplification
matrix splits into a background and a perturbed contribu-
tion as $\mathcal{A} = \mathcal{A}^{(0)} + \mathcal{A}^{(1)}$. The background part reduces to
$\mathcal{A}^{(0)} = \mathbf{1}$ and the first-order contribution is

$$
(\mathcal{A}^{(1)})_{ab}(\eta, \mathbf{e}) = -2 \int_{\eta_0}^{\eta} d\eta' \frac{\eta'(\eta - \eta')}{\eta} \partial_a \partial_b \Psi, \quad (90)
$$

where the integral runs along the unperturbed geodesic
**x** = **x**_o + **e**(*η*₀ − *η*). The trace of the amplification matrix is
the convergence **κ**, while its tracefree part defines the
cosmic shear (*γ*₁, *γ*₂). The convergence is given by

$$
\kappa^{(1)}(\eta, \mathbf{e}) = -2 \int_{\eta_0}^{\eta} d\eta' \frac{\eta'(\eta - \eta')}{\eta} \Delta\Psi \quad (91)
$$

for a source located at $\eta$. The total convergence is obtained
by integrating this contribution weighted by the number of
sources along the line of sight,

$$
\kappa^{(1)}(\boldsymbol{e}) = -2 \int d\eta p_{\eta}(\eta) \int_{\eta_0}^{\eta} d\eta' \frac{\eta'(\eta - \eta')}{\eta} \Delta\Psi, \quad (92)
$$

where $p_\eta$ is a function that takes into account the redshift distribution of sources, $d\eta p_\eta(\eta) = dz p_z(z)$. Equation (92) can be rewritten as

$$
\kappa^{(1)}(\boldsymbol{e}) = -2 \int d\eta g(\eta) \eta \Delta\Psi, \quad (93)
$$

with
$$g(\eta) \equiv \int_{\eta}^{\eta_H} d\eta' p_{\eta}(\eta') \frac{\eta' - \eta}{\eta'}, \quad (94)$$

and $\eta_H$ is the value of $\eta$ associated with the size of the observable Universe. In this analysis, we just consider the total convergence, but our result can be trivially extended to redshift bins if tomographic information exists.

The convergence (92) depends on the lensing potential $2\Psi$ and on the direction of observation. Since the GW energy density also depends on the perturbation variables, both quantities are correlated. We introduce the cross-correlation

$$B(\theta, \nu_o) \equiv \langle \mathcal{E}(\eta_o, \nu_o, \mathbf{x}_o, \mathbf{e}_1, \nu_o) \kappa(\eta_o, \mathbf{x}_o, \mathbf{e}_2) \rangle, \quad (95)$$

where $\mathbf{e}_1 \cdot \mathbf{e}_2 = \cos \theta$. This function can be expanded in Legendre polynomials as

$$B(\theta, \nu_o) = \sum_{\ell} \frac{2\ell + 1}{2\pi} B_{\ell}(\nu_o) P_{\ell}(\mathbf{e}_1 \cdot \mathbf{e}_2), \quad (96)$$

with $\mathbf{e}_1 \cdot \mathbf{e}_2 = \cos \theta$. Following the same steps as in Sec. IV, we expand $\kappa$ in multipoles as

$$\kappa(\eta_o, \mathbf{x}_o, \mathbf{e}) = \sum_{\ell m} \kappa_{\ell m}(\eta_o, \mathbf{x}_o) Y_{\ell m}(\mathbf{e}), \quad (97)$$

with

$$\kappa_{\ell m}(\eta_o, \mathbf{x}_o) = \int d^2\mathbf{e} \kappa(\eta_o, \mathbf{x}_o, \mathbf{e}) Y_{\ell m}^*(\mathbf{e}). \quad (98)$$

Inserting these decompositions together with Eqs. (75)–(76) and (97)–(98) for the multipolar decomposition of $\mathcal{E}$ into Eq. (95), we get

$$(2\ell + 1)B_{\ell}(\nu_o) = \sum_{m} \langle a_{\ell m}(\nu_o) \kappa_{\ell m}^* \rangle. \quad (99)$$

The calculation of this quantity follows the same line as in Sec. IV: one needs to expand $\mathcal{E}(\eta_o, \mathbf{x}_o, \nu_o, \mathbf{e})$ and $\kappa(\eta_o, \mathbf{x}_o, \mathbf{e})$ in Fourier modes, hence defining

$$\hat{\kappa}_k(\eta_o, \mathbf{e}) = 2k^2 \int d\eta g(\eta) \eta \hat{\Psi}_k(\eta') e^{ik\mu\Delta\eta'}. \quad (100)$$

Decomposing the Bardeen potential in terms of a transfer function $\Psi_k(\eta)$ and the unit random variable $\hat{a}_k$, one gets

$$\hat{\kappa}_k(\eta_o, \mathbf{e}) = \hat{\kappa}_k(\eta_o, k, \mu) \hat{a}_k. \quad (101)$$

After decomposing the exponential in spherical harmonics, we get

$$\hat{\kappa}_k(\eta_o, \mu) \equiv 4\pi \sum_{\ell m} i^{\ell} \hat{\kappa}_{\ell}(k) Y_{\ell m}^{*}(\hat{\mathbf{k}}) Y_{\ell m}(\mathbf{e}), \quad (102)$$

with

$$\hat{\kappa}_{\ell}(k) = 2k^2 \int d\eta g(\eta) \eta j_{\ell}(k\Delta\eta) \hat{\Psi}(\eta, k). \quad (103)$$

To conclude, we just need to substitute the expressions of the Fourier decompositions of $\mathcal{E}$ and $\kappa$, i.e. Eqs. (84)–(85) and (101)–(102), into the definitions (76) and (98) of the multipoles $a_{\ell m}$ and $\kappa_{\ell m}$ to get

$$B_{\ell}(\nu_o) = \frac{2}{\pi} \int dk k^2 \mathcal{E}_{\ell}^{*}(k, \nu_o) \kappa_{\ell}(k). \quad (104)$$

## VI. CORRELATION WITH GALAXY NUMBER COUNTS

Galaxy surveys provide catalogs with information about galaxy redshift $z_G$, angular coordinates, fluxes, and shapes. Redshift space coordinates need to be converted into physical coordinates, with which theoretical predictions are drawn.

The galaxy number counts [76,77], defined as the perturbation of the number of galaxies as a function of direction and redshift,<fnref n="17" /> are also correlated with the GW background. We consider the number of galaxies in a direction $\mathbf{e}$ at redshift $z_G$, which we denote as

$$N_G(\mathbf{e}, z_G) d^2\Omega d\mathbf{z}_G. \quad (105)$$

The redshift distribution $\langle N_G \rangle(z_G) d\mathbf{z}_G$ is deduced from this quantity taking an average over the sky. The galaxy density perturbation at a fixed redshift $z_G$ observed in a direction $\mathbf{e}$ is given by

$$\begin{aligned} \delta_z(\mathbf{e}, z_G) &\equiv \frac{n_G(\mathbf{e}, z_G) - \langle n_G \rangle(z_G)}{\langle n_G \rangle(z_G)} \\ &\equiv \frac{N_G(\mathbf{e}, z_G) - \langle N_G \rangle(z_G)}{\langle N_G \rangle(z_G)} - \frac{\delta\nu(\mathbf{e}, z_G)}{\nu(z_G)}, \end{aligned} \quad (106)$$

where we have introduced the galaxy number density

$$n_G(\mathbf{e}, z_G) \equiv \frac{N_G(\mathbf{e}, z_G)}{\nu(\mathbf{e}, z_G)}, \quad (107)$$

and the physical survey volume density, $\nu$, per redshift bin and solid angle. $\nu$ is also a perturbed quantity since the solid angle of observation as well as the redshift bin are distorted between the source and the observer. Expressing the physical volume as $V(\mathbf{e}, z_G) \equiv \nu(\mathbf{e}, z_G) d\mathbf{z}_G d^2\Omega$, it follows that

$$\nu(\mathbf{e}, z_G) \equiv \nu(z_G) + \delta\nu(\mathbf{e}, z_G). \quad (108)$$

<footnote n="17">¹⁷Being an observable quantity and therefore gauge invariant, we have the freedom to derive it in an arbitrary gauge.</footnote>
We can therefore define the galaxy number counts as

$$ \Delta(\mathbf{e}, z_G) \equiv \frac{N_G(\mathbf{e}, z_G) - \langle N_G \rangle(z_G)}{\langle N_G \rangle(z_G)} \quad (109) $$

$$ \equiv \delta_z(\mathbf{e}, z_G) + \frac{\delta\nu(\mathbf{e}, z_G)}{\nu(z_G)} \quad (110) $$

in terms of redshift perturbation $\delta_z$ and volume density perturbation $\delta\nu/\nu$. All these quantities are gauge invariant. An explicit calculation of $\Delta$ in terms of cosmological quantities can be found in Ref. [76].

In the following, we consider the most dominant term of the galaxy number counts (beside the intrinsic density contrast), following the approximation first physically motivated by Kaiser [78]. Since then it has been considered to estimate the effect of peculiar velocities on the observed redshift. In the following, we refer to the FLRW comoving coordinate frame as real space and to the observed coordinates as redshift space, indicated with subscript obs. Galaxy number conservation must bring the same result in real space and redshift space,

$$ n_{\text{G,obs}}(\mathbf{x}_{\text{obs}})d^3\mathbf{x}_{\text{obs}} = n_{\text{G}}(\mathbf{x})d^3\mathbf{x}. \quad (111) $$

We now define

$$ n_{\text{G,obs}}(\mathbf{x}_{\text{obs}}) \equiv \bar{n}_{\text{G}}[1 + \delta_{\text{obs}}(\mathbf{x}_{\text{obs}})], \quad (112) $$

$$ n_{\text{G}}(\mathbf{x}) \equiv \bar{n}_{\text{G}}[1 + \delta(\mathbf{x})], \quad (113) $$

where we have used that in the absence of perturbations, we have no space distortion. From Eq. (111), it follows

$$ [1 + \delta_{\text{G,obs}}(\mathbf{x}_{\text{obs}})]d^3\mathbf{x}_{\text{obs}} = [1 + \delta_{\text{G}}(\mathbf{x})]d^3\mathbf{x}. \quad (114) $$

To isolate $\delta_{\text{G,obs}}(\mathbf{x}_{\text{obs}})$, we need to calculate the Jacobian of the transformation from redshift space coordinates to real space. Since only the radial coordinate is affected by redshift perturbations, we have

$$ \chi_{\text{obs}} \equiv \chi(z_G) \approx \chi(\bar{z}) + \frac{\partial\chi}{\partial\bar{z}} \delta z, \quad (115) $$

where we have used $z_G = \bar{z} + \delta z$ and $\partial\chi/\partial\bar{z} = \mathcal{H}^{-1}(1 + \bar{z})^{-1}$. Keeping only the Doppler contribution in Eq. (55), we find

$$ \chi_{\text{obs}} \approx \chi + \frac{\mathbf{e} \cdot \mathbf{v}}{\mathcal{H}}. \quad (116) $$

The Jacobian,

$$ J \equiv \left| \frac{d^3 \mathbf{x}}{d^3 \mathbf{x}_{\text{obs}}} \right| = \frac{d\chi}{d\chi_{\text{obs}}} \frac{\chi^2}{\chi_{\text{obs}}^2}, \quad (117) $$

reduces at linear order to

$$ J \approx 1 - \frac{\partial}{\partial \chi} \left[ \frac{\mathbf{e} \cdot \mathbf{v}}{\mathcal{H}} \right] - 2 \frac{\mathbf{e} \cdot \mathbf{v}}{\mathcal{H} \chi}, \quad (118) $$

where we have also kept the dominant term in the limit $|\mathbf{v}| \to 0$. Going back to Eq. (114), it follows

$$ \delta_{\text{G,obs}}(\mathbf{x}) = \delta_{\text{G}}(\mathbf{x}) + \frac{\partial}{\partial \eta} \left[ \frac{\mathbf{e} \cdot \mathbf{v}}{\mathcal{H}} \right] - 2 \frac{\mathbf{e} \cdot \mathbf{v}}{\mathcal{H} \chi}. \quad (119) $$

The last term is usually neglected, since the derivative of perturbations are expected to dominate<fnref n="18" /> together with the dominant intrinsic clustering term [83]. Identifying the observed number counts with the redshift-space overdensity, we obtain the so-called Kaiser approximation for redshift-space distortions

$$ \Delta(\mathbf{x}) \approx b\delta_{\text{CDM}}(\mathbf{x}) + \frac{\partial}{\partial \eta} \left[ \frac{\mathbf{e} \cdot \mathbf{v}}{\mathcal{H}} \right], \quad (120) $$

where $b$ is the bias (which we consider to be a function of time) and $\mathbf{v}$ the cold dark matter velocity flow. We can integrate along the line of sight to get the effective number count, function of the direction of observation,

$$ \Delta(\eta_o, \mathbf{x}_o, \mathbf{e}) \equiv \int d\eta \Delta(\eta_o, \mathbf{x}_o, \mathbf{e}, \eta), \quad (121) $$

where we have explicitly indicated the dependence on the observer spacetime coordinates. Explicitly, in the Kaiser approximation

$$ \Delta(\eta_o, \mathbf{x}_o, \mathbf{e}) \approx \int d\eta \left[ b\delta_{\text{CDM}}(\mathbf{x}) + \frac{\partial}{\partial\eta} \left( \frac{\mathbf{e} \cdot \mathbf{v}}{\mathcal{H}} \right) \right]. \quad (122) $$

Let us now turn to the correlation between $\Delta$ and the GW energy density defined in Eq. (71). It is given by

$$ D(\theta, \nu_o) \equiv \langle \mathcal{E}(\eta_o, \mathbf{x}_o, \mathbf{e}_1, \nu_o) \Delta(\eta_o, \mathbf{x}_o, \mathbf{e}_2) \rangle, \quad (123) $$

with $\mathbf{e}_1 \cdot \mathbf{e}_2 = \cos\theta$. Again, it can be expanded in Legendre polynomials as

$$ D(\theta, \nu_o) = \sum_{\ell} \frac{2\ell + 1}{2\pi} D_{\ell}(\nu_o) P_{\ell}(\mathbf{e}_1 \cdot \mathbf{e}_2). \quad (124) $$

Following the same procedure as in Sec. IV, we expand $\Delta$ in multipoles as

<footnote n="18">¹⁸Only on very large scales comparable to the Hubble radius $\mathcal{H}^{-1}$, the contribution of the last term is non-negligible and should be carefully considered when interested in wide-angle surveys [79–82].</footnote>
$$ \Delta(\eta_o, \mathbf{x}_o, \mathbf{e}) = \sum_{\ell m} d_{\ell m}(\eta_o, \mathbf{x}_o) Y_{\ell m}(\mathbf{e}), \quad (125) $$

with

$$ d_{\ell m}(\eta_o, \mathbf{x}_o) = \int d^2\mathbf{e} \Delta(\eta_o, \mathbf{x}_o, \mathbf{e}) Y_{\ell m}^*(\mathbf{e}), \quad (126) $$

so that

$$ (2\ell + 1)D_{\ell}(\nu_o) = \sum_m \langle a_{\ell m}(\nu_o) d_{\ell m}^* \rangle. \quad (127) $$

In order to calculate $D_\ell$, we expand $\Delta(\eta_o, \mathbf{x}_o, \mathbf{e})$ in Fourier space as

$$ \Delta(\eta_o, \mathbf{x}_o, \mathbf{e}) = \int \frac{d^3\mathbf{k}}{(2\pi)^3} \hat{\Delta}_k(\eta_o, \mathbf{e}), \quad (128) $$

the phase $\exp(i\mathbf{k} \cdot \mathbf{x}_o)$ being absorbed in the definition of the Fourier coefficient. We find

$$ \hat{\Delta}_k(\eta_o, \mathbf{e}) = \int d\eta \left[ b\hat{\delta}_{\text{CDM}k}(\eta) + ik\mu\partial_\eta\left(\frac{\hat{v}_k(\eta)}{\mathcal{H}}\right) \right] e^{ik\mu\Delta\eta}, \quad (129) $$

with $\mathbf{e} \cdot \mathbf{k} = k\mu$. Decomposing $\hat{v}_k(\eta)$ and $\hat{\delta}_{\text{CDM},k}(\eta)$ as the product of a random variable $\hat{a}_k$ and transfer function $\hat{v}_k(\eta)$ and $\hat{\delta}_{\text{CDM},k}$, we get

$$ \hat{\Delta}_k(\eta_o, \mathbf{e}) = \hat{a}(\mathbf{k}) \hat{\Delta}_k(\eta_o, \mu). \quad (130) $$

Expanding the exponential functions in spherical harmonics by using Eq. (83), we can define

$$ \hat{\Delta}_k(\eta_o, \mu) = 4\pi \sum_{\ell m} i^\ell \hat{\Delta}_\ell(k) Y_{\ell m}^*(\hat{\mathbf{k}}) Y_{\ell m}(\mathbf{e}), \quad (131) $$

with

$$ \hat{\Delta}_{\ell}(k) = \int d\eta \left[ b \hat{\delta}_{\text{CDM},k}(\eta) j_{\ell}(k\Delta\eta) + k \partial_{\eta} \left( \frac{\hat{v}_k}{\mathcal{H}} \right) j'_{\ell}(k\Delta\eta) \right]. \quad (132) $$

The cross-correlation is then obtained exactly in the same way as in the previous section, so

$$ D_{\ell}(\nu_o) = \frac{2}{\pi} \int dk k^2 \mathcal{E}_{\ell}^*(k, \nu_o) \Delta_{\ell}(k). \quad (133) $$

## VII. CONCLUSIONS

We have developed an analytic approach for the computation of the anisotropies of the background of gravitational radiation from unresolved galactic and extragalactic

sources. We have distinguished three different scales in the problem: a cosmological, a galactic, and an astrophysical scale. We have explained how to relate quantities measured by the observer (i.e. the directional flux of gravitational radiation per unit of frequency) to effective and local quantities characterizing the GW emission. At each scale, different mathematical tools need to be exploited: from cosmological perturbation theory at the cosmological scale to statistical mechanics and astrophysical calculation in order to characterize the local process of GW emission. This coarse graining approach enables us to obtain a parametrization for the observed flux of GW in a given direction, as a function of local quantities characterizing the process of GW emission and already presented in the literature, for various source types.

Our final results for the anisotropies of the AGWB energy density contains a cosmological component, which depends on the distribution of structures at large scales, and a local one, related to the local process of GW emission and to the intrinsic properties of the sources; see Eq. (36) for the general expression for AGWB anisotropies and Eq. (71) for the expression specialized to a FLRW universe with scalar perturbations. We have derived the expression of the AGWB angular power spectrum; see Eq. (86) as well as its correlation with lensing and galaxy number counts, Eqs. (104) and (133), respectively.

In the next few decades the LISA and ET experiments will provide new data on an unexplored frequency band, where many galactic and extragalactic astrophysical processes are expected to contribute. In principle, observations will allow us to reconstruct an angular map of the AGWB energy density, similar to the CMB map for electromagnetic radiation. On the other hand, modeling the astrophysical background as precisely as possible to extract information on its strength, frequency range, and statistical properties is crucial to distinguish it from the cosmological signal, foreground, or instrumental noise.

Even though this observational landscape may seem far away, our analysis (and its companion in [53]) will enable us to quantify the relative influences of the astrophysical and cosmological sources. It will either open a window on astrophysics (e.g. measurement of binaries rates, ...) if the former dominate, or provide a new cosmological probe if the result is robust to the unknown astrophysics. We observe that, since the angular power spectrum of the AGWB energy density has a frequency dependence, the relative importance of cosmological and astrophysical effects could depend on the frequency band chosen, hence offering the possibility to distinguish different astrophysical processes. The cross-correlation with lensing and galaxy distribution opens a window on the understanding of the distribution of GW sources, in analogy, e.g. with Ref. [57], albeit in a different context.

Our treatment assumed that all the GW sources are located inside galaxies. This can easily be generalized to
include more exotic scenarios, e.g. models in which dark matter is made of primordial back holes with masses of the order of 20–100 solar masses. Cross-correlations can indicate to what extent black-hole mergers are correlated with galaxies and hence study this hypothesis.

To conclude, we stress that in order to get quantitative predictions for the angular power spectrum of the AGWB energy density, the local physics processes generating GWs have to be parametrized in detail. A numerical study of the AGWB anisotropies, specifying local physics models, will be detailed in our companion article [53].

## ACKNOWLEDGMENTS

G. C. thanks IAP for warm hospitality during this work. We thank Chiara Caprini, Irina Dvorkin, Guillaume Faye, and Elisabeth Vangioni for discussions and comments on the early version of this text. We thank Enrico Barausse, Carlo Contaldi, Ruth Durrer, Pierre Fleury, and Atsushi Nishizawa for their comments on the first version of the paper. The work of G. C. is supported by the Swiss National Science Foundation. The work of C. P. and J. P. U. is made in the ILP LABEX (under Reference No. ANR-10-LABX-63) and was supported by French state funds managed by the ANR within the Investissements d'Avenir programme under Reference No. ANR-11-IDEX-0004-02.

## APPENDIX: REFERENCE FRAMES

From now on we focus on a single source $i$ emitting GWs in a galaxy. To simplify the notation, we omit the superscript ($i$). We introduce the following reference frames:

(1) $S: (E_A^\mu)_S \equiv (u^\mu, E_i^\mu)_S$ reference frame of the source, with the $z$ axis parallel to the axis of the source $A_S$.

(2) $K: (E_A^\mu)_K \equiv (u^\mu, E_i^\mu)_K$ reference frame of the galaxy, with the $z$ axis parallel to the boosted axis of the source $A_K$. It is related to the first one by a boost with velocity $\Gamma$.

(3) $G: (E_A^\mu)_G \equiv (u^\mu, E_i^\mu)_G$ reference frame of the galaxy, with the $z$ axis aligned with the axis of the galaxy $B_G$. It is related to the reference frame $(E_A^\mu)_K$ through a rotation with Euler angles $\alpha, \beta, \gamma$.

We write explicitly the relation among the various frames. $K$ and $G$ are related by a rotation

$$ (E_i)_K = R_j^i(\alpha, \beta, \gamma)(E^j)_G \\ = R_{E_G^\alpha}(\gamma)_j R_{E_G^\beta}(\beta)_k R_{E_G^\gamma}(\alpha)_\ell (E^\ell)_G. \quad (A1) $$

$S$ and $K$ are related by a boost with velocity $\Gamma$. Explicitly,

$$ (E_\mu^A)_K = \Lambda_B^A (E_\mu^B)_S, \quad (A2) $$

$$ (E_A^\mu)_K = (E_A^\mu)_S (\Lambda^{-1})_A^B = \Lambda_A^B (E_B^\mu)_S, \quad (A3) $$

with

$$ \begin{gathered} \Lambda^0_0 = \gamma, \quad \Lambda^0_i = \Lambda^i_0 = -\gamma\Gamma_i, \\ \Lambda^i_j = \delta^i_j + \frac{\gamma^2}{1+\gamma}\Gamma^i\Gamma_j, \quad (A4) \end{gathered} $$

and $\gamma^{-2} = 1 - \Gamma_i\Gamma^i$ and $\beta^2 \equiv \Gamma_i\Gamma^i$. A schematic representation of these reference frames is presented in Fig. 2.

### 1. Situation in S

The energy spectrum emitted in the source reference frame $S$ depends on the angle between the line of sight and the axis of the source, i.e. $e_S \cdot A_S \equiv \mu_S$, on the frequency $\nu_S$,

<figure><img src="image_3.png" /><figcaption>FIG. 2. Schematic representation of the change of reference frame we are considering. In the second line we represent in blue a generic source of GWs and in yellow the galaxy in which the source is contained.</figcaption></figure>
on the usual set of parameters $\theta_G$ and on the parameters $\theta^{(i)}$ characterizing the source under examination. To simplify the notation, we omit in the following this last dependence. We therefore have

$$ \frac{d^3 E_S}{d\nu_S d^2\Omega_S} (\nu_S, \mathbf{e}_S \cdot \mathbf{A}_S, \theta_G) = \frac{1}{2\pi} \frac{d^2 E_S}{d\nu_S d\mu_S} (\nu_S, \mu_S, \theta_G). \quad (A5) $$

## 2. Situation in K

We need to boost the previous spectrum defined in S. The 4-vector $k^\mu$ of a graviton is decomposed in the two frames S and K as

$$ k^{\mu} = \nu_{S}(u_{S}^{\mu} - \mathbf{e}_{S}^{\mu}) = \nu_{K}(u_{K}^{\mu} - \mathbf{e}_{K}^{\mu}). \quad (A6) $$

The relation between the frequency and the line of sight observed in the two frames is derived in the Appendix of Ref. [84]. At linear order in the boost velocity $\Gamma$,

$$ \nu_S = (1 + \Gamma \cdot \mathbf{e}_K)\nu_K, \quad (A7) $$

$$ \mathbf{e}_S^i = (1 - \Gamma \cdot \mathbf{e}_K) \mathbf{e}_K^i + \Gamma^i, \quad (A8) $$

where $\mathbf{e}_S^i \equiv E_\mu^i \mathbf{e}_S^\mu$ and $\mathbf{e}_K^i \equiv E_\mu^i \mathbf{e}_K^\mu$. Analogously, introducing a vector $A^\mu$ defining the axis of the source, we introduce the decomposition

$$ A^{\mu} = A_{S}(u_{S}^{\mu} + A_{S}^{\mu}) = A_{K}(u_{K}^{\mu} + A_{K}^{\mu}). \quad (A9) $$

We find up to linear order in $\Gamma$

$$ A_S = (1 - \Gamma \cdot A_K)A_K, \quad (A10) $$

$$ A_S^i = (1 + \Gamma \cdot A_K) A_K^i - \Gamma^i, \quad (A11) $$

where $A_S^i \equiv E_\mu^i A_S^\mu$ and $A_K^i \equiv E_\mu^i A_K^\mu$. Defining $\mu_S \equiv \cos \alpha_S = A_S \cdot e_S$ and $\mu_K \equiv \cos \alpha_K = A_K \cdot e_K$, we get at linear order in $\Gamma$

$$ \mu_S = \mu_K - (1 + \mu_K)(\mathbf{e}_K - \mathbf{A}_K) \cdot \Gamma. \quad (A12) $$

This expression relates the angle identifying the orientation of a source in the source rest frame to the angles defining the orientation of the source in the galaxy rest frame. If the source spectrum in the source rest frame has a rotational symmetry such that it depends only on $\mu_S$, then this is not the case in the galaxy frame due to the aberration of the boost relating both frames, as can be seen by the dependence in $\mathbf{e}_K \cdot \Gamma$. Indeed, the spherical coordinates of $\mathbf{e}_K$ in the frame K are obtained from

$$ e_K^i = (E_\mu^i)_K e_K^\mu = (\cos \alpha_K \cos \beta_K, \cos \alpha_K \sin \beta_K, \sin \alpha_K). \quad (A13) $$

## 3. Situation in G

We need to write the axis of the source in the reference frame G, using (A1). We get

$$ A_G^i \equiv R_j^i(0, \beta, \gamma)A_K^j, \quad (A14) $$

$$ A_G = (\sin \beta \cos \gamma, \sin \beta \sin \gamma, \cos \beta). \quad (A15) $$

The line of sight components will transform with (A1). Explicitly,

$$ \begin{aligned} e_G^i &\equiv (E_\mu^i)_G e_G^\mu = (\cos \alpha_G \cos \beta_G, \cos \alpha_G \sin \beta_G, \sin \alpha_G) \\ &= R_j^i(0, \beta, \gamma) e_K^j. \end{aligned} \quad (A16) $$

We now recall the relation between the reference frames G and K and we make use of (A12) to relate the angle between the axis of the source and the line of sight in the source frame to the spherical angles defining the orientation of the source, the line of sight, and the boost velocity in the galaxy frame. We indicate these angles and the magnitude of the boost velocity collectively as $\lambda_G$, with $\lambda_G \equiv (A_G, \Gamma, \mathbf{e}_G)$ and the resulting relation is of the type $\mu_S(\lambda_G)$. We then relate the spectrum in the frame G to the one in the source frame thanks to

$$ \frac{d^3 E_G}{d\nu_G d^2\Omega_G} (\nu_G, \lambda_G, \theta_G) = \left| \frac{d^2\Omega_S}{d^2\Omega_G} \right| \frac{d^3 E_S}{d\nu_S d^2\Omega_S} (\nu_S, \mu_S, \theta_G). \quad (A17) $$

## 4. Spectrum transformation

The proper time of the observer in K (and hence in G) is related to the proper time in S by

$$ dt_K = \frac{1}{\gamma} dt_S = dt_S + \mathcal{O}(\Gamma^2). \quad (A18) $$

It follows that at linear order in $\Gamma$ the two times can be identified.

The spectrum (A17) in the galaxy rest frame can be written in terms of quantities defined in the source rest frame as

$$ \begin{aligned} & \frac{d^3 E_G}{d^2 \Omega_G d \nu_G} (\nu_G, \lambda_G, \theta_G) \\ &= \frac{1}{2\pi} \frac{d^2 E_S}{d \mu_S d \nu_S} (\nu_S, \mu_S, \theta_G) [1 + \mathcal{O}(\Gamma)], \end{aligned} \quad (A19) $$

where we have used that the energy spectrum transforms as a scalar under boosts and substituted Eq. (A12) to assess that $d^2\Omega_K/(2\pi d\mu_S)$ is unity up to linear terms in $\Gamma$.
We insert these results into Eqs. (26) and (29) and make use of Eq. (19). By doing this, we get expressions for the effective luminosity of a galaxy in terms of local quantities defined in the source local frame. Since linear contributions in $\Gamma$ vanish when averaging over $\Gamma$ with the distribution function $f(\theta_G, \Gamma)$, in the final result for the effective luminosity Eqs. (26) and (29) energy, frequency, and time in the source and galaxy frame can be simply identified, i.e. $\nu_G = \nu_S$, $E_G = E_S$, and $t_G = t_S$, and the energy spectrum depends on directions only through $\mu_G = \mu_S$ (up to linear order in $\Gamma$). We find, explicitly up to first order in $\Gamma$,

$$ \mathcal{L}_G^I(\nu_G, \theta_G) = \sum_{(i)} \int d\theta^{(i)} \mathcal{N}^{(i)}(\theta^{(i)}, \theta_G) \\ \times \int d\mu_S \frac{d^3 E_S^{(i)}}{d\mu_S dt_S d\nu_S} (\nu_S, \mu_S, \theta^{(i)}, \theta_G), \quad (A20) $$

and

$$ \mathcal{L}_G^M(\nu_G, \theta_G) = \sum_{(i)} \int d\theta^{(i)} \frac{d\mathcal{N}^{(i)}}{dt_S} (\theta^{(i)}, \theta_G) \\ \times \int d\mu_S \frac{d^2 E_S^{(i)}}{d\mu_S d\nu_S} (\nu_S, \mu_S, \theta^{(i)}, \theta_G). \quad (A21) $$

Summarizing, to get how the energy spectrum and the power emitted by a GW source transform under a boost with velocity $\Gamma$, three effects have to be considered: the Doppler shift of frequencies, aberration effect, and transformation of time. We have computed these effects and we have showed that up to linear order in $\Gamma$, we can identify the energy spectrum and the power in the reference frames of the source and of the galaxy since linear terms in $\Gamma$ coming from the Lorentz transformation cancel once averaged with the velocity distribution. We emphasize that no other effect comes into play. On the largest scale of our problem (cosmological scale), galaxies are considered as pointlike sources, and at this scale the effect of a relativistic beaming is hidden in the distance-duality relation, see e.g. Sec. III.2.4. of [85].

[1] B. P. Abbott et al. (VIRGO and LIGO Scientific Collaborations), Observation of Gravitational Waves from a Binary Black Hole Merger, Phys. Rev. Lett. 116, 061102 (2016).

[2] B. P. Abbott et al. (VIRGO and LIGO Scientific Collaborations), GW151226: Observation of Gravitational Waves from a 22-Solar-Mass Binary Black Hole Coalescence, Phys. Rev. Lett. 116, 241103 (2016).

[3] B. P. Abbott et al. (VIRGO and LIGO Scientific Collaborations), GW170104: Observation of a 50-Solar-Mass Binary Black Hole Coalescence at Redshift 0.2, Phys. Rev. Lett. 118, 221101 (2017).

[4] B. P. Abbott et al. (VIRGO and LIGO Scientific Collaborations), GW170814: A three-detector observation of gravitational waves from a binary black hole coalescence, arXiv:1709.09660.

[5] B. P. Abbott et al. (VIRGO and LIGO Scientific Collaborations), Upper Limits on the Stochastic Gravitational-Wave Background from Advanced LIGO First Observing Run, Phys. Rev. Lett. 118, 121101 (2017); Erratum, Phys. Rev. Lett. 119, 029901(E) (2017).

[6] J. Aasi et al. (VIRGO and LIGO Scientific Collaborations), Improved Upper Limits on the Stochastic Gravitational-Wave Background from 2009–2010 LIGO and Virgo Data, Phys. Rev. Lett. 113, 231101 (2014).

[7] M. Maggiore, Gravitational wave experiments, and early universe cosmology, Phys. Rep. 331, 283 (2000).

[8] B. Allen, The Stochastic gravity wave background: Sources and detection, in Relativistic Gravitation and Gravitational Radiation, Proceedings, School of Physics, Les Houches,

France, September 26–October 6, 1995 (Cambridge Contemporary Astrophysics, Cambridge, 1996), p. 373.

[9] T. L. Smith, E. Pierpaoli, and M. Kamionkowski, A new cosmic microwave background constraint to primordial gravitational waves, Phys. Rev. Lett. 97, 021301 (2006).

[10] S. Henrot-Versille et al., Improved constraint on the primordial gravitational-wave density using recent cosmological data and its impact on cosmic string models, Classical Quantum Gravity 32, 045003 (2015).

[11] R. M. Shannon et al., Gravitational-wave limits from pulsar timing constrain supermassive black hole evolution, Science 342, 334 (2013).

[12] B. Allen and A. C. Ottewill, Detection of anisotropies in the gravitational wave stochastic background, Phys. Rev. D 56, 545 (1997).

[13] N. J. Cornish, Mapping the gravitational wave background, Classical Quantum Gravity 18, 4277 (2001).

[14] S. Mitra, S. Dhurandhar, T. Souradeep, A. Lazzarini, V. Mandic, S. Bose, and S. Ballmer, Gravitational wave radiometry: Mapping a stochastic gravitational wave background, Phys. Rev. D 77, 042002 (2008).

[15] E. Thrane, S. Ballmer, J. D. Romano, S. Mitra, D. Talukder, S. Bose, and V. Mandic, Probing the anisotropies of a stochastic gravitational-wave background using a network of ground-based laser interferometers, Phys. Rev. D 80, 122002 (2009).

[16] J. D. Romano, S. R. Taylor, N. J. Cornish, J. Gair, C. M. F. Mingarelli, and R. van Haasteren, Phase-coherent mapping of gravitational-wave backgrounds using ground-based laser interferometers, Phys. Rev. D 92, 042003 (2015).
[17] J. D. Romano and N. J. Cornish, Detection methods for stochastic gravitational-wave backgrounds: A unified treatment, Living Rev. Relativity **20**, 2 (2017).

[18] B. P. Abbott et al. (VIRGO and LIGO Scientific Collaborations), Directional Limits on Persistent Gravitational Waves from Advanced LIGO's First Observing Run, Phys. Rev. Lett. **118**, 121102 (2017).

[19] C. M. F. Mingarelli, T. Sidery, I. Mandel, and A. Vecchio, Characterizing gravitational wave stochastic background anisotropy with pulsar timing arrays, Phys. Rev. D **88**, 062005 (2013).

[20] S. R. Taylor and J. R. Gair, Searching for anisotropic gravitational-wave backgrounds using pulsar timing arrays, Phys. Rev. D **88**, 084001 (2013).

[21] J. Gair, J. D. Romano, S. Taylor, and C. M. F. Mingarelli, Mapping gravitational-wave backgrounds using methods from CMB analysis: Application to pulsar timing arrays, Phys. Rev. D **90**, 082001 (2014).

[22] T. Regimbau, The astrophysical gravitational wave stochastic background, Res. Astron. Astrophys. **11**, 369 (2011).

[23] P. Peter and J.-P. Uzan, *Primordial Cosmology*, Oxford Graduate Texts (Oxford University Press, New York, 2005).

[24] J. F. Dufaux, A. Bergman, G. N. Felder, L. Kofman, and J.-P. Uzan, Theory and numerics of gravitational waves from preheating after inflation, Phys. Rev. D **76**, 123517 (2007).

[25] A. Vilenkin, Gravitational radiation from cosmic strings, Phys. Lett. **107B**, 47 (1981).

[26] C. J. Hogan and M. J. Rees, Gravitational interactions of cosmic strings, Nature **311**, 109 (1984).

[27] T. Vachaspati and A. Vilenkin, Gravitational radiation from cosmic strings, Phys. Rev. D **31** (1985) 3052.

[28] R. R. Caldwell and B. Allen, Cosmological constraints on cosmic string gravitational radiation, Phys. Rev. D **45**, 3447 (1992).

[29] S. Kuroyanagi, K. Takahashi, N. Yonemaru, and H. Kumamoto, Anisotropies in the gravitational wave background as a probe of the cosmic string network, Phys. Rev. D **95**, 043531 (2017).

[30] C. Caprini, R. Durrer, T. Konstandin, and G. Servant, General properties of the gravitational wave spectrum from phase transitions, Phys. Rev. D **79**, 083519 (2009).

[31] C. Caprini et al., Science with the space-based interferometer eLISA. II: Gravitational waves from cosmological phase transitions, J. Cosmol. Astropart. Phys. 04 (2016) 001.

[32] C. Caprini and R. Durrer, Gravitational wave production: A Strong constraint on primordial magnetic fields, Phys. Rev. D **65**, 023517 (2001).

[33] P. Binetruy, A. Bohe, C. Caprini, and J.-F. Dufaux, Cosmological backgrounds of gravitational waves and eLISA/NGO: phase transitions, cosmic strings and other sources, J. Cosmol. Astropart. Phys. 06 (2012) 027.

[34] A. Buonanno and B. S. Sathyaprakash, Sources of Gravitational Waves: Theory and Observations, arXiv:1410.7832.

[35] B. P. Abbott et al. (VIRGO and LIGO Scientific Collaborations), GW150914: Implications for the Stochastic Gravitational Wave Background from Binary Black Holes, Phys. Rev. Lett. **116**, 131102 (2016).

[36] T. Regimbau, M. Evans, N. Christensen, E. Katsavounidis, B. Sathyaprakash, and S. Vitale, Digging Deeper: Observing

Primordial Gravitational Waves Below the Binary Black Hole Produced Stochastic Background, Phys. Rev. Lett. **118**, 151105 (2017).

[37] V. Mandic, S. Bird, and I. Cholis, Stochastic Gravitational-Wave Background due to Primordial Binary Black Hole Mergers, Phys. Rev. Lett. **117**, 201102 (2016).

[38] I. Dvorkin, J.-P. Uzan, E. Vangioni, and J. Silk, Synthetic model of the gravitational wave background from evolving binary compact objects, Phys. Rev. D **94**, 103011 (2016).

[39] K. Nakazato, Y. Niino, and N. Sago, Gravitational-wave background from binary mergers and metallicity evolution of galaxies, Astrophys. J. **832**, 146 (2016).

[40] I. Dvorkin, E. Vangioni, J. Silk, J.-P. Uzan, and K. A. Olive, Metallicity-constrained merger rates of binary black holes and the stochastic gravitational wave background, Mon. Not. R. Astron. Soc. **461**, 3877 (2016).

[41] E. F. D. Evangelista and J. C. N. Araujo, The gravitational wave background from coalescing compact binaries: A new method, Braz. J. Phys. **44**, 824 (2014).

[42] L. Z. Kelley, L. Blecha, L. Hernquist, and A. Sesana, The gravitational wave background from massive black hole binaries in illustris: Spectral features and time to detection with pulsar timing arrays, Mon. Not. R. Astron. Soc. **471**, 4508 (2017).

[43] K. Crocker, V. Mandic, T. Regimbau, K. Belczynski, W. Gladysz, K. Olive, T. Prestegard, and E. Vangioni, Model of the stochastic gravitational-wave background due to core collapse to black holes, Phys. Rev. D **92**, 063005 (2015).

[44] M. Surace, K. D. Kokkotas, and P. Pnigouras, The stochastic background of gravitational waves due to the f-mode instability in neutron stars, Astron. Astrophys. **586**, A86 (2016).

[45] D. Talukder, E. Thrane, S. Bose, and T. Regimbau, Measuring neutron-star ellipticity with measurements of the stochastic gravitational-wave background, Phys. Rev. D **89**, 123008 (2014).

[46] P. D. Lasky, M. F. Bennett, and A. Melatos, Stochastic gravitational wave background from hydrodynamic turbulence in differentially rotating neutron stars, Phys. Rev. D **87**, 063004 (2013).

[47] K. Crocker, T. Prestegard, V. Mandic, T. Regimbau, K. Olive, and E. Vangioni, Systematic study of the stochastic gravitational-wave background due to stellar core collapse, Phys. Rev. D **95**, 063015 (2017).

[48] I. Kowalska, T. Bulik, and K. Belczynski, Gravitational wave background from population III binaries, Astron. Astrophys. **541**, A120 (2012).

[49] C. J. Moore, R. H. Cole, and C. P. L. Berry, Gravitational-wave sensitivity curves, Classical Quantum Gravity **32**, 015014 (2015).

[50] B. P. Abbott et al. (LIGO Scientific Collaboration), Exploring the sensitivity of next generation gravitational wave detectors, Classical Quantum Gravity **34**, 044001 (2017).

[51] S. R. Taylor et al., Limits on anisotropy in the nanohertz stochastic gravitational-wave background, Phys. Rev. Lett. **115**, 041101 (2015).

[52] A. Sesana, A. Vecchio, and C. N. Colacino, The stochastic gravitational-wave background from massive black hole
binary systems: implications for observations with pulsar timing arrays, Mon. Not. R Astron. Soc. **390**, 192 (2008).

[53] G. Cusin, J.-P. Uzan, and C. Pitrou, Anisotropy of the astrophysical gravitational waves background I: Predictions for the angular power spectrum (unpublished).

[54] B. Carr, F. Kuhnel, and M. Sandstad, Primordial black holes as dark matter, Phys. Rev. D **94**, 083504 (2016).

[55] S. Bird, I. Cholis, J. B. Muoz, Y. Ali-Hamoud, M. Kamionkowski, E. D. Kovetz, A. Raccanelli, and A. G. Riess, Did LIGO Detect Dark Matter?, Phys. Rev. Lett. **116**, 201301 (2016).

[56] S. Clesse and J. Garca-Bellido, Detecting the gravitational wave background from primordial black hole dark matter, arXiv:1610.08479.

[57] A. Raccanelli, E. D. Kovetz, S. Bird, I. Cholis, and J. B. Munoz, Determining the progenitors of merging black-hole binaries, Phys. Rev. D **94**, 023516 (2016).

[58] T. Namikawa, A. Nishizawa, and A. Taruya, Detecting black-hole binary clustering via the second-generation gravitational-wave detectors, Phys. Rev. D **94**, 024013 (2016).

[59] C. R. Contaldi, Anisotropies of gravitational wave backgrounds: A line of sight approach, Phys. Lett. B **771**, 9 (2017).

[60] T. Namikawa, A. Nishizawa, and A. Taruya, Anisotropies of Gravitational-Wave Standard Sirens as a New Cosmological Probe without Redshift Information, Phys. Rev. Lett. **116**, 121302 (2016).

[61] D. Bertacca, A. Raccanelli, N. Bartolo, and S. Matarrese, Cosmological perturbation effects on gravitational-wave luminosity distance estimates, arXiv:1702.01750.

[62] S. Olmez, V. Mandic, and X. Siemens, Anisotropies in the Gravitational-Wave Stochastic Background, J. Cosmol. Astropart. Phys. 07 (2012) 009.

[63] M. Oguri, Measuring the distance-redshift relation with the cross-correlation of gravitational wave standard sirens and galaxies, Phys. Rev. D **93**, 083511 (2016).

[64] J.-P. Uzan, N. Aghanim, and Y. Mellier, Distance duality relation from x-ray and SZ observations of clusters, Phys. Rev. D **70**, 083533 (2004).

[65] I. Etherington, LX. On the definition of distance in general relativity, Philos. Mag. **15**, 761 (1933).

[66] M. Maggiore, *Gravitational Waves, Vol. 1: Theory and Experiments*, Oxford Master Series in Physics (Oxford University Press, New York, 2007).

[67] N. Deruelle and J.-P. Uzan, *Théories de la relativité*, Collection Echelles (Belin, Paris, 2014).

[68] R. K. Sachs, Gravitational waves in general relativity. 6. The outgoing radiation condition, Proc. R. Soc. A **264**, 309 (1961).

[69] R. A. Isaacson, Gravitational radiation in the limit of high frequency. I. The linear approximation and geometrical optics, Phys. Rev. **166**, 1263 (1968).

[70] B. A. Bassett and M. Kunz, Cosmic distance-duality as a probe of exotic physics and acceleration, Phys. Rev. D **69**, 101305 (2004).

[71] G. F. R. Ellis, R. Poltis, J.-P. Uzan, and A. Weltman, Blackness of the cosmic microwave background spectrum as a probe of the distance-duality relation, Phys. Rev. D **87**, 103530 (2013).

[72] C. Bonvin, C. Caprini, R. Sturani, and N. Tamanini, Effect of matter structure on the gravitational waveform, Phys. Rev. D **95**, 044029 (2017).

[73] P. Schneider, J. Ehlers, and E. E. Falco, *Gravitational Lenses* (Springer-Verlag, Berlin, 1992).

[74] M. Bartelmann and P. Schneider, Weak gravitational lensing, Phys. Rep. **340**, 291 (2001).

[75] P. Fleury, J. Larena, and J.-P. Uzan, The theory of stochastic cosmological lensing, J. Cosmol. Astropart. Phys. **11** (2015) 022.

[76] C. Bonvin and R. Durrer, What galaxy surveys really measure, Phys. Rev. D **84**, 063505 (2011).

[77] A. Challinor and A. Lewis, The linear power spectrum of observed source number counts, Phys. Rev. D **84**, 043516 (2011).

[78] N. Kaiser, Clustering in real space and in redshift space, Mon. Not. R Astron. Soc. **227**, 1 (1987).

[79] T. Matsubara, Correlation function in deep redshift space as a cosmological probe, Astrophys. J. **615**, 573 (2004).

[80] F. Montanari and R. Durrer, Measuring the lensing potential with tomographic galaxy number counts, J. Cosmol. Astropart. Phys. **10** (2015) 070.

[81] T. Matsubara, The Correlation function in redshift space: General formula with wide angle effects and cosmological distortions, Astrophys. J. **535**, 1 (2000).

[82] P. H. F. Reimberg, F. Bernardeau, and C. Pitrou, Redshift-space distortions with wide angular separations, J. Cosmol. Astropart. Phys. **01** (2016) 048.

[83] S. R. Green and R. M. Wald, How well is our universe described by an FLRW model?, Classical Quantum Gravity **31**, 234003 (2014).

[84] G. Cusin, C. Pitrou, and J.-P. Uzan, Are we living near the center of a local void?, J. Cosmol. Astropart. Phys. **03** (2017) 038.

[85] P. Fleury, Ph.D. thesis, Paris Institute of Astrophysics, 2015, arXiv:1511.03702.