Constraining extended cosmologies with GW×LSS cross-correlations

The rapid development of gravitational wave astronomy provides the unique opportunity of exploring the dynamics of the Universe using clustering properties of coalescing binary black hole mergers. Gravitational wave data, along with information coming from future galaxy surveys, have the potential of shedding light about many open questions in Cosmology, including those regarding the nature of dark matter and dark energy. In this work we explore which combination of gravitational wave and galaxy survey datasets are able to provide the best constraints both on modified gravity theories and on the nature of the very same binary black hole events. In particular, by using the public Boltzmann code Multi_CLASS, we compare cosmological constraints on popular ΛCDM extensions coming from gravitational waves alone and in conjunction with either deep and localized or wide and shallow galaxy surveys. We show that constraints on extensions of General Relativity will be at the same level of existing limits from gravitational waves alone or one order of magnitude better when galaxy surveys are included. Furthermore, cross-correlating both kind of galaxy survey with gravitational waves datasets will allow to confidently rule in or out primordial black holes as dark matter candidate in the majority of the allowed parameter space.


Introduction
The ΛCDM, i.e., the cosmological Standard Model, passed numerous tests and provided a convincing explanation for the physics behind many cosmological observables, both from the early [1] and late Universe [2][3][4].However, some tensions are still present [5], and a deeper investigation might in fact reveal hints of so-called new physics.
One potential avenue to investigate these tensions is to include in the analysis new and independent datasets.In this respect, gravitational waves (GWs) present the unprecedented opportunity to explore the properties and evolution of the Universe in a completely new fashion, complementary to that of more traditional cosmological probes.Thanks to the extraordinary effort of the GW community, we have already detected almost a hundred GW events [6][7][8], and the number is expected to rapidly grow for every future LIGO-Virgo-KAGRA Collaboration run.Moreover, once next generation GW observatories like Einstein Telescope [9] (ET) and Cosmic Explorer [10] (CE) will start operating during the next decade, the number of detected events per year will be of order O(10 4 −10 5 ), allowing for the unprecedented opportunity to map the entire sky using GWs.
Similarly to electromagnetic radiation, GWs carry information about the sources generating them, the environment where these sources reside and about the Universe the gravitational radiation travelled across.In this work we show how combining GW and Large Scale Structure (LSS) datasets we can extract both astrophysical and cosmological information, shedding light both on the nature of dark matter and gravity.
In recent years, Primordial Black Holes (PBHs) emerged as a possible candidate to make up a non-negligible fraction of the dark matter.PBHs are (still hypothetical) black holes that -in the most standard scenario -formed in the early Universe before big bang nucleosynthesis, deep in the radiation-dominated era [11][12][13][14], and since they are not generated through a stellar process, they can span a wide range of masses.In particular, if PBH masses are larger than M PBH ≳ 10 −18 M ⊙ , their present day abundance can account for a fraction or even the totality of dark matter [15,16].Moreover, similarly to astrophysical BHs, PBHs can form binaries and be detected by present and future groundbased GW observatories if they have masses of order O(1 − 10 2 ) M ⊙ .However, PBH binaries (and their relative GW signal) are expected to trace LSS differently from binaries of astrophysical BHs, providing a potential way to infer the nature of the binary constituents [17][18][19].
In this work we study if galaxy surveys, deep and localized (DL), or shallow and wide (SW), are suited for extracting cosmological information in combination with GW datasets.On the GW side, we consider binaries made of either astrophysical or primordial BHs detected by second and third generation detectors.The GW events are modeled using the external modules of CLASS GWB [51], while the galaxy and GW clustering signal is computed by extending the Multi CLASS [52,53] Boltzmann code 1 to account for modified gravity models.
We show that both kinds of survey we consider will be effective in constraining different extensions of the ΛCDM model.Regarding the modified gravity/dark energy sector, we show that common extensions of General Relativity like scalar-tensor theories or nDGP models can be constrained by GW alone using only linear scales, obtaining limits comparable to current ones, in a completely independent way and without relying on non-linearity modeling.When cross-correlating with galaxy surveys, constraints improve by almost an order of magnitude, in similar ways for the dark energy models we considered.On the other hand, when investigating potential dark matter candidates, we find that the cross-correlation of galaxies and GW can definitively rule in or out with high degree of confidence solar mass PBHs as the main component of dark matter.
This article is structured as follows: in section 2 we present our formalism and the statistical tools we use to forecast the cosmological constraints.In section 3 we characterize the tracer we are considering, i.e., galaxies and gravitational waves.In section 4 we quantify which kind of survey gives the best constraints.Finally we conclude in section 5. Appendices provide additional information regarding the number count fluctuation in General Relativity (A), modifications implemented in Multi CLASS (B), second and third generation GW detector network specifications (C), PBH binary formation channels (D), optimistic versions of Modified Gravity parameters (E) and astrophysical-primordial BHs second generation detected events (F).

Formalism and methodology
In this section we briefly summarize the most important aspects of the galaxy and gravitational wave clustering, along with the basic tools we use to assess the constraining power of future galaxy surveys and GW observatory datasets.The interested reader can find a broader and richer discussion of these topics in refs.[52,54].

Clustering statistics
Anisotropies of the Universe LSS contain a wealth of information concerning the origin and the growth of cosmological perturbations.In this work we are interested in the fluctuation of the number density of objects ∆ X (z, n) at redshift z in the direction n [55][56][57], where X labels the tracer, i.e., the objects we observe (galaxies or GWs).Given the spherical symmetry of the sky, it is convenient to expand the number density fluctuation in spherical harmonics as where the a X,z ℓm are the spherical harmonic coefficients and Y ℓm are the spherical harmonics.The two-point statistics in harmonic space, i.e., the angular power spectrum C ℓ , is given by [58,59] a X,zi ℓm a where * indicate the complex conjugate and δ K is the Kronecker delta.The angular power spectrum is defined as [60,61] where P(k) = k 3 P (k)/2π 2 is the almost scale-invariant primordial power spectrum.The number count fluctuation is where dN X /dz is the redshift distribution of the tracer X and W (z, z i , ∆z i ) is a window function (normalized to unity) centered at redshift z i with half-width ∆z.The number density fluctuation transfer functions ∆ X ℓ (k, z) is given by the sum of different contribution commonly called density, velocity, lensing and gravity effects [55][56][57]: and the explicit form of these terms is reported in appendix A. The angular power spectra are computed using Multi CLASS [52,53], a public extension of the Boltzmann code CLASS [60] that allows to compute the angular power spectra for different combinations of tracers including all the projection effects included in equation (2.5).Any tracer population is characterized by four functions of redshift and, possibly, scale: the number density redshift distribution d 2 N X /dzdΩ; the clustering bias b X , which connects the tracer overdensity to the matter overdensity by δ X = b X δ m ; the magnification bias [62][63][64][65] that describes how a flux-limited observable is affected by cosmic lensing effects; and the evolution bias [56,57,66] which describes the formation rate of new tracers.Moreover, for the purpose of this analysis, we introduce a new parameter called effective bias where z min and z max are the minimum and maximum redshift of the tracer X survey, respectively.
Finally, real measurements are affected by multiple sources of noise, which we describe via the noise angular power spectrum N XY ℓ (z i , z j ).In other words, the observed angular power spectrum Cℓ has the form CXY The first source of noise we consider is the shot noise due to the discrete nature of the tracer, both galaxies and GWs, which is characterized by a white noise of the form Errors due to poor galaxy localisation, both in terms of angular position and rdshift determination, are expected to be subdominant in our analysis.However, contrarily to galaxies, GWs suffer from a limited spatial resolution due to instrumental noise.We take an approach similar to that used in weak lensing analysis, see e.g., refs.[67][68][69], to treat the error due to poor determination of the GW redshift.First we define a probability p obs for a GW emitted at redshift z GW to be observed at redshift z obs due to instrumental noise, where, for simplicity, we assume (2.10) In the case of an infinitely precise instrument the dispersion σ z tends to zero and the probability p obs becomes a Dirac delta centered at the true value z GW , as we naively expect to be the case.Then we introduce an observed angular number of objects dN obs GW /dΩ and we use this quantity in the shot noise angular power spectrum instead of the true angular number of objects, as done in equation (2.9).More explicitly, for each redshift bin we have where, also in this case, in the limit of no instrumental error σ z → 0 the function in square brackets becomes 2Θ H (z GW − z min )Θ H (z max − z GW ) and we recover the expected form of the shot noise.On the other hand, we describe the effect of poor angular resolution via a Gaussian beam enhancing factor of the noise angular power spectrum, as suggested in ref. [70].Therefore the total GW angular power spectrum Cℓ has the form where θ avg res is the average resolution of a GW detector network and "log" refers to the natural logarithm.The maximum multipole ℓ max = 180 • /θ max res we include in our analysis is determined by the maximum resolution of a GW detector network θ max res .

