Properties of the stochastic astrophysical gravitational wave background: astrophysical sources dependencies

TL;DR

This paper investigates how different astrophysical models influence the anisotropies of the gravitational wave background, revealing that certain parameters significantly affect the angular power spectrum and that anisotropies offer unique astrophysical insights.

Contribution

It demonstrates the sensitivity of AGWB anisotropies to specific astrophysical parameters, highlighting their potential as a new observable for astrophysical research.

Findings

Metallicity and merger time delay significantly affect anisotropies.

Some parameters have negligible impact on the anisotropy spectrum.

AGWB anisotropies can provide unique astrophysical information.

Abstract

This article explores the properties (amplitude and shape) of the angular power spectrum of the anisotropies of the astrophysical gravitational wave background (AGWB) focusing on the signatures of the astrophysical models describing sub-galactic physics. It demonstrates that while some parameters have negligible impact others, and in particular the stellar evolution models, the metallicity and the merger time delay distribution can result in relative differences of order 40% in the angular power spectrum of anisotropies in both the LIGO/Virgo and LISA frequency bands. It is also shown that the monopole and the anisotropic components of the AGWB are complementary and sensitive to different astrophysical parameters. It follows that AGWB anisotropies are a new observable with the potential to provide new astrophysical information that can not be accessed otherwise.

Tables

Table 1. Table 1: Values of the normalization parameter β 𝛽 \beta for the different models described in §  IV.3 and  IV.4 . The normalization was determined in order to obtain 10 10 10 detected events during the aLIGO O1+O2 timespan for all of the models.

Model	Ref	Limongi	imf-low	imf-hi	dMco	uMco	aconst
$\beta$	$0.01$	$0.013$	$0.004$	$0.032$	$0.006$	$0.023$	$0.17$

Equations

$${\Omega }_{\mathrm{GW}}(\mathbf{e},f)$$

$$=\frac{{\bar{\Omega }}_{\mathrm{GW}}(f)}{4\pi }+\delta {\Omega }_{\mathrm{GW}}(\mathbf{e},f),$$

$${\Omega }_{\mathrm{GW}}(\mathbf{e},f)=\frac{c^2}{4G{\rho }_c}\frac{1}{T_{{}_{\mathrm{O}}}}{f}^3\sum _{A=+,\times }|{\tilde{h}}_A(f,\mathbf{e}){|}^2,$$

$${\bar{\Omega }}_{\mathrm{GW}}(f)=$$

$${\partial }_{\eta }{\bar{\Omega }}_{\mathrm{GW}}=\frac{f}{{\rho }_c}\mathcal{A}(f,\eta ),$$

$$\mathcal{A}(\eta ,f)\equiv a^4\int \text{d}{\theta }_{{}_{\mathrm{G}}}{\bar{n}}_{{}_{\mathrm{G}}}(\eta ,{\theta }_{{}_{\mathrm{G}}}){\mathcal{L}}_{{}_{\mathrm{G}}}(\eta ,f_{{}_{\mathrm{G}}},{\theta }_{{}_{\mathrm{G}}}),$$

$$f_{{}_{\mathrm{G}}}=(1+z_{{}_{\mathrm{G}}})f$$

$$\delta {\Omega }_{\mathrm{GW}}(\mathbf{e},f)=$$

$$=\frac{f}{4\pi {\rho }_c}{\int }_{{\eta }_{\ast }}^{{\eta }_{{}_{\mathrm{O}}}}\text{d}\eta \mathcal{A}\left(\eta ,f\right)\left[{\delta }_G+4\Psi -2\mathbf{e}\cdot \nabla v+6{\int }_{\eta }^{{\eta }_{{}_{\mathrm{O}}}}\text{d}{\eta }^{\prime }\dot{\Psi }\right]$$

$$+\frac{f}{4\pi {\rho }_c}{\int }_{{\eta }_{\ast }}^{{\eta }_{{}_{\mathrm{O}}}}\text{d}\eta \mathcal{B}(\eta ,f)\left[\mathbf{e}\cdot \nabla v-\Psi -2{\int }_{\eta }^{{\eta }_{{}_{\mathrm{O}}}}d{\eta }^{\prime }\dot{\Psi }\right],$$

$$\displaystyle\mathcal{B}(\eta,f)\equiv f\,a^{3}\bar{n}_{{}_{\rm G}}(\eta)\int d\theta_{{}_{\rm G}}\frac{\partial\mathcal{L}_{{}_{\rm G}}}{\partial f_{{}_{\rm G}}}\Big{|}_{{\bar{f}}_{{}_{\rm G}}}(\eta,f_{{}_{\rm G}},\theta_{{}_{\rm G}})\,,$$

$${\delta }_{\mathrm{G}}^{\mathrm{c}}=b{\delta }_{\mathrm{m}}^{\mathrm{c}}\Rightarrow {\delta }_{\mathrm{G}}+3\mathcal{H}v=b\left({\delta }_{\mathrm{m}}+3\mathcal{H}v\right),$$

$$C_{\ell }(f)=\frac{2}{\pi }\int \text{d}k{k}^2|\delta {\Omega }_{\ell }(k,f){|}^2,$$

$$\displaystyle{\delta\Omega}_{\ell}(k,f)=\frac{f}{4\pi\rho_{c}}\Bigg{\{}\int_{\eta_{*}}^{\eta_{{}_{\rm O}}}\text{d}\eta\,\mathcal{A}(\eta,f)\times$$

$$[\left(4{\Phi }_k(\eta )+b{\delta }_{\mathrm{m},\mathrm{k}}(\eta )+(b-1)3\mathcal{H}{v}_k(\eta )\right)j_{\ell }(k\Delta \eta )$$

$$-2k{v}_k(\eta )j_{\ell }^{\prime }(k\Delta \eta )]$$

$$+{\int }_{{\eta }_{\ast }}^{{\eta }_{{}_{\mathrm{O}}}}\text{d}\eta \mathcal{B}(\eta ,f)\left[-{\Phi }_k(\eta )j_{\ell }(k\Delta \eta )+k{v}_k(\eta )j_{\ell }^{\prime }(k\Delta \eta )\right]$$

$$\displaystyle+\int_{\eta_{*}}^{\eta_{{}_{\rm O}}}\text{d}\eta\left[6\mathcal{A}(\eta,f)-2\mathcal{B}(\eta,f)\right]\int_{\eta}^{\eta_{{}_{\rm O}}}\text{d}\tilde{\eta}\,{\Phi}^{\prime}_{k}(\tilde{\eta})j_{\ell}(k\Delta\tilde{\eta})\Bigg{\}}\,,{}$$

$$C_{\ell }^{\text{rel}}\equiv C_{\ell }\frac{(4\pi {)}^2}{{\bar{\Omega }}_{\mathrm{GW}}^2}.$$

$$B_{\ell }^X(f)\equiv \frac{2}{\pi }\int \text{d}k{k}^2\frac{4\pi }{{\bar{\Omega }}_{\mathrm{GW}}(f)}\delta {\Omega }_{\ell }^{\ast }(k,f)X_{\ell }(k).$$

$${\kappa }_{\ell }=-\frac{\ell (\ell +1)}{2}{\int }_0^{{\chi }_H}d\chi g(\chi ){\hat{\Psi }}_k(\chi )j_{\ell }(k\chi ),$$

$$B_{\ell }(f)$$

$$g(\chi )\equiv \frac{1}{\chi }{\int }_{\chi }^{{\chi }_H}d{\chi }^{\prime }{p}_{\chi }({\chi }^{\prime })\frac{({\chi }^{\prime }-\chi )}{{\chi }^{\prime }},$$

$$D_{\ell }(f,z)$$

$$b(z)=b_0\sqrt{1+z}{b}_0=1.5.$$

$$\delta {\Omega }_{\mathrm{GW}}(\mathbf{e},f)={\int }_{{\eta }_{\ast }}^{{\eta }_{{}_{\mathrm{O}}}}\text{d}\eta {\partial }_{\eta }\left(\frac{{\bar{\Omega }}_{\mathrm{GW}}}{4\pi }\right){\delta }_G(\mathbf{e},\eta ).$$

$$C_{\ell }\simeq \frac{2}{\pi }\int k^2\text{d}k{P}_{\mathrm{Gal}}(k){\left|\int \text{d}\eta {\partial }_{\eta }\left(\frac{{\bar{\Omega }}_{\mathrm{GW}}}{4\pi }\right)j_{\ell }(k\Delta \eta )\right|}^2,$$

$$\int \text{d}x{j}_{\ell }(x)f(x)\simeq \sqrt{\frac{\pi }{2\ell +1}}f(\ell +1/2).$$

$$C_{\ell }^{\mathrm{Limber}}\simeq$$

$$k\Delta \eta =\ell +\frac{1}{2}.$$

$${\ell }_{\max }\equiv k[{\eta }_0-\eta (z_{\max })].$$

$$C_{\ell }^{\text{Limber+static}}\propto \frac{1}{\ell +\frac{1}{2}}{\left|{\partial }_{\eta }\left(\frac{{\bar{\Omega }}_{\mathrm{GW}}}{4\pi }\right)\right|}_{\eta ={\eta }_0}^2\times$$

$${\int }_{k_{\min }(\ell )}{P}_{\mathrm{Gal}}(k)\text{d}k,$$

$$C_{\ell }^{\mathrm{Limber}}\simeq {\left(\ell +\frac{1}{2}\right)}^{-1}\int \text{d}\log rk{P}_{\mathrm{Gal}}(k){\left|{\partial }_r\left(\frac{{\bar{\Omega }}_{\mathrm{GW}}(r)}{4\pi }\right)\right|}^2,$$

Full text

1. I Introduction

2. II Anisotropies of AGWB

3. II.1 Isotropic contribution

4. II.2 Anisotropies

