# Anisotropies in the astrophysical gravitational-wave background: Predictions for the detection of compact binaries by LIGO and Virgo

Alexander C. Jenkins, Mairi Sakellariadou, Tania Regimbau, Eric Slezak

► To cite this version:

Alexander C. Jenkins, Mairi Sakellariadou, Tania Regimbau, Eric Slezak. Anisotropies in the astrophysical gravitational-wave background: Predictions for the detection of compact binaries by LIGO and Virgo. Physical Review D, 2018, 98 (6), pp.063501. 10.1103/PhysRevD.98.063501. hal-01817918

HAL Id: hal-01817918

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

Submitted on 2 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
# Anisotropies in the astrophysical gravitational-wave background: Predictions for the detection of compact binaries by LIGO and Virgo

Alexander C. Jenkins* and Mairi Sakellariadou†

*Theoretical Particle Physics and Cosmology Group, Physics Department, King's College London, University of London, Strand, London WC2R 2LS, United Kingdom*

Tania Regimbau‡

Université Savoie Mont Blanc, CNRS/IN2P3, Laboratoire d'Annecy-le-Vieux de Physique des Particules (LAPP), 74941 Annecy, France and Université Côte d'Azur, Observatoire de la Côte d'Azur, CNRS, Laboratoire Artemis, CS 34229, 06304 Nice Cedex 4, France

Eric Slezak§

Université Côte d'Azur, Observatoire de la Côte d'Azur, CNRS, Laboratoire Lagrange, CS 34229, 06304 Nice Cedex 4, France

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

(Received 6 June 2018; published 5 September 2018)

We develop a detailed anisotropic model for the astrophysical gravitational-wave background, including binary mergers of two stellar-mass black holes, two neutron stars, or one of each, which are expected to be the strongest contributions in the LIGO-Virgo frequency band. The angular spectrum of the anisotropies, quantified by the $C_\ell$ components, is calculated using two complementary approaches: (i) a simple, closed-form analytical expression, and (ii) a detailed numerical study using an all-sky mock light cone galaxy catalogue from the Millennium simulation. The two approaches are in excellent agreement at large angular scales, and differ by a factor of order unity at smaller scales. These anisotropies are considerably larger in amplitude than e.g., those in the temperature of the cosmic microwave background, confirming that it is important to model these anisotropies, and indicating that this is a promising avenue for future theoretical and observational work.

DOI: 10.1103/PhysRevD.98.063501

## I. INTRODUCTION

The first direct detections of gravitational waves (GW) by the Advanced LIGO [1,2] and Advanced Virgo [3] detectors, from the inspiral and merger of pairs of black holes [4–7], and the recent first detection of GWs from the inspiral and merger of a pair of neutron stars [8], have opened a new window to the Universe, its physical processes, and the astrophysical sources that populate it. Each of these detections are associated with individual loud events, but it is expected that there are many more quiet compact binary mergers in the Universe that are too faint or too distant to be individually resolved. Signals from these quiet events superimpose to create a stochastic gravitational-wave background (SGWB), which is discernible from instrumental noise by cross-correlating the output from multiple detectors, and which is expected to be

detected by Advanced LIGO and Advanced Virgo after a few years of operation at design sensitivity [4]. Many other sources, both astrophysical (e.g., rotating neutron stars, supernovae, core collapse formation of black holes or neutron stars) and cosmological (e.g., phase transitions, inflation, cosmic strings), are expected to contribute to the SGWB. The cosmological background can provide unique information about the early stages of our Universe, while the astrophysical background probes the Universe’s evolution since the beginning of star formation.

Gravitational-wave sources with an inhomogeneous spatial distribution lead to a SGWB characterized by preferred directions, and therefore anisotropies. This is analogous to the measured temperature anisotropies in the cosmic microwave background (CMB) radiation. Methods to detect these SGWB anisotropies using radiometer and spherical harmonic techniques have been developed [9–16] and applied to data from Advanced LIGO’s first observational run to set upper limits on the energy flux and energy density of the SGWB as a function of sky position [17]. Recently, Refs. [18–20] have investigated a general framework to model the anisotropies in the SGWB induced by

<footnote>* Alexander.Jenkins@kcl.ac.uk</footnote>

<footnote>† Mairi.Sakellariadou@kcl.ac.uk</footnote>

<footnote>‡ Tania.Regimbau@oca.eu</footnote>

<footnote>§ Eric.Slezak@oca.eu</footnote>
various types of sources, which was applied to the case of cosmic strings in Ref. [20] and to the case of binary black hole (BBH) mergers in Ref. [21].

In what follows we use a result [Eq. (26)] of Ref. [20] to study anisotropies induced in the SGWB by astrophysical sources, namely from compact object mergers within anisotropically distributed galaxies. We assume that the background is stationary, which implies that the statistical properties of the SGWB do not vary over the time scale of the observation. Our analysis follows two distinct approaches. In the former (analytical approach), we use simple analytical functions for the galaxy number density and the galaxy-galaxy two-point correlation function to derive a simple, closed-form expression for the two-point correlation function and its multipole components $C_\ell$. Let us point out that while Ref. [21] also used a two-point correlation function, our study differs in several ways: (i) we include all three different types of merger: binaries consisting of two neutron stars, two stellar-mass black holes, or one of each; (ii) we calculate the kinematic dipole, which we then subtract since it interferes with the anisotropy statistics (as demonstrated in Ref. [20]); (iii) we follow the fiducial model of the LIGO/Virgo collaboration paper [22] rather than the astrophysical model of Refs. [23–25] to estimate the rate of mergers; (iv) we use a nonlinear power law expression Eq. (64) to model the two-point galaxy clustering, rather than basing this on linear transfer functions of the matter overdensity as in Ref. [21]. In the latter (catalogue approach), we make use—for the first time in the SGWB literature—of a realistic mock catalogue of galaxies and employ recipes to infer the production and merger rates of compact objects.

In both approaches, we restrict our attention to the most important sources for LIGO, Virgo, and other planned ground-based detectors. Future space-based detectors such as the Laser Interferometer Space Antenna (LISA) [26] will probe a lower-frequency window of the GW spectrum, and will therefore be sensitive to different populations of astrophysical systems, which will produce a different stochastic background signal. However, much of Sec. II is still valid in this case, and we intend to repeat our analysis for the LISA frequency band in a future work.

Our present study is organized as follows. In Sec. II we derive a general expression for the astrophysical gravitational wave background (AGWB) induced by an anisotropic distribution of galaxies hosting compact binary mergers. We first write the SGWB density parameter $\Omega_{gw}$ in terms of the average number density of galaxies per comoving volume, the galaxy number overdensity, the rate of binary mergers per galaxy, the observer’s peculiar velocity, and the GW strain spectrum of each binary merger. These quantities depend on the star formation rates (SFRs) and metallicities of the galaxies, and the spins and masses of the two binary components. We perform our analysis using the fiducial model of Ref. [22], and thus

consider the merger rate, the distribution of binary parameters, and the gravitational waveforms emitted by the binary mergers as those of Ref. [22]. We then decompose the density parameter as the average (isotropic) part $\Omega_{gw}$ and the GW density contrast $\delta_{gw}$, which encodes the anisotropies. In Sec. III we use analytical functions to describe the average galaxy number density and the galaxy-galaxy two-point correlation function and derive a simple closed-form expression for each of the $C_\ell$ coefficients, up to a frequency-dependent factor which we evaluate numerically. In Sec. IV we follow a more accurate approach and use a realistic mock catalogue of galaxies, hence relaxing some of the assumptions made in the analytical approach. In Sec. V we compare the analytical and numerical approaches, and comment on the imprint to the SGWB from the anisotropic distribution of compact object mergers as compared to those inferred by a cosmic string network [20]. Some technical details are given explicitly in the appendixes.

## II. GENERAL EXPRESSION FOR THE AGWB

The dimensionless density parameter expressing the intensity of a SGWB, with observed frequency between $\nu_o$ and $\nu_o + d\nu_o$ and arriving from an infinitesimal solid angle $d^2\sigma_o$ centered on the direction $\hat{e}_o$, is defined as

$$ \Omega_{gw}(\nu_o, \hat{e}_o) = \frac{1}{\rho_c} \frac{d^3 \rho_{gw}}{d(\ln \nu_o) d^2 \sigma_o} = \frac{8\pi G \nu_o}{3H_o^2} \frac{d^3 \rho_{gw}}{d\nu_o d^2 \sigma_o}, \quad (1) $$

where we have used the customary normalization with respect to the critical density $\rho_c = 3H_o^2/(8\pi G)$. To study anisotropies in the SGWB induced by astrophysical sources, we consider a Friedmann-Lemaître-Robertson-Walker spacetime, and neglect cosmological perturbations, keeping only the anisotropy due to the source density contrast and the dipole induced by the peculiar motion of the observer. Following Ref. [20], we thus have

$$ \Omega_{gw}(\nu_o, \hat{e}_o) = \frac{\pi}{3} (t_H \nu_o)^3 \int_0^\infty dz \frac{1+z}{E(z)} \times \int d\xi \bar{n} R (1 + \delta_n + \hat{e}_o \cdot \mathbf{v}_o) \int_{S^2} d^2 \sigma_s r_s^2 \tilde{h}^2, \quad (2) $$

where $t_H \equiv 1/H_o$ is the Hubble time, $z$ is the redshift, $\xi$ represents the set of source parameters, $\bar{n}$ is the average (homogeneous) source number density per comoving unit volume, $R$ is the rate of gravitational-wave bursts per source, $\delta_n \equiv (n - \bar{n})/\bar{n}$ is the source number overdensity, $\hat{e}_o$ is the observation direction, $\mathbf{v}_o$ is the observer’s peculiar velocity, $\tilde{h}$ is the GW strain spectrum of a burst, and the final integral is over a sphere centered on the source. Note that Eq. (2) has been modified with respect to Eq. (26) of Ref. [20] by integrating over redshift rather than conformal time, and by using the number density per comoving
volume, rather than per physical volume. Assuming the standard flat ΛCDM cosmology, we have

$$E(z) = \frac{H(z)}{H_o} = \sqrt{\Omega_m (1+z)^3 + \Omega_\Lambda}, \quad (3)$$

with $\Omega_m = 0.3065$, $\Omega_\Lambda = 0.6935$, and $H_o = 67.9 \text{ km s}^{-1} \text{ Mpc}^{-1}$.

Note that, while Eq. (26) of Ref. [20] imposed a cutoff time on the line-of-sight integral to remove nearby, loud, infrequent, resolvable sources and isolate the stochastic part of the signal (see Sec. II C of Ref. [20]), this is not done in Eq. (2) above. This is because the background from binary mergers consists of far fewer GW-emitting events than the cosmic string background considered in Ref. [20]. As a result, the problem of distinguishing different events from each other is greatly reduced, and the detector noise becomes the main barrier to resolving events individually. Indeed, a network of future “third-generation” detectors is expected to be able to resolve the vast majority of binary mergers in this frequency window, including more than 99.9% of BBH mergers [22,27,28]. (For a discussion of how such a large ensemble of individually resolved mergers may be used as a cosmological probe, see Refs. [29,30].) The astrophysical background from unresolvable binary mergers therefore depends on the sensitivity of the detector network, and future detectors will be able to greatly reduce the level of the background by resolving binaries out to much larger redshifts. In order to phrase our results in a detector-independent way, we have opted to include all binary merger events, no matter how close they are to the observer.