Fisher matrix analysis
Despite its simplicity, a Fisher matrix analysis is able to quickly forecast the constraining power of future experiments.In particular, by estimating the curvature of the log-likelihood around its maximum, the Fisher analysis returns the best error an experiment can achieve, i.e., the Cramér-Rao bound [71][72][73][74].Given a pair of cosmological parameters {θ α , θ β }, their respective Fisher matrix element reads as where to compute the RHS we assume a Gaussian likelihood L for the spherical harmonic coefficients, see, e.g., ref. [52].The parameter f sky represents the observed fraction of sky in common between the different (galaxy or GW) surveys considered in the forecast.Assuming the tracers are binned into n redshift bins z 1 , ..., z n , the (symmetric) covariance matrix C ℓ reads as The matrix ∂ α C ℓ ≡ ∂C ℓ /∂θ α contains derivatives of the elements of the covariance matrix with respect a single cosmological parameter keeping all the others fixed.
Errors on cosmological parameters are obtained from the error covariance matrix Σ = F −1 , in particular the marginalized errors are given by σ θα = √ Σ αα .Moreover, using the marginalized errors we can perform a sort of null hypothesis testing comparing two models, the "fiducial" and the "alternative", each one characterized by the value of a parameter, in our case the GW effective bias.By computing a Signal-to-Noise ratio (SNR) of the form we can assess whether the alternative model, labelled by "alt", is statistically different from the fiducial one, labelled by "fid".For values of SNR ≲ 1 we consider the two models as statistically indistinguishable.

Cosmological model
In this work we investigate both the nature of dark matter and that of dark energy using GWs, however, given the peculiarity of each scenario, we have to assume different cosmological fiducial models for the two scenarios.In the case where we study the origin of coalescing BHs, we assume a standard ΛCDM model, in which the gravitational section is described by General Relativity and a cosmological constant is responsible for the accelerated expansion.On the other hand, in the second case, where we investigate the nature of dark energy, we allow for deviations from General Relativity.Many different theories have been proposed to extend General Relativity by introducing some dynamical mechanism that explains the late time acceleration [37,38], however in this article we focus only on two classes, scalar-tensor and extra-dimension theories.
The most common class of Modified Gravity models are scalar-tensor theories where a scalar field dominates the dynamics of the Universe at late times, in particular we highlight Horndeski and beyond-Horndeski models [75][76][77][78][79].Despite having been constrained by early and late time cosmological probes, see, e.g., ref. [80], and by the remarkable binary neutron star event [81][82][83], the parameter space in this class of models can still allow for interesting deviations from General Relativity at late times.Following the spirit of the Effective Field Theory of Dark Energy [84][85][86][87][88][89][90], we choose to conveniently parametrize deviations from General Relativity instead of starting from a well defined Lagrangian density.In particular, at background level we impose a ΛCDM expansion history, while at perturbative level we allow for deviations by introducing two free functions of redshift z and scale k. Between many different parametrization choices [35], we choose the µ − η one given by [91] where the Bardeen potentials Φ, Ψ are defined by the line element the conformal time and spatial coordinates are τ and x j , respectively, a is the scale factor, D = δ +3H(1+w)θ/k 2 is the gauge-invariant energy density fluctuation, δ is the energy density fluctuation, θ is the velocity divergence, H = a ′ /a is the Hubble expansion rate in conformal time, w is the equation of state of the Universe, G is Newton constant.Regarding the µ − η Modified Gravity functions, the former affects the growth of structure by effectively changing the strength of Newton constant, while the latter changes the relation between the metric potentials.
Regarding the form of µ and η, several authors have considered a scale-independent form for these functions; however, in total generality, Modified Gravity theories often have additional scales appearing both at linear [90] and nonlinear [92][93][94][95][96][97] level.We build on the redshift-dependent parametrization introduced in ref. [98] and we complement it with an additional scale dependence to account phenomenologically for the presence of an additional scale k mg in the theory as [99,100] where µ 0 , η 0 , z th , ∆z th and k mg are constant parameters.The parameters z th and ∆z th control when and how fast we transition from the General Relativity to the Modified Gravity regime.The General Relativity limit is recovered for µ 0 = η 0 = 1, hence when µ = η = 1.Moreover, this parametrization allows to recover the ΛCDM limit in the Early Universe (z ≫ z th ) or at small scales (k ≫ k mg ) independently from the value of µ 0 and η 0 , while, at the same time, allowing for deviations at late time and large scales.Following refs.[30,101] we choose z th = 6, ∆z th = 0.05.We choose as benchmark scale k mg = 10 −3 , 10 −2 and 10 −1 Mpc −1 to study the dependence of the constraints on the extra scale of the theory.The different parametrizations have been implemented in Multi CLASS, and the code has been validated against the results of refs.[102,103] obtained using the MGCLASS II code (more details are included in appendix B).
In alternative, other popular Modified Gravity models involve the existence of extra dimensions, as the Dvali-Gabadadze-Porrati (DGP) model [104]: this 1-parameter extension of ΛCDM describes the Universe as a 4D brane embedded in a 5D Minkowski space.While matter is confined on the brane, gravity propagates in the extra dimension above the scale r c , which represents the cross-over scale between the 4D and 5D behaviour.In this work we focus on only one of the two branches of the theory [105], the so-called "normal branch" (hereafter denoted as nDGP), since the "selfaccelerating branch" has been showed to be unstable [106][107][108].This simple model has an incredibly rich phenomenology, allowing for potential deviations from GR at the background, linear and nonlinear level [109][110][111][112][113][114].However, in this work we focus on the model proposed in ref. [115], which leaves unaltered the expansion history of the Universe but changes the growth of structure.Despite having a substantially different theoretical framework from the case of scalar-tensor theories, also this model can be described at perturbative level using the parametrization given in equation (2.16) by imposing [116,117] where and H = ȧ/a and " ˙" represents a derivative with respect to cosmic time.In this case we do not include any additional scale dependence in the Modified Gravity functions.
In conclusion, in the case where we investigate dark matter properties, we consider as cosmological parameters for the forecast where ω cdm is the cold dark matter physical density, h is the Hubble parameter, n s is the spectral index of scalar perturbations and b eff g , b eff GW are the effective bias parameters for galaxies and GWs.The fiducial value of GW effective bias depends on binary formation channels, as explained in section 3. On the other hand, in the case where we forecast sensitivity for dark energy models, we consider and for the Modified Gravity and nDGP cases, respectively.In the latter two cases, the GW effective bias refers to the case of astrophysical BBHs. 2 The fiducial value of the standard cosmological parameters are {ω cdm , h, n s } = {0.12038,0.67556, 0.9619}, while for the Modified Gravity models we use {µ 0 , η 0 } = {0.87,1.3} and for nDGP we choose {Ω rc } = {0.2}.