5. II.3 Cross-correlation

6. II.4 Coarse graining approach

7. III General properties of the angular power spectrum

8. III.1 Analytic approximation of the angular power spectrum

9. III.2 General treatment of shot noise

10. III.3 Shot-noise and cross-correlation

11. III.4 Frequency-direction factorization

12. IV Astrophysical models

13. IV.1 Model framework

14. IV.2 Normalization to the number of detected events

15. IV.3 Reference model

16. IV.4 Modified models

17. IV.4.1 BH formation model [limongi model]

18. IV.4.2 BH mass cutoff due to PISN [dMco and uMco models]

19. IV.4.3 Stellar initial mass function [imf-high and imf-low models]

20. IV.4.4 Distribution of initial separations [aconst model]

21. IV.4.5 Metallicity of progenitor stars

22. V Results

23. V.1 Monopole

24. V.2 Angular power spectrum

25. V.3 Cross-correlations

26. VI Contribution of binary neutron stars

27. VII Lower frequencies and the LISA band

28. VIII Conclusion

29. A Computation of shot noise

I Introduction

Diffuse stochastic backgrounds arise from the incoherent superposition of signals from resolved and unresolved sources. Many such backgrounds for different kinds of radiation have been observed in astronomy. Electromagnetic backgrounds of radiation include the cosmic microwave background (CMB) with its black body spectrum Penzias and Wilson (1965), the cosmic infrared background (CIB) from stellar dust Hauser and Dwek (2001) and the extragalactic background made up of all the electromagnetic radiation emitted by stars, galaxies, galaxy clusters etc. since their formation Partridge and Peebles (1967); Shanks et al. (1991). Similarly, there should exist a neutrino background Hannestad (2006) and a background of gravitational waves (GW).

The GW background can be split into a stochastic background of gravitational radiation of cosmological origin, e.g. produced during inflation, and one of astrophysical origin (AGWB). The latter results from the superposition of a large number of resolved and unresolved sources from the onset of stellar activity until today. The nature of the AGWB is expected to be significantly different from its cosmological counterpart, which is expected to be, at least for inflation, stationary, unpolarized, almost statistically Gaussian and isotropic, by analogy with the cosmic microwave background. Many different astrophysical sources contribute to the AGWB, including merging stellar-mass black hole (BH) and neutron star (NS) binaries Abbott et al. (2016); Regimbau et al. (2016); Mandic et al. (2016); Dvorkin et al. (2016a); Nakazato et al. (2016); Dvorkin et al. (2016b); Evangelista and Araujo (2014), merging supermassive black hole binaries Kelley et al. (2017), rotating neutron stars Surace et al. (2016); Talukder et al. (2014); Lasky et al. (2013), stellar core collapse Crocker et al. (2017, 2015) and population III binaries Kowalska et al. (2012).

The observational landscape is growing and covers a large range of frequencies; see e.g. Ref. Moore et al. (2015) for a review. At extremely low frequencies $\sim 10^{-16}$ Hz observational bounds come mainly from the analysis of CMB B-modes while frequencies in the range $10^{-10}-10^{-6}$ Hz are covered by pulsar timing arrays: the Parkes Pulsar Timing Array111http://www.atnf.csiro.au/research/pulsar/ppta/ (PPTA), the Large European Array for Pulsar Timing222http://www.leap.eu.org (LEPTA), and the North American Nanohertz Observatory for Gravitational Waves (NANOGrav), all of which form the International Pulsar Timing Array Consortium333http://www.ipta4gw.org (IPTA). Frequencies in the range $10^{-4}-10^{-1}$ Hz will be probed with the space-based Laser Interferometer Space Antenna444www.lisamission.org (LISA) scheduled to be launched in 2034. Higher frequencies ($1-10^{3}$ Hz) are accessible with ground-based interferometers, including Advanced LIGO (aLIGO) LIGO Scientific Collaboration (2015) and Advanced Virgo (aVirgo) Acernese et al. (2015) which already conducted two observational runs (O1 and O2) during 2014-2017 and are currently entering into the third observational run, KAGRA interferometer which is expected to become operational in 2018-2019 and LIGO India which is currently under construction. A third generation of ground-based interferometers, the Einstein Telescope555http://www.et-gw.eu (ET) and the Cosmic Explorer (CE) Abbott et al. (2017a) are in their design stages.

The latest upper bounds obtained in Ref. The LIGO Scientific Collaboration and the Virgo Collaboration (2019a) using the first and second aLIGO observing runs are $\Omega_{\rm GW}(f=25\text{Hz})<4.8\times 10^{-8}$, assuming a population of compact binary sources, and $\Omega_{\rm GW}(f=25\text{Hz})<6\times 10^{-8}$ for a frequency-independent background for the frequency ranges $20-92$ Hz and $20-80$ Hz, respectively, where $\Omega_{\rm GW}=\text{d}\rho_{\rm GW}/\text{d}\ln f/\rho_{c}$ is the energy density in GW per logarithmic frequency interval in units of the critical density of the Universe. This improves bounds on the stochastic background obtained from the analysis of big-bang nucleosynthesis Maggiore (2000); Allen (1996), and of the cosmic microwave background Smith et al. (2006); Henrot-Versille et al. (2015) at 100 Hz. LIGO-Virgo upper bounds are well above current theoretical predictions, for example the population model derived from the O1+O2 source catalog The LIGO Scientific Collaboration and the Virgo Collaboration (2018a, 2019a) predicts an amplitude of the total background (binary black holes and binary neutron stars) of $\Omega_{\rm GW}(f=25\text{Hz})=8.9_{-5.6}^{+12.6}\times 10^{-10}$ and $\Omega_{\rm GW}(f=25\text{Hz})=5.3_{-2.5}^{+4.2}\times 10^{-10}$ from binary black holes alone, where the uncertainties are due to $90\%$ confidence limits on the merger rates. Assuming the most probable rate for compact binary mergers, Ref. Abbott et al. (2018) concludes that the total background may be detectable with a signal-to-noise-ratio of 3 after 40 months of total observation time. At low frequencies, Pulsar Timing Arrays give a bound $\Omega_{\rm GW}<1.3\times 10^{-9}$ for $f=2.8\times 10^{-9}$ Hz Shannon et al. (2013). The possibility of measuring and mapping the gravitational wave background is discussed in Refs. Allen and Ottewill (1997); Cornish (2001); Mitra et al. (2008); Thrane et al. (2009); Romano et al. (2015); Romano and Cornish (2016); Renzini and Contaldi (2018) while different methods employed by LIGO and LISA to reconstruct an angular resolved map of the sky are presented in Ref. Abbott et al. (2017b). An analogous discussion for Pulsar Timing Arrays can be found in Refs. Mingarelli et al. (2013); Taylor and Gair (2013); Gair et al. (2014).

The latest observational constraints from the first and second aLIGO runs The LIGO Scientific Collaboration and the Virgo Collaboration (2019b) are derived for multipoles up to $\ell=4$ with upper limits on the amplitude in the range $\Omega_{\rm GW}(f=25\text{Hz},\Theta)<0.64-2.47\times 10^{-8}$ sr*-1* for a population of binary compact objects and assuming that the angular and frequency dependencies factorize, an assumption we shall investigate below for realistic astrophysical models.

From a theoretical perspective, as any background of radiation, the AGWB is fully characterized in terms of Stokes parameters, intensity and polarization, as a function of direction of frequency; see Refs. Romano and Cornish (2016); Cusin et al. (2019a) for a definition of Stokes parameters for a background of spin-2 radiation. The first prediction of the AGWB angular power spectrum was presented in our analysis Cusin et al. (2018a) following our seminal formalism666Note that a first attempt to describe anisotropies of the AGWB with a Boltzmann approach was proposed by Ref. Contaldi (2017) while Ref. Cusin et al. (2019a) refines it by introducing an emissivity function that realistically describes GW emission at the galactic scale. introduced in Refs. Cusin et al. (2018b, 2017). This formalism is very flexible and splits the cosmological, large-scale structure and sub-galactic scales so that it can be applied to any source contributions and any frequency band. It significantly differs from the simpler computation of the monopole Regimbau (2011); Dvorkin et al. (2016a) in which the sources were assumed to be homogeneously and isotropically distributed. The anisotropies were further studied in Refs. Jenkins et al. (2018a, b) using a different set of astrophysical models (see Ref. Cusin et al. (2018c) for a critical analysis of these works). The first study of the generation of polarization induced by the diffusion by massive structures is presented in Ref. Cusin et al. (2019a).

The goal of this article is to extend our previous analysis Cusin et al. (2019a) for the contribution of BH mergers in the LIGO/Virgo frequency band, in the following directions:

describe in details the properties of the angular power spectrum of AGWB; 2. 2.

explore the astrophysical dependencies of the angular power spectrum; 3. 3.

study the contribution of binary NS mergers; 4. 4.

extend the analysis to lower frequencies (LISA band and lower).

In the LIGO band where the background is dominated by mergers of compact objects, we explore different stellar models for the evolution of BH stellar progenitors and we study the dependence of the result on the distribution function of orbital parameters of binary objects, mass distribution and initial mass function. We also analyze how different BH populations contribute to the background. We explicitly show that anisotropies are very sensitive to changes in the astrophysical model used to describe the sub galactic process of formation and evolution of GW sources and the sub-galactic process of GW emission. We present a detailed discussion of which are the astrophysical parameters and functions the angular power spectrum is most sensitive to and hence that we will be able to constrain the first, in both the LIGO and the LISA frequency bands.