In what follows, we study the imprint of an anisotropic distribution of binary mergers, within galaxies, on the SGWB. Since any anisotropies on scales smaller than the typical size of a galaxy are inaccessible to us, we treat galaxies as point sources. We therefore interpret $n$ as the number density of galaxies, and $R$ as the rate of binary mergers per galaxy. Let us write

$$\zeta = (\zeta_b, \zeta_g), \quad (4)$$

where $\zeta_g$ are the parameters of the galaxy (mass, age, luminosity, metallicity, etc.) and $\zeta_b$ are the parameters of the compact binary (masses and spins of the components). We then have

$$n = n(z, \hat{e}_o, \zeta_g), \quad \tilde{h} = \tilde{h}(\nu_s, \hat{e}_s, \zeta_b), \quad R = R(z, \zeta_g, \zeta_b), \quad (5)$$

since the galaxy number density is independent of the source parameters, the strain is independent of the galaxy parameters, and the rate per galaxy depends on both. Note that $\tilde{h}$ is a function of the source-frame frequency $\nu_s$, which is related to the observed frequency $\nu_o$ by

$$\nu_s = \nu_o(1+z)[1 + \hat{e}_o \cdot (\nu_g - \nu_o)], \quad (6)$$

where $\nu_g$ and $\nu_o$ are the peculiar velocities of the galaxy and the observer, respectively. We allow each source to emit GWs anisotropically, so that $\tilde{h}$ is also a function of the direction of propagation away from the source’s position, $\hat{e}_s$.

We can capture much of the relevant astrophysical information with just six parameters: the SFR $\psi$ and metallicity $Z$ of the galaxies, and the masses and spins of the two binary components. The compact objects we consider are all the end products of stellar evolution, so we can use the SFR history of a galaxy to calculate its population of compact binaries, and therefore parametrize the rate of mergers of these objects. The formation of massive black holes from stars is inhibited by stellar winds in high-metallicity environments, so we also include the galaxy metallicity as a parameter to account for this. All other parameters describing the galactic environment are neglected. We further assume that the compact binary orbits are circular and have nonprecessing spins (this is expected to be the case for mergers resulting from isolated binary evolution—such binaries are expected to circularize long before merger due to gravitational backreaction).

Let $\psi(z)$ be the SFR of a galaxy at redshift $z$ in units of $M_{\odot}yr^{-1}$, and $Z$ be the metallicity of the galaxy (i.e., the fraction of the galaxy’s mass that is in elements heavier than helium). In order to use the SFR to parametrize the merger rate, we must take into account the time delay between stars being formed, evolving to become compact objects, and eventually merging with each other. These delay times are typically much longer than the timescales over which the SFR of a galaxy changes, so cannot be neglected. We therefore define the delayed SFR,

$$\psi_d(z) = \int_0^\infty dt_d p(t_d) \psi(z_f), \quad (7)$$

where $z_f(z, t_d)$ is the redshift at which the stars are formed, at a time $t_d$ before the merger occurs at redshift $z$. This is the convolution of the SFR with the probability distribution for the delay times, which varies depending on the objects we consider. Certainly, not all stars become merging compact objects, so $\psi_d$ is not equal to the merger rate—but it is proportional to it. It is therefore much simpler to use $\psi_d$ as a galaxy parameter, rather than $\psi$. Similarly, rather than using the metallicity $Z$, it is more convenient to use a logarithmic scaling, and to normalize relative to the solar metallicity $Z_{\odot} \approx 0.02$. So we define

$$Z \equiv \log_{10} \frac{Z}{Z_{\odot}}, \quad Z \in (-\infty, Z_{\max}], \quad (8)$$

where $Z_{\max} \equiv -\log_{10}Z_{\odot} \approx 1.70$.

Hence we write the galaxy and binary parameter vectors as
$$ \xi_g = (\psi_d, \mathcal{Z}), \quad \xi_b = (m_1, m_2, \chi_1, \chi_2), \quad (9) $$

where $\psi_d$ and $\mathcal{Z}$ are the delayed SFR and log-normalized metallicity introduced above, $m_1$ and $m_2$ are the masses of the two compact objects in a binary (measured in units of $M_\odot$), and

$$ \chi_1 = \frac{S_1}{m_1^2}, \quad \chi_2 = \frac{S_2}{m_2^2}, \quad \chi_1, \chi_2 \in [-1, +1], \quad (10) $$

are the dimensionless spin parameters of the compact objects (where $S_i$ is the spin angular momentum). Using these parameters, we consider three different types of merger: binaries consisting of two neutron stars, two stellar-mass black holes, or one of each. We call these BNS, BBH, and BHNS, respectively. Each type of binary has a different average merger rate, a different distribution over the parameters $\xi_b$, and a different distribution for the delay time $t_d$ (giving a different delayed SFR $\psi_d$).

Equation (2) therefore becomes

$$ \Omega_{\text{gw}}(\nu_o, \hat{e}_o) = \sum_i \frac{\pi}{3} (t_H \nu_o)^3 \int_0^{z_{\text{max}}} dz \frac{1+z}{E(z)} \\ \times \int d\xi_g \bar{n}(1+\delta_n + \hat{e}_o \cdot \mathbf{v}_o) \int d\xi_b R_i S_i, \quad (11) $$

where index $i$ runs over BBH, BNS, and BHNS, and for brevity we have introduced the variable

$$ S_i(\nu_s, \xi_b) \equiv \int_{S^2} d^2 \sigma_s r_s^2 \tilde{h}_i^2. \quad (12) $$

We have set an upper limit $z_{\text{max}} = 10$ on the redshift integral, reflecting the fact that the merger rate is essentially 0 at redshifts greater than this. The galaxy and binary parameter integration measures are

$$ \int d\xi_g = \int_0^{+\infty} d\psi_{d,i} \int_{-\infty}^{Z_{\text{max}}} d\mathcal{Z}, \quad (13) $$

$$ \int d\xi_b = \int_0^\infty dm_1 \int_0^\infty dm_2 \int_{-1}^{+1} \chi_1 \int_{-1}^{+1} d\chi_2, \quad (14) $$

where we have allowed the delayed SFR to be different for different types of binary, by allowing a different distribution of delay times for each case,

$$ \psi_{d,i}(z) \equiv \int_0^\infty dt_d p_i(t_d) \psi(z_f). \quad (15) $$

For the $\Lambda$CDM cosmology we consider, the formation redshift $z_f$ of an object that takes a time $t_d$ to eventually merge at redshift $z$ is given by

$$ 1 + z_f(z, t_d) = (1 + z) \left[ \cosh\left(\frac{3\Omega_\Lambda^{1/2}t_d}{2t_H}\right) - \frac{E(z)}{\Omega_\Lambda^{1/2}}\sinh\left(\frac{3\Omega_\Lambda^{1/2}t_d}{2t_H}\right) \right]^{-2/3}. \quad (16) $$

## A. Intragalactic details of the model

In what follows we describe in detail how we model intragalactic quantities (i.e., those that are relevant within each galaxy)—the merger rate, the distribution of binary parameters, and the gravitational waveforms emitted by the binary mergers. This is almost identical to the fiducial model in Ref. [22], to which we refer the reader for a more thorough discussion and justification of the various modeling choices. The only differences from Ref. [22] are (i) the inclusion of BHNS mergers, for completeness, and (ii) the use of the merger rate per galaxy, rather than per comoving volume (this simplifies our later analysis regarding the anisotropies in the background). The modeling of intergalactic quantities—the galaxy number density $n$ and its clustering statistics—is different in our two approaches, and is discussed in Secs. III and IV.

First, in order to calculate $S_i$, we use the hybrid wave-form models of Refs. [31,32]. These are valid for BBH, and when integrated over the sphere give

$$ S_{\text{BBH}} = \int_{S^2} d^2 \sigma_s r_s^2 \tilde{h}_{\text{BBH}}^2 = \frac{5(GM)^{5/3}}{6\pi^{1/3}} \times \begin{cases} \nu_s^{-7/3} [1 + \sum_{i=2}^3 \alpha_i (\pi GM \nu_s)^{i/3}]^2, & \nu_s < \nu_1 \\ c_1 \nu_s^{-4/3} [1 + \sum_{i=1}^2 \epsilon_i (\pi GM \nu_s)^{i/3}]^2, & \nu_1 \le \nu_s < \nu_2 \\ c_2 [1 + (\frac{\nu_s - \nu_2}{\nu_3})^2]^{-2}, & \nu_2 \le \nu_s < \nu_4. \end{cases} \quad (17) $$

Here we have defined the total mass and chirp mass,

$$ M \equiv m_1 + m_2, \quad \mathcal{M} \equiv \frac{(m_1 m_2)^{3/5}}{M^{1/5}}. \quad (18) $$

Recall that the source-frame frequency $\nu_s$ is a function of the redshift and the peculiar velocities [as given by Eq. (6)], which means that $S_i$ implicitly depends on these as well. The frequencies $\nu_1, \nu_2, \nu_3, \nu_4$ are numerical constants for each system, given by

$$ \nu_i = \frac{1}{\pi GM} \left[ \nu_i^{(0)}(\chi) + \sum_{j=1}^{3} \sum_{k=0}^{3-j} y_i^{(jk)} \left(\frac{M}{M}\right)^{5j/3} \chi^k \right], \quad (19) $$

where the spin parameter $\chi$ is defined as

$$ \chi \equiv \frac{m_1}{M} \chi_1 + \frac{m_2}{M} \chi_2. \quad (20) $$

For the BNS and BHNS cases, we truncate the waveform at the merger frequency $\nu_1$ and ignore any higher frequencies, since the merger and ringdown phases of the above expression are only valid for BBH. (BNS mergers occur
at frequencies above the LIGO-Virgo band anyway, so for observational purposes this has little effect. However, it should be possible to extend our model to higher frequencies by using waveforms that account for neutron-star (NS) matter effects. This may be desirable as future detectors improve our sensitivity to GWs at these frequencies). The spectrum is therefore determined for all types of binaries by the numerical constants $\alpha_i$, $\epsilon_i$, $y_i^{(jk)}$ (which are derived from fits between numerical simulations and post-Newtonian expansions), and $c_i$ (which are chosen to ensure the waveform is continuous). The values used for each of these match those in Refs. [31,32], and are given in Appendix B.

For the binary merger rate per galaxy, $R_i$, we have

$$R_i(z, \xi_g, \xi_b) = \epsilon f_Z p_i(\xi_b) \psi_{d,i}, \quad (21)$$