Tracers
In this section we present the specifics of the two kinds of galaxy survey (SW or DL, § 3.1) and the details on the different astrophysical ( § 3.2) or primordial ( § 3.3, 3.4) BH binaries.Regarding PBHs, we consider binaries that form either at late or early times, since they trace LSS in well defined ways.
Including other formation channels in specific environment, see, e.g., ref. [118], or formation channels that include also neutron stars [119] is left for future work.Finally, in § 3.5 we show how to combine GW events coming from different binary models in a single framework.Binary models have been studied using the external modules of CLASS GWB [51], and the detectability of these sources has been done both for second and third generation GW detector networks, as detailed in appendix C. In particular, we consider as detector networks: • HLVIK, a second generation network made by Hanford and Livingston LIGO, Virgo, KAGRA and LIGO India L-shaped interferometers; • ET2CE, a third generation network made by Einstein Telescope (triangle-shaped interferometer) and two Cosmic Explorer (L-shaped interferometers).
We assume all the detector to work at their design sensitivity, and individual events are considered "detected" if their Signal-to-Noise ratio for the network is ρ ≥ 12, see, e.g., ref. [51] and references therein.Errors on the redshift determination are of the order of δz/z ≃ δd L /d L ≃ 3/ρ ≲ 0.3 [70,120] and are discussed further in appendix C along with typical average and maximum resolution of the two detector networks.For any effective purpose, both second and third generation detector networks cover the entire sky (even if with different sensitivities at different times), therefore for GWs we assume f sky = 1.0.In the next subsections we provide fitting formulas to allow interested readers to reproduce our results, without claiming that they provide any particular insight on the properties of the tracers. 3

Galaxy surveys
First we consider a DL survey that observes up to 326 × 10 6 galaxies over 8000 deg 2 (f sky = 0.2), covering the [0.5, 4.1] redshift range.For the sake of this analysis we assume data is grouped into 18 redshift bins of half-width ∆z = 0.1.On the other hand, for the SW survey we consider a mission that covers the entire sky (f sky = 0.7) in the redshift range [0.2, 4.2].In this case we bin galaxies into 10 redshift bins with half-width ∆z = 0.2.These specifics are inspired by those of SPHEREx [43][44][45]. 4n both cases we consider tophat window functions when binning the galaxy catalogs.
The redshift distributions of the galaxy populations for both surveys can be parametrized as where for the DL survey A = 194575 gal/deg 2 , z 0 = 0.07, α = 0.34 and β = 0.42; while for SW one we consider A = 25509 gal/deg 2 , z 0 = 0.09, α = 1.75 and β = 0.69, which take into account the contribution of each sample.
The bias for the DL survey is approximated as b g = 0.84/D(z), where D(z) is the linear growth factor normalized to 1 at redshift z = 0, while for SW one we consider a parametric form where b 0 = 0.53, b 1 = 1.59 and b 2 = −0.08.Therefore, the effective galaxy bias for these two surveys reads as b eff g = 1.76 and b eff g = 1.69 for DL and SW, respectively.Regarding the magnification bias s g , since this is a preliminary study and no accurate data is available, we choose s g = 0.6 for both surveys.Finally, the evolution bias f evo g is inferred from the galaxy redshift distribution given in equation (3.1).

Astrophysical black hole binaries
First of all we consider GWs emitted by coalescing BHs generated as a result of stellar evolution.These BHs are expected to live in galaxies, where the star formation history is more intense, hence they are expected to trace LSS in a similar (yet not identical) way to what galaxies do.We create catalogs 5 of GW events following the model and the procedure presented in ref. [51], hence, for the sake of conciseness, we will not repeat all the details in this article.Suffice to say that every GW event is described by its time of arrival at Earth, merger redshift, binary formation redshift, dark matter halo mass both at formation and at merger, masses and spins of the compact objects, sky localization of the event, inclination angle of the binary and polarization of the signal.Assuming that all the BH binaries are astrophysical in origin, the LIGO-Virgo-KAGRA Collaboration measured a local merger rate of RLVK 0 ≃ 20 Gpc −3 yr −1 [121].Given the catalogs, we compute the observed GW redshift distribution both for second and third generation detector networks and we normalized it to T obs = 10 years of observation time.We parametrize the GW redshift distribution as in equation (3.1), finding {A 2G , z 2G 0 , α 2G , β 2G } = {84.26GW/sr, 0.04, 2.79, 0.68} and {A 3G , z 3G 0 , α 3G , β 3G } = {1.35×10−7 GW/sr, 2.02×10 −3 , 6.12, 0.41} for second and third generation detectors, respectively.The bias of each individual GW event is computed starting from the dark matter halo mass at merger and by using the halo bias model, as proposed in refs.[51,54].The biases are then averaged in every redshift bin and the overall bias function is parametrized as where we confirm the result first presented in ref. [18] that s 2G GW ≃ [0.0, 0.6] and s 3G GW ≃ 0 for second and third generation of detectors, respectively (see also figures 7 and 1).The GW evolution bias, as for galaxies, is computed from the GW redshift distribution.