The article is organized as follows. Section II summarizes the main results of our formalism Cusin et al. (2018b, 2017) describing the computation of the anisotropies of the AGWB (power spectrum and cross-correlations with other cosmological probes). We stress that three main building blocks are present: choice of a cosmological framework, description of large scale structures and astrophysical modeling of the GW sources. While the first two are standard lore in cosmology, we focus on the sub-galactic astrophysics which is less constrained today. Section III describes the general properties of the angular power spectrum of anisotropies. Using the Limber approximation, we provide an analytic approximation of the angular power spectrum of anisotropies, useful to derive order of magnitude estimates and to understand our numerical results. We then propose a general derivation of shot-noise and we explain that cross-correlation with galaxy number counts can help to extract a map of AGWB anisotropies, even for shot-noise dominated background maps. We then discuss the frequency-direction factorization hypothesis used in current directional searches (e.g. by LIGO-Virgo). We demonstrate that this approximation fails in capturing the physics in the upper part of the LIGO-Virgo frequency band. Section IV defines the building blocks of the astrophysical modeling. It describes the normalisation of the model, which is necessary to have a meaningful model comparison. Section V explores the signature of these models on the monopole, the angular power spectrum and the various cross-correlations. In Section VI this analysis is extended to include the contribution of NS mergers and in Section VII to the LISA frequency band.

II Anisotropies of AGWB

The dimensionless energy density of the GW background per unit of solid angle, $\text{d}^{2}{\bf{e}}$, and logarithmic frequency, $\text{d}f/f$, can be split into an homogenous and isotropic component and a directional dependent one as

[TABLE]

where $\rho_{c}$ is the critical energy density of the universe and $\rho_{\rm GW}$ is the energy density of GW. Using the standard expression for the energy density in terms of the wave amplitude, see e.g. Ref. Maggiore (2007), and recalling that the definition of energy requires an average over several periods of the wave, we find

[TABLE]

where $T_{{}_{\rm O}}$ comes from the time average and represents the period of observation of the detector.777The dependence on $T_{{}_{\rm O}}$ may appear strange at first, but note that for a continuous background signal centered around a frequency $f$ with width $\Delta f$ one has for $T_{{}_{\rm O}}\gg 1/f$, $T_{{}_{\rm O}}\propto\int_{0}^{T_{{}_{\rm O}}}|h_{A}(t)|^{2}\text{d}t\simeq\int|h_{A}(f^{\prime})|^{2}\text{d}f^{\prime}\simeq\Delta f|h_{A}(f)|^{2}$ and hence $|h_{A}(f)|^{2}\propto T_{{}_{\rm O}}$.

In this section we describe the generic properties of monopole and anisotropies. To that purpose, we rely on the reference model, fully defined in § IV.3 below. We first recall the computation of the isotropic contribution (§ II.1), of the power spectrum of anisotropies (§ II.2) and of the cross-correlations with other cosmological probes (§ II.3). We conclude in § II.4 with a schematic illustration of the general computation strategy and numerical implementation, focusing on the different scales that enter in the discussion.

II.1 Isotropic contribution

The background component in Eq. (II) is given by the line of sight integration on our past lightcone

[TABLE]

where $\eta_{*}$ stands for a maximal distance (or, equivalently, maximal redshift) above which there are no astrophysical sources. For future convenience, we introduce the redefinition

[TABLE]

with

[TABLE]

where $\mathcal{L}_{{}_{\rm G}}$ is the effective GW luminosity of a galaxy per unit of emitted frequency, $f_{{}_{\rm G}}$, characterized by the set of parameters $\theta_{{}_{\rm G}}$ (mass, metallicity…). This effective luminosity has been introduced in Ref. Cusin et al. (2017). In the galaxy rest-frame, it represents the sum of the luminosity of the GW emitted by all astrophysical sources contained in that galaxy, averaged on the distribution function of their peculiar velocity. As shown in Ref. Cusin et al. (2017), at linear order, the effects of the peculiar motion of a source in its host galaxy can be neglected on average. The relation between the effective luminosity and the emitted strain is established in Eqs. (79) and (80) of Ref. Cusin et al. (2018b). In Eq. (5), $n_{{}_{\rm G}}$ is the comoving number density of galaxies, $a$ is the scale factor of the Friedman-Lemaître spacetime normalized to 1 today and $\eta$ its comoving time. The frequencies at emission and observation are related by

[TABLE]

where $z_{{}_{\rm G}}$ is the redshift.

Figure 1 depicts the background contribution for the reference astrophysical model. It shows that the low frequency part is well-fitted by a power law $\propto f^{2/3}$. Figure 2 presents the astrophysical kernel $\mathcal{A}(z,f)$ defined in Eq. (4) as a function of redshift (left panel) and frequency (right panel).

II.2 Anisotropies

The anisotropic component of the background energy density of Eq. (II) is given by

[TABLE]

where the astrophysical kernel $\mathcal{A}$ is defined in Eq. (5) while

[TABLE]

in which the relation between the frequencies at emission and observation needs only to be evaluated at lowest order so that Eq. (6) reduces to $f_{{}_{\rm G}}=f/a$.

In Eq. (II.2), $n_{{}_{\rm G}}$ stands for the comoving number density of galaxies, ${\Psi}$ for the gravitational potential, $v$ for the comoving velocity field and $\delta_{\rm G}$ is the galaxy over-density. This latter is related to the dark matter over-density, $\delta_{\rm m}$, by the bias $b$, defined in comoving gauge, such that

[TABLE]

where ${\cal H}\equiv\text{d}\ln a/\text{d}\eta$ is the comoving Hubble parameter.

The statistical properties of $\delta\Omega_{\rm GW}({\bf{e}},f)$ are first encoded in its angular correlation function or, equivalently, in its angular power spectrum

[TABLE]

where ${\delta\Omega}_{\ell}(k,f)$, derived in Ref. Cusin et al. (2017), is given by

[TABLE]

where $\Delta\eta=\eta_{0}-\eta$, $j_{\ell}$ stands for the spherical Bessel function, $k$ is the wavenumber and $X_{k}(\eta)$ stands for the Fourier modes of $X$ and it has been assumed that the bias has no scale dependence.

To conclude, we introduce the reduced angular power spectrum by normalizing over the monopole as

[TABLE]

From now on we will refer to this dimensionless quantity as the angular power spectrum of relative anisotropies.

II.3 Cross-correlation

Since the AGWB anisotropy depends on cosmological perturbations, see Eq. (II.2), it correlates with any other cosmological probe, such as galaxy number counts and weak lensing convergence. The cross-correlation power spectra have been presented in Ref. Cusin et al. (2017),

[TABLE]

For weak lensing, $X_{\ell}$ is the cosmic convergence

[TABLE]

so that $B_{\ell}\equiv B_{\ell}^{\kappa}$ is given by

[TABLE]

In these expressions, $\chi_{H}$ corresponds to the maximal depth of a given survey and

[TABLE]

where the function $p_{\chi}$ is the sources distribution function. For the cross-correlation with galaxy, $X_{\ell}$ is number counts $\Delta_{\ell}(k,z)$, defined e.g. in Eq. (44) of Bonvin and Durrer (2011) so that $D_{\ell}\equiv B_{\ell}^{\Delta}$

[TABLE]

For our reference astrophysical model, the cross-correlations with weak lensing convergence (for both SKA and Euclid source distributions)888For SKA, we use the source distribution of Ref. Andrianomena et al. (2014), and we verified that the results for the cross-correlation are quantitatively similar to what one would obtained using the SKA2 source distribution Ref. Harrison et al. (2016). and with galaxy number counts are presented in Fig. 3. Note that cross-correlating with galaxy number counts at different redshifts is equivalent to filtering the astrophysical kernel $\mathcal{A}$ with a window function that selects different redshift bins. The shift of the peak of $D_{\ell}$ to higher multipoles as we consider higher redshift bins can be understood from the Limber relation between redshift and multipoles, as detailed in § III. Interestingly, the cross-correlation with galaxy number counts can help to reconstruct the astrophysical kernel $\mathcal{A}$ as a function of redshifts. Furthermore, the study of the cross-correlation with galaxy number counts is useful to distinguish in observations a AGWB from cosmological backgrounds, which are not expected to be correlated with the galaxy distribution.

II.4 Coarse graining approach

As can be seen for the previous sections, three main ingredients enter the computation of the angular power spectrum. They are related to the three main building blocks describing the astrophysics of GW sources (${\cal A}$ and ${\cal B}$), the large scale distribution of the sources (through the cosmological variables) and the properties of the cosmological model. We can thus distinguish three scales in the problem: 999We stress that we refer here to different physical processes that take place on different scales. These processes impact the shape of the power spectrum on all angular scales.

cosmological scale: the large scale structure can be effectively described by cosmological scalar perturbations in the metric and the matter distribution; structures are then assumed to move with the cosmic flow; 2. 2.

galactic scale: each galaxy is characterized by a set of parameters $\theta_{{}_{\rm G}}$ (mass, metallicity…) and an effective GW luminosity resulting from the contributions of the various sources it contains. The galaxy number density is computed from the halo mass function; 3. 3.

sub-galactic scale: different classes of GW sources (binary compact objects, rotating neutron stars, etc.) are characterized by parameters (masses, orbital parameters…). Each source emits an energy spectrum which depends on these parameters, in a typical range of frequencies.

This approach is schematically summarized in Fig. 4. It illustrates how the effective GW luminosity associated with each galaxy is first obtained from the properties of the galaxy by integrating over all its GW sources and then integrated over the halo mass distribution to give the astrophysical functions $\mathcal{A}$ and $\mathcal{B}$. These two functions then enter the equation for the background energy density (II) and the master equation (II.2) for the anisotropies. To finish, anisotropies depend also on the cosmological perturbation variables which need to be evolved during the cosmic history, hence depending on both the initial primordial power spectrum and a set of transfer functions. This last step is nowaday part of the standard lore of cosmology. Putting all these ingredients together, the angular power spectrum can be computed along with the cross-correlation with the various cosmological probes.