where $p_i(\xi_b) = p_i(m_1, m_2, \chi_1, \chi_2)$ is the joint probability distribution for the masses and spins for a compact binary of type i. The fact that not all stars ultimately end up in compact binaries is accounted for by multiplying the delayed SFR by the efficiency factor $\epsilon$ (which is some unknown constant), as well as the factor $f_Z$, which accounts for the suppression of massive BH formation in high-metallicity environments. The factor $f_Z$ enforces the assumption that black holes with mass greater than 30 $M_\odot$ can only be formed in galaxies with

$$Z \le \frac{1}{2} Z_{\odot} \Rightarrow Z \le \log_{10} \frac{1}{2} \approx -0.301. \quad (22)$$

Hence one considers a rate correction factor,

$$f_Z(Z, m_1, m_2) = \begin{cases} 1, & m_1, m_2 < 30 M_{\odot} \\ \Theta(\log_{10} \frac{1}{2} - Z), & \text{else} \end{cases} \quad (23)$$

$$p_{\text{BNS}} \propto \delta(\chi_1)\delta(\chi_2),$$

$$p_{\text{BBH}} \propto m_1^{-2.35}/(m_1 - 5 M_{\odot}),$$

$$p_{\text{BHNS}} \propto m_1^{-2.35}\delta(\chi_2),$$

$$1 M_{\odot} \le m_1 \le 2 M_{\odot},$$

$$5 M_{\odot} \le m_2 \le m_1 \le 95 M_{\odot},$$

$$5 M_{\odot} \le m_1 \le 95 M_{\odot},$$

$$1 M_{\odot} \le m_2 \le 2 M_{\odot},$$

$$m_1 + m_2 \le 100 M_{\odot},$$

$$1 M_{\odot} \le m_2 \le 2 M_{\odot},$$

with a normalizing constant chosen appropriately in each case such that $\int d\xi_b p_i(\xi_b) = 1$.

The rate $R_i$ is not yet fully determined, as the constant $\epsilon$ is unknown. However, we can eliminate this constant using the local (i.e., redshift zero) rates inferred by LIGO/Virgo. We do this by calculating the total rate of mergers of type i per unit comoving volume at redshift $z$,

$$\begin{aligned} \mathcal{R}_i(z) &= \int d\xi_g \bar{n} \int d\xi_b R_i(z, \xi_g, \xi_b) \\ &= \epsilon \int d\xi_g \bar{n} \int d\xi_b f_Z p_i(\xi_b) \psi_{d,i}. \end{aligned} \quad (27)$$

where $\Theta(x)$ is the Heaviside step function. The delay time distribution is modeled as $p(t_d) \propto 1/t_d$ between some minimum delay time $t_{\text{min},i}$ (20 Myr for BNS, 50 Myr for BBH and BHNS) and the maximum $t_{\text{max}}$ equal to the age of the Universe at redshift $z$, namely,

$$t(z) = \frac{2t_H}{3\Omega_{\Lambda}^{1/2}} \operatorname{arcsinh}\left(\sqrt{\frac{\Omega_{\Lambda}}{\Omega_{\text{m}}(1+z)^3}}\right). \quad (24)$$

At redshift zero this reduces to $t_0 \approx 0.958t_H$, the current age of the Universe. The delayed SFR is therefore

$$\psi_{d,i} = \frac{1}{\ln(t(z)/t_{\text{min},i})} \int_{t_{\text{min},i}}^{t(z)} d(\ln t_d) \psi(z_f), \quad (25)$$

where the prefactor gives the appropriate normalization of the probability distribution.

It remains to specify the source parameter distributions $p_i(\xi_b)$. We take neutron stars as having masses uniformly distributed between 1 $M_\odot$ and 2 $M_\odot$, with zero spin. The latter is motivated by pulsar observations and by parameter estimation on GW170817, which both indicate that BNS systems typically have low spins ($\chi \lesssim 0.05$) by the time they merge [33] (in any case, we find that spin has little effect on the final results). For black holes, the primary mass is given by a Salpeter initial mass function $\propto m_1^{-2.35}$ between 5 $M_\odot$ and 95 $M_\odot$. If the secondary mass is a black hole, then it is distributed $\propto 1/(m_1 - 5 M_\odot)$. We require the sum of the masses to be less than or equal to 100 $M_\odot$, and take the black hole spins as uniform between -1 and +1. This information can be summarized by writing

We require that in the limit $z \to 0$ this matches the local rates inferred by LIGO/Virgo, so $\mathcal{R}_i(0) = \mathcal{R}_i^{(\text{local})}$. The most up-to-date values are [22]

$$\mathcal{R}_{\text{BNS}}^{(\text{local})} = 1.54 \times 10^{-6} \text{ Mpc}^{-3} \text{ yr}^{-1},$$

$$\mathcal{R}_{\text{BBH}}^{(\text{local})} = 1.03 \times 10^{-7} \text{ Mpc}^{-3} \text{ yr}^{-1}. \quad (28)$$

LIGO and Virgo have not yet observed any BHNS events, so it is only possible to place an upper limit of [34]

$$\mathcal{R}_{\text{BHNS}}^{(\text{local})} \le 3.60 \times 10^{-6} \text{ Mpc}^{-3} \text{ yr}^{-1}. \quad (29)$$
Based on previous population synthesis studies, the true BHNS rate is expected to be somewhere between the BNS and BBH rates (e.g., see Ref. [35] for a review), and therefore roughly an order magnitude less than this upper limit. In order to remain agnostic about the true value of the BHNS rate, we consider two cases: one where we include BHNS mergers at the maximal rate Eq. (29), and one where we set the BHNS rate to 0.

Matching to these local rates allows us to eliminate $\epsilon$ to find

$$R_i = \frac{\mathcal{R}_i^{(\text{local})}}{\mathcal{I}_i} p_i(\zeta_b) f_Z \psi_{d,i}, \quad (30)$$

where we have defined the normalizing constants

$$\mathcal{I}_i = \int d\zeta_g (\bar{n}\psi_{d,i})|_{z=0} \int d\zeta_b p_i(\zeta_b) f_Z. \quad (31)$$

This is particularly simple in the BNS case, since $f_Z = 1$, so $\mathcal{I}_{\text{BNS}} = \int d\zeta_g (\bar{n}\psi_{d,i})|_{z=0}$.

Our expression for the SGWB density parameter therefore becomes

$$\Omega_{\text{gw}} = \sum_i \frac{\pi \mathcal{R}_i^{(\text{local})}}{3\mathcal{I}_i} (t_H \nu_o)^3 \int_0^{z_{\text{max}}} dz \frac{1+z}{E(z)} \\ \times \int d\zeta_g \psi_{d,i} \bar{n} (1 + \delta_n + \hat{\mathbf{e}}_o \cdot \mathbf{v}_o) \int d\zeta_b p_i(\zeta_b) f_Z S_i. \quad (32)$$

## B. Decomposing the background

We can decompose the density parameter $\Omega_{\text{gw}}$ as

$$\Omega_{\text{gw}} \equiv \bar{\Omega}_{\text{gw}}(1 + \delta_{\text{gw}}), \quad (33)$$

where $\bar{\Omega}_{\text{gw}}$ is the average (isotropic) value of the density parameter over the sky, and $\delta_{\text{gw}}$ is the GW density contrast, which encodes the anisotropies. The latter can itself be decomposed as

$$\delta_{\text{gw}} \equiv \delta_{\text{gw}}^{(s)} + \mathcal{D} \hat{\mathbf{e}}_o \cdot \hat{\mathbf{v}}_o, \quad (34)$$

where $\delta_{\text{gw}}^{(s)}$ is the density contrast due to the true cosmological anisotropies, and the latter term is the kinematic dipole, with direction $\hat{\mathbf{v}}_o \equiv \mathbf{v}_o/|\mathbf{v}_o|$ and magnitude given by the dipole factor $\mathcal{D}$. This factor can be calculated by performing a Taylor expansion in $x \equiv 1 + \hat{\mathbf{e}}_o \cdot \mathbf{v}_o$ around $x=1$, giving

$$\mathcal{D}(\nu_o) \equiv v_o \bar{\Omega}_{\text{gw}}^{-1} \left. \frac{\partial \Omega_{\text{gw}}}{\partial x} \right|_{x=1, \delta_n=0}, \quad (35)$$

where $v_o \equiv |\mathbf{v}_o|$. For our particular case, we therefore have

$$\bar{\Omega}_{\text{gw}} = \sum_i \frac{\pi \mathcal{R}_i^{(\text{local})}}{3\mathcal{I}_i} (t_H \nu_o)^3 \int_0^{z_{\text{max}}} dz \frac{1+z}{E(z)} \int d\zeta_g \psi_{d,i} \bar{n} \int d\zeta_b p_i(\zeta_b) f_Z S_i, \quad (36)$$

$$\delta_{\text{gw}}^{(s)} = \bar{\Omega}_{\text{gw}}^{-1} \sum_i \frac{\pi \mathcal{R}_i^{(\text{local})}}{3\mathcal{I}_i} (t_H \nu_o)^3 \int_0^{z_{\text{max}}} dz \frac{1+z}{E(z)} \int d\zeta_g \psi_{d,i} \bar{n} \delta_n \int d\zeta_b p_i(\zeta_b) f_Z S_i, \quad (37)$$

$$\mathcal{D} = v_o \bar{\Omega}_{\text{gw}}^{-1} \sum_i \frac{\pi \mathcal{R}_i^{(\text{local})}}{3\mathcal{I}_i} (t_H \nu_o)^3 \int_0^{z_{\text{max}}} dz \frac{1+z}{E(z)} \int d\zeta_g \psi_{d,i} \bar{n} \int d\zeta_b p_i(\zeta_b) f_Z \frac{\partial}{\partial x}(x S_i). \quad (38)$$

The anisotropies are then characterized by the overdensity field $\delta_{\text{gw}}^{(s)}$, either directly or in terms of its statistics. One particularly useful statistical descriptor is the two-point correlation function (2PCF), defined as the second moment of the overdensity field,