Late time primordial black hole binaries
PBHs inside virialized dark matter halos can form binaries at late times through a direct capture process of unbound compact objects [122,123], similarly to standard astrophysical BHs.In this article we follow the approach of ref. [122]  ≃ 2 Gpc −3 yr −1 to be compared to the (overestimated) local merger rate RLPBH m,0 ≃ 9.5 Gpc −3 yr −1 obtained without including the time delay effect.As for astrophysical BHs in the previous section, we use the catalog approach to derive the late time PBH merging binary redshift distribution and bias.Parametrizing again the GW redshift distribution with equation (3.1), we find {A 2G , z 2G 0 , α 2G , β 2G } = {119.29GW/sr, 0.20, 2.40, 1.05} and {A 3G , z 3G 0 , α 3G , β 3G } = {0.2GW/sr, 5 × 10 −3 , 2.42, 0.33} for second and third generation detectors, respectively.For the bias we follow the same methodology used for astrophysical BHs and, parametrizing the bias as in equation ( 3 for second and third generation detector networks, respectively.The effective GW bias is given by b 2G,eff GW = 0.59 and b 3G,eff GW = 0.61 in both the redshift ranges probed by the DL and SW surveys.Also in this case, we find that s 2G GW ≃ [0.0, 0.6] and s 3G GW ≃ 0 and the GW evolution bias is calculated from the GW redshift distribution.

Early time primordial black hole binaries
Alternatively, PBH binaries can form before matter-radiation equality [124][125][126].Thanks to tidal forces generated by other PBHs, the binary components avoid an head-to-head collision and start their inspiralling motion, eventually merging at later times.Numerous authors analysed both the binary creation stage and the binary evolution history, finding that overall the merger rate of such PBH binaries can be modelled as [127][128][129] REPBH where t 0 = t(0) is the age of the Universe, and A m is an amplitude factor that accounts for effects connected to the distribution of initial binaries [129][130][131], cosmological perturbations [132], interaction with surrounding environment [133], and many more.Typical values of the amplitude factor are of order O(10 5 ) Gpc −3 yr −1 , which, if taken at face value, would constraint PBH abundance to be of order f PBH ∼ O(10 −3 ) [134][135][136].However, numerical simulations seem to suggest that binary disruption due to three body encounters might be more frequent than what was previously thought [137,138], and that the real value of the amplitude factor is of order A m ∼ O(10) Gpc −3 yr −1 in the tens of solar masses range, making PBH a viable candidate for dark matter.In the following we fix as benchmark value A m = 18 Gpc −3 yr −1 : in this fashion, when f PBH = 1, we have that late and early time PBH binaries completely saturate the measured local merger rate RLVK 0 .Also for early time PBH binaries we create catalogs of events and assess how many are detected by different GW detector networks.Parametrizing the GW redshift distribution according to equation (3.1), we find {A 2G , z 2G 0 , α 2G , β 2G } = {133.75GW/sr, 0.09, 3.41, 0.84} and {A 3G , z 3G 0 , α 3G , β 3G } = {9.40GW/sr, 0.03, 2.85, 0.41} for second and third generation detectors, respectively.Since early time PBH binaries formed during the radiation dominated era and they follow the subsequent dark matter evolution, we assume they trace very well the underlying matter distribution, hence we assume their bias to be equal to b GW = 1, and therefore also the effective bias is b eff GW = 1.As for astrophysical BH and late time PBH binaries, we find s 2G GW ≃ [0.0, 0.6] and s 3G GW ≃ 0 for second and third generation detectors and we calculate the GW evolution bias from the GW redshift distribution.

Complete black hole binary population
So far we consider cases where only a single BH binary population contributes to observations, however in reality all these populations contribute to the totality of detected events.Hence the true functions characterizing the two-point statistics of GW are the total number density and, by definition, total bias, magnification bias and evolution bias The relative importance of PBH with respect to astrophysical BHs depends on the f PBH and A m parameters.We investigate the PBH parameter space given by A m ∈ [0, 18] and f PBH ∈ [0, 1], where, at every point, the local merger rate of astrophysical BHs is defined by in order to saturate by definition the measured local merger rate.We show in figure 1 the total redshift distribution, bias, magnification bias and evolution bias for different models spanning the f PBH − A m parameter space for a third generation detector network, along with three reference scenarios corresponding to the astrophysical BHs only (f PBH = A m = 0), mixed astrophysical and late time PBH (f PBH = 1, A m = 0), and PBHs only (f PBH = 1, A m = 18) cases.In the astrophysical BH only scenario we observe a redshift distribution peaked z ≃ 1.5, when also the cosmic star formation rate peaks.Increasing the abundance of PBHs slightly alters the distribution, if only late time PBH binaries are included, or significantly shifts it at higher redshift of order z ≃ 3.5, when early time PBH binaries are included.At the same time, the inclusion of PBH significantly lowers the total bias of the population, especially at high redshift, while affecting only partially the total magnification and evolution bias.The analogue figure for GW events detected by a second generation detector network is reported in appendix F.
Finally, we note that, even keeping fixed the total local merger rate to the measured value, the total number of detected GW events changes due to the different merger rate redshift dependence of different models.Since the average number of detected events is N ABH ≃ 33000 GW yr −1 when f PBH = A m = 0 (astrophysical BHs only scenario) or N PBH ≃ 43000 GW yr −1 when f PBH = 1, A m = 18 (PBH only scenario, 41000 early time and 2000 late time events), over the span of T obs = 10 yr we would observe approximately between 2.9 × 10 5 and 4.3 × 10 5 GW events, depending on how much each population contributes to the total.These estimates are considerably influenced by the shape of the merger rate at high redshift, which is currently unknown, however we can already use them to infer how much the shot noise changes in our forecast.In particular we can see that it would change at most by a factor 1.5 over the entire f PBH − A m parameter space.The astrophysical BH model we adopt can be seen as conservative compared to existing estimates that predict up to N ABH ≃ 7 − 8 × 10 4 observed events per year [139], hence also our errors and detectability prospects can be interpreted as conservative.In this plot we consider only events detected by a third generation GW detector network.The light blue, blue and dark blue lines correspond to the astrophysical BHs only (fPBH = Am = 0), mixed astrophysical and late time PBH (fPBH = 1, Am = 0), and PBHs only (fPBH = 1, Am = 18) scenarios, respectively.In the mixed astrophysical and late time PBH scenario, late time PBHs contribute to approximately 6% of the total number of observed events.The green lines refer to mixed models where different populations contribute to the total number of detected events.

Testing extended cosmologies
In this section we present the capability of GWs to test cosmological model that deviates from the standard ΛCDM either in the dark matter or in the dark energy sectors.We present results for GW alone and in synergy with galaxy surveys, in particular we discuss whether a DL or a SW survey are more effective in constraining our set of extended cosmologies.In order to draw fair conclusions regarding both kind of surveys, first we report in table 1 errors on standard cosmological parameters: our results suggest that the DL survey is equal to or even a factor two more efficient than the SW survey in constraining them.(orange, red and dark red lines, respectively).We show the effect only for galaxies, however also the GW angular power spectrum and the galaxy-GW angular power spectra share a similar phenomenology.

Constraints on modified gravity theories
First of all, we consider the effect of the additional Modified Gravity scale dependence in the µ(k, z) and η(k, z) functions parametrized by k mg .We present in figure 2 the relative difference between the angular power spectrum in two different redshift bins in Modified Gravity models with k mg = 10 −1 , 10 −2 and 10 −3 Mpc −1 with respect to the ΛCDM case.We also show in the same figure the case in which no extra scale dependence is included as a useful comparison obtained: this scenario is obtained by taking the limit k mg → ∞, in which µ and η become functions only of redshift.The suppression effect showed in figure 2 can be understood by comparing the Modified Gravity scale k mg to the typical scale associated to each multipole k ≃ (ℓ + 1/2)/r(z) [140,141], where r(z) is the comoving distance to the redshift bin centered at z.In the cases showed in the figure we have r(0.4)≃ 1.6 Gpc and r(3.6) ≃ 7 Gpc, therefore at multipoles such that ℓ ≫ ℓ mg (z) = k mg r(z) − 1/2 we expect deviations from ΛCDM to be exponentially suppressed and to recover the ΛCDM limit.In this example the critical multipole at which the two scales are approximately equal is ℓ mg (0.4) ≃ 160, 16, 1 (ℓ mg (3.6) ≃ 700, 70, 7) for k mg = 10 −1 , 10 −2 and 10 −3 Mpc −1 , respectively, as it could have been estimated by eye in the figures.
In other words, Modified Gravity effects are highly suppressed for small values of k mg , therefore in those scenarios we expect to find weaker constraints on Modified Gravity parameters, as we show in table 2 both for GW alone and by cross-correlating them with galaxy surveys for our conservative astrophysical BH model.The numbers outside parenthesis refer to the errors obtained by the full Fisher matrix, while numbers in parenthesis refer to the case where standard cosmological models are perfectly known, and we quote them to show how tighten the constraints can become when a strong prior on those parameters is imposed. 6In particular, the difference in errors between the k mg = 10 −3 Mpc −1 and the k mg = 10 −2 Mpc −1 scenarios is about one order of magnitude, both for µ 0 and η 0 , while the difference between the k mg = 10 −2 Mpc −1 and k mg = 10 −1 Mpc −1 cases is a factor few, both with and without galaxy surveys.Note that for large values of k mg GWs alone can provide Tracers  independent and complementary constraints that are comparable to existing ones [1,35].Moreover, we find that errors on the standard cosmological parameters of table 1 are almost unaffected by the exact choice of k mg .In general, we conclude that SW and DL surveys will be equally competitive in constraining Modified Gravity parameters when cross-correlated to GWs, despite the latter being intrinsically better at constraining standard cosmological parameters, and they will provide constraints that are around one order of magnitude tighter than existing ones.
Moving to the nDGP model, we report errors on the Ω rc parameter also in table 2. Also in this case we notice that both galaxy surveys have the same constraining power in terms of ruling out the nDGP model.Moreover, we notice that errors in this case are comparable to existing ones, however in this case the constraints have been obtained considering exclusively linear scales.In other words, our constraints can be considered more robust in the sense that does not depend on including the effect of non-linearities in the nDGP model.
We also provide in table 5 of appendix E the constraints on Modified Gravity and nDGP parameters for the optimistic case corresponding to observe twice as many GW events per year, as some estimates suggest.While improving by a factor two the estimated errors on the Modified Gravity and nDGP parameters in the GW alone case, the optimistic astrophysical BH model shows only a marginal improvement in the magnitude of the errors when galaxy surveys are included.The general considerations regarding the importance of the Modified Gravity additional scale dependence we made for the conservative case apply also for the optimistic one.

Constraints on black hole origin
For the purpose of estimating the SNR in equation (2.15), we consider the astrophysical BH population as fiducial model, and the astrophysical-primordial BH population described by the parameters f PBH − A m as alternative model.We report the marginalized errors on the GW effective bias both for second and third generation detector networks in table 3 for both DL and SW surveys.While for the second generation case a SW galaxy survey performs significantly better than a DL one, due to the wider coverage of the sky, for the third generation case the two surveys tighten the constraint up to one order of magnitude and have substantially equivalent constraining power, since now GW events can also explore the Universe at higher redshift similarly to the DL survey.
Even though we are not interested in investigating models where Modified Gravity and PBHs describe simultaneously the dark energy and dark matter sectors, respectively, we can briefly comment  on the magnitude of errors on GW effective bias in table 2 and 3. We notice that errors on the GW effective bias appear to be minimally affected by the choice of dark energy model, whether it is a pure cosmological constant, a Modified Gravity or a nDGP model, even if there is some residual correlation that make the constraints looser when deviations from General Relativity are stronger.Therefore, we can rather safely conclude that our prospect of disentangling the origin of BH binary progenitors are only minimally dependent on our description of the dark energy sector.
In figure 3 we show the GW effective bias for the mixed BH populations (top panels) and the corresponding SNR for second and third generation GW detector networks (bottom panels).The SNRs refer to the case where GWs are cross-correlated with the SW galaxy survey since in the second generation case it outperforms the DL one, while in the third they are almost equivalent.We note that second generation GW detector networks are inefficient in discriminating the origin of the BHs because of (i) the low number of detected events, which implies a more noisy measurements of the GW clustering statistics, and (ii) the fact that detected sources are mainly at low redshift, where relative differences between the models are less accentuated, as the interested reader can see in figure 7 of appendix F. On the other hand, the third generation case is radically different.In particular, we will be able to statistically infer the presence of a primordial component of sources at more than 3σ confidence in a large portion of the f PBH − A m .Most notably, if PBHs are the dark matter, i.e., if f PBH = 1, we will be able to assess it at the 20σ level, making GW-galaxy cross-correlation one of the most important probes to determine PBH existence in the O(10) M ⊙ mass range.If we consider an optimistic astrophysical BH model with twice the number of events of the conservative one, the improvement in the SNR in the second generation detector network case is of the order of 30%, hence it is still insufficient in effectively discriminating the nature of the compact objects.On the other hand, the same improvement in the third generation detector network case further expands the parameter space that can be effectively probed by the cross-correlation technique.
In parallel to this work, another article investigated PBH clustering [142].The two approaches are comparable yet different, since we focus our attention to produce realistic catalogs of detected events and on the existence of an intrinsic correlation between different PBH binary formation channels given by the abundance parameter f PBH .The treatment of astrophysical and binary PBH formation uncertainties is also somewhat different but conservative in both cases.Nevertheless, constraints are comparable once we rescale appropriately the early time PBH binary local merger rate and the PBH abundance.

Conclusions
The emerging field of GW astronomy offers outstanding opportunities to address multiple open questions on the nature of our Universe.Future GW detectors, thanks to their increased sensitivity, are expected to observe tens of thousand of GW events per year, shedding light on the physics behind GW sources.Along with this unprecedented opportunity to explore the Universe through GW radiation, the Cosmology community has to face a new challenge: developing new techniques and ideas to fully maximize the scientific outcome of this new kind of dataset.In this work we explore this avenue by investigating whether GW clustering alone or in combination with galaxies can improve our understanding of both the dark matter and dark energy sectors.In fact GWs, since they are emitted by sources which live in a variety of environments (galaxies, dark matter halos, and so on), represent a legitimate tracers of the large-scale structure of the Universe.
Since a priori it is not known which kind of galaxy survey is more suitable to be correlated with GWs, in this work we consider both a shallow and wide (SW) and a deep and localized (DL) galaxy surveys.Moreover, since the intrinsic origin of BHs is currently unknown, we analyze both astrophysical and primordial BHs as potential sources for detected GW events.We created realistic catalogs of GW events following specific state-of-the-art binary models and we determine the clustering properties of GW events detected by second and third generation detector networks using the external modules of CLASS GWB.We generalized the public code Multi CLASS to include extension of General Relativity described by the standard µ − η parametrization, and we use it to compute the GW and galaxy clustering statistics.
Using astrophysical BHs alone, we could constrain extensions of General Relativity approximately at the same level of existing constraints using exclusively linear scales.In other words, GWs represent a tracer of the large-scale structure of the Universe that is independent from, for instance, galaxies, and that can be used to place constraints that are subject to a different set of systematic errors.Moreover, these constraints do not depend on the modelling of non-linear scales, which is not a well explored territory in Modified Gravity theories.By including cross-correlations with galaxies we find competitive constraints even for our conservative astrophysical benchmark model.We also showed how the strength of the constraints explicitly depend on the presence of extra scales in the Modified Gravity theories, providing a fair representation (which is sometimes lacking in the literature) of the constraining power of this new probe.
Regarding the possibility of assessing whether merging binary BHs might have astrophysical or primordial origin, we explored two popular PBH binary formation channels and we compared them to the standard astrophysical one.Even though the number of events detected by a second generation GW detector network at design sensitivity appear to be insufficient to establish the origin of the binary due to volume selection effects that smears out the relative differences between models, we find that third generation detectors will be able to either rule in or out PBHs as a major component of the dark matter at more than 3σ confidence level in the tens of solar masses range.Therefore GW clustering can be effectively used to robustly constraint the abundance of PBH in a mass range that is heavily probed [143][144][145][146][147][148][149][150][151][152][153][154][155][156][157][158][159][160], but still suffers from large theoretical uncertainties due to the theoretical modelling of the observables and the unknown PBH mass distribution, see, e.g., refs.[160,161].Finally, we note that it is more difficult to probe PBH abundances of f PBH ≲ 0.1, independently on the relative abundance of late time and early time PBH binaries.
In this work we considered only dark sirens as a GW tracer to provide a conservative result.However other types of merging events, such as BH-neutron star or neutron star-neutron star mergers, can be reasonably included in this kind of analysis.On top of a straightforward increase in the number of objects used to create GW maps of the sky, which certainly helps in diminishing the noisiness of the maps, in certain instances these kind of events might have a electromagnetic counterpart, which also allows for a better localization.However, when comparing the total number of events detected both by second and third generation GW detector networks, we still find that the number of detected binary BH mergers dominates over BH-neutron star and neutron star-neutron star events.Therefore, even in the (very) optimistic case where every single event involving a neutron star has an electromagnetic counterpart, we do not expect the improvement to be substantial for the kind of analysis performed here since the resolution of the maps would still be limited by binary BH events.
Given the potential of the cross-correlation technique to test beyond General Relativity models, it would be interesting to perform this kind of analysis using a Modified Gravity model for which we know its strong field regime predictions, and to account in a self-consistent way for all aspects of these theories ranging from differences in the emitted GW waveform to differences in the GW propagation in a perturbed Universe.Further investigations on different GW sources also appear to be a promising avenue, including in different frequency bands such as the mHz (probed by LISA) or the nHz (probed by PTA), especially when the properties of these sources heavily rely on properties of their environment.

A Number count contributions
The contributions to ∆ X ℓ (k, z) introduced in equation (2.5) read as (A.1)As we already mentioned, b X is the bias parameter of the tracer X, s X is the magnification bias parameter and f evo X is the evolution bias parameter.Then, r(z) is the radial comoving distance at redshift z, τ z = τ 0 − r(z) and τ = τ 0 − r are conformal times of perturbations that we observe at the comoving distances r(z) and r, with τ 0 being the conformal age of the Universe.Besides, H = a ′ /a is the Hubble expansion rate and the prime ′ denotes a derivative with respect to conformal time, D is the gauge invariant density perturbation, V is the gauge invariant velocity perturbation, Φ and Ψ are Bardeen potentials.Finally, j ℓ = j ℓ (kr(z)) are Bessel functions and

B Implementation of Modified Gravity models in Multi CLASS
The implementation of the Modified Gravity and nDGP models in Multi CLASS follows the approach described in ref. [102].In particular, we modify the evolution equations for the metric potentials, which in the µ − η parametrization now read as where ρtot and ptot are the total background energy density and pressure, respectively, σ tot is the total anisotropic stress, Ω m (z) is the fractional abundance of matter at redshift z and θ m is the matter velocity divergence.The ΛCDM limit is recovered for constant µ = η = 1, as discussed in ref. [102].
It is trivial to take derivatives of the µ − η functions in equations (2.18) in our parametrization of the Modified Gravity scenario, therefore we do not report them here.On the other hand, in the nDGP model, we have Deriving equation (2.20) with respect to the conformal time we obtain where the explicit expression for H ′′ can be calculated from the Friedmann equation 2 Ḣ + 3H  [162,163].Notation follows that of ref. [51].
finding that with pΛ and pm being the cosmological constant and matter background pressure, respectively.The RHS of equation (B.5) is expected to be very suppressed during the matter and dark energy dominated eras.

C Gravitational wave detector networks
In this work we consider both second and third generation GW observatories to compare their performances.Regarding the former, we consider a network made by LIGO Hanford and Livingston, Virgo, KAGRA and LIGO India L-shaped interferometers.On the other hand, for the latter, we choose a network made by Einstein Telescope (composed by three V-shaped interferometer) and Cosmic Explorer (two L-shaped interferometer, with arm length of 40 and 20 km for the main and secondary site, respectively).Existing and future detector location (latitude and longitude) and orientation angles of the arms are reported in table 4.
Regarding ET, we consider the "ET-D" base sensitivity for the one-sided noise spectral density S n .The nominal curve is given for a 10 km L-shaped interferometer hence it has to be rescaled: following the ET community guidelines we use . The frequency range considered for the three detectors is [1,10000] Hz.Noise is assumed to be uncorrelated in the three interferometers.Regarding CE, we consider two L-shaped interferometers, one with 40-km arms, the other with 20km arms.The comparison between second and third generation noise spectral densities is presented on the left panel of figure 4. We choose ρ det = 12 as SNR threshold for both second and third generation detector networks.We show on the right panel of figure 4 the average detection efficiency as a function of redshift for both detector network and for all binary populations.For each binary population the displayed average efficiency is an average of the detection efficiency of five different catalogs of 10 5 events.The noisy behaviour of the curves at high redshift is due to the low number of events in that epoch.Given the catalogs of detected events for different populations, we can estimate what are the typical errors in determining the redshift of the source.We show in top panels of figure 5 the redshift error probability distribution function for different detector networks and binary populations.The errors are estimated for each event using δz event = 3z event /ρ event [70,120] and should be compared to the redshift bin width 2∆z = 0.2 and 2∆z = 0.4 we use for DL and SW galaxy surveys, respectively.We find that the median redshift error is δz ≃ 0.1 for all populations and GW detector networks, making our binning choice broadly compatible with typical error on the source distance.Moreover we explicitly check the magnitude of the average localization error ⟨δz⟩ in each redshift bin, and we show it on the bottom panels of figure 5.These figures further support our binning choice, and they also show the redshift-dependent dispersion σ z = ⟨δz⟩ we later use in the computation of the GW noise in equation (2.12).While for second generation networks the three classes of binaries have similar ρ 2 distributions, thus a similar ⟨δz⟩ value.On the other hand, in the case of third generation detectors, our model of PBH describes binaries that are detected with a higher signal-to-noise ratio compared to astrophysical ones, hence the better localization error in redshift.Finally, for third generation detectors we find a high δz tail due to events located at high redshift, suggesting that the optimal binning strategy might be one that has narrow bins at low redshift and wide bins at high redshift, however we leave the investigating of an optimal binning strategy for future work.
Regarding the determination of the average and maximum resolution, we use the results in figure 3 of ref. [70].In particular, we find that on average the number of GW detected by second and third generation detector network peaks at redshift z 2G peak ≃ 0.3, 0.4, 0.5 and z 3G peak ≃ 1.4, 2.2, 3.4 for astrophysical binaries, late-time and early-time primordial binaries, respectively.Therefore the average resolution for the two classes of detectors is given by θ avg,2G res The maximum resolution is given by θ max,2G res = 0.6 • , 0.8 • , 1.0 • and θ max,3G res = 0.1 • , 0.2 • , 0.4 • , which would correspond to ℓ 2G max = 300, 225, 180 and ℓ 3G max = 1200, 900, 450; however we consider as maximum multipole ℓ max = 200 to avoid modelling the non-linear regime, see, e.g., ref. [52].We note that these estimates are consistent also with more recent ones, see, e.g., ref. [139].

D Late time primordial black hole binary formation
In total generality, any PBH population is described by the fractional abundance function [161] where the abundance parameter f PBH = ρPBH /ρ dm describes the fraction of dark matter in form of PBHs, and dΦ PBH /dM PBH describes the shape of the PBH mass distribution.The latter function is normalized to unity by construction.For the rest of this appendix we consider two benchmark cases: a monochromatic mass distribution (MMD) where all compact object have mass where M ⋆ PBH = 30 M ⊙ is our benchmark value, and a lognormal extended mass distribution [161] (EMD) characterized by two parameters, mean µ and standard deviation σ, with (µ, σ) = (40, 0.5) being our benchmark value.With this choice of values we have an EMD peaked around 30 M ⊙ , mimicking a MMD but with some width.Additionally, we assume PBH to be spinless, i.e., S PBH = 0, hence the (dimensionless) spin parameter reads as χ PBH = cS PBH /GM 2 PBH = 0.This assumption is consistent with the predictions of the popular models where PBH form from the collapse of large overdensities during the radiation-dominated era [164][165][166].
Independently from the origin of the BHs, binaries are described by the total and reduced masses respectively, where M 1 , M 2 are the two compact object masses.

D.1 Late time binaries
The binary formation rate density appearing in equation (3.5) is given by [122,123] Rbf = where M h is the dark matter halo mass, dn h /dM h is the halo number density given in ref. [167] and Rbf,h is the binary formation rate per halo with ⟨σ dc v rel ⟩ is the velocity-averaged direct capture cross section.Accordingly with ref. [122], we consider M min h = 400 M ⊙ as minimum halo mass and, as maximum halo mass, the redshift-dependent maximum mass M max h (z) a virial structure can have at redshift z [168], even though the latter does not influence the total merger rate because it is dominated by low mass halos.

D.1.1 Dark matter halos
Assuming that dark matter halos made of PBH have a similar phenomenology to halo in cold matter scenarios, 7 we describe the spatial distribution of matter inside dark matter halos by a Navarro-Frenk-White density profile [169]  where ρ s and r s are the characteristic density and radius.Halos are commonly described by virial radius r h , average density ρ c ∆ c , where ρ c = 3H 2 (z)/8πG is the critical density of the Universe and ∆ c = 200 is the most common choice for the virial overdensity threshold, and mass where we define the concentration parameter C = r h /r s and the function From equation (D.8) we have that the characteristic density and radius read in terms of the concentration parameter as therefore the integration over the halo profile in equation (D.6) factorizes out and is given by The concentration parameter is a function both of halo mass and redshift, and its value can be inferred and fitted from numerical simulations [170,171].In this work we use the fitting formulas of ref. [171], which reads as where ν = δ c /σ(M, z) is the "peak height", δ c = 1.686 is the linearly extrapolated critical overdensity for spherical collapse in Einstein-de Sitter, σ is the variance of the matter field, the fitting parameters are given by and D(z) is the linear growth factor.For the matter variance and linear growth factor we use the fitting formulas of refs.[167] and [172,173], respectively (see also appendix D of ref. [51]).
Regarding the distribution of velocities inside dark matter halos, it is helpful to define two quantities starting from the expression for the circular velocity v c = GM (r)/r for a halo in virial equilibrium: v h and v dm .The former is called "virial" velocity and it is the circular velocity at the edge of the halo, i.e., while the latter is defined as the escape velocity at radius r max = C max r s , with C max = 2.1626, which represents the position of the maximum of the circular velocity, and reads as Numerical N-body simulations suggest that the relative velocity v rel probability distribution function of pair of objects within a dark matter halo can be described by a truncated Maxwell-Boltzmann distribution, i.e., as [122,174,175] p

D.1.2 Direct capture process
The direct capture cross section reads as [176] σ Since the cross section strongly depends on the inverse of the relative velocity, the process will be more effective in low-mass dark matter halos.Using equation (D.16) we find the velocity-averaged cross section where a useful reference value is given in the MMD case by M MMD PBH = 4 5/7 , independently from where the mass distribution is peaked.
We follow ref.[175] to model the initial properties of the binaries: they form with an initial semi-major axis and eccentricity where the energy of the binary at the beginning of the inspiral phase We show in figure 6 the initial eccentricity and time delay probability distribution function obtained for our choice of MMD and EMD as example of our study of the binary properties.First of all we note by comparing the MMD and EMD histograms that binaries properties are quite independent from the details of the PBH mass function.As expected all the binaries are created with very high eccentricities.Second, we prove that the time delay probability distribution function scales approximately as p(t d ) ∝ t −1 d , with minimum delay of t d,min ≃ 1 s and maximum one that can exceed the lifetime of the Universe, corresponding to binaries that never merge.Therefore, we conclude that the maximum time delay should be set to be t d,max = t(z), with t(z) being the age of the Universe at redshift z.

E Optimistic cross-correlation constraints
Table 5 shows the marginalized errors on Modified Gravity and nDGP parameters with and without considering perfect knowledge of the standard cosmological parameters and assuming an optimistic version of our benchmark model.The optimistic model assumes that the number of observed GW events per year is twice that of the conservative case (hence the shot noise is roughly half of that of the conservative case), in line with current estimates of certain astrophysical BH models [139].) parameters for single (GW only) and multi-tracer analysis.SW and DL refer to "shallow and wide" and "deep and localized" galaxy surveys, respectively.Errors refer to our optimistic benchmark model.Numbers in parenthesis are the marginalized errors obtained assuming perfect knowledge of the standard cosmological parameters.In all the scenarios we observe a redshift distribution peaked z ≃ 0.3 − 0.5, with differences between models suppressed by the detector volume selection effects.The dispersion of bias, magnification bias and evolution bias values around the peak of the redshift distribution is rather small, lowering the possibility of disentangling how much each component contributes to the total population of observed events.In this plot we consider only events detected by a second generation GW detector network.The light blue, blue and dark blue lines correspond to the astrophysical BHs only (fPBH = Am = 0), mixed astrophysical and late time PBH (fPBH = 1, Am = 0), and PBHs only (fPBH = 1, Am = 18) scenarios, respectively.The green lines refer to mixed models where different populations contribute to the total number of detected events.