In the present work, we use the standard $\Lambda$CDM cosmological model in which the universe is described by a Friedmann-Lemaître spacetime with perturbations that describe the large scale structure. In the theory of cosmological perturbations, any variable, $X(\eta,x^{i})$ say, is a stochastic field. It can be decomposed in Fourier modes, $X(\eta,{\mathbf{k}})$, which can be expressed as the product of a transfer function and of the initial metric perturbation: $X(\eta,{\mathbf{k}})=X_{k}(\eta)\Phi^{P}({\mathbf{k}})$. The power spectrum of $\Phi^{P}({\mathbf{k}})$ is predicted e.g. from inflation and constrained from CMB analysis. We use Planck satellite Ade et al. (2016) cosmological parameters. Linear transfer functions are obtained from CMBquick CMB and we use Halofit Smith et al. (2003) to account for the non-linearities in the matter power spectrum.

Figure 6 presents the galaxy correlation function used in this work, at a redshift $z=0.6$. It assumes that galaxies follow the evolution of the underlying dark matter field so that they are related by a bias function. We assume a scale-independent bias model scaling as $\propto\sqrt{1+z}$ Marin et al. (2013); Rassat et al. (2008). Explicitly,

[TABLE]

At $z=0.6$ the bias is $1.8$, and one can check that our correlation function is consistent with the one of SDSS VIPER, see Figs. 3 and 4 of Ref. Anderson et al. (2012). Since very small scales do not contribute to the final power spectrum, as we will demonstrate in § III, the use of a more refined scale-dependent model for the bias would not significantly affect our results. As already mentioned, the non-linear evolution of the density growth has been taken into account by using the Halofit approximation Smith et al. (2003).

A detailed description of how this work treats the sub-galactic physics is postponed to § IV.

III General properties of the angular power spectrum

Before turning to the definition and investigation of the imprint of various astrophysical models, we take some time to describe some general features of the angular power spectrum of the AGWB anisotropies. Again, in order to illustrate our purpose, we use our reference astrophysical model described below in § IV.3. The goal of this section is three-fold. First in § III.1 we describe an analytic approximation that will allow us to understand our numerics, then in § III.2 we discuss the shot noise contribution that will inevitably limit our predictions. To conclude, in § III.4 we investigate to which extent the direction-frequency factorization is a good hypothesis.

III.1 Analytic approximation of the angular power spectrum

Keeping only the dominant contribution in Eq. (II.2) and using the expression (4) of the astrophyical kernel, we get

[TABLE]

Hence the angular power spectrum simplifies to

[TABLE]

where $P_{\rm Gal}(k)$ is the galaxy power spectrum and $\Delta\eta=\eta_{0}-\eta$. The Limber approximation LoVerde and Afshordi (2008); Bernardeau et al. (2011) can then be used to derive the slope of the angular power spectrum. One method consists in noticing that for any test function $f(x)$

[TABLE]

It follows that

[TABLE]

where $k$ and $\Delta\eta$ must satisfy the Limber constraint

[TABLE]

Note that the integrand function can be thought as the product of $kP_{\rm Gal}$ and of the window function $\left|\partial_{\eta}\bar{\Omega}_{\rm GW}(\eta)\right|^{2}$ evaluated at the value of $\eta$ satisfying the constraint (23). In Fig. 7 we plot the galaxy power spectrum today together with the window functions at different values of $\ell$. It clearly shows that the window function selects different areas below the power spectrum. For sufficiently small $\ell$, the dominant contribution to the integral comes from the peak of the function $kP_{\rm Gal}$, around $k\sim 0.05$ Mpc*-1* corresponding to length scale of $\sim 120$ Mpc. For larger $\ell$, the peak of the power spectrum is cut out and the integral is dominated by large modes.

Formally, this is similar to the computation of CMB angular power spectrum. However for the CMB the visibility function is sharply peaked around the recombination, hence selecting (for different multipoles $\ell$) a very narrow region of the power spectrum. For the GW background, the visibility function extends typically from $z\simeq 4$ (corresponding to a comoving distance of order $7000$ Mpc) down to $z=0$ as depicted in Fig. 2. It follows that there is no direct relation between a wavemode $k$ and an angular mode $\ell$. For a fixed observed frequency, a given wavemode contributes to all multipoles such that $\ell\lesssim\ell_{\rm max}(k)$, where $\ell_{\rm max}(k)$ is related to $z_{\rm max}$ (the typical maximum redshift for sources observed at that frequency) by

[TABLE]

For a given multipole $\ell$ and a given wavemode $k$, if the corresponding distance is too large, the number of GW sources is suppressed and so is the $C_{\ell}$. Fig. 7 illustrates the contributions of various bins of $k$ to the total signal. It is clear that for each bin in $k$ the contribution scales as $1/\ell$ on large scales, and drops beyond $\ell_{\rm max}$.

It is useful to introduce a further assumption, namely that the emission depends mildly on redshift. More precisely, using Eq. (22) and assuming that we can ignore the time variation of $\partial_{\eta}\bar{\Omega}_{\rm GW}$, we find that the multipoles are approximated by

[TABLE]

where $k_{\rm min}(\ell)$ is set by the fact that there is a maximum distance $r_{\rm max}$ at which we can find GW sources and thus a minimum Fourier mode set by the Limber constraint (23). We note that with our galaxy power spectrum which describes correctly the large scales, and thus the small Fourier modes, the proportionality relation (III.1) is insensitive to $k_{\rm min}$ and one can replace it by [math]. Indeed for $k<k_{\rm eq}$ (with $k_{\rm eq}\simeq 0.01\,{\rm Mpc}^{-1}$ the Fourier mode entering the horizon at matter-radiation equivalence), $P_{\rm Gal}(k)\propto k^{\alpha}$ with $\alpha\simeq 1$. In particular we find that $C_{\ell}\propto\ell^{-1}$.

An equivalent formulation of the Limber approximation is given by

[TABLE]

where we have used the fact that the comoving distance and conformal time are related by $r=\Delta\eta$, hence the Limber constraint can be rewritten as $k\,r=\ell+1/2$. In this case, the function $\left|\partial_{r}\bar{\Omega}_{\rm GW}(r)\right|^{2}$ is fixed while the power spectrum has to be evaluated at $k$ satisfying the Limber constraint, for different $\ell$. The left panel of Fig. 8 shows this fixed window function and the power spectrum evaluated at the Limber constraint for $\ell=1,10,100$.

Figure 8 shows that a given bin in $r$ contributes more efficiently to multipoles such that $\ell\sim 1/(rk_{\rm peak})$. This can also be understood by separating the contributions coming from various bins of distance, or equivalently various bins of redshift, as in the right panel of Fig. 8. As expected, the lowest redshifts, corresponding to the shortest comoving distances, contribute the most on large scales (small $\ell$).

III.2 General treatment of shot noise

When working with galaxy data (e.g. with a galaxy catalogue) it is important to keep in mind that the background angular power spectrum suffers a shot noise component which adds to the theoretical predictions,

[TABLE]

where “exp” indicates the angular power spectrum computed from data, $C_{\ell}$ is the theoretical power spectrum (i.e. what we compute in this work) and $S_{n}$ denotes the shot noise contribution (see e.g. Ref. Jenkins and Sakellariadou (2019)). This contribution is flat in $\ell$ space and gives a constant offset to the angular power spectrum

[TABLE]

The derivation of this result can be found in appendix A. There are two important points to keep in mind: (1) this is an offset. Since the angular power spectrum decreases with $\ell$, it will affect more large $\ell$ and (2) the prediction (28) diverges for $r=0$, hence the contribution of the Poisson noise depends on the cut-off used to regularize this integral. From an observational point of view, the physical quantity on which the cut-off has to be set is the observed flux: sources with a flux bigger than a given threshold can be resolved and are therefore filtered out.

Using the fact that the flux $\Phi$ received per unit of frequency from a source at redshift $z$ is related to the luminosity per unit of emitted frequency by $\Phi(f)=1/(4\pi)\mathcal{L}_{{}_{\rm G}}/(1+z)$, we see that an upper bound on $\Phi$ defines the region of integration in the plane $(z,\mathcal{L}_{{}_{\rm G}})$. In other terms, we introduce a selection function which is 1 for $\Phi<\Phi_{\text{cut}}$, and 0 otherwise, i.e. we multiply the integrand function in Eq. (28) by the selection function

[TABLE]

Of course, if one can assume that all galaxies have the same luminosity, then the cut-off in the flux translates directly into a lower cut-off in redshift (or analogously in $r$). The same selection function has to be applied to the integral defining the theoretical curve of the angular power spectrum. A detailed derivation of the shot noise component and comments on the regularization procedure can be found in appendix A.

We conclude this section with a few remarks. First of all, we observe that when one derives a prediction for the angular power spectrum using a galaxy catalogue, this predictions contains both the ”theoretical part” and the shot noise components. Since the angular power spectrum decreases with $\ell$ and the shot noise component is an offset, the high-$\ell$ part of the angular power spectrum will be dominated by the shot-noise part of the result. This explains the shape of the curve of Ref. Jenkins et al. (2018a), where a simulated galaxy catalogue is used: at high $\ell$, $C_{\ell}\sim\text{cnst}$ indicating the fact that shot noise dominates over the signal for those multipoles. Second, the part of the angular power spectrum dominated by Poisson noise also contains astrophysical information so that an adequate understanding and modeling of both the ”theoretical” and shot noise components is necessary to extract astrophysical quantities out of future observations.

We emphasize that the shot noise we explore in this section is due to the discreteness of the spatial distribution of galaxies. This shot noise component (that we will refer to as spatial shot noise) is always present, and has to be added to the theoretical predictions for the angular power spectrum in any frequency band. It is due to the fact that when deriving theoretical predictions, we treat the galaxy number density as a continuous field. The same type of shot noise is present for the case of the CIB, see e.g. Ref. Ade et al. (2014).