$$C_{\text{gw}}(\theta_o, \nu_o) \equiv \langle \delta_{\text{gw}}^{(s)}(\nu_o, \hat{\mathbf{e}}_o) \delta_{\text{gw}}^{(s)}(\nu_o, \hat{\mathbf{e}}'_o) \rangle, \quad (39)$$

where $\theta_o \equiv \cos^{-1}(\hat{\mathbf{e}}_o \cdot \hat{\mathbf{e}}'_o)$, and the angle brackets denote an averaging over all pairs of directions $\hat{\mathbf{e}}_o, \hat{\mathbf{e}}'_o$ whose angle of separation is $\theta_o$, as well as an ensemble averaging over possible random realizations of the SGWB. The first

moment (i.e., mean) vanishes by definition, and if the background is a Gaussian random field (GRF) then all higher moments either vanish or are expressed in terms of the second moment by Wick's theorem. The 2PCF therefore uniquely characterizes the anisotropies in the Gaussian part of the background. It is common practice (particularly in the CMB literature) to perform a multipole expansion of the 2PCF,

$$C_{\text{gw}}(\theta_o, \nu_o) = \sum_{\ell=0}^{\infty} \frac{2\ell + 1}{4\pi} C_{\ell}(\nu_o) P_{\ell}(\cos \theta_o), \quad (40)$$
where $P_\ell(x)$ denotes the $\ell$th Legendre polynomial. The anisotropies are then described in terms of the $C_\ell$ components, which are given by

$$C_\ell(\nu_o) = 2\pi \int_{-1}^{+1} d(\cos \theta_o) C_{gw}(\theta_o, \nu_o) P_\ell(\cos \theta_o). \quad (41)$$

The quantity $\ell(\ell+1)C_\ell/2\pi$ is roughly the contribution to the anisotropic variance of $\delta_{gw}^{(s)}$ per logarithmic bin in $\ell$, as can be seen by considering

$$\begin{aligned} \operatorname{var}(\delta_{gw}^{(s)}) &= C_{gw}(\theta_o = 0) = \sum_{\ell} \frac{2\ell + 1}{4\pi} C_{\ell} \\ &\approx \int d(\ln \ell) \frac{\ell(\ell + 1)C_{\ell}}{2\pi}. \end{aligned} \quad (42)$$

While the overdensity $\delta_{gw}^{(s)}$ is a random field, the $C_\ell$'s are treated as deterministic quantities, averaged over some ensemble of realizations of the SGWB (e.g., as observed at different locations in the Universe). As a result, any measurement of the $C_\ell$'s has some cosmological uncertainty associated with it, due to the fact that we can only access a single realization of the SGWB. This is called the cosmic variance, and is given by

$$\operatorname{var}(C_\ell) = \frac{2}{2\ell + 1} C_\ell^2. \quad (43)$$

### C. Non-Gaussianity in the AGWB

As mentioned above, the analysis of the SGWB is simplified if it is a GRF, as this eliminates the need for anything other than the second moment of the overdensity (i.e., the 2PCF). In Ref. [20], it was found that a GW background composed of independent sources and discretized into $N_{pix}$ pixels on the sky is a GRF at frequency $\nu$ if it satisfies

$$\nu T \gg \frac{N_{\text{pix}}}{\Lambda} + \frac{N_{\text{pix}}^2}{\Lambda^2}, \quad (44)$$

where $T$ is the observation time and $\Lambda$ is the average number of signals in-band at any given moment (this is referred to as the duty cycle). However, this was derived by assuming that the duration of a signal at frequency $\nu$ is roughly $1/\nu$, which is a good approximation for burstlike signals, but is inaccurate for the chirp signals emitted by coalescing compact binaries. As mentioned in Ref. [20], the duration of a compact binary signal in some small frequency interval $[\nu, \nu + \delta\nu]$ is given by

$$\Delta t \approx \left( \frac{96}{5} \pi^{8/3} M^{5/3} \nu^{11/3} \right)^{-1} \delta \nu, \quad (45)$$

where $M$ is the chirp mass of the binary, as defined in Eq. (18). Taking this into account, the appropriate limit on the observing time becomes

$$T \gg \frac{N_{\text{pix}}}{R} + \frac{N_{\text{pix}}^2}{R^2 \Delta t}, \quad (46)$$

where $R$ is the average rate of arrival of GW signals.

If we focus on the AGWB at a single frequency, then we are choosing a frequency interval $\delta\nu$ equal to the frequency resolution $1/T$. The above implies that the background is only a GRF at this frequency if

$$\frac{96}{5} \pi^{8/3} M^{5/3} \nu^{11/3} \frac{N_{\text{pix}}^2}{R^2} \ll 1, \quad (47)$$

which is impossible for the sources considered here. It is therefore inconsistent to assume Gaussianity when considering the AGWB at a single frequency. (We note in passing that due to the large exponent on the frequency in the equation above, this requirement is much easier to fulfil in the LISA frequency band.)

There are two possible ways of addressing this: (i) we can choose a larger frequency interval, integrate the signal over this interval, and treat the result as a GRF; (ii) we can characterize the anisotropies at a single frequency by computing higher-order correlators as well as the 2PCF (e.g., the AGWB bispectrum and trispectrum). The former is computationally very expensive, due to the required sampling of many points in GW frequency space (however, it is worth noting that LIGO/Virgo stochastic searches typically integrate the data over a frequency bin that is significantly larger than $1/T$ anyway, so it may desirable to compute predictions for this approach regardless of considerations about Gaussianity). The latter requires a more detailed study, but has the potential to reveal more detailed astrophysical and cosmological information. This will be investigated in a future work.

For the purposes of this work, we focus on the 2PCF computed at a single frequency, as this still contains a great deal of important information, and is the first step in either of the approaches described above.

### III. ANALYTICAL APPROACH

In this section, we use simple analytical functions for the galaxy number density $\bar{n}$ and galaxy-galaxy 2PCF $\langle \delta_n \delta_n \rangle$. This allows us to derive a simple closed-form expression for each of the $C_\ell$'s, up to a frequency-dependent factor $A_{gw}$ which can be integrated numerically.

#### A. Galaxy distribution

In Sec. II we have not yet specified the mean number density distribution $\bar{n}$. For this initial analytical approach, we assume that the two galaxy parameters are
independent—i.e., that there is no correlation between the metallicity of a galaxy and its delayed SFR. This allows us to write

$$ \bar{n}(z, \psi_d, \mathcal{Z}) = \bar{N}(z) p(\psi_d|z) p(\mathcal{Z}|z), \quad (48) $$

where $p(\psi_d|z)$ and $p(\mathcal{Z}|z)$ are redshift-dependent probability distributions for the two parameters, and

$$ \bar{N}(z) = \int d\zeta_g \bar{n} \quad (49) $$

is the total number density of galaxies per comoving volume. With this assumption, the density parameter $\Omega_{gw}$ reads

$$ \begin{aligned} \Omega_{gw} = & \sum_i \frac{\pi R_i^{(\text{local})}}{3\mathcal{I}_i} (t_H \nu_o)^3 \\ & \times \int_0^{z_{\text{max}}} dz \frac{1+z}{E(z)} \bar{N}(z) \left( \int_0^{\infty} d\psi_{d,i} p(\psi_{d,i}|z) \psi_{d,i} \right) \\ & \times \int_{-\infty}^{z_{\text{max}}} d\mathcal{Z} p(\mathcal{Z}|z) \int d\zeta_b p_i(\zeta_b) f_Z S_i(1 + \delta_n + \hat{e}_o \cdot \mathbf{v}_o). \end{aligned} \quad (50) $$

Notice that since the integrand is proportional to $\psi_d$, we only need the first moment (i.e., the mean) of the corresponding probability distribution,

$$ \bar{\psi}_d = \int_0^\infty d\psi_d p(\psi_d|z)\psi_d, \quad (51) $$

and there is no need to specify anything else about the distribution. The quantity in Eq. (51) is just the mean delayed SFR, averaged across all galaxies at the appropriate redshifts. Equivalently, we can think of this as the delayed form of the cosmic mean SFR per galaxy $\bar{\psi}$,

$$ \bar{\psi}_d = \frac{1}{\ln(t(z)/t_{\text{min}})} \int_{t_{\text{min}}}^{t(z)} d(\ln t_d) \bar{\psi}(z_f). \quad (52) $$

Expressions for the cosmic mean SFR are more appropriately given in units per unit comoving volume rather than per galaxy, so we define

$$ \bar{\psi}^{(V)}(z) = \bar{N}(z) \bar{\psi}(z), \quad (53) $$

with the mean number of galaxies per comoving volume $\bar{N}$ converting between the two. This function is given in Ref. [36] by

$$ \bar{\psi}^{(V)}(z) = \bar{\psi}_{\text{peak}}^{(V)} \frac{\alpha \exp[\beta(z - z_{\text{peak}})]}{\alpha - \beta + \beta \exp[\alpha(z - z_{\text{peak}})]}, \quad (54) $$

with dimensionless constants $\alpha = 2.80$, $\beta = 2.62$, and with the normalization relative to the peak value

$\bar{\psi}_{\text{peak}}^{(V)} = 0.145 M_{\odot}\text{yr}^{-1}\text{Mpc}^{-3}$ at redshift $z_{\text{peak}} = 1.86$. We thus write the mean of the delayed SFR distribution as

$$ \bar{\psi}_d = \frac{1}{\ln(t(z)/t_{\text{min}})} \int_{t_{\text{min}}}^{t(z)} d(\ln t_d) \frac{\bar{\psi}^{(V)}(z_f)}{\bar{N}(z_f)}. \quad (55) $$

No further information or assumptions about the SFR distribution are needed. We do, however, need an expression for the total galaxy number density $\bar{N}(z)$, in order to convert the SFR per comoving volume $\bar{\psi}^{(V)}$ back into SFR per galaxy. We use Ref. [37], where a fit to observational data over a redshift range $0 \le z \le 8$ has been performed, to get

$$ \log_{10} \left( \frac{\bar{N}(z)}{1 \text{ Mpc}^{-3}} \right) = -0.26 - 1.08 \log_{10} \left( \frac{t(z)}{1 \text{ Gyr}} \right), \quad (56) $$

where $t(z)$ is the age of the Universe at redshift $z$, as given by Eq. (24).

The dependence of the integrand in Eq. (50) on the metallicity is more complicated due to the factor $f_Z$; we must therefore specify the full distribution $p(\mathcal{Z}|z)$, and not just its mean. Following Ref. [22], we take the metallicity distribution as a Gaussian with variance 1/4 centered on the cosmic mean metallicity at that redshift, namely,

$$ p(\mathcal{Z}|z) \propto \exp[-2(\mathcal{Z} - \bar{\mathcal{Z}}(z))^2]. \quad (57) $$

Note that this is not strictly a Gaussian, as $\mathcal{Z}$ has a finite maximum value of $Z_{\text{max}}$. The distribution Eq. (57) must be reweighted accordingly to give

$$ p(\mathcal{Z}|z) = \frac{\sqrt{\frac{8}{\pi}} \exp[-2(\mathcal{Z} - \bar{\mathcal{Z}}(z))^2]}{1 - \text{erf}[\sqrt{2}(\bar{\mathcal{Z}}(z) - Z_{\text{max}})]}. \quad (58) $$

The cosmic mean metallicity is well modeled by [38]

$$ \begin{aligned} \bar{\mathcal{Z}}(z) &= \log_{10} \frac{\bar{\mathcal{Z}}(z)}{Z_{\odot}} \\ &= Z_{\text{max}} + \log_{10} \left( \frac{1}{\psi_{\text{norm}}^{(V)}} \int_{z}^{z_{\text{max}}} dz' \frac{\bar{\psi}^{(V)}(z')}{(1+z')E(z')} \right), \end{aligned} \quad (59) $$

with $\bar{\psi}^{(V)}$ being the mean SFR from Eq. (54) and $\psi_{\text{norm}}^{(V)}$ a normalizing constant with units of $M_{\odot}\text{yr}^{-1}\text{Mpc}^{-3}$, given by [38]

$$ \psi_{\text{norm}}^{(V)} = \frac{\rho_{\text{bary}} H_o}{\sqrt{10y(1-R)}} = 8.90 M_{\odot}\text{yr}^{-1}\text{Mpc}^{-3}, \quad (60) $$