F Mixed black hole population for second generation detectors
65, 1.45, −0.14, 0.01} for second and third generation detector networks, respectively.The effective GW bias is given by b 2G,eff GW = 1.29 and b 3G,eff GW = 2.86, in the redshift range probed by the DL survey, and by b 2G,eff GW = 1.14 and b 3G,eff GW = 2.84 in the SW one.Regarding the GW magnification bias s GW (z) = − d log 10 d2NGW dzdΩ dρ ρ=12 (3.4)

Figure 1 .
Figure 1.GW mixed population redshift distribution (top left panel ), bias (top right panel ), magnification bias (bottom left panel ) and evolution bias (bottom right panel ) for different values of fPBH and Am.In this plot we consider only events detected by a third generation GW detector network.The light blue, blue and dark blue lines correspond to the astrophysical BHs only (fPBH = Am = 0), mixed astrophysical and late time PBH (fPBH = 1, Am = 0), and PBHs only (fPBH = 1, Am = 18) scenarios, respectively.In the mixed astrophysical and late time PBH scenario, late time PBHs contribute to approximately 6% of the total number of observed events.The green lines refer to mixed models where different populations contribute to the total number of detected events.

Figure 2 .
Figure 2. Galaxy angular power spectrum relative difference between Modified Gravity (MG) scenarios with respect to the ΛCDM case for redshift bins centered at z = 0.4 (left panel ) and z = 3.6 (right panel ).The Modified Gravity extra scale dependence is set at kmg → ∞ (yellow line) and kmg = 10 −1 , 10 −2 and 10 −3 Mpc −1(orange, red and dark red lines, respectively).We show the effect only for galaxies, however also the GW angular power spectrum and the galaxy-GW angular power spectra share a similar phenomenology.