When the signal is dominated by popcorn-like events (i.e. events with a short duration with respect to the time of integration, and no time overlap), there is a second contribution to shot noise coming from the “discreetness” of events in time, which has to be added to the spatial shot noise. This second component of the shot noise is present in the LIGO frequency band where BH and NS mergers give the dominant contribution. It is absent in the LISA frequency band, where the background comes from radiation emitted by binary systems of compact objects in the inspiraling phase, and can be treated as a continuous (almost) stationary background. To take into account in an effective way this time-like shot noise component, one can multiply the factor $n_{G}$ at the denominator in Eq. (28) by the fraction of galaxies which contain a merger in the observation time.

III.3 Shot-noise and cross-correlation

As observed in Ref. Jenkins and Sakellariadou (2019), in the LIGO-Virgo frequency band the contribution of temporal shot noise to the angular power spectrum of anisotropies is dominating the signal, for observation times of the order of years. We suggest here that using cross-correlation with galaxy number counts can help in overcoming the problem. By considering the cross-correlation with galaxy number counts it may be possible to extract anisotropies of the AGWB even if the map of anisotropies is shot-noise dominated. We sketch here the derivation. A future work will be dedicated to a detailed study of the properties of the cross-correlation map.

We schematically write galaxy number counts and AGWB anisotropies as

[TABLE]

where we denote as $N_{G}$ the number of galaxies in a given position in the sky and as $\bar{N}_{G}$ its spatial average. $n_{G}$ is the number of galaxies which contain a merger in the time $T_{o}$, and $\bar{n}_{G}=f\bar{N}_{G}$ where the fraction $f\leq 1$. The quantity $L$ is a typical galaxy luminosity. We assume that the variable number of galaxies $N_{G}$ follows a Poisson distribution with average $\bar{N}_{G}$ and variance $\bar{N}^{2}_{G}$. Using basic properties of Poisson distributions it is easy to verify that

[TABLE]

where $P$ stays for Poisson. 101010To derive the expression for the cross-correlation we used the fact that $N_{G}=n_{G}+m_{G}$ where $m_{G}$ is the number of galaxies which do not contain a merger over the time of observation. The variables $n_{G}$ and $m_{G}$ follow independent Poisson distributions with average $\bar{n}_{G}=f\bar{N}_{G}$ and $\bar{m}_{G}=(1-f)\bar{N}_{G}$. Then $\langle(n_{G}-\bar{n}_{G})(N_{G}-\bar{N}_{G})\rangle=\langle n_{G}N_{G}\rangle-f\bar{N}_{G}^{2}=\langle n_{G}n_{G}\rangle+\langle n_{G}\rangle\langle m_{G}\rangle-f\bar{N}_{G}^{2}=Var(n_{G})=f\bar{N}_{G}$.

We can compute now the signal part of these cross-correlations. We simply use the fact that we can rewrite Eqs. (30) and (31) as $\Delta_{N}=\delta_{G}$ and $\delta\Omega_{GW}=Lf\bar{N}_{G}\delta_{G}$ and we use that the galaxy over density $\delta_{G}$ is a stochastic variable. Then we have

[TABLE]

where $S$ stays for signal. Then comparing Eqs. (32) and (35) we can have a rough estimate of the ratio signal over Poisson-noise for the autocorrelation $\langle\delta\Omega_{GW}\delta\Omega_{GW}\rangle$ and cross correlation $\langle\delta\Omega_{GW}\Delta_{N}\rangle$ We have

[TABLE]

in other words the signal to noise of cross-correlation is boosted with respect to the one of the auto-correlation of a factor $1/f\gg 1$. We stress that this is a simplistic derivation to illustrate the idea of using cross-correlations to overcome the shot-noise problem. A full and realistic study of cross-correlation map will be presented in a future work.

III.4 Frequency-direction factorization

Searches for anisotropies (e.g. at LIGO-Virgo) usually rely on the assumption that the frequency and the direction dependencies of the background energy density can be factorized LIGO Scientific Collaboration and Virgo Collaboration (2017); The LIGO Scientific Collaboration and the Virgo Collaboration (2019b), i.e. that

[TABLE]

where

[TABLE]

with $f_{\text{ref}}$ some reference frequency and $\alpha=2/3$ for a background from merging compact binaries. The factorization assumption used in the LIGO-Virgo analysis The LIGO Scientific Collaboration and the Virgo Collaboration (2019b) relies on the models in Ref.Jenkins et al. (2018a, b). Note that since $H(f)$ is dimensionless in natural units, the quantity $P({\bf{e}})$ in Eq. (40) has dimensions of a time. The angular power spectrum of the direction-dependent factor in Eq. (40) is then defined using the decomposition

[TABLE]

as

[TABLE]

Observational constraints are then set on the dimensionless spectrum

[TABLE]

It is easy to verify that the relation between this spectrum and the angular power spectrum defined above in Eq. (10) is

[TABLE]

Since the angular spectrum on the r.h.s. of Eq. (45) does not depend on frequency while the one on the l.h.s. does, a convenient way to check the consistency of the factorization assumption (40) is to verify that the frequency dependence cancels on the r.h.s., i.e. that the angular power spectrum for mergers predicted by our model scales with frequency as

[TABLE]

with $\alpha=2/3$. If the factorization hypothesis is valid, the relative fluctuation should not depend on frequency. Figure 10 shows that the scaling with frequency is not exactly a power law. Analogously, for a fixed value of frequency, for larger $\ell$ the assumption is slightly worse than for small $\ell$. At a frequency of 80 Hz the error one makes assuming a power law is of the order of 20%.

IV Astrophysical models

We are now in the position to fully investigate different astrophysical scenarios. In § IV.1 we start by describing the computation of the astrophysical kernel and then discuss the normalisation procedure in § IV.2. We then describe the reference model (§ IV.3) that was used so far in this article. In § IV.4 we present a series of modified models to investigate the parameters that have an effect on the angular power spectrum of anisotropies.

IV.1 Model framework

The astrophysical sources of GW - such as merging binary BHs and NSs, spinning NSs and supernovae - reside in galaxies and therefore reflect the processes of galactic evolution and star formation. In this work we concentrate on the background from merging binary BHs.

In order to compute the astrophysical kernels $\mathcal{A}$ and $\mathcal{B}$ defined in Eqs. (5) and (8) we need to compute the GW luminosity $\mathcal{L}_{{}_{\rm G}}(z,f_{{}_{\rm G}},M_{{}_{\rm G}})$ for each galaxy as a function of its halo mass $M_{{}_{\rm G}}$ and redshift $z$. We can then sum over the entire galactic population, where the number densities are provided by the halo mass function $\text{d}n/\text{d}M_{{}_{\rm G}}(M_{{}_{\rm G}},z)$. Our computation follows the formalism we developed in Refs. Dvorkin et al. (2016a, b); Dvorkin et al. (2018), which we now briefly describe.

The first step is to calculate the star formation rate (SFR) $\psi(M_{{}_{\rm G}},t)$, given in units of $M_{\odot}/{\rm yr}$ and the stellar-to-halo mass ratio of a given galaxy using a modified version of the abundance-matching relations of Ref. Behroozi et al. (2013). We use a Salpeter-like initial mass function (IMF) Salpeter (1955) to describe the number of stars per unit total stellar mass formed,

[TABLE]

where $M_{*}$ is the mass of the star at birth.

Having described the total mass of stars formed in each galaxy as a function of time, we also need to model the evolution of massive stars and the nature and mass of their remnants. We assume that the latter depends only on the mass of the progenitor star $M_{*}$, and on its metallicity $Z$ and is encoded in the function $m=g_{s}(M_{*},Z)$, to be specified for each model. Typically, massive stars ($M_{*}\gtrsim 8M_{\odot}$) explode as supernovae or collapse to form BH on a timescale of a few Myr. If we assume such short stellar lifetimes, the stellar metallicity tracks the metallicity of the interstellar medium (ISM) given by $Z=Z(M_{{}_{\rm G}},z)$. We adopt the observational relation of Ref. Ma et al. (2016) for the ISM metallicity as a function of galaxy stellar mass and redshift. We also introduce a cut-off at high BH masses $M_{\rm co}$ which can arise due to pair-instability supernovae, as explained below.

Under these assumptions, the instantaneous BH formation rate at a given cosmic time $t$ (or, equivalently, redshift $z$) for a galaxy with halo mass $M_{{}_{\rm G}}$, in units of events per unit BH mass $m$, is given by

[TABLE]

where $M_{*}(m)$ and $\text{d}M_{*}/\text{d}m$ are deduced from the relation $m=g_{s}(M_{*},Z)$ (we assume negligible stellar lifetimes). We then assume that only a fraction $\beta$ of these BHs resides in binary systems that merge within the age of the Universe, so that the rate of formation of the latter is

[TABLE]

As we will show later, this overall factor $\beta$ is used to normalize our model with respect to the number of events observed by aLIGO/aVirgo. Following Ref. Dvorkin et al. (2018), the birth rate of binaries with component masses $(m,m^{\prime}\leq m)$ is

[TABLE]

where the distribution function of binary masses $P(m,m^{\prime})$ is normalized so that $\int{\cal R}_{2}(m){\cal R}_{2}(m^{\prime})P(m,m^{\prime})\text{d}m\text{d}m^{\prime}=0.5\int{\cal R}_{2}(m)\text{d}m$.