where $\rho_{\text{bary}}$ is the present-day baryon density, $y = 0.019$ is the net metal yield created by each generation of stars, and $R = 0.27$ is the fraction of metals returned to the interstellar medium by stars.
Given these distributions, the isotropic GW energy density is therefore

$$ \bar{\Omega}_{\text{gw}} = \sum_{i} \frac{\pi R_{i}^{(\text{local})}}{3\mathcal{I}_{i}} (t_{H}\nu_{o})^{3} \int_{0}^{z_{\text{max}}} \frac{dz(1+z)}{E(z) \ln(t(z)/t_{\text{min},i})} \\ \times \left( \int_{t_{\text{min},i}}^{t(z)} d(\ln t_{d}) \bar{\psi}^{(V)}(z_{f}) \frac{\bar{N}(z)}{\bar{N}(z_{f})} \right) \int_{-\infty}^{z_{\text{max}}} d\mathcal{Z} p(\mathcal{Z}|z) \int d\xi_{b} p_{i}(\xi_{b}) f_{\mathcal{Z}} S_{i}, \quad (61) $$

with the rate normalizing factor given by

$$ \mathcal{I}_i = \left[ \frac{1}{\ln(t(z)/t_{\text{min},i})} \left( \int_{t_{\text{min},i}}^{t(z)} d(\ln t_d) \bar{\psi}^{(V)}(z_f) \frac{\bar{N}(z)}{\bar{N}(z_f)} \right) \int_{-\infty}^{z_{\text{max}}} d\mathcal{Z} p(\mathcal{Z}|z) \int d\xi_b p_i(\xi_b) f_{\mathcal{Z}} \right] \Bigg|_{z=0} . \quad (62) $$

The result is shown in Fig. 1.

## B. Anisotropies

We now consider anisotropies in the AGWB. Since we do not know the form of $\delta_n$ as a function of sky position, we resort to a statistical description. If the signal is Gaussian, then a complete and simple description is given by the 2PCF. Even if the signal is non-Gaussian (as suggested by Sec. II C), this is still a convenient and informative descriptor of the anisotropies. We write the 2PCF of the stochastic background as