SNRFigure 3 .
Figure 3.Total effective bias of the mixed BH binary population (top panels) and SNR obtained assuming astrophysical BHs are the fiducial model (bottom panels), both for second (left panels) and third generation (right panels) GW detector network.The fiducial model is given by the conservative astrophysical BH benchmark model.The effective bias of the alternative model is defined in equation (3.8).For the HLVIK case we consider the cross-correlation with the SW survey since it is the most promising case, while for the ET2CE case the results are equivalent for both kind of galaxy surveys.

Figure 4 .
Figure 4. Left panel : One-sided noise spectral density for individual detectors belonging to second and third generation detector networks.Right panel : Average detection efficiencies for different BBH sources and detector networks.

Figure 5 .
Figure 5. Top panels: redshift error probability distribution function for second (left panel ) and third (right panel ) generation GW detector network for all binary populations considered in this work.Bottom panels: average redshift localization error in different redshift bins for second (left panel ) and third (right panel ) generation GW detector network for all binary populations considered in this work.Solid and dashed lines refer to the redshift binning chosen for the SW and DL surveys, respectively.For comparison, the width of redshift bins in our analysis is 2∆z = 0.4 and 2∆z = 0.2 for SW and DL surveys, respectively.

2/ 7 G 2 M 12 /7 µ 2/ 7 c
10/7 v 11/7 dm5Γ[5/7, 0] − 5Γ[5/7, V 2 ] − 7V 10/7 e −V 2 3 √ πErf[V ] − (6V + 4V 3 )e −V 2 , (D.18)whereV = v h /v dm and Γ[n, x] = ∞ x dte −t t n−1 is the incomplete Gamma function.Interestingly, in equation (D.6), the dependence on the compact objects mass factorizes out of the merger rate, i.e., as pointed out in ref. [18], we can rescale the binary formation rate for different PBH mass distributions and abundances just by recomputing the quantity M PBH = dM 1 dM 2 dΦ PBH dM 1 impact parameter b.The impact parameter is assumed to be uniformly distributed between the minimum and maximum impact parameters time delay between binary formation and merge is given by statistically sampling the initial properties of the binaries in a similar fashion to what presented in ref.[51], we can determine the time delay probability distribution function p(t d ) entering in equation(3.5).