The merger rate depends on the time to coalescence of the binaries which can be expressed as a function of the distribution of the orbital parameters $P(a_{\rm f},e_{\rm f})$ at the time of formation. Since BH binaries are expected to circularize due to gravitational wave radiation reaction before reaching the LIGO/Virgo band (see e.g. Ref. Lower et al. (2018)), we assume circular orbits in what follows, so that the only distribution left to determine is $f(a_{\rm f})$. Hence, the birth rate of BH binaries (per unit mass squared per unit time and per unit $a_{\rm f}$) is

[TABLE]

from which we deduce that the merger rate at time $t$ is

[TABLE]

with $\tau_{\rm m}(m,m^{\prime},a_{\rm f})$ the merger time of the system $(m,m^{\prime},a_{\rm f})$.

The GW luminosity of the galaxy is then

[TABLE]

We then need to sum over the entire galactic population as in Eq. (5), where the integral over the galactic properties $\theta_{{}_{\rm G}}$ reduces to an integration over $M_{{}_{\rm G}}$ weighted by the halo mass function given in Ref. Tinker et al. (2008). The result is the quantity $\mathcal{A}$ presented on Fig. 2 (for the reference model defined in § IV.3 below).

IV.2 Normalization to the number of detected events

The overall normalization parameter $\beta$ is adjusted so as to match the number of detections by aLIGO/aVirgo during the O1+O2 observing runs The LIGO Scientific Collaboration and the Virgo Collaboration (2018a). We therefore require that all our models result in $10$ detectable events over the span of the O1+O2 observation time. In our estimate of the detection rates we follow Ref. Dvorkin et al. (2018), namely we calculate the signal-to-noise rate (SNR) for each binary BH merger produced in the model:

[TABLE]

where the index $i$ refers to either Hanford or Livingston detector, $h(f)$ is the GW strain in the observed frequency domain and $S_{n,i}(f)$ are the O2 noise power spectral densities of the corresponding detectors Collaboration (a, b) and we use the correction factor in Ref. Finn and Chernoff (1993) to account for different source orientations. The strain $h(f)$ is calculated using the PhenomB template Ajith et al. (2011) and assuming zero spins. We define observed events as those with $\rho>8$. The number of sources detectable during O1+O2 is given by multiplying the detection rate by the total observation time $T_{\rm obs}=169.7$ days.

The resulting values of $\beta$ for each of the models discussed below is shown in Table 1. This already shows that while all the astrophysical models are adjusted to predict $10$ detected events during the aLIGO O1+O2 timespan, they predict that the fraction of BH in binaries can range from 0.6% to 17%. It can be anticipated that this differences will imprint the AGWB. Also note that while individual mergers are resolved only at low $z$, the AGWB is affected also by higher redshifts. Therefore, even though all models are calibrated to the same number of resolved sources, the resulting AGWB may vary if the high-redshift population of sources differs among the models.

IV.3 Reference model

We start by describing the parameter choices we made for our Reference model, used so far in this article.

The IMF slope for this model is set to $p=2.35$. 2. 2.

The BH formation model was chosen as the ‘delayed’ model in Ref. Fryer et al. (2012). Specifically, we used the functional form provided by Ref. Fryer et al. (2012) to calculate the function $m=g_{s}(M_{*},Z)$. This choice affects the distribution of BH masses and will be discussed in Sec. IV.4.1. We also introduced a cutoff mass of $M_{\rm co}=45M_{\odot}$ in the BH mass distribution. The recent analysis The LIGO Scientific Collaboration and the Virgo Collaboration (2018b) of the population of BHs detected by aLIGO/aVirgo suggests a cutoff at this value, although further observations are needed to confirm it. We will discuss the possible causes of this cutoff in Sec. IV.4.2. 3. 3.

We assume $P(m,m^{\prime})=$ cnst for the distribution of masses in the binaries. 4. 4.

We assume that the distribution of the semi-major axis at formation is $f(a_{\rm f})\propto a_{\rm f}^{-1}$ with cut-off at $a_{\rm f,min}=0.014$ AU and $a_{\rm f,max}=4000$ AU. The lower bound was chosen so as to ensure that the lightest BH binaries in our model ($5M_{\odot}-5M_{\odot}$) merge within a Hubble time.

This reference model predicts that the merger rate observed by LIGO/Virgo can be explained provided $\sim 1\%$ of BHs reside in binaries that merge within the Hubble time. We stress that this estimate relies on the assumption of isolated stellar evolution, i.e. that BH formation is not influenced by binary interactions. This may not be the case of close stellar binaries, where mass exchange processes and in particular the co-evolution of the primary compact object and its stellar companion during the common envelope phase may play a major role in the later stages of the secondary evolution Belczynski et al. (2016a); Dominik et al. (2012); Belczynski et al. (2016b). These processes and their effects on the stochastic background will be further explored in future work.

IV.4 Modified models

In order to explore the astrophysical dependencies of the AGWB anisotropies, we consider several models, each varying from the reference model described above in one key aspect, keeping the others fixed. In addition to varying the corresponding parameter, we also need to change the overall efficiency factor $\beta$, as explained above, so that all of the models discussed here results in the same total number of detectable events.

IV.4.1 BH formation model [limongi model]

BH masses depend on the properties of their stellar progenitors, in particular the mass prior to core collapse and chemical composition, as well as other parameters such as the rotation velocity, see e.g. Refs. O’Connor and Ott (2011); Mandel and de Mink (2016); Marchant et al. (2016); Limongi (2017). The formation of binaries may further depend on such processes as common envelope evolution, dynamical processes in stellar clusters and evolution of hierarchical triple systems, see e.g. Refs. Belczynski et al. (2016a); Dominik et al. (2012); Belczynski et al. (2016b); Rodriguez et al. (2016).

In this article we assume that binary formation process is encoded in the efficiency parameter $\beta$ and the distribution of merger time delays. Furthermore, we explore only one aspect of this complex problem, namely the evolution of isolated massive stars. Recent studies Fryer and Kalogera (2001); Fryer et al. (2012); Blondin et al. (2003) suggest that the explosion is powered by neutrinos stored behind the shock and that the explosion energy depends on neutrino heat transport mechanisms, the nature of the hydrodynamic instabilities that convert neutrino thermal energy into kinetic energy that can power the supernova, and the resulting time delay between shock bounce and explosion.

The Reference model uses the description by Ref. Fryer et al. (2012) which provides an analytic model for a neutrino-driven explosion and calculate the explosion energy, as well as the remnant mass, using numerical pre-collapse stellar models from Ref. Woosley et al. (2002).

Another set of stellar evolution models is provided in Ref. Limongi (2017). These models differ from the ones in Ref. Fryer et al. (2012) in two aspects. First, Ref. Limongi (2017) uses a different set of pre-collapse stellar models which vary from Woosley et al. (2002) in their treatment of convection, mass-loss rate and angular momentum transport. Second, Limongi (2017) assumed a constant explosion energy in the calculation of the remnant mass, contrary to Fryer et al. (2012). As was shown in Dvorkin et al. (2018), these models predict different mass distributions of detectable BHs. In the following, the model limongi uses the model described in Ref. Limongi (2017) without stellar rotation.

The parameter $\beta$ derived for this model is very similar to the one in the Reference model, but the mass distribution of the BHs is different, and will affect the resulting AGWB, as we will show below.

IV.4.2 BH mass cutoff due to PISN [dMco and uMco models]

Very massive stars (typically in the range $[130-260]M_{\odot}$) are unstable to electron-positron pair creation which may lead to pair-instability supernova (PISN) that disrupt the entire star. In this case, no BH is formed Kasen et al. (2011). Although direct observational evidence is lacking, the absence of BHs in the mass range $[60-260]M_{\odot}$ may provide an indirect confirmation of this effect. A cutoff in BH mass may be present at even lower masses due to pulsational PISN, where the instability causes short episodes of mass ejection followed by periods of quiescent evolution Woosley et al. (2007); Woosley and Heger (2015). As a result, the stellar mass is reduced below the limit of the onset of the instability, and it was suggested in Ref. Talbot and Thrane (2018) that this process may lead to an excess of BHs around $\sim 40M_{\odot}$. Recent analysis of the aLIGO/aVirgo events detected during O1+O2 observational runs The LIGO Scientific Collaboration and the Virgo Collaboration (2018b) provides a tentative measurement of the BH mass cutoff at $M_{\rm co}=45M_{\odot}$ which can be due to PISN . As we will show in what follows, this mass cutoff has an important effect both on the isotropic and anisotorpic stochastic background, suggesting an alternative way to measure this effect.

In order to explore the sensitivity of the stochastic background to the PISN-induced mass cutoff we varied to $M_{\rm co}=40M_{\odot}$ [dMco model] and $M_{\rm co}=50M_{\odot}$ [uMco model].

IV.4.3 Stellar initial mass function [imf-high and imf-low models]

The reference model assumes a Salpeter-like IMF with slope $p=2.35$. Interestingly, some studies show that the IMF slope may not be universal (see e.g. the discussion in Ref. Krumholz (2014)), for example a recent a hint to a more shallow IMF in the Large Magellanic Cloud Schneider et al. (2018). In order to estimate the influence of the IMF we explored two models where the slope was taken to be $p=2.6$ [imf-high model] and $p=2.1$ [imf-low model].

IV.4.4 Distribution of initial separations [aconst model]

The reference model assumes the initial separation of the BHs is distributed like $P(a)\propto a^{-1}$. This separation then translates into a distribution of merger delay times, favoring short delay times. We consider the extreme scenario [aconst model] of a flat distribution of the initial separations $P(a)\sim const$, which results in longer delay times.

IV.4.5 Metallicity of progenitor stars

Metallicity plays an important role in the evolution of massive stars, in particular, high-metallicity stars experience strong winds throughout their lives. As a result, the remnant mass is reduced relative to the low-metallicity case Belczynski et al. (2010); Fryer et al. (2012); Spera et al. (2015). In the reference model the metallicity is evolved with the stellar mass following the observational relations of Ref. Ma et al. (2016). To test the effect of metallicity on the stochastic background, we used a model with constant metallicity of $Z=10^{-3}Z_{\odot}$, which leads to the formation of heavier BHs.