$$ C_{\text{gw}}(\theta_o, \nu_o) = \frac{\pi^2 (t_H \nu_o)^6}{9 \bar{\Omega}_{\text{gw}}^2} \int_0^{z_{\text{max}}} dz \frac{1+z}{E(z)} \int_0^{z_{\text{max}}} dz' \frac{1+z'}{E(z')} \left( \sum_i \int d\xi_g \bar{n}(z, \xi_g) \int d\xi_b R_i(z, \xi_g, \xi_b) S_i(\nu_s, \xi_b) \right) \\ \times \left( \sum_j \int d\xi'_g \bar{n}(z', \xi'_g) \int d\xi'_b R_j(z', \xi'_g, \xi'_b) S_j(\nu'_s, \xi'_b) \right) \langle \delta_n(z, \hat{e}_o, \xi_g) \delta_n(z', \hat{e}'_o, \xi'_g) \rangle, \quad (63) $$

which is determined by the final factor on the rhs, namely, the galaxy-galaxy 2PCF, $\xi_{gg} \equiv \langle \delta_n \delta_n \rangle$. There exists a simple analytical estimate for this,

<figure><img src="image_3.png" /><figcaption>FIG. 1. The total GW flux spectrum $4\pi\bar{\Omega}_{\text{gw}}$ given by numerically integrating Eq. (61). (The factor of $4\pi$ is included to give the total flux, rather than the average flux per solid angle.) The total flux with (without) BHNS mergers is given in black (blue), while the flux coming from BNS, BBH, and BHNS mergers is indicated in red, green, and magenta, respectively. The curves including the BHNS contribution are dashed to indicate that this is an upper bound only.</figcaption></figure>

$$ \xi_{gg}(d) \approx \left(\frac{d}{d_1}\right)^{-\gamma} \delta(z - z'), \quad (64) $$

where $d(z, \hat{e}_o, \hat{e}'_o)$ is the comoving distance between the two points, $d_1$ is the comoving scale at which the correlation is unity, and $\gamma$ is some positive real number. We assume for simplicity that all galaxies cluster in the same way, regardless of their SFR and metallicity, and that $d_1$ and $\gamma$ are constant, and take the values<fnref n="1" /> $d_1 = (4.29 \pm 0.19)h^{-1}$ Mpc and $\gamma = 1.63 \pm 0.04$. Writing the comoving distance $d$ as

$$ d = 2r \tan \frac{\theta_o}{2}, \quad \theta_o \equiv \cos^{-1}(\hat{e}_o \cdot \hat{e}'_o), \quad (65) $$

where $r = \int dz/H$ is the comoving distance between the observer and the galaxy, the 2PCF can be simply written as

$$ C_{\text{gw}}(\theta_o, \nu_o) = A_{\text{gw}}(\nu_o) \left(\tan \frac{\theta_o}{2}\right)^{-\gamma}, \quad (66) $$

where the “amplitude” of the correlation $A_{\text{gw}}$ is a function of the GW frequency alone, and has no angular dependence,

<footnote n="1">¹The quoted values are taken from the VIMOS Public Extragalactic Redshift Survey (VIPERS) [39], and are valid for the brightest galaxies in a redshift bin centered on $z=1$, which are precisely the galaxies that contribute most strongly to $\Omega_{\text{gw}}$.</footnote>
$$A_{\text{gw}}(\nu_o) = \bar{\Omega}_{\text{gw}}^{-2} \left(\frac{d_1}{2}\right)^{\gamma} \int_{0}^{z_{\text{max}}} \frac{dz}{r^{\gamma}} \left(\frac{\partial \bar{\Omega}_{\text{gw}}}{\partial z}\right)^2 . \quad (67)$$

The multipole components $C_{\ell}$ are then given by

$$C_{\ell}(\nu_o) = 2\pi A_{\text{gw}} \int_{-1}^{+1} d(\cos\theta_o) P_{\ell}(\cos\theta_o) \left(\tan\frac{\theta_o}{2}\right)^{-\gamma}, \quad (68)$$

which can be evaluated explicitly (see Appendix A) to give

$$C_{\ell} = 4\pi A_{\text{gw}} \frac{F_2(-\ell, \ell + 1, 1 - \frac{\gamma}{2}; 1, 2; 1)}{\text{sinc}(\pi\gamma/2)}, \quad (69)$$

where $\text{sinc}(x) = \sin(x)/x$, and $F_2$ is a generalized hypergeometric function, which is simple to evaluate numerically.

We emphasize that Eq. (69) is only a simple approximation of the true angular spectrum of the anisotropies. In reality, the galaxy clustering is more complicated than we have assumed, since it evolves with redshift and varies between different populations of galaxies with different physical attributes. We also note that we have extended Eq. (64) beyond its realm of validity by assuming that it holds for all distances $d$. The galaxy 2PCF drops below this power law for distances smaller than $\approx 0.1 h^{-1}$ Mpc and tends to some finite value as $d \to 0$, and also drops below the power law for distances larger than $\approx 10 h^{-1}$ Mpc, eventually becoming negative for large separations. There are two resulting inaccuracies in Eqs. (66) and (69) that are immediately evident. First, the correlation $C_{\text{gw}}$ diverges as $\theta_o \to 0$, even though this should give a well-defined finite value equal to the variance of the field at each point, $\text{var}(\delta_{\text{gw}}^{(s)})$. Secondly, the $\ell = 0$ moment of the 2PCF evaluates to $C_0 = 4\pi A_{\text{gw}} / \text{sinc}(\pi\gamma/2) > 0$, even though $C_0$ should be exactly equal to 0 (as shown in Appendix A). These inaccuracies are due to the extrapolation of Eq. (64) down to $d = 0$ and out to $d \to \infty$, respectively. We therefore turn to a more detailed and accurate approach in the following section, to address the deficiencies of this simple model.

## IV. CATALOGUE APPROACH

The study of anisotropies in the SGWB induced by the large-scale distribution of astrophysical sources can be also performed using a large enough realistic mock catalogue of galaxies with recipes to infer the production and merger rates of compact sources. In such catalogues, the growth of dark matter structure is first simulated and the history of the hierarchical mass assembly is then recorded by means of merger trees of haloes and subhaloes within snapshots stored at different time steps. To simulate the visible galaxies without performing a cosmological hydrodynamic simulation, one uses a semianalytic model aiming to

reproduce the observed properties of these galaxies (e.g., clustering, counts, scaling relations) at different redshifts. Such models (see, e.g., Refs. [40–42]) differ by the number of astrophysical processes included to describe the fate of baryons within the dark matter haloes, the phenomenological description of these processes, and their consistency with the various observational data sets. In that respect, quantitative results strongly depend on the treatment and implementation of feedback mechanisms. Mock light cones can then be generated from the postprocessed snapshots of the dynamical collisionless N-body simulation, taking care to shuffle them in order to suppress replication effects.

In what follows, we make use of a mock light cone catalogue [43] constructed by applying the “L-galaxies” model [40] to the Millennium simulation [44,45]. Very relevant for our study are the all-sky coverage and depth of this mock catalogue, even if the random tiling process used to generate a mock light cone (the size of which is larger than the simulated volume) introduces discontinuities in the density field and hence a loss of clustering information. Hence, the clustering signal suffers from a small negative bias of less than 10% from 1 to $10 h^{-1}$ Mpc around the two-point correlation length. In addition, it exhibits finite volume effects due to the limited size of the simulated box which for instance erases the two-point angular correlation signal at scales larger than a tenth of the simulation box size. However, both effects are not important for our study since only a fraction of the spatial information is lost affecting scales that are not very relevant for the anisotropies we are looking at. The comoving box of size $500 h^{-1}$ Mpc on a side of the Millennium simulation has indeed been chosen that large in order to encompass typical scales related to the large-scale features of the spatial distribution of galaxies.

The all-sky mock light cone catalogue [40,43] we use contains 5,715,694 galaxies, all limited at apparent AB magnitude of 18 in the r filter from SDSS. It has been built using a random tiling technique applied to 64 postprocessed snapshots saved during the Millennium simulation with a time step of about 100 Myr. We queried the database and retrieved for each galaxy the following data: its sky location, cosmological redshift, metallicity, and peculiar velocity. In order to calculate the delayed star formation rate of each galaxy on the light cone, it was necessary to access information about its star formation rate at earlier snapshots by querying the full Millennium simulation. This was made more complicated by the fact that the light cone galaxies are the result of a sequence of mergers of smaller galaxies, each with their own independent star formation history. In order to account for this, we queried the Millennium simulation to extract the full star formation history of each light cone galaxy, given by the star formation rates and redshifts of its progenitors at earlier snapshots. This included a total of 973,224,532 redshift and SFR measurements from the progenitor galaxies.
We note that due to the magnitude-limited sample used to construct the light cone catalogue, it only extends out to a redshift of $z \approx 0.78$. While this is sufficient to study anisotropies on the angular scales we are interested in, it means that the total energy density $\Omega_{gw}$ given by the catalogue will be significantly less than the true value, due to the missing contribution from redshifts $z > 0.78$. However, since we describe the anisotropies in terms of the GW overdensity $\delta_{gw}^{(s)}$, this will have minimal effect on the $C_\ell$ spectrum. The main advantage of the mock catalogue lies in its accurate nonlinear modeling of the galaxy clustering statistics, which remains valid on the scales of interest even when the redshift limit is introduced.

In performing our analysis of the mock catalogue, we adhere to the (now outdated) WMAP values of the cosmological parameters that were used in the Millennium simulation: $H_0 = 73 \text{ kms}^{-1} \text{ Mpc}^{-1}$, $\Omega_m = 0.25$, $\Omega_\Lambda = 0.75$. For comparison with the analytical prediction Eq. (69), we use values of $\gamma$ and $d_1$ that match those of the simulation, which are themselves consistent with the 2-degree Field Galaxy Redshift Survey [44,46]: $\gamma = 1.67 \pm 0.03$, $d_1 = 5.05 \pm 0.26 \text{h}^{-1} \text{ Mpc}$. This ensures consistency between our results and the underlying simulation. We confirm that replacing these with the more up-to-date values mentioned previously ($H_0 = 67.9 \text{ kms}^{-1} \text{ Mpc}^{-1}$, $\Omega_m = 0.3065$, $\Omega_\Lambda = 0.6935$, $d_1 = 4.29 \text{h}^{-1} \text{ Mpc}$, $\gamma = 1.63$) has little impact on the results.

### A. Reconstructing the SGWB from pointlike sources

In order to use the information extracted from the catalogue, we must first explicitly rewrite our equations for $\Omega_{gw}$ in terms of the available data for each galaxy. We do this by expressing the galaxy number density $n$ as a weighted sum of Dirac delta functions, reflecting the fact that we treat each galaxy as a point source.

Let us consider a catalogue of $N$ galaxies, indexed by a label $k$. We write their redshift, sky location, delayed star formation rate, and log-normalized metallicity as $z_k$, $\hat{e}_k$, $\bar{\psi}_{d,k}(z)$, and $Z_k$, respectively. We also allow each galaxy to have a peculiar velocity $v_k$ along the line of sight. This means that galaxy $k$ has a source-frame frequency given by

$$ \nu_{s,k} = \nu_o(1 + z_k)(1 + v_k - \hat{e}_k \cdot \mathbf{v}_o). \quad (70) $$

By integrating the number density per comoving volume $n$ over redshift, SFR, and metallicity, we must have

$$ \begin{aligned} \mathcal{N} &= \int_{z=0}^{z_{\text{max}}} \mathrm{d}^3 V(z) \int \mathrm{d}\zeta_g n \\ &= \int_{S^2} \mathrm{d}^2 \sigma_o \int_0^{z_{\text{max}}} \mathrm{d} \frac{r^2}{H} \int_0^\infty \mathrm{d}\bar{\psi}_d \int_{-\infty}^{Z_{\text{max}}} \mathrm{d}Z n. \end{aligned} \quad (71) $$

This implies that

$$ n = H_0 \sum_k \frac{E(z_k)}{r^2(z_k)} \delta^{(2)}(\hat{e}_o, \hat{e}_k) \delta(z - z_k) \delta(\bar{\psi}_d - \bar{\psi}_{d,k}(z_k)) \delta(Z - Z_k). \quad (72) $$

The SGWB is then given by

$$ \begin{aligned} \Omega_{gw} &= \sum_i \frac{\pi}{3} (t_H \nu_o)^3 \int_0^{z_{\text{max}}} \mathrm{d}z \frac{1+z}{E(z)} \int \mathrm{d}\zeta_g n (1 + \hat{e}_o \cdot \mathbf{v}_o) \int \mathrm{d}\zeta_b R_i(z, \mathcal{Z}, \zeta_b) S_i(\nu_s, \zeta_b) \\ &= \sum_k \sum_i \frac{\pi H_o}{3} (t_H \nu_o)^3 \frac{1+z_k}{r^2(z_k)} (1 + \hat{e}_k \cdot \mathbf{v}_o) \int \mathrm{d}\zeta_b R_i(z_k, \mathcal{Z}_k, \zeta_b) S_i(\nu_{s,k}, \zeta_b) \delta^{(2)}(\hat{e}_o, \hat{e}_k). \end{aligned} \quad (73) $$

Let us remind the reader that $i \in \{\text{BNS, BBH, BHNS}\}$ indexes a sum over the different types of binary merger in each galaxy, while $k \in \{1, 2, \dots, N\}$ indexes a sum over the galaxies in the catalogue. The above expression for $\Omega_{gw}$ becomes much simpler if we define the weight

$$ w_k(\nu_o) = \sum_i \frac{\pi H_o}{3} (t_H \nu_o)^3 \frac{1+z_k}{r^2(z_k)} (1+\hat{e}_k \cdot \mathbf{v}_o) \int \mathrm{d}\zeta_b R_i(z_k, \mathcal{Z}_k, \zeta_b) S_i(\nu_{s,k}, \zeta_b) \quad (74) $$

for each galaxy, so that

$$ \Omega_{gw} = \sum_k w_k \delta^{(2)}(\hat{e}_o, \hat{e}_k). \quad (75) $$

In the analytic approach we took in the previous section, we were limited by the fact that $\delta_n$ was not known as a function of sky location; all we knew was the (approximate) 2PCF of $\delta_n$, so we were forced to use the 2PCF of $\Omega_{gw}$. Using a catalogue we now have an explicit expression for $\Omega_{gw}$ as a function of sky location, and we can use whatever statistics we like to describe it.
One simple choice is to calculate the spherical multipole components of $\Omega_{gw}$,

$$ \Omega_{\ell m}(\nu_o) = \int_{S^2} d^2\sigma_o \Omega_{gw} Y_{\ell m}^*(\hat{e}_o) = \sum_k w_k Y_{\ell m}^*(\hat{e}_k). \quad (76) $$

However, these quantities become expensive to compute if our catalogue contains a large number of galaxies, as each component requires $\mathcal{N}$ evaluations of a spherical harmonic function. This work focuses on the 2PCF of the GW energy overdensity $C_{gw} = \langle \delta_{gw}^{(s)} \delta_{gw}^{(s)} \rangle$, which captures much of the information contained in the anisotropies. It is possible to write down an expression for the $C_\ell$ components of this function in terms of the weights $w_k$ and position vectors $\hat{e}_k$ of each galaxy; however, these are even more expensive to compute, as they require a sum over all pairs of galaxies which scales as $\mathcal{N}^2 - \mathcal{N}$, or more than $3 \times 10^{13}$ operations for each $C_\ell$.

In practice, it is much more efficient to compute the statistics of the background using HEALPix<fnref n="2" /> [47], a powerful software package for manipulating pixelized maps of the sphere. We construct a HEALPix map by computing the weight $w_k$ of each galaxy [given by Eq. (74)] and adding this to the pixel corresponding to the galaxy's sky position $\hat{e}_k$. The $C_\ell$'s (or other relevant quantities) can then be calculated directly in a matter of seconds using HEALPix routines.

## V. RESULTS AND DISCUSSION

We have numerically integrated Eq. (61) to find the frequency spectrum of the isotropic component of the AGWB as shown in Fig. 1, which matches the corresponding figure in Ref. [22]. At frequencies up to ≈20 Hz, all sources are included and are in the inspiral regime, giving a monopole $\bar{\Omega}_{gw}$ that scales as $\nu_o^{2/3}$. The number of emitting sources decreases at higher frequencies, causing the spectrum to taper off and eventually vanish at ≈3000 Hz. Using this spectrum, we are able to compute the factor $A_{gw}$ given in Eq. (67), and hence obtain a simple analytical prediction for the $C_\ell$ spectrum of the background using Eq. (69). The value of this factor as a function of frequency is shown in Fig. 2. In the inspiral-dominated regime below ≈20 Hz, the value of $A_{gw}$ is essentially constant at ≈9 × 10⁻⁵. As the frequency increases above this regime, the falloff in the number of sources increases the granularity of the background, and therefore the size of the anisotropies, causing $A_{gw}$ to increase by many orders of magnitude, and eventually diverge as $\bar{\Omega}_{gw} \to 0$.

Using the above, we can compute the $C_\ell$ at any frequency in the LIGO-Virgo band. For concreteness, we focus on a frequency of 65.75 Hz, as this has been used as

<figure><img src="image_4.png" /><figcaption>FIG. 2. The amplitude factor $A_{gw}$ defined in Eq. (67), as a function of GW frequency $\nu_o$. The lines in red, green, and magenta show what this factor would be if the background included only BNS, only BBH, or only BHNS events, respectively. The blue line shows the combined result, based on the spectra in Fig. 1 (note that the combined curve is not simply the sum of the BBH, BNS and BHNS curves, due to the normalization with respect to $\bar{\Omega}_{gw}^2$, which is different in all four cases). The factor diverges at large frequencies because it is normalized with respect to the monopole $\bar{\Omega}_{gw}^2$, which vanishes at these frequencies. The enclosed plot shows $\bar{\Omega}_{gw}^2 A_{gw}$, which encapsulates the absolute size of the anisotropies rather than their relative size compared to the monopole, and therefore does not diverge.</figcaption></figure>

the reference frequency in previous directional searches for an inspiral-dominated astrophysical background conducted by LIGO/Virgo [17]. We find

$$ A_{gw}(\nu_o = 65.75 \text{ Hz}) \approx 1.284 \times 10^{-4}. \quad (77) $$

This includes BHNS mergers at the maximal rate given in Eq. (29); setting the BHNS rate to 0 gives a 6% decrease in $A_{gw}$.

We have also computed the $C_\ell$ spectrum of the AGWB from an all-sky mock light cone catalogue based on the Millennium simulation, giving a much more accurate description of the anisotropies by relaxing many of the assumptions about the galaxy-galaxy 2PCF that were used to derive Eq. (61). This was done by constructing a GW overdensity map as described in Sec. IV, which is shown in Fig. 3. The $C_\ell$ spectrum computed from this map using the HEALPix software package is shown in Fig. 4, along with the analytical prediction. Despite the simplicity of the analytical model, and the fact that it contains no free parameters, it is in very strong agreement with the corresponding result from the catalogue at the largest angular scales (i.e., the lowest $\ell$-modes). It also captures the scaling of $C_\ell$ with $\ell$ at small scales (i.e., large $\ell$), although it overestimates the anisotropic variance at these scales by a

<footnote n="2">²http://healpix.sourceforge.net.</footnote>
<figure><img src="image_5.png" /><figcaption>FIG. 3. A HEALPix [47] map of the GW energy overdensity $\delta_{gw}^{(s)}$ constructed from the all-sky mock light cone catalogue [40,43–45], as described in Sec. IV. This was generated with the HEALPix $N_{side}$ parameter set to 256, corresponding to an angular resolution of 13.7 arcminutes, and an average of 7.3 galaxies per pixel.</figcaption></figure>

constant factor of ≈4. It is worth mentioning that such high values of $\ell$ are inaccessible to the current detector network due to their poor angular resolution—current directional searches only probe the first few $\ell$ [17], although this is expected to improve as more detectors are added to the network. Note that in order to ensure a fair comparison in Fig. 4, the analytical curve shown is for the cosmological parameters $H_0$, $\Omega_m$, $\Omega_\Lambda$ and galaxy 2PCF parameters $\gamma$, $d_1$ corresponding to the simulation. Replacing these with the

<figure><img src="image_6.png" /><figcaption>FIG. 4. The quantity $\ell(\ell + 1)C_\ell/2\pi$, which is the approximate contribution to the anisotropic variance of $\delta_{gw}^{(s)}$ as a function of $\ln \ell$. The red curve shows the simple analytical prediction given in Eq. (69), while the blue curve shows the spectrum computed from the map in Fig. 3 using HEALPix. Both curves include error regions from cosmic variance, while the blue curve includes Poisson errors associated with the finite number of galaxies per pixel in the HEALPix map.</figcaption></figure>

more up-to-date values from Planck [48] and VIPERS [39] has negligible effect.

The results in Fig. 4 can be compared directly with Fig. 2 of Ref. [21] [with the caveat that Fig. 4 includes BNS and BHNS mergers as well as BBH, whereas Ref. [21] includes only BBH—however, Fig. 2 indicates that this should only introduce a factor of $\mathcal{O}(1)$ between the two]. We note that the $C_\ell$'s in Fig. 4 are larger in amplitude than in Ref. [21], and do not fall off as quickly with $\ell$. This perhaps indicates that the linear perturbation-theory approach adopted in Ref. [21] to describe the galaxy distribution underestimates the clustering at small scales, where the perturbations to the density field become nonlinear. A more detailed comparison of the two sets of results is needed to fully assess the accuracy and applicability of each approach.

Figure 4 can also be compared with Fig. 5 of Ref. [20], which shows the corresponding $C_\ell$ spectra for cosmic string networks with different string tensions $G\mu$. The anisotropies in the astrophysical background are clearly much larger in amplitude, and become nonlinear for moderate values of $\ell$. This is unsurprising, given that the bulk of the GW flux from astrophysical sources is emitted at much lower redshifts, so that angular separations on the sky correspond to much smaller physical separations than in the cosmic string case, and therefore to scales where the departure from homogeneity is more evident.

We have also computed the kinematic dipole factor, which in the case shown in Figs. 3 and 4 has a value of

$$ \mathcal{D}(\nu_o = 65.75 \text{ Hz}) \approx 7.099 \times 10^{-4}. \qquad (78) $$

This is an order of magnitude smaller than in the cosmic string case, which is intuitively sensible given that the galaxies we consider are generally at much lower redshifts and therefore their velocities due to the Hubble flow are much smaller, giving a weaker Doppler effect. The fact that $\mathcal{D}$ is smaller in this case, while $\delta_{gw}^{(s)}$ is simultaneously larger, means that the effects of the kinematic dipole are unimportant for the astrophysical background.

## VI. CONCLUSION

We have developed a detailed anisotropic model for the AGWB, including the most important sources for the LIGO-Virgo frequency band (BBH, BNS, and BHNS). The angular spectrum of the anisotropies, quantified by the $C_\ell$ components, has been calculated through two complementary approaches: a simple, closed-form analytical expression Eq. (69), and a detailed numerical study using an all-sky mock light cone galaxy catalogue from the Millennium simulation [40,43–45]. The two approaches are in excellent agreement at large angular scales ($\ell \lesssim 10$), and differ only by a factor of order unity at smaller scales, following the variation of $\ell(\ell + 1)C_\ell/2\pi$ over many orders of magnitude despite the simplicity of the analytical model and the lack of free parameters. These anisotropies
are considerably larger in amplitude than those in the temperature of the CMB, or those in the SGWB due to cosmic strings [20], and become nonlinear at higher multipoles $ℓ$. This shows that modeling the AGWB in a purely isotropic manner neglects a great deal of astrophysical and cosmological information, thereby motivating future theoretical and observational work. We expect that in the near future, anisotropies in the AGWB (and the SGWB more generally) will become an important probe of the large-scale structure of the Universe, and of the astrophysical processes that occur within it.

We have highlighted several key avenues for future work to explore. One of the key differences between the SGWB and the CMB is the very broad frequency spectrum over which the SGWB can be investigated. This motivates the development of models that are valid at frequencies outside of the LIGO-Virgo band; the study of AGWB anisotropies in the LISA band will be of particular interest in the coming decades. As shown in Sec. II C, the use of narrow frequency bins makes it impossible to guarantee that the AGWB is Gaussian, and a full characterization of the anisotropies in this bin requires higher-order correlators of the field (such as the bispectrum and trispectrum). This is a daunting task, both theoretically and observationally, but it holds the promise of incredibly rich new astrophysical and cosmological information.

There will also be a need to develop increasingly accurate models of the AGWB anisotropies. This will be addressed in future work by improving upon the fiducial astrophysical model of Ref. [22], using past and future LIGO-Virgo observing runs to enhance our knowledge of the relevant sources. The catalogue approach described in Sec. IV will also be improved upon by using more accurate and complete galaxy catalogues. It will also be interesting to explore different ways in which the statistical properties of the AGWB might differ from those of the model described here, e.g., the extension to models which are statistically nonstationary, or which are no longer statistically isotropic due to the inclusion of sources in the galactic plane of the Milky Way.

Finally, we plan to produce realistic mock data from this model and from future models, and use them to test and optimize data analysis and parameter estimation methods, leading the way to future detections and measurement of AGWB anisotropies.

## ACKNOWLEDGMENTS

We thank Joe Romano and Andrew Matas for reading the manuscript carefully and providing us with valuable comments. The Millennium Simulation databases used in this paper and the web application providing online access to them were constructed as part of the activities of the German Astrophysical Virtual Observatory. Some of the results in this paper have been derived using the HEALPix package [47]. A. C. J. is supported by King's College

London through a Graduate Teaching Scholarship. M. S. is supported in part by the Science and Technology Facility Council (STFC), United Kingdom, under Grant No. ST/P000258/1.

## APPENDIX A: EXPLICITLY EVALUATING THE ANALYTICAL $C_ℓ$ SPECTRUM

The $C_ℓ$'s are given by

$$C_{\ell}(\nu_{o}) = 2\pi A_{\text{gw}} \int_{-1}^{+1} d(\cos \theta_{o}) P_{\ell}(\cos \theta_{o}) \left(\tan \frac{\theta_{o}}{2}\right)^{-\gamma} \quad (\text{A1})$$

Applying simple trigonometry, we have

$$\tan \frac{\theta_o}{2} = \sqrt{\frac{1 - \cos \theta_o}{1 + \cos \theta_o}}$$

and therefore

$$C_{\ell} = 2\pi A_{\text{gw}} \int_{-1}^{+1} dx P_{\ell}(x) \left(\frac{1+x}{1-x}\right)^{\gamma/2}.$$

This can be evaluated by writing the Legendre polynomial as a sum,

$$P_{\ell}(x) = \sum_{k=0}^{\ell} \frac{(\ell + k)!}{(k!)^2 (\ell - k)!} \left(\frac{x-1}{2}\right)^k,$$

so that we have

$$C_{\ell} = 2\pi A_{\text{gw}} \sum_{k=0}^{\ell} \frac{(\ell + k)!}{(k!)^2 (\ell - k)!} \int_{-1}^{+1} dx \left(\frac{1+x}{1-x}\right)^{\gamma/2} \left(\frac{x-1}{2}\right)^k.$$

With a change of variables to $y \equiv \frac{1-x}{2}$, the integral becomes

$$\int_{-1}^{+1} dx \left(\frac{1+x}{1-x}\right)^{\gamma/2} \left(\frac{x-1}{2}\right)^k = 2(-1)^k \int_{0}^{+1} dy (1-y)^{\gamma/2} y^{k-\frac{\gamma}{2}}.$$

This is a Euler integral of the first kind, and evaluates to

$$\int_{0}^{+1} dy (1-y)^{\gamma/2} y^{k-\frac{\gamma}{2}} = \frac{\Gamma(1+k-\frac{\gamma}{2})\Gamma(1+\frac{\gamma}{2})}{\Gamma(k+2)}$$

where we assume $0 < \gamma < 2$. We therefore have
$$ 
\begin{aligned}
C_{\ell} &= 4\pi A_{\text{gw}} \sum_{k=0}^{\ell} \frac{(-1)^k (\ell + k)! \Gamma(1 + k - \frac{\ell}{2}) \Gamma(1 + \frac{\ell}{2})}{(k!)^2 (\ell - k)! \Gamma(k + 2)} \\
&= 4\pi A_{\text{gw}} \sum_{k=0}^{\ell} (-1)^k \frac{\Gamma(\ell + k + 1)}{\Gamma(k + 1)^2 \Gamma(\ell - k + 1)} \frac{\Gamma(1 + k - \frac{\ell}{2}) \Gamma(1 + \frac{\ell}{2})}{\Gamma(k + 2)} \\
&= 4\pi A_{\text{gw}} \sum_{k=0}^{\infty} (-1)^k \frac{\Gamma(\ell + k + 1)}{\Gamma(k + 1)^2 \Gamma(\ell - k + 1)} \frac{\Gamma(1 + k - \frac{\ell}{2}) \Gamma(1 + \frac{\ell}{2})}{\Gamma(k + 2)}
\end{aligned}
$$

where the final line follows because $\Gamma(\ell + k + 1)/\Gamma(\ell - k + 1)$ vanishes for $k > \ell$. We can simplify further by defining the rising Pochammer symbol,

$$ (x)_k \equiv x(x+1)\cdots(x+k-1), $$

and by using Euler's reflection formula,

$$ \mathrm{sinc}(\pi x) \equiv \frac{\sin \pi x}{\pi x} = \frac{1}{\Gamma(1+x)\Gamma(1-x)}. $$

This gives

$$ C_{\ell} = \frac{4\pi A_{\text{gw}}}{\mathrm{sinc}(\pi\gamma/2)} \sum_{k=0}^{\infty} \frac{(-\ell)_k (\ell + 1)_k (1 - \frac{\ell}{2})_k}{(1)_k (2)_k k!}. $$

The series above defines a generalized hypergeometric function ${}_3F_2$, so that we get our final expression,

$$ C_{\ell} = 4\pi A_{\text{gw}} \frac{{}_3F_2(-\ell, \ell+1, 1-\frac{\ell}{2}; 1, 2; 1)}{\mathrm{sinc}(\pi\gamma/2)}. \quad (A2) $$

Note that when $\ell = 0$, the hypergeometric function evaluates to unity, and we have

$$ C_0 = \frac{4\pi A_{\text{gw}}}{\mathrm{sinc}(\pi\gamma/2)} > 0. \quad (A3) $$

However, the $C_0$ component is just an average of $C_{gw}$ over the sphere,

$$ 
\begin{aligned}
C_0 &= 2\pi \int_{-1}^{+1} d(\cos \theta_o) C_{gw} = \int_0^{2\pi} d\phi_o \int_0^{\pi} d\theta_o \sin \theta_o C_{gw} \\
&= \int_{S^2} d^2 \sigma_o C_{gw},
\end{aligned}
$$

and therefore must be equal to 0, since the spherical integral and the ensemble averaging process commute, and since the average of $\delta_{gw}^{(s)}$ over the sphere is 0 by definition,

$$ C_0 = \left\langle \delta_{gw}^{(s)} \int_{S^2} d^2 \sigma_o \delta_{gw}^{(s)} \right\rangle = 0. \quad (A4) $$

This is indicative of the inaccuracies inherent to Eq. (69).

## APPENDIX B: WAVEFORM NUMERICAL CONSTANTS

The following expressions and numerical values are needed to fully specify the hybrid waveforms given in Eq. (17). They exactly match those given in Ref. [32].

<table><thead><tr><th>i</th><th>ν<sub>i</sub><sup>0</sup></th><th>y<sub>i</sub><sup>(10)</sup></th><th>y<sub>i</sub><sup>(11)</sup></th><th>y<sub>i</sub><sup>(12)</sup></th><th>y<sub>i</sub><sup>(20)</sup></th><th>y<sub>i</sub><sup>(21)</sup></th><th>y<sub>i</sub><sup>(30)</sup></th></tr></thead><tbody><tr><td>1</td><td>1 - 4.455(1 - χ)<sup>0.217</sup> + 3.521(1 - χ)<sup>0.26</sup></td><td>+0.6437</td><td>+0.827</td><td>-0.2706</td><td>-0.05822</td><td>-3.935</td><td>-7.092</td></tr><tr><td>2</td><td>[1 - 0.63(1 - χ)<sup>0.3</sup>]/2</td><td>+0.1469</td><td>-0.1228</td><td>-0.02609</td><td>-0.0249</td><td>+0.1701</td><td>+2.325</td></tr><tr><td>3</td><td>[1 - 0.63(1 - χ)<sup>0.3</sup>](1 - χ)<sup>0.45</sup>/4</td><td>-4.098</td><td>-0.03523</td><td>+0.1008</td><td>+1.829</td><td>-0.02017</td><td>-2.87</td></tr><tr><td>4</td><td>0.3236 + 0.04894χ + 0.01346χ<sup>2</sup></td><td>-0.1331</td><td>-0.08172</td><td>+0.1451</td><td>-0.2714</td><td>+0.1279</td><td>+4.922</td></tr></tbody></table>

$$ \alpha_2 = -\frac{323}{224} + \frac{451}{168} \left(\frac{M}{M}\right)^{5/3}, $$

$$ \epsilon_1 = 1.4547\chi - 1.8897, $$

$$ c_1 = \nu_1^{-1} \left[ \frac{1 + \sum_{i=2}^{3} \alpha_i (\pi GM\nu_1)^{i/3}}{1 + \sum_{i=1}^{2} \epsilon_i (\pi GM\nu_1)^{i/3}} \right]^2, $$

$$ \alpha_3 = \left[ \frac{27}{8} - \frac{11}{6} \left(\frac{M}{M}\right)^{5/3} \right] \chi, $$

$$ \epsilon_2 = -1.8153\chi + 1.6557 $$

$$ c_2 = c_1 \nu_2^{-4/3} \left[ 1 + \sum_{i=1}^{2} \epsilon_i (\pi GM\nu_1)^{i/3} \right]^2. $$
[1] G. M. Harry (LIGO Scientific Collaboration), *Classical Quantum Gravity* **27**, 084006 (2010).

[2] J. Aasi *et al.* (LIGO Scientific Collaboration), *Classical Quantum Gravity* **32**, 115012 (2015).

[3] F. Acernese *et al.* (VIRGO Collaboration), *Classical Quantum Gravity* **32**, 024001 (2015).

[4] B. P. Abbott *et al.* (Virgo, LIGO Scientific Collaboration), *Phys. Rev. Lett.* **116**, 241103 (2016).

[5] B. P. Abbott *et al.* (VIRGO, LIGO Scientific Collaboration), *Phys. Rev. Lett.* **118**, 221101 (2017).

[6] B. P. Abbott *et al.* (Virgo, LIGO Scientific Collaboration), *Astrophys. J.* **851**, L35 (2017).

[7] B. P. Abbott *et al.* (Virgo, LIGO Scientific Collaboration), *Phys. Rev. Lett.* **119**, 141101 (2017).

[8] B. Abbott *et al.* (Virgo, LIGO Scientific Collaboration), *Phys. Rev. Lett.* **119**, 161101 (2017).

[9] S. W. Ballmer, *Classical Quantum Gravity* **23**, S179 (2006).

[10] B. Abbott *et al.* (LIGO Scientific Collaboration), *Phys. Rev. D* **76**, 082003 (2007).

[11] E. Thrane, S. Ballmer, J. D. Romano, S. Mitra, D. Talukder, S. Bose, and V. Mandic, *Phys. Rev. D* **80**, 122002 (2009).

[12] J. Abadie *et al.* (LIGO Scientific Collaboration), *Phys. Rev. Lett.* **107**, 271102 (2011).

[13] C. Messenger, H. J. Bulten, S. G. Crowder, V. Dergachev, D. K. Galloway, E. Goetz, R. J. G. Jonker, P. D. Lasky, G. D. Meadors, A. Melatos, S. Premachandra, K. Riles, L. Sammut, E. H. Thrane, J. T. Whelan, and Y. Zhang, *Phys. Rev. D* **92**, 023006 (2015).

[14] C. Chung, A. Melatos, B. Krishnan, and J. T. Whelan, *Mon. Not. R. Astron. Soc.* **414**, 2650 (2011).

[15] L. Sun, A. Melatos, P. D. Lasky, C. T. Y. Chung, and N. S. Darman, *Phys. Rev. D* **94**, 082004 (2016).

[16] J. Aasi *et al.* (VIRGO, LIGO Scientific Collaboration), *Phys. Rev. D* **88**, 122004 (2013).

[17] B. P. Abbott *et al.* (Virgo, LIGO Scientific Collaboration), *Phys. Rev. Lett.* **118**, 121102 (2017).

[18] G. Cusin, C. Pitrou, and J.-P. Uzan, *Phys. Rev. D* **96**, 103019 (2017).

[19] G. Cusin, C. Pitrou, and J.-P. Uzan, *Phys. Rev. D* **97**, 123527 (2018).

[20] A. C. Jenkins and M. Sakellariadou, arXiv:1802.06046.

[21] G. Cusin, I. Dvorkin, C. Pitrou, and J.-P. Uzan, *Phys. Rev. Lett.* **120**, 231101 (2018).

[22] B. P. Abbott *et al.* (Virgo, LIGO Scientific Collaboration), *Phys. Rev. Lett.* **120**, 091101 (2018).

[23] I. Dvorkin, J.-P. Uzan, E. Vangioni, and J. Silk, *Mon. Not. R. Astron. Soc.* **479**, 121 (2018).

[24] I. Dvorkin, J.-P. Uzan, E. Vangioni, and J. Silk, *Phys. Rev. D* **94**, 103011 (2016).

[25] I. Dvorkin, E. Vangioni, J. Silk, J.-P. Uzan, and K. A. Olive, *Mon. Not. R. Astron. Soc.* **461**, 3877 (2016).

[26] P. Amaro-Seoane *et al.* (LISA Collaboration), *arXiv: 1702.00786*.

[27] T. Regimbau, M. Evans, N. Christensen, E. Katsavounidis, B. Sathyaprakash, and S. Vitale, *Phys. Rev. Lett.* **118**, 151105 (2017).

[28] R. Smith and E. Thrane, *Phys. Rev. X* **8**, 021019 (2018).

[29] T. Namikawa, A. Nishizawa, and A. Taruya, *Phys. Rev. Lett.* **116**, 121302 (2016).

[30] T. Namikawa, A. Nishizawa, and A. Taruya, *Phys. Rev. D* **94**, 024013 (2016).

[31] P. Ajith, S. Babak, Y. Chen, M. Hewitson, B. Krishnan, A. M. Sintes, J. T. Whelan, B. Brugmann, P. Diener, N. Dorband, J. Gonzalez, M. Hannam, S. Husa, D. Pollney, L. Rezzolla, L. Santamaria, U. Sperhake, and J. Thornburg, *Phys. Rev. D* **77**, 104017 (2008); **79**, 129901(E) (2009).

[32] P. Ajith, M. Hannam, S. Husa, Y. Chen, B. Brugmann, N. Dorband, D. Muller, F. Ohme, D. Pollney, C. Reisswig, L. Santamaria, and J. Seiler, *Phys. Rev. Lett.* **106**, 241101 (2011).

[33] B. P. Abbott *et al.* (Virgo, LIGO Scientific Collaboration), *arXiv:1805.11579*.

[34] B. P. Abbott *et al.* (Virgo, LIGO Scientific Collaboration), *Astrophys. J.* **832**, L21 (2016).

[35] J. Abadie *et al.* (VIRGO, LIGO Scientific Collaboration), *Classical Quantum Gravity* **27**, 173001 (2010).

[36] E. Vangioni, K. A. Olive, T. Prestegard, J. Silk, P. Petitjean, and V. Mandic, *Mon. Not. R. Astron. Soc.* **447**, 2575 (2015).

[37] C. J. Conselice, A. Wilkinson, K. Duncan, and A. Mortlock, *Astrophys. J.* **830**, 83 (2016).

[38] K. Belczynski, D. E. Holz, T. Bulik, and R. O'Shaughnessy, *Nature (London)* **534**, 512 (2016).

[39] F. Marulli *et al.*, *Astron. Astrophys.* **557**, A17 (2013).

[40] G. De Lucia and J. Blaizot, *Mon. Not. R. Astron. Soc.* **375**, 2 (2007).

[41] S. Hatton, J. E. G. Devriendt, S. Ninin, F. R. Bouchet, B. Guiderdoni, and D. Vibert, *Mon. Not. R. Astron. Soc.* **343**, 75 (2003).

[42] B. M. B. Henriques, S. White, P. Thomas, R. Angulo, Q. Guo, G. Lemson, V. Springel, and R. Overzier, *Mon. Not. R. Astron. Soc.* **451**, 2663 (2015).

[43] J. Blaizot, Y. Wadadekar, B. Guiderdoni, S. T. Colombi, E. Bertin, F. R. Bouchet, J. E. G. Devriendt, and S. Hatton, *Mon. Not. R. Astron. Soc.* **360**, 159 (2005).

[44] V. Springel *et al.*, *Nature (London)* **435**, 629 (2005).

[45] G. Lemson (VIRGO Collaboration), *arXiv:astro-ph/ 0608019*.

[46] E. Hawkins *et al.*, *Mon. Not. R. Astron. Soc.* **346**, 78 (2003).

[47] K. M. Gorski, E. Hivon, A. J. Banday, B. D. Wandelt, F. K. Hansen, M. Reinecke, and M. Bartelmann, *Astrophys. J.* **622**, 759 (2005).

[48] P. A. R. Ade *et al.* (Planck Collaboration), *Astron. Astrophys.* **594**, A13 (2016).