Figure 6 .
Figure 6.Initial eccentricity (left panel ) and time delay (right panel ) probability distribution function for our choice of benchmark MMD (blue lines) and EMD (orange lines).

Figure 7
Figure7shows the total redshift distribution, bias, magnification bias and evolution bias for different models spanning the f PBH − A m parameter space for a second generation detector network, along with three reference scenarios corresponding to the astrophysical BHs only (f PBH = A m = 0), mixed astrophysical and late time PBH (f PBH = 1, A m = 0), and PBHs only (f PBH = 1, A m = 18) cases.In all the scenarios we observe a redshift distribution peaked z ≃ 0.3 − 0.5, with differences between models suppressed by the detector volume selection effects.The dispersion of bias, magnification bias and evolution bias values around the peak of the redshift distribution is rather small, lowering the possibility of disentangling how much each component contributes to the total population of observed events.

Figure 7 .
Figure 7. GW mixed population redshift distribution (top left panel ), bias (top right panel ), magnification bias (bottom left panel ) and evolution bias (bottom right panel ) for different values of fPBH and Am.In this plot we consider only events detected by a second generation GW detector network.The light blue, blue and dark blue lines correspond to the astrophysical BHs only (fPBH = Am = 0), mixed astrophysical and late time PBH (fPBH = 1, Am = 0), and PBHs only (fPBH = 1, Am = 18) scenarios, respectively.The green lines refer to mixed models where different populations contribute to the total number of detected events.
to compute the intrinsic merger rate density for PBH late time binaries as RLPBH is the time delay between merger and binary formation, p(t d ) is the time delay probability density function, z f is the binary formation redshift and Rbf is the intrinsic binary formation rate density (see appendix D for details on its computation).By simulating catalogs of binaries in a similar fashion to what we do to create GW events, we derived for the first time the expected time delay probability distribution function and, as we show in appendix D, we can approximate it as p(t d ) ∝ t −1 d .However, since t d,max ≃ t(z), where t(z) is the age of the Universe at redshift z, the time delay is non-negligible.The local merger rate of RLPBH

Table 1 .
Constraints on standard cosmological parameters from the DL and SW galaxy surveys.Errors are calculated for the Modified Gravity model with kmg = 10 −3 Mpc −1 , however they are consistent also with the errors found for other choices of kmg and for the nDGP model.

Table 2 .
Marginalized errors on Modified Gravity (upper table) and nDGP (lower table) parameters for single (GW only) and multi-tracer analysis.SW and DL refer to "shallow and wide" and "deep and localized" galaxy surveys, respectively.Errors refer to our conservative bechmark model.Numbers in parenthesis are the marginalized errors obtained assuming perfect knowledge of the standard cosmological parameters.

Table 3 .
Constraints on the GW effective bias parameter using second (left table) and third generation (right table) detector networks.SW and DL refer to "shallow and wide" and "deep and localized" galaxy surveys, respectively.

Table 5 .
Marginalized errors on Modified Gravity (upper table) and nDGP (lower table