V Results

We have now defined all the quantities that are required to go through the general computation described in Fig. 4 to compute the monopole (§ V.1), the power spectrum (§ V.2), and the cross-correlations (§ V.3) for the different astrophysical models.

V.1 Monopole

Figure 11 shows the AGWB monopole as a function of frequency for the various models described in the previous section, for both the LISA and LIGO/Virgo frequency bands. The left panel depicts the results for the 9 individual models discussed above (the AGWB is multiplied by a power-law $f^{-2/3}$ to accentuate the differences in amplitude). The amplitude varies by a factor of $\sim 2$ between the different models, reflecting the differences in the normalization $\beta$ and source mass distribution.

The right panel of Fig. 11 compares the reference model (black) to the contributions of low-mass $<25M_{\odot}$ and high-mass $>25M_{\odot}$ BH populations in red. Both populations are derived from the reference model, where we took into account only the low-mass (high-mass) BHs, respectively, and neither of these models includes ’mixed’ (low-mass/high-mass) binaries. As a result, the two red curves do not sum up to the reference model. It can be seen that the AGWB signal is dominated by the high-mass binaries, even though they are outnumbered by the low-mass ones as a result of the power-law stellar mass function.

V.2 Angular power spectrum

We now focus on the LIGO/Virgo frequency band. In Fig. 12 (left panel) we present the fractional difference between the power spectra of the various models described and the angular power spectrum of the reference model, defined as

[TABLE]

where $C_{\ell}^{\rm ref}$ and $C_{\ell}^{\rm mod}$ are the angular power spectra of the reference and modified model, respectively.

We can easily conclude that changing the stellar evolution model (and hence the mass distribution of the BH population) and changing the cut-off in mass with respect to the reference model gives a relative variation of the order of 50%.

The right panel of Fig. 12 shows the same quantity but for the angular power spectra normalized over the monopole (of each model). It can be concluded from these results that changing the distribution of the initial semi-major axis leads to a variation of up to 80% at low multipoles. We stress that the modified distribution considered here, namely flat in the semi-major axis, has not a strong astrophysical motivation and is taken here for illustrative purposes. The effect of varying the orbital semi-major axis is enhanced when plotting the fractional difference of relative anisotropies (i.e. normalized over the monopole). This can be understood considering that the amplitude of the green dotted curve on the right panel can be (roughly) obtained from the corresponding one in the left panel multiplying it by the ratio between the monopole of the reference model and of the model with new distribution of the semi-major axis. This gives a multiplicative factor of order $1.3$, see Fig. 11, which shifts the curve of fractional differences when going from the left to the right panel. A similar reasoning holds for the other models.

The difference between the models in the left panel of Fig. 12 can be explained by noting that the angular power spectrum is sensitive to the astrophysical kernel $\mathcal{A}$, which we plot as a function of redshift in the left panel of Fig. 13. For a fixed frequency, the monopole of the energy density is sensitive to the integral over redshift of $\mathcal{A}$ while the amplitude of the anisotropies at a given multipole $\ell$ is sensitive to the amplitude of the kernel $\mathcal{A}$ at the redshift $z$ corresponding (through Limber) to the multipole $\ell$ considered. Similarly, the amplitude of relative anisotropies (i.e. normalized over the monopole) at a given $\ell$ is sensitive to the amplitude of $\mathcal{A}$ normalized over the monopole, see right panel of Fig. 13. This also explains why the green dotted curve in Fig. 12 (representing the fractional difference between the anisotropies of the model with modified distribution of initial semi-major axis and the reference one) decreases with multipoles: the slope of $\mathcal{A}$ for this model does not increase between redshift 0.1 and 1, while it does for the reference model.

In Fig. 14 we present the angular power spectra for the models testing different black hole populations and the fractional difference between the power spectrum of the models with only high and only low black hole mass (described in section IV) and the angular power spectrum of the reference model. We observe that the fact that the sum of two angular power spectra (for low and high black hole masses) does not give back the total angular power spectrum is due to the fact that in the division between the sub-population low-high mass binaries we do not consider binary systems with one high-mass and one low-mass black hole.

V.3 Cross-correlations

Figure 15 shows the cross-correlation between anisotropies of the AGWB and weak lensing convergence, for Euclid and SKA source distributions and for a frequency $f=32$ Hz while Figure 16 shows the cross-correlation with galaxy number counts for two different redshift bins and integrated over frequency in the LIGO/Virgo band ($10$ Hz $<f<100$ Hz). For both cross-correlations, the effect of changing metallicity and distribution of semi axis gives a typical variation of order 30%.

Notice that for each model, the amplitude of the cross-correlation with lensing as a function of multipoles is a biased tracer of the amplitude of the astrophysical kernel $\mathcal{A}$ as a function of redshift, normalized over the monopole (we recall that we are defining the cross-correlations normalizing over the monopole of the AGWB, see eq. (13)). On the other side, the cross-correlation with galaxy number counts selects different bins in redshift in the astrophysical kernel $\mathcal{A}$. This can be seen in Fig. 16: as we shift the window in redshift, the peak of the cross-correlation shifts towards higher multipoles, as expected from the Limber relation between multipoles and redshifts. In a given redshift bin, the amplitude of the cross-correlation is bigger for models with the bigger value of $\mathcal{A}/\bar{\Omega}_{\rm GW}$ in the redshift bins considered. Hence, interestingly, the cross-correlation with galaxy number counts can help to reconstruct the astrophysical kernel $\mathcal{A}$ as a function of redshifts.

VI Contribution of binary neutron stars

The detection of the binary NS merger by the LIGO/Virgo network The LIGO Scientific Collaboration and the Virgo Collaboration (2017); Abbott et al. (2017) and the estimated rate of mergers in the local Universe of $R=920^{+2220}_{-790}$ Gpc*-3* yr*-1* The LIGO Scientific Collaboration and the Virgo Collaboration (2018a) led to the conclusion that these sources may have a comparable contribution to the AGWB relative to binary BHs Abbott et al. (2018); The LIGO Scientific Collaboration and the Virgo Collaboration (2019a). We may therefore expect that their contribution to the anisotropies of the AGWB will also be important.

While it will be difficult to disentangle the relative contributions of binary BHs and NSs to the overall AGWB, especially in view of the large modeling uncertainty in the binary NS merger rates Chruslinska et al. (2018); Giacobbo and Mapelli (2019), it is interesting to note that their host galaxies are expected to have different properties. In the isolated BH formation scenario discussed here, BHs masses are heavily influenced by the metallicity of their progenitor stars, as discussed above. Specifically, metal-poor stars retain most of their mass throughout their evolution and collapse to form heavier BHs. As a consequence, these BHs (that also produce stronger GW signal when they merge) form preferentially in high-redshift and/or low-mass galaxies Lamberts et al. (2016); Cao et al. (2018); Mapelli et al. (2018); Artale et al. (2019). In contrast, NSs can also form in metal-rich environments. In view of the different clustering properties of the host galaxy populations, binary BHs and binary NSs can in principle give rise to very different anisotropic components of the AGWB.

In order to estimate the contribution from binary NSs we calculate the formation rate of NSs in our astrophysical model. For simplicity we assume that NSs form from stars in the mass range $(8,11)M_{\odot}$ with a constant mass of $1.3M_{\odot}$. We follow the formalism described above to calculate the number of detectable sources in the local Universe and normalize the fraction of NSs that reside in binaries that merge within the Hubble time $\beta$ (similarly to the case of BHs) to result in $1$ detection during O1+O2 observing period.

Our results for the anisotropies from binary NS mergers are shown in Fig. 17. We present the results in the LIGO-Virgo frequency band. The results in the LISA band look exactly the same. The value of the background energy density in both the LISA and LIGO-Virgo frequency band is given by $\bar{\Omega}_{GW}(f)\sim 2.15\times 10^{-11}f^{2/3}$.

Fig. 17 shows that the contribution from binary NSs can be dominant relative to the one from binary BHs. We stress however that this result is model-dependent, and that moreover depends on the (highly uncertain) merger rate of binary NSs.

VII Lower frequencies and the LISA band

The early inspiral phase of merging stellar-mass binary BH that may be observable with LISA space-borne interferometer is another probe of the astrophysical and cosmological processes discussed in this article. The AGWB from binary BH is expected to be dominant in the LISA band Dvorkin et al. (2016) and below, and may become a source of confusion noise for some of the other types of sources. Observations with LISA will allow one to study some aspects of resolved and unresolved stellar-mass BH binaries that are difficult to observe with ground-based interferometers.

For example, at the mHz frequencies accessible to LISA, some of the binaries may not be fully circularized, and their residual eccentricities may provide an indication to their formation channel. In particular, binaries formed through dynamical processes in dense stellar cluster can have measurable eccentricities. These can be constrained for the subset of resolved merger, and in addition the distribution of eccentricities of the entire population may also affect the resulting AGWB.

LISA will also allow one to study the AGWB from other types of sources such as close white dwarf binaries (see e.g. Ref. Vecchio (2002)), which may also produce anisotropies Ungarelli and Vecchio (2001); Kudoh and Taruya (2005). Moreover, these anisotropies can potentially be used to distinguish the astrophysical source of stochastic background from the cosmological ones in the early Universe Geller et al. (2018).

In Fig. 3 of Ref. Cusin et al. (2019b) we presented the first prediction of the angular power spectrum of anisotropies in the LISA bands. We note that the relative anisotropies have the same frequency dependence as the monopole, i.e. $C_{\ell}\propto f^{4/3}$, as expected. It follows that the frequency factorization discussed in § III.4 is a solid approximation in the LISA frequency band. A plot for relative anisotropies is presented in Fig. 17, for our reference astrophysical model. The value of the background energy density in both the LISA frequency band is given by $\bar{\Omega}_{GW}(f)\sim 1.29\times 10^{-11}f^{2/3}$.

We compare the anisotropies in the LISA band produced by stellar-mass BH for the same set of astrophysical models in Fig. 18. As already explained for the LIGO band, for a given model and a given frequency, the amplitude of anisotropies for a given $\ell$ is a tracer of the amplitude of the astrophysical kernel $\mathcal{A}$ for that frequency and at a redshift $z$ related by Limber to the multipole under consideration. Similarly, the amplitude of relative anisotropies as a function of multipoles is sensitive to the amplitude of $\mathcal{A}$ normalized over the monopole, as a function of redshift. In Fig. 19 we represent the function $\mathcal{A}$ for the various astrophysical models and this function normalized over the monopole (left and right panel respectively), as a function of redshift.111111We observe that there is a (small) difference between the results obtained in the LIGO and LISA band; Figs. 12 and 18 respectively. This is due to the fact that in the LIGO band anisotropies have a non-trivial frequency dependence, see § III.4.

VIII Conclusion

This article explored the astrophysical dependencies of the anisotropies of the AGWB and analyzed the properties of its angular power spectrum. It focused on the contribution to the background coming from BH mergers in both the LIGO-Virgo and LISA frequency bands. In particular, the properties of anisotropies in the LISA band are studied here for the first time.

After a summary of the framework developed in Ref. Cusin et al. (2018b, 2017) to study AGWB anisotropies and cross-correlations, we analyzed in details the general properties of the angular power spectrum. In particular, using an analytic approach we showed how different redshift bins contribute to the various multipoles of the correlation function. We also showed that low multipoles are dominated by the contribution of large scales (i.e. small $k$ Fourier modes). It follows that a semi-analytic description of clustering, like the one used in this work, has to be preferred to a catalogue approach to describe large angular separations (which are the first one we will hopefully observationally access). We also tested the standard factorization assumption used in most of present searches (e.g. by LIGO-Virgo), showing that it breaks down in the upper-part of the LIGO spectrum. A general analytic derivation of the shot noise contribution in different frequency bands has also been presented. In particular, we showed that the spatial component of shot noise, due to the discreteness of the galaxy distribution, is present in both LISA and LIGO frequency bands and it dominates the spectrum at high multipoles. Using a toy model, we have also illustrated the possibility to use cross-correlation with galaxy number counts to extract AGWB anisotropies from a shot-noise dominated map. A future work will be dedicated to studying the possibility of filtering out shot noise using cross-correlations and masking.

In the second part of the article, we explored the signatures of 9 different astrophysical models, differing from a reference model in one physical ingredient (e.g. distribution function of orbital parameters, cut-off in the BH mass distribution, stellar evolution model for BH generation, IMF, metallicity…). First, from the requirement that the models are all compatible with the LIGO O1+O2 detection rate, we concluded that the fraction of BH in binaries varies from 0.06% to 17% across the models. This highlights the actual degeneracy of the models. We then demonstrated that the anisotropies are very sensitive to some astrophysical parameters such as the distribution of the initial semi-major axis of the binary systems, the choice of the metallicity profile and of the astrophysical model to describe BH generation out of stellar progenitors. More quantitatively, we concluded that changing the distribution of the initial semi axis $a_{\rm f}$ of binary systems from a distribution $\propto 1/a_{\rm f}$ to a flat one induces relative variation in anisotropies (normalized over the monopole) of order 80% at small multipoles ($\ell<10$). Changing the stellar evolution model or keeping a low constant metallicity leads to a variation of order 40%. These results hold for both LIGO and LISA frequency bands. This result contradicts the claim of Ref. Jenkins et al. (2018b), where the authors conclude that the anisotropies do not depend on astrophysical properties of the BH populations.

The present analysis concludes that the prediction of the angular power spectrum is very sensitive to the astrophysical modeling: astrophysical models differing for the description of stellar evolution and for the choice of the distribution of the orbital semi-major axis give predictions differing of $>40$% in both the LIGO and LISA frequency bands. For those models, we also compared the results for the cross-correlation with weak lensing convergence (for both SKA and Euclid source distributions) and galaxy number counts, at different redshifts. At low multipoles, variations are up to order 30%. Results for the angular power spectrum in the LISA band are shown to scale with frequency as the monopole. i.e. $C_{\ell}\propto f^{4/3}$. The spectrum of relative anisotropies in the LISA band is therefore frequency independent.

This result has interesting astrophysical implications. The AGWB angular power spectrum is an observable we will hopefully have access to observationally in the next few years. From the comparison between future observed intensity maps of the background (for different frequencies) and our theoretical predictions, we will be in the position to set constrains on astrophysical quantities, such as stellar evolution model, metallicity distribution and distribution of orbital parameters, that cannot be accessed otherwise. This is the first stage of a long term research program. We analyzed which ingredients of the astrophysical parametrization give sizable differences in the anisotropies and which do not play a role. It will allow us to refine the astrophysical modeling. We also provide a theoretical understanding of relation between the shape of the angular power spectrum and the shape of the astrophysical kernel.

We emphasize that the monopole and the anisotropies of the AGWB give access to complementary information. While the former is sensitive to the integral over redshift of the astrophysical kernel describing sub-galactic physics (the function $\mathcal{A}$ in this work), anisotropies at different angular separations are sensitive to the amplitude of this kernel at different redshifts. Anisotropies open a new window on the reconstruction of the redshift dependence of the astrophysical kernel: this new observable will help us to test properties of source populations at different redshifts and eventually reconstruct their evolution in time.

Acknowledgements.

It is a pleasure to acknowledge D. Alonso for a careful reading of this work and very interesting discussions on shot noise. We thank G. Lavaux, F. Lacasa, B. Whiting, I. Michaloliakos, P. Ferreira, Ruth Durrer and E. Vangioni for stimulating discussions. This work was done within the Labex ILP (reference ANR-10-LABX-63), part of the Idex SUPER, and received financial state aid managed by the ANR, as part of the programme Investissements d’avenir under the reference ANR-11-IDEX-0004-02. GC acknowledges financial support from ERC Grant No: 693024 and Beecroft Trust.

Appendix A Computation of shot noise

We assume a binned galaxy survey. The quantity $\Delta_{i}$ describes the galaxy overdensity in the bin $i$. Then the quantity $\langle\Delta_{i}\Delta_{j}\rangle$ contains the theoretical correlation function and a shot noise contribution

[TABLE]

where $V_{p}\bar{n}_{G}$ is the mean number of galaxies per pixel and $\delta_{ij}$ is the Kronecker symbol. We can take the continuous limit of a binned survey as

[TABLE]

where a hat denotes observed quantities and $V_{p}$ is the pixel volume. Using Eq. (19), we can rewrite the measured background anisotropies as

[TABLE]

or equivalently, using that $\partial_{z}\bar{\Omega}_{\rm GW}(z,f)=H\partial_{r}\bar{\Omega}_{\rm GW}(r,f)$

[TABLE]

We use this result when computing the correlation function of the background anisotropies

[TABLE]

The first part in the square brackets is the theoretical result for the correlation function while the second term is the Poisson noise contribution. Writing

[TABLE]

after standard manipulations one finds

[TABLE]

where the first contribution is the theoretical one, the second the shot noise component. Using the standard definition

[TABLE]

we immediately find for the angular power spectrum

[TABLE]

where

[TABLE]

We see that this integral diverges in $r=0$, hence the contribution of Poisson noise depends on the cut-off used to regularize this integral. From an observational point of view, the physical quantity on which the cut-off has to be set is the observed flux: sources with a flux bigger than a given threshold can be resolved and are therefore filtered out. To see where in Eq. (65) the cut-off in flux appears, we use Eq. (4), i.e.

[TABLE]

and we rewrite the astrophysical kernel Eq. (5) as an integral over emitted luminosity as

[TABLE]

where $\bar{n}_{{}_{\rm G}}(\mathcal{L}_{{}_{\rm G}},r)$ is the average number of galaxies at distance $r$ with luminosity $\mathcal{L}_{{}_{\rm G}}$. We replace Eq. (67) in Eq. (65) and change the variable of integration to redshift. Using that the flux $\Phi$ received per units of frequency from a source in $z$ is related to the luminosity per units of emitted frequency by $\Phi(f)=1/(4\pi)\mathcal{L}_{{}_{\rm G}}/(1+z)$, we see that an upper bound on $\Phi$ is translated into a lower bound in redshift and upper bound in luminosity (more precisely, it defines the region of integration in the plane $(z,\mathcal{L}_{{}_{\rm G}})$. In other terms, we introduce a selection function which is 1 for $\Phi<\Phi_{\text{cut}}$, and 0 otherwise. Explicitly in Eq. (58) we introduce the selection function

[TABLE]

Of course, if one can assume that all galaxies have associated the same luminosity, then the cut-off on flux translates directly into a lower cut-off in redshift (or analogously in $r$). Following the steps described in this section, we find that in Eq. (64) both the theoretical prediction for the angular power spectrum and the Poisson noise involve an integration with the selection function. The selection function is regularizing the integral defining the shot noise part of the result, see Fig.9 for an illustration. Choosing a cut-off in flux correspond to filtering out resolvable sources: for a given detector, one can compute the signal to noise ratio as a function of flux and define $\Phi_{\rm cut}$ as the flux corresponding to a signal to noise bigger than a five threshold. For example choosing 5 as such a threshold one can find the cut in flux from

[TABLE]

The fact that the level of Poisson noise depends on the choice of cut-off in flux, is not surprising and the same result holds for the cosmic infrared background case, see e.g. Ref. Ade et al. (2014).