Number count of gravitational waves and supernovae in luminosity distance space for ΛCDM and scalar-tensor theories

The clustering of gravitational waves in luminosity distance space is emerging as a promising probe of the growth of structure. Just like for galaxies, its observation is subject to a number of relativistic corrections that affect the measured signal and need to be accounted for when fitting theoretical models to the data. We derive the full expression for the number count of gravitational waves in luminosity distance space, including all relativistic corrections, in ΛCDM and in scalar-tensor theories with luminal propagation of tensors. We investigate the importance of each relativistic effect and the detectability of the total signal by current and planned GW detectors. We consider also supernovae in luminosity distance space, highlighting the differences with gravitational waves in the case of scalar-tensor theories. We carry out a thorough comparison among the number count of gravitational waves and supernovae in luminosity distance space, and that of galaxies in redshift space. We show how the relativistic corrections contain useful complementary information on the growth of perturbations and on the underlying theory of gravity, highlighting the synergy with other cosmological probes.


Introduction
The direct detection of gravitational waves (GW) [1] has opened a new window onto the Universe, that can be highly complementary to EM observations.The use of GW as probes of cosmology was first advanced in the '80s by Schutz [2], who proposed the standard sirens method to measure the present value of the Hubble constant.This became a reality in 2017, when LIGO-Virgo-Kagra (LVK) detected the first merger of two binary neutron stars, GW170817 [3], while a network of EM detectors was able to measure the coincident EM emission from the merger and remnant [4], placing the value of H 0 at 70.0 +12.0 −8.0 [5].However, events such as GW170817 are a lucky rarity, while most of the GW detected in LVK are not accompanied by a EM counterpart, nor they are expected to be in the future.Moreover, standard sirens also present a big limitation, in that they are only sensitive to the subset of cosmological parameters that determines the background expansion.At the same time, they are totally unaffected by other important parameters that drive the growth of structure.To efficiently test the full extent of the cosmological model then, other methodologies must be investigated.
One possibility that took hold in recent years is to use GWs as tracers of the large scale structure (LSS), analogously to what is normally done with galaxies.This builds on the reasonable assumption that, with the possible exception of exotic sources such as primordial black holes, the binaries that generate GWs are located within or in the proximity of galaxies, and share similar clustering properties.On top of that, the GWs that reach the detector have travelled through a Universe rich of small, local perturbations, which have imprinted on the wave by slightly altering its shape [6][7][8][9].There are two complementary roads to approaching LSS studies with GW: the weak lensing convergence field can be probed by measuring perturbations induced on the amplitude of the incoming wave, as explored e.g. in [10][11][12][13][14][15]; or, the number count, or clustering, of compact objects -defined as the local abundance in the number of detected GW sources in the sky -which can by itself constitute the observable.This second possibility has been explored before in [16][17][18][19][20][21][22], and it is also the focus of this work.
It is well known that, when performing clustering studies with galaxies, the signal is affected by relativistic effects due to the fact that the number density of galaxies is measured in redshift bins, whose actual central value and volume are distorted by perturbations locally and along the line of sight.These effects have first been derived in the linear regime in [23][24][25][26][27][28][29].Similar caution must be applied in the case of GWs, with the crucial difference that the number density is now measured in luminosity distance space (or D-space).In fact, except for few cases in which the redshift of the source can be determined, GWs are dark standard sirens.The inclusion of the relativistic effects for the clustering of GWs in luminosity distance space has been explored in [17][18][19]30], where lensing and luminosity distance space distortions where taken into account.More recently, [22] provided a full theoretical expression of the number count of GWs in luminosity distance space, including all corrections for a ΛCDM Universe.This work retraced the calculation performed in redshift space (or z-space) by [28], connecting the GWs number count to the luminosity distance perturbations of photons and GWs derived in [8,[31][32][33][34].In this paper, we build on [8,9] and [22] by further extending the calculation to the broader class of scalar-tensor theories of gravity with luminal speed of sound for the propagating tensors.
The paper is organized as follows.In Sec. 2 we re-derive the result of [22] in a more general form, without ever making any assumption about the background luminosity distanceredshift relation, so that it can be readily generalized to theories beyond GR.In Sec. 3 we specialize scalar-tensor (ST) theories, writing the full expression of the luminosity distance number count for the sub-set of Horndeski theories with c T = 1.We caution that, in the GR limit, our result reproduces that of [22], up to a minus sign in the coefficients of six relativistic corrections.We comment more about this in the text.In Sec. 4, we extend our result to supernovae (SN).The latter also live naturally in luminosity distance space, however, in the case of ST, they do not receive any explicit corrections from the modifications of the theory (see e.g.[9]).Thus, the number count of SN in luminosity distance space will have a different form than that of GWs.We provide both expressions and comment on the implications for testing gravity theories in Sec. 4. In Sec. 4, we also build the 2point angular correlation functions of the luminosity distance number count for GWs and SN, and compare them with the case of galaxies (in redshift space), exploring the potentialities of a joint cosmological analysis.We investigate the weight that the individual relativistic components have on the total signal, indicating under which circumstances they can be safely neglected.Finally, in Sec. 5 we assess the potentiality for observation of the number count correlation of compact objects, by plotting its cumulative signal-to-noise ratio for a set of different benchmark scenarios, considering 3rd and 4th generation GW detectors.We draw our conclusions in Sec. 6.

Number Count in Luminosity Distance Space
Throughout our derivation for the number count in luminosity distance space, we will distinguish between two different reference frames: the smooth, background, frame and the observed, or perturbed, frame.This distinction is necessary because the local inhomogeneities in the Universe perturb the propagation of photons and gravitational waves [8,31].The direct consequence is that all astrophysical parameters inferred by measuring the perturbed waves are no longer the true, physical ones, but will differ by a small amount.
More rigorously, the perturbed frame contains all measurable quantities.In this frame, the past comoving lightcone of GWs (or photons) is not flat, as the wave propagation is perturbed by relativistic effects induced by the large scale structures.We use un-barred comoving coordinates (η, x), and we describe the past lightcone in terms of the corresponding comoving distance χ.In the observed frame, every celestial object has a redshift z(χ) and a corresponding luminosity distance D(χ).Both quantities are the ones that we can infer directly from data, modulus the measurement uncertainty.
At the same time, each object also exists in the background frame.Here, it has a comoving position identified by barred coordinates (η, x) and its past lightcone, parametrized in terms of the true comoving distance χ, is flat, so that the background comoving geodesic has coordinates xµ = (η 0 − χ, n χ) . (2.1) Here, n is the direction of arrival of the incoming wave, i.e. the unit vector is oriented from the observer to the source, ni ≡ xi / χ.This means that the total derivative along the wave geodesic can be computed as In this frame, the object has a background redshift z( χ), with corresponding luminosity distance D( χ).These can be thought of as the true redshift and luminosity distance of the object or, analogously, the redshift and distance that would be inferred through measurements of photons and gravitational waves if the Universe was perfectly homogeneous on all scales.We now define the quantity ∆D as the difference between the luminosity distance in the perturbed frame and the one in the background frame, i.e.

∆D(χ, n)
where we have further defined δD(z, n) ≡ D(z(χ), n) − D(z( χ)).In other words, we can split the luminosity distance perturbations in two contributions: δD is the effect resulting from the perturbation of the past lightcone, which introduces a difference between the observed and physical luminosity distance even when these are evaluated at the same value of the affine parameter; the term proportional to δχ accounts instead for the variation of the affine parameter itself in going from one frame to the other.Note that for few terms in Eq. (2.3) we introduced a dependency on the specific line of sight n.This is because the perturbed luminosity distance observed for each object depends on the specific configuration of the inhomogeneities inducing the wave perturbations, both at the source location (as is the case for perturbations induced by the local potentials and peculiar velocities) and along the wave path (as is the case for integrated effects, such as lensing).Eq. ( 2.3) can also be re-written in terms of the observed redshift z(χ) and the true redshift z( χ) as The mapping of the observed universe onto the smooth background Friedmann universe is necessary, if one wishes to exploit all standard cosmological relations connecting the background quantities such as D, z, χ, η etc.However, in general it introduces a gauge dependence for the perturbations, as they will now depend on the specific choice made for the background, which is known as the cosmological gauge problem [35,36].This gauge dependence will cancel out in the final expression of the observed number count defined below (see Eq. (2.5)), as ∆ is built fully from observable quantities [28].
Having made this distinction between frames, we can now give a formal definition for the number count, following the standard definition used for galaxies, translated in luminosity distance space [16,17,19,22] where N (D, n) is the number of objects detected within a small volume at observed luminosity distance D and in direction n, and N (D) is its angular average over the whole sky.As it is built from observational data, N lives in the perturbed frame.On the contrary, N (D) is an ambiguous object, as its explicit dependence on the specific perturbations of the wave lightcone (i.e. on n) has been removed by the average, but it is still evaluated at the perturbed luminosity distance.Thus, to fully work out the theoretical expression for ∆, we need to further expand N and refer it to background quantities.Before doing so, we follow [28] and use that N ≡ ρV , where ρ is the number density of events, and V is the small sky-volume considered.In terms of averaged quantities, we also have N (D) = ρ(D) V (D), and we can relate averaged and un-averaged as analogously to what we did in Eq. (2.3).Here δρ and δV are assumed to be small.Inserting Eqs.(2.6) in (2.5) and keeping only linear order terms in the perturbations, we obtain i.e. we can split the number count into the sum of two separate perturbations, one on the number density and one on the volume.We can now further expand the barred quantities at the RHS of Eq. (2.7), and refer them to the background distance, obtaining, at linear order The first term in Eq. (2.8) is the standard definition of the density contrast, or clustering, of gravitational waves δ N gw in Newtonian gauge.The same quantity in synchronous gauge can be linked to the density contrast of the underlying dark matter distribution through a bias parameter b gw so that δ S gw = b gw δ S m .The calculation of the three remaining terms can be pushed forward explicitly, keeping into account the dependence of the background density and volume on the background distance.We do so in App.A.1 for the density term and in App.A.3 for the volume terms.In doing so, we follow closely the calculation made in ΛCDM by [28], that was recently adapted to the luminosity distance case by [22].However, we maintain an extra degree of generality in our calculations, so that we can then proceed to extend the result to scalar-tensor theories of gravity.In particular, we first compute the RHS of Eq. (2.8) without ever substituting i) the explicit expression of the luminosity distance perturbations ∆D and ii) the explicit dependence of the background luminosity distance D and its first and second derivatives on the background redshift z.
We also work in Newtonian gauge, for which the metric perturbations in the perturbed frame are our result for the luminosity distance number count in its general form is2 (2.11) In the first line of Eq. (2.11) we can already recognise a number of perturbative terms that have similar origin to those affecting the number count of galaxies.In order, these are the density contrast term, the time delay term, a lensing term, the Doppler perturbation and the perturbations due to the local potentials at the source position.More perturbative terms are hidden in the second line of the equation, and will become explicit once we write the full expression for the luminosity distance perturbations in the desired MG theories, and compute their derivatives along geodesics.

GW number count in scalar-tensor theories of gravity
We consider the class of scalar-tensor theories described by the following action: where φ is the additional scalar degree of freedom, X its kinetic term, X = −∂ µ φ∂ µ φ/2, F (φ) represents a non-minimal coupling and G, K are general functions of φ and X which can lead to kinetic braiding and non-trivial kineticity of the scalar field.Action (3.1) corresponds to the subset of the broader Horndeski class of theories [37] that leads to tensors propagating at the speed of light.This requirement is inspired by the LIGO observations of the binary neutron star system GW170817 [38], which placed significant limitations to the class of surviving MG theories, as investigated in [39][40][41] (see also [42] for a discussion of caveats).More concretely, the c 2 T = 1 restriction is dictated by the fact that relativistic corrections to the luminosity distance of GWs have been calculated only for Horndeski theories with c 2 T = 1 [9].In [43], some of the authors generalized the calculation to DHOST [44] theories with c 2 T = 1.It would be straightforward to apply our calculation to this broader class, however the difference is negligible in terms of impacts on the number count and, for simplicity, in this work we focus on the Horndeski case.
We use the following convention for splitting the scalar field into a background and perturbation component: where again, δ represents variations of the field evaluated at the same affine parameter, while ∆ describes the full perturbation.We stress here that the perturbations of the scalar field, similarly to Ψ and Φ for the perturbed metric in (2.9), are considered to be long wavelength perturbations, clearly distinguished from the high frequency gravitational waves.The coefficient functions in the EFT action are perturbed accordingly.In particular, the function F (φ) that introduces the non-minimal coupling can be expanded as where we have defined F φ( φ) ≡ d F /dφ| φ.It was shown in [9] that, for this set of theories, the background luminosity distance inferred through the measurement of a GW event is itself modified as where DEM can be thought of as the corresponding background luminosity distance that would be probed with an electromagnetic source in the same location, such as a supernova. 3e can use this relation to calculate explicit expressions for the first and second derivative of the background luminosity distance w.r.t. the background redshift in Eq. (2.11).For the first derivative we have where in the last step we have used the definition α M ≡ Ḟ/( F H).For the subset of theories that we are considering in this study, this function is indeed the same α M parameter that was first proposed in [45] as part of a minimal parametrization of Horndeski theories.This quantity is a measurement of the running of the Planck mass, more properly defined as α M ≡ H −1 d log M 2 P /dη.This definition in general implies that α M gets contributions from different operators in the most general Horndeski lagrangian.However, considering Eq. (3.1), the only surviving function that still contributes to α M in our case is F (φ), so that our definition matches the standard one for the theories explored in this work.In Eq. (3.5) we have introduced α M for two reasons: first, it will lighten significantly the notation for the calculations that follow; additionally, it is the function most commonly used in the literature to capture deviations from ΛCDM through GW at the background level, as it modifies the GW propagation in a simple way, introducing an additional friction term in the propagation equation ( [46]).However, this does not mean that from here on we will fully adopt the parametrization of [45], and indeed, we will clarify this more in detail in the next section.For the time being, the introduction of α M should be regarded primarily as a convenient way to simplify the notation.
Recalling the definition of γ (2.10), Eq. (3.5) leads to For the second derivative of the background luminosity distance, we obtain The last quantity we need to make explicit is the total perturbation of the luminosity distance ∆D.This was calculated in full in the ΛCDM case in [8,31,33], and subsequently extended to the class of ST theories considered in this work by [9] (hereafter, G20) where the modifications due to the scalar field manifest explicitly in the second line.
One important consideration about Eq. (3.8) is that ∆D G20 ≡ D[z(χ), n] − D[z( χ)] was derived under the condition 1 + z(χ) = 1 + z( χ) 4 .This was done in the spirit of working primarily in terms of observable quantities, where the only relevant redshift is z, the observed one.When observed, this redshift can be interpreted either as the perturbed redshift of a wave propagating along a perturbed geodesics, or the unperturbed redshift of an unperturbed wave originally emitted at redshift z.In simpler words, ∆D G20 can be expressed in terms of redshift as ∆D G20 = D(z, n) − D(z), where D(z) is the true distance of a source at true redshift z.Therefore, ∆D G20 accounts for the full perturbations of the past lightcone, as can be seen combining Eqs.(5.39) and (5.63) in [9], and to connect it to the ∆D of Eq. (2.4) required for the number count, it is sufficient to perform a further expansion in terms of the background redshift z, so that In the last equation, we have used the definition of γ Eq. (2.10), while ∆z is the difference between the perturbed and unperturbed redshift of the source, and ∆ ln a is the corresponding perturbation of the scale factor, that can be expressed in Newtonian gauge as (see [9]) Several coefficients have simplified, and now α M enters only in the coefficient for the scalar field contribution.We compute the derivative of ∆D/ D along the null geodesics as (3.12) Finally, we substitute Eqs.(3.6), (3.7) and (3.12) into Eq.(2.11) to obtain where we have further defined Eq. (3.13) represents the final theoretical expression, at linear order, for the number count of gravitational waves in the set of scalar-tensor theories of gravity considered.We can recognise a number of perturbative terms contributing to the number count, analogous to the ones affecting galaxies (see e.g.[28]).All of them, except from the GW clustering δ N GW , that constitutes the core of our signal, are induced by relativistic effects responsible for the GW perturbed propagation.In order, these terms are: the density contrast of GW events in Newtonian gauge, the lensing term, the luminosity distance space distortions (DSD) termanalogous to the redshift space distortions of galaxies -the Doppler term, the Shapiro time delay term, the ISW term, four local potential terms containing the metric perturbations Ψ and Φ and their derivatives evaluated at the source position, and two additional terms, that vanish in the ΛCDM limit, which are the relativistic effects linked to variations in the scalar field.We stress that these last two terms are not the only explicit contribution to the GW number count introduced by the scalar field, as the coefficient of all terms contain the functions γ and β, which in turn explicitly contain α M .Moreover, each of the terms will be implicitly affected by the modifications of gravity, as the perturbations grow and evolve differently in different theories.As such, the GW number count has the potential to become a valuable probe to test gravity models on cosmological scales.We will investigate this potential in the next sections.
We note that our calculation for the number count does not immediately apply to other EM sources living in luminosity distance space, such as supernovae.Indeed, the propagation of photons is affected differently by modifications of gravity.At the background level, the luminosity distance of an EM source is defined in Eq. (3.4), and this will affect the following Eqs.(3.5), (3.6) and (3.7),where the α M contribution vanishes.Eq. (3.11) is also different in the case of photons: here α M = 0 implies that there is no perturbative term due to the fluctuations of the scalar field δφ.The net effect of these differences on the final result however is quite easy to work out.Indeed, the number count of EM sources in ST theories will just be Eq.(3.13), with the coefficients of all remaining terms taken in the ΛCDM limit, i.e. it will be These differences open the interesting possibility of adding constraining power over MG models not only by combining the GW signal with galaxies, but also with other EM sources living in luminosity distance space, similarly to what was done in [43] to isolate the scalar field fluctuations.
Let us conclude this part, noticing that the ΛCDM limit of Eq. (3.13) does not match the results of [22] (Hereafter, F23) -note that their definition of β is such that The discrepancy is due to the sign difference pointed out in Eq. (2.11), propagated to the final result, and affects all terms involving β, i.e. lensing, Doppler, time delay, ISW and two potential terms.Though, as we will show, most of these terms are negligible corrections to the number count, the differences induced on the lensing term are important.We comment more on this in Sec.4.5.

Magnification and evolution biases
The observation of any object in the sky is limited by a detector sensitivity.For both supernovae and gravitational waves, this takes the form of a minimum flux, below which the object goes undetected.As flux is linked to the luminosity distance of the source, a variation in the observed luminosity distance induced by perturbations can change the observed flux, pushing sources above or below the detection limit, and in turn affecting the number count.This effect is widely known even in the case of galaxy surveys, and goes by the name of magnification bias.The magnification bias for the luminosity distance number count was derived for a ΛCDM universe in [22].Here, we retrace their derivation, highlighting, if any, the additional contributions coming from the presence of the scalar field.
For SN, the flux limit of the survey is translated to an upper limit on their apparent magnitude m, which in turn is connected to the luminosity distance of the source through M = m − 5 log 10 (D/Mpc) + 25 . (3.16) Being SN standard candles, any perturbation in the luminosity distance is translated into a corresponding perturbation to the apparent magnitude This in turn implies a correction to the SN number density, i.e.Eq (A.1) in App.A.1 becomes Remembering that ρ = a −3 n, with n the comoving density of sources, the net perturbation to the number density can be written as Any MG effect will enter explicitly in the above expression through the luminosity distance perturbations Eq. (3.11).
Gravitational waves have also a detection threshold, which is typically chosen to be a minimum value for the signal-to-noise of the GW in the detector.This SNR Θ is defined as where S n (f ) is the noise spectral density of the specific detector, and h(f ) is the GW strain in frequency domain, which in general reads [47] h(f Here, A(f ) and ψ(f ) are respectively the amplitude and phase of the GW wave, ι is the inclination of the orbital plane of the binary with respect to the line of sight and F + and F × are the frequency-dependent pattern functions that model the detector response to each polarization of the wave.Then, the signal spectral power becomes i.e., for the SNR we can write where we have used that, at any post-newtonian order, the amplitude of the wave can always be written as A(f ) = Q(f )/D, with Q(f ) a function including in its definition the intrinsic parameters of the binary.From Eq. (3.23), we now appreciate that a perturbation on the luminosity distance induces a fluctuation in the GW SNR that can be written as from which we obtain In total analogy with SN surveys, the corresponding perturbation produced on ∆GW (D, n) then reads This expression carries the imprint of deviations from ΛCDM at two levels: first, through the luminosity distance perturbations Eq. (3.11); then, because of the implicit dependence of Θ( D) on the non-minimal coupling, through Eq. (3.4), which can in turn impact the SNR evolution of the comoving density of sources, ∂ ln n/∂ ln Θ.
The magnification bias discussed in this section affects the luminosity distance number count in addition to the evolution bias ∂ ln n/∂ ln a presented in Sec. 2. The proper modeling of both quantities for GW requires the accurate characterization of the comoving density of events, which in turns needs to rely heavily on astrophysical models for the redshifts evolution of the binary rates.This task is complex and delicate, and we do not address it in the present study, though we stress its investigation is an important ingredient of any complete cosmological forecast analysis.See [19,48] for some recent work in this direction.For the rest of the paper, we will simply set 4 Tomographic observables in Scalar-Tensor theories The number count ∆(D, n) is not yet the actual observable.This can be obtained by integrating Eq. (3.13) over an observed luminosity distance window function, so that we get where we are taking a tomographic approach and W i (D) becomes then the window function identifying the i th tomographic luminosity distance bin, centered around D i .Binning has the unfortunate effect of partially smoothing out the signal, decreasing its strength, but is necessary as this window accounts for the measurement uncertainty on the observed luminosity distances, which depends on the specific detector considered.The smaller this error, the finer the binning can be, the strongest the observed number count signal.The last equality in Eq. (4.1) holds as we are doing a linear calculation and ∆ is already a first order perturbation.In what follows, as there is no more room for confusion, we drop the bar to indicate background quantities and always implicitly intend all objects evaluated on the background.Eq. (4.1) is correctly integrated in luminosity distance space, where ∆ can be measured.However, we are interested in exploring the potential of this observable as a cosmological probe, by itself and in correlation with other tracers, such as galaxies and SN.To this extent we re-parametrize this integral in terms of conformal time with D = (η 0 − η) √ F /a, η 0 being the conformal time at observation and These expressions can then be implemented in an Einstein-Boltzmann solver; in our case, we employ EFTCAMB [49,50], an effective field theory extension of the publicly available CAMB [51,52].
The standard cosmological observables associated with the number count are the 2-point angular correlation functions.To extract these correlations from Eq. (4.2), we manipulate ∆ obs going to Fourier-space, and expanding each mode in spherical harmonics such that The coefficients ∆ obs ℓm (D i ) are given by the term in squared parenthesis.The 2-point correlation functions are then defined in the standard way as where the bracket on the LHS indicates an average over large volumes.Computing the average and exploiting the normalization properties of the spherical harmonics, it follows easily that where P δ (k, η, η ′ ) is the a-dimensional power spectrum defined as connected to the primordial power spectrum where ∆ T is the transfer function encoding the time evolution of ∆ from the primordial inflationary perturbation.

Scalar-Tensor models in the EFT formalism
We work with EFTCAMB [49,50] (a patch to the publicly available Einstein-Boltzmann solver CAMB [51,52]), which uses the framework of effective field theory of dark energy (EFTofDE) [53][54][55][56][57][58][59] to explore ST theories.We adopt the convention for the EFTofDE action defined in [49,50], where background and linear dynamics for our action (3.1) are equivalently described by the following quadratic action: Model  where we have adopted the notation used in EFTCAMB [49,50] and m 0 = (8πG) −1 .The corresponding running Planck mass is which is consistent with the fact that the mapping of Eq. (3.1) into Eq.(4.9) gives The free functions of time multiplying the different operators, {Ω, c, Λ, γ 1 , γ 2 } are dubbed EFT functions.The first three impact both the background and linear dynamics, while the γs affect only perturbations.We work in the so-called designer approach, where we make a choice for the expansion history and Ω, and rely on the Friedmann equations to correspondingly fix Λ and c.Consequently, the full phenomenology of action Eq.(3.1) on linear scales can be reproduced by varying four free functions of time: {Ω, Λ, γ 1 , γ 2 }.
We explore different ST models belonging to Eq. (3.1), via a parametrization of the expansion history and {Ω(η), γ 1 (η), γ 2 (η)}.For the background, we parametrize the equation of state of dark energy following [60,61] Here, ρ DE is the background dark energy density, and ρ 0 DE its value today.For the EFT function, we choose: With this set of definitions, any MG model is completely specified with the choice of five parameters: w 0 , w a , Ω 0 , γ 1,0 and γ 2,0 .Additionally, we must consider the six additional parameters that determine the standard background cosmology.These are: the present value of the Hubble parameter H 0 , the density contrast of matter Ω m,0 , baryons Ω b,0 and curvature Ω K,0 today, the primordial power spectrum spectral index n s and the amplitude σ 8 of the linear matter power spectrum at 8Mpc/h.We keep these parameters fixed at the best fit of Planck (2018) [62].
We explore the number count correlations in five different models within action (3.1), namely: the ΛCDM cosmological model; a wCDM model, i.e. a model with a non-trivial dark energy equation of state, but all EFT functions set to zero; a Generalized Brans-Dicke model (GBD, [63]), characterised by non-minimal coupling and a standard kinetic term; a Kinetic Gravity Braiding model (KGB, [64]); and an f (R) model ( [65]) on a ΛCDM background.For the latter, we adopt the built-in full mapping module of EFTCAMB [66], that implements a family of designer f (R) models on ΛCDM background, labeled by the parameter B 0 .This parameter is representative of the mass scale of the scalar degree of freedom today, in units of horizon scale, and can be related to Ω 0 [66].We summarize the fiducial values of the DE and EFT parameters adopted for each model in Tab. 1.For GBD and KGB we chose as fiducial model the best fit values of [67], while the fiducial for f (R) is chosen to be compatible with the bounds of [50].
4.2 GW, SN and galaxies: angular correlations in D-space vs z-space In this section and the following, we plot different quantities related to the theoretical correlations of Eq. (4.6).In doing so, we make a number of choices.First, as discussed in Sec.3.1 we neglect the evolution and magnification biases, i.e. we fix For the bias factor b gw linking the GW density contrast to the underlying dark matter density contrast, δ S gw = b gw δ S m (and similar for SN), for simplicity, in this work we limit to b gw = b sn = 1.Then, we need to fix the shape of the window function in Eq. (4.3).We choose the adopt the redshift gaussian window function already implemented in CAMB, i.e. we use with each gaussian bin having center z i and width σ i z .Since we are concerned with a theoretical analysis, we can parametrize the window function in terms of the redshift of the source exploiting the background relation D i (z i ) of Eq. (3.4).This way, we can study how our expected signal changes as a function of redshift, providing a direct comparison with the number count of galaxies.Moreover, using the exact same window to calculate angular correlations for GW/SN and galaxies will facilitate the comparison of the two.In practice, of course, the GW redshift will not be available, and in a forecast analysis we should use a window parametrized in terms of D, as in Eq. (4.3).For the main part of our analysis, we will use σ i z = 0.2, unless specified otherwise.We will analyze further the implications of this choice in Sec 4.3.
In Fig. 1, we plot the auto-correlations of the number counts signal for GW, SN, and galaxies, at different redshifts.The dark solid line in the top of each plot shows the total signal, while the other lines correspond to the auto-correlations of each perturbative component of the signal, i.e. the auto-correlations of the number counts obtained switching off all perturbative terms in Eq. (3.13) except the one under consideration.We choose to show the results directly in GBD instead of ΛCDM, because in this model we have a non-vanishing α M , which potentially allows us to show a difference between the number count of GW and SN, and additionally sources the scalar field clustering contribution to the total signal of GW.We find, however, that the deviations from the ΛCDM curves are small, as we will show in more detail below.So small in fact that , this plot can be considered also representative of the ΛCDM scenario (with the only difference that there would be no scalar field curve in the left column panels).
From the Figure, one can see that the total signal for GW and SN at low and intermediate redshifts (z = 0.3 − 1.0) is comparable to the galaxy signal, with the galaxy one being only the number count signal obtained switching off all other perturbations.The dark blue solid curve shows instead the total signal.We use GBD to show the possible difference between GW and SN, and the individual contribution of the scalar field clustering.Comparing with the ΛCDM case, we see no difference in these curves appreciable by eye, except the absence of the scalar field curve in left column panels, thus this plot can also be used as representative of ΛCDM.
slightly higher at ℓ < 25.For these redshifts, the density contrast dominates the signal, except for a small contribution from the lensing and velocity terms at very low multipoles.For all sources, the lensing auto-correlation becomes more relevant as the redshift increases; this is the case also for the other integrated effects, i.e. the time delay for which the auto-correlation changes significantly at different redshifts, and the ISW that grows mildly in redshift.We notice that, comparing our results with [22], in our case the lensing auto-correlation grows way slower in redshift, remaining about 1 − 2 orders of magnitude lower than the density contrast correlation even at z = 3.0.We find that that is imputable to the sign difference we spotted in Eq. (2.11), propagated on the final amplitude of the lensing convergence term.We caution that this also impacts the total number count correlation in possibly a non-negligible way, as consequently the lensing convergence will carry a higher weight on the total signal.As for the local effects, the Doppler auto-correlation term decreases with increasing redshifts, consistent with the impact of the peculiar velocity becoming less and less relevant, while the potential terms auto-correlation remain essentially the same at different redshifts.The DSD correlation of GW and SN, similarly to the galaxy case, increase at increasing redshifts.Finally, we see that the scalar field clustering contributes very little to the GW number counts, being several orders of magnitude lower than the signal from any other perturbation.Though rich in information, Fig. 1 still makes it difficult to quantify how much the different perturbative terms contribute to the total signal, as it shows only the auto-correlation of the individual perturbations.However, also their cross-correlation with all other perturbations can carry significant weight.We quantify more accurately the contribution of each perturbation in Fig. 2, where we plot it as a fraction of the total signal, at different redshifts from 0.5 to 3.0.Specifically, for each contribution to the number count, we plot the quantity where C ∆ ℓ is the total auto-correlation of the number count, while Ĉℓ is the same auto-correlation, but obtained removing from the number count the perturbation under analysis.For example, in the case of DSD Ĉℓ would be the correlation of We group together the scalar field, ISW, time delay and potentials terms since their contributions are significantly smaller than the others, as can be anticipated from Fig. 1.
We explore the fractional contributions in ΛCDM, and compare the GW and galaxies case.We can see that at low redshifts, the number count signal is strongly dominated by the density contrast.At higher redshifts, the other components start to become important, especially on larger angular scales.In the case of GW, at ℓ ≥ 100 the second most important component is always lensing, which contributes about 6% of the signal at low redshifts, growing up to 10 − 20% at z = 1.5 and z = 3.0.The other most important contributors are the DSD, which however are only relevant on larger scales (ℓ ≤ 100).At low z, the DSD contribution remains below ∼ 5%, while it weights up to 50−70% on the total signal at higher z and low multipoles.Finally, the Doppler term always remains subdominant with respect to lensing and DSD, contributing no more than 3% of the signal, at all redshifts.Similarly, all other perturbations jointly only contribute, at best, a few % of the signal, and only for ℓ ≤ 20.We note that Fig. 2, aside from giving a visualization of the dominant components of the number count, can also be used to quantify the error committed in neglecting some of the less relevant perturbative terms in a cosmological analysis.Neglecting terms can sometimes be useful, as it can speed up the numerical calculation of the angular correlations, but it comes at the price of loosing accuracy.Whether or not this is acceptable, needs to be decided on a case-by-case basis, depending on the specific GW detector targeted, which also determines the noise associated to the correlation, as we will discuss in Sec. 5. Generally speaking, we conclude from Fig. 2 that in any accurate analysis, the density, DSD and lensing perturbations should all always be included.The other perturbative terms can sometimes be safely neglected, depending on the required accuracy.While Fig. 2 refers to ΛCDM, we have performed the same analysis for all the MG models listed in Tab. 1, and found negligible differences.
Comparing the number counts in D-space (GWs) with those in z-space (galaxies), we observe similar trends, with some crucial differences involving primarily lensing.We find that the net lensing convergence contribution is always positive in D-space at the scales considered, while in z-space the lensing curve at z = 1.5 presents the well known V-shaped (divergent) feature, which happens when the fractional contribution flips sign.The DSD (RSD) contribution, on the contrary, is always positive in both spaces, and at all redshifts has a similar impact on the signal for both GWs and galaxies; the remaining perturbations weights remain small even for the galaxies, rarely raising above 1%, with another relevant difference in the Doppler term, whose impact at intermediate and high redshifts is slightly larger in D-space.Relative contributions to the total signal (same as Fig. 2) of the density contrast, lensing and DSD.We plot these contributions at different redshifts and for different choices of the tomographic bin width σ i z .We notice that the binning choice has an important impact in determining the relative weight of each perturbation within the number count, with density contrast and DSD dominating thinner bins, while lensing becomes more and more relevant for wider bins, as both the density contrast and the DSD average out faster.

Impact of the binning choice
Fixing the free parameters in Eq. (4.14) corresponds to making a tomographic binning choice, with z i being the center of the bin, and σ i z its width.In a forecast analysis these parameters can be chosen with a certain degree of arbitrariness, though with the caution that the lower value of σ i z (or the corresponding σ i D ) is limited by the luminosity distance measurement uncertainty.This will, in general, depend on the specific GW/SN survey considered, though it can impact the results significantly.In general, wider bins have the effect of smoothing out the number count correlation.However, they can also impact differently the individual components of the signal.
We illustrate the effects of the binning choice in Fig. 3, where we plot again the fractional contributions to the total correlation, this time focusing on the dominant terms, density, lensing and DSD, and for different choices of the bin width.We can appreciate that, for any binning choice, the density contrast is always the dominant ingredient of the clustering signal, with the only exception already pointed out for Fig. 2 of high redshifts and low multipoles.However, in the case of fine binning, the impact of the DSD on the total correlation becomes more and more important, coming to dominate over the lensing contribution at all scales at high z, and at ℓ ≤ 50, 300 for z = 0.5, 1.5.On the contrary, for wider bins the DSD smooth out extremely fast, letting lensing become the second most important component at all scales above ℓ ≥ 10.For this binning choice, lensing becomes particularly relevant also at high redshifts and very low scales, where its relative weight is comparable to that of the DSD and of the density contrast.
As a final remark, we stress that the quantities plotted in Fig. 3 are relative weights.Thus, the fact that the lensing curves increase for larger bins must not be deceiving: all correlations, and lensing in particular, are partially averaged out over large bins, i.e. all absolute signals are decreasing.What is changing is their relative importance within the total number count correlation, with lensing becoming more and more important as DSD and density contrast average out faster.
We conclude that, in general, the binning choice has a huge impact on the individual components, and this should be taken into account carefully in forecast analyses.We note that the reference choice of σ i z = 0.2 that we used for most of the paper is compatible with a luminosity distance measurement uncertainty of ∼ 10% or less.This is consistent with several future observational scenarios, but it might be too narrow for some of the more conservative setups (see discussion in Sec. 5).

GW number count in ST theories: deviations from ΛCDM
To better highlight the effects of MG on the GW number count, in Fig. 4 we plot the deviations from GR of the angular correlations in units of cosmic variance (referred to GR), i.e.
for the models described in sub-section 4.1 .Any deviation higher than 1 σ GR ℓ could potentially be observed, provided that the right multipoles are accessible for observation, and that the measurement error is small enough.We highlight with black horizontal lines the corresponding limits, i.e. ∆C ℓ /σ GR ℓ = ±1, in the Figure .For wCDM models, we keep w a fixed at −0.05, while we vary w 0 in the range [−0.825, −0.95].These values are all compatible within 1σ with the latest DES Year 3 bounds and DES+SN joint constraints, and within 2σ of the joint DES+BAO+BBN and DES+Planck bounds [68].We see that a simple modification of the dark energy equation of state does not induce big deviations in the correlation signal with respect to the GR case.Both at low and at high redshifts, these deviations almost always remain within the cosmic variance limit.Only deviations induced by w 0 ≤ −0.85 could be detectable at low redshifts, though the measurement requires a high level of accuracy and the access to multipoles higher than 400.Exploring the GBD case, we find that this is one class of models that, in general, induces the highest deviations in the correlation curves.We keep w 0 and w a fixed to their fiducial values, while we vary the value of the non minimal coupling today, Ω 0 .We observe that for any value of Ω 0 > 0.01, deviations from ΛCDM are potentially detectable, provided that the highest multipoles are available.In particular, a value of Ω 0 ≥ 0.02 generates deviations from GR that can be detectable both at z = 0.5 and z = 1.5, for ℓ ≥ 300 and ℓ ≥ 400 respectively, while the deviation signal itself can reach up to 5 σ GR ℓ .All values of Ω 0 explored here are consistent within 3σ with the bounds of [67] which are, to the best of our knowledge, the most recent constraints put on the specific parametrization of the EFT functions we adopted in this work.In particular, Ω 0 = 0.02 is closest to their best fit value for the GBD model, Ω 0 = 0.018.
For KGB models, we keep fixed w 0 , w a , Ω 0 and γ 0 1 .Past works (see e.g.[67]) have shown that LSS and CMB depend very mildly on γ 1 (η) and we have no reason to believe that GW would be more sensitive to this EFT function5 , (in fact, we checked that that is the case).We therefore vary only on γ 2,0 , in the interval [−0.6, −0.1] which is compatible with the bounds of [67].We find small deviations from the GR case.The deviations are particularly difficult to capture at low redshifts, with only values of γ 2,0 ≤ −0.5(0.6)give a detectable signal for multipoles higher than 600(800).They become slightly more significant at z = 1.5, giving a detectable signal for γ 2,0 ≤ −0.4,provided that ℓ ≥ 500 are observable.1.
Finally, for the f (R) family, we yet again fix the background and vary only B 0 .Like GBD, f (R) models give potentially detectable deviations from GR in the GW number counts, some of which already higher than 1σ ℓ for ℓ ≥ 400.However, the relevance of these deviations is limited to low redshifts, while at z = 1.5 the signal is essentially indistinguishable from ΛCDM as the MG effects quickly fade.We also notice that f (R) and two of the KGB scenarios are the only cases in which the modifications of gravity actually increase the GW clustering signal.In all other cases explored here, ∆C ℓ is negative, which means that the MG signal is lower than the fiducial GR one.
Let us note that the deviations shown in Fig. 4 correspond to values of the EFT parameters that are still loosely compatible with the observational bounds placed in [50,67,69] using CMB data and state of the art galaxy surveys.As such, they are models that in general depart very little from GR.If we were to explore the GW signal separately, as an independent cosmological probe, in principle we could loosen the priors on these parameters and possibly explore higher departures from GR.

Synergy between GW and EM tracers
As already discussed, photons and GWs respond differently to the underlying theory of gravity, with GW geodesics being explicitly affected by the non-minimal coupling and perturbations in the scalar field.This introduces an interesting difference between the number count in luminosity distance space of SN and that of GWs in the case of ST theories.This comes in addition to the differences due to working in luminosity distance space as opposed to the redshift space, commonly used for galaxies.In this Section we perform a preliminary exploration of the joint constraining power of these three tracers.
In Fig. 5, we plot the difference between the GW number count angular correlations and the analogous EM ones, where EM stands for either galaxies (GAL) or supernovae.We focus on ΛCDM and the KGB model, which among the different ST scenarios that we explore is the one allowing for larger values of Ω 0 (and α M ).This is particularly useful for the specific comparison we perform in this Figure, since a bigger α M implies a bigger ∆C GW−SN ℓ .We can appreciate the potentialities of a joint galaxies-GW analysis, as ∆C GW−GAL ℓ reaches up to 20% at lower multipoles, plateauing at high ℓ around a fixed value between 5 and 20% for all redshifts below z = 3.0, and for both models.On the contrary, the difference between GW and SN auto-correlations remains very small, reaching peaks of up to 0.3% at lower multipoles, and settling around 0.05−0.1% for ℓ ≥ 100.Whether these differences are big enough to add valuable cosmological information is better shown in the bottom row of Fig. 5, where we plot the same quantity ∆C GW−EM ℓ , but in units of cosmic variance σ EM ℓ .We see that in the case of SN, ∆C ℓ is only a small fraction of the cosmic variance, indicating that the cross-correlations between GW and SN is unlikely to add any insight about the specific MG model -though we note that a joint analysis would still be worthwhile, as their cross-correlation adds more strength to the clustering signal they both probe independently.On the contrary, the galaxy signal is clearly distinguished from the GW one for all redshifts below z = 3.0, differing of up to 5 − 6 σ GAL ℓ at high redshifts and high multipoles, hinting that their cross correlation is likely rich of extra cosmological information.In general, we conclude again that the access to small scales is key, as it is easier to distinguish between GW and EM signals at low z and high multipoles, where the cosmic variance itself is much smaller.We also note, however, that at z ∼ 1.0 it is already possible to distinguish the two signals for ℓ ≥ 30.
The fact that ∆C GW−SN ℓ is always lower than the C SN ℓ cosmic variance means that, as far as the number count cross-correlations are concerned, these two sources give, essentially, the same signal.However, for the sake of completeness, we mention that there is another possibility that concerns these number counts, i.e. the definition of a new observable We plot the difference between the GW and the galaxy signal in ΛCDM (left panel) and in KGB (central panel), and the difference between GW and SN auto-correlations in KGB (right panel).Bottom row: same as above, but this time the difference between the auto-correlations ∆C ℓ is plotted in units of cosmic variance σ EM ℓ and the black horizontal lines mark ∆C ℓ = ± σ EM ℓ .This shows that, even if the percentage deviations between the GW and the EM signal is slightly higher at lower multipoles, it is actually easier to capture it at higher ℓ, where the cosmic variance noise is smaller.

reads
and This quantity essentially isolates the explicit MG contribution to the GW number count.
In ΛCDM, ∆GW−SN (D, n) is exactly zero, thus the sheer measurement of a non-vanishing value can in principle serve as a smoking gun of physics beyond the standard cosmological model.The practical detection of this quantity might, however, be subject to significant observational challenges.Among these, the request that the number counts of GW and SN need to be observable in the same sky areas down to sufficiently small scales, which might demand unrealistic requirements in terms of numbers, locations and overlap of the two sources (see e.g.[43]).Further explorations in this sense are beyond the scopes of the present work.

Observability of the number count in GW detectors
It is interesting to build on the modeling of the previous Sections, to establish what are the observational requirements to carry out number count studies.In this Section, we perform a first step towards this, investigating the possibility of detecting correlations in the number count of compact objects above some noise, for different benchmark observational configurations.
There are several limiting factors affecting the actual detectability of this signal.The experimental error on the luminosity distance determination from the GW amplitude forces a lower limit on the width of the tomographic bins.It thus increases the noise of the measured correlations, while weakening the signal itself, as the perturbative effects partially average out over larger bins (see discussion in Sec.4.3).On top of this, the total number of detected sources plays a key role, with denser catalogues allowing for a stronger measured correlation.Finally, a severely limiting factor is the uncertainty in the sky-localization of the sources, which reduces the accessible scales cutting out the higher multipoles, where most of the cosmological information is contained.
To determine whether the number count of Eq. (3.13) will be observable in future GW detectors, we define a cumulative signal-to-noise ratio (SNR) as The noise N i ℓ takes the form of a shot noise dependent on the luminosity distance of the source, and can be written as (see e.g.[16]) Table 2. Specifications for the benchmark observational scenarios explored in this work.In particular, we set the average fractional uncertainty on the luminosity distance measurement σ GW (D i )/D i to be used for the shot noise calculation, assumed constant at all redshifts; the number of compact objects detected in each tomographic bin N i GW ; the maximum multipole accessible ℓ max ; and the width σ bin D of the luminosity distance tomographic bin.Specifications are chosen to be loosely compatible with expected observational setups for 3rd generation ground based detectors, a LISA-like space based interferometer and a space based Deci-Hz observatory, in their pessimistic and optimistic configurations.
Here, N i GW is the total number of sources in the i th tomographic bin, while f sky is the sky fraction covered by the survey; for future GW surveys we can set f sky = 1, as the proposed detectors are designed to have essentially no blind spots.σ i GW instead is rigorously defined through where W i (D) is the bin window function and σ GW (D) is the average observational uncertainty on the luminosity distance measurement for a GW source at distance D, such that D obs = D + σ GW (D).The last equality in Eq. (5.3) follows from approximating the fractional error σ GW (D)/D constant over the bin.
The values for the observational uncertainties in Eq. (5.1) need to be chosen relying on forecasts for the next generation detector performances.We consider three different classes of detectors: the 3rd generation (3G) ground based interferometers, a LISA-like space based interferometer and the deci-Hertz space based observatories.The next generation of ground based detectors, such as the Einstein Telescope (ET, [70]) and Cosmic Explorer (CE, [71]), will be most sensitive to the high frequency band (∼ 1 − 10 4 Hz), where there is a richness of inspiraling and merging stellar mass binaries (SMB).Based on the current rate of events measured by the LIGO-Virgo-Kagra network [72,73], forecasts indicate that the number of sources that could be detected in 10 years of data taking by ET is as high as O(10 5 ) for binary black holes (BBH) and neutron star-black hole binaries (NSBH), and O(10 6 ) for binary neutron stars (BNS) [70,74,75], while the number and quality of the measurements improve hugely if the observation is carried out jointly by a network of ET and one or two CE [75][76][77][78][79].We select two scenarios for the 3G detectors: a conservative scenario in which we assume a measurement error on the luminosity distance of 30% and require at least 10 4 sources in each tomographic bin; and an optimistic scenario, with σ GW (D)/D = 5% and N i GW = 5 • 10 4 .We stress that in both cases, the number of sources that we put in each redshift bin is 1 or even 2 orders of magnitude lower than the total number of expected detections, i.e. we are implicitly considering the (realistic) case in which the requested accuracy on the luminosity distance is met only by a fraction of the detected sources.
Next generation space-based detectors such as the Laser Interferometer Space Antenna (LISA, [80]), Taiji [81] and TianQin [82] will be sensitive to low frequency GW (∼ 10 −4 − 0.1 Hz).Their extra-galactic targets will primarily be massive black holes binaries (MBHB) and extreme mass-ratio inspirals (EMRIs) systems.Though the merger rates for these sources are still highly uncertain, even in the most promising astrophysical scenarios LISA is not expected to measure more than a few hundreds MBHB over the whole 4 years mission duration [83].On the contrary, EMRIs are potentially much more numerous, with rates ranging between ∼ O(10 2−3 ) detections per year, depending on the astrophysical model [83,84].LISA will provide excellent luminosity distance measurements for these sources, thus for LISA-like detectors we consider an optimistic scenario, where we choose σ GW (D)/D = 1% and a more conservative scenario, with σ GW (D)/D = 10%.In both cases, we impose N i GW = 10 3 .As such, both LISA-like scenarios are to be regarded as upper limits, as they require, to begin with, a high number of detections.We note that if these numbers are not met, it is unlikely that LISA detections can be used efficiently for number counts studies.We further note that the LISA cases can also be representative of two corresponding 3G scenarios, in which a few thousands of golden SBH are detected with extremely accurate luminosity distances.
Further away in time, the frequecy gap between LISA and 3G detectors will be covered by space based deci-Hertz interferometers.This band contains, among other signals, the inspirals of millions of SBH.Proposed detectors operating in these frequency range are the Advanced Laser Interferometer Antenna (ALIA, [85,86]), the Big Bang Observer (BBO, [87]) and the DECi-hertz Interferometer Gravitational wave Observatory (DECIGO [88,89]).While the actual technical feasibility of these missions is under investigation, a slightly less ambitious family of detectors in the deci-Hertz band, generically named Deci-hertz Observatories (DO) has been proposed in [90].Given the early planning stages of these experiments, solid forecasts on the measurement accuracy are unavailable, but based on preliminary estimates (see e.g.[87]), they are likely to outperform both LISA and 3G observatories.Tentatively, for this family of detectors we again consider a conservative and an optimistic scenario, requiring 5 • 10 4 and 10 5 detections in each bin, and σ GW (D)/D = 1% and 0.5% respectively.
We summarize the assumed specifications for each family of detectors in Tab. 2. While we vary the accuracy on the luminosity distance measurement to be used in Eq. (5.2) for the shot noise calculation, we keep fixed the relative gaussian bin width to σ bin D /D i = 10%, with the only exception of the 3G conservative scenario, where the higher uncertainty on the distance measurements imposes wider bins.This choice is made consistently with the requirement of a high number of sources within each bin, while assuming constant ratios for σ i GW /D i and σ bin D /D i automatically accounts for a redshift evolution in the luminosity distance uncertainty, and consequently of the bin width.
In Fig. 6 we plot the cumulative SNR for the different scenarios described above, and for different redshifts, clearly marking the cumulative SNR of 1 and 5.We plot all 3G and LISA curves up to ℓ max = 180.This value corresponds, roughly, to an average sky-localization of 1 deg 2 for each source.Better localizations are extremely unlikely to be achieved with these detectors, especially for a high number of sources, and indeed though the 1 deg 2 limit might be met with BBH observations [78], BNS detection are expected to be way more uncertain, with only ∼ 20% of detected source localizaed within 10 deg 2 [75], roughly corresponding to ℓ max = 50.The DO family instead is way more likely to achieve optimal sky-localization.Expectations depend massively on the number of interferometers that can be put in a network in space, as more detectors allow a better triangulation of the GW.However, it is expected that a DO-like detector alone could measure signals up to ℓ max ∼ 300 [90], while a network -25 - of DO could localize each source within a few arcmin 2 [87], i.e. ℓ max ≥ 1000.From Fig. 6, we can see that for both LISA and 3G configurations, only the optimistic case produces a number count correlation signal with SNR ≥ 5.For the conservative scenarios, a correlation is detectable only at low redshifts (z ≤ 0.5) and with very low significance (SNR ≤ 1.5).The optimistic configurations instead produce angular correlations that are detected almost at all multipoles, even with a 5σ significance for the smallest redshifts.However, for redshifts z = 2 and higher, detecting a correlation signal requires access to the smallest scales, ℓ ≥ 70 for 3G detectors and ℓ ≥ 60 for LISA.Going up to ℓ = 100, however, could produce a cumulative SNR of 50 − 70 for the smallest redshifts.This is once again proof that optimal source localization will be key to determine the feasibility of number count studies with the nest generation of detectors.For the DO case instead, both explored scenarios give optimal SNR: the signal is detectable with 5σ significance at all multipoles, and even at multipoles around ℓ ∼ 100, the cumulative SNR is higher than 100 at essentially all redshifts.A network of space based DO could probe scales up to ℓ ∼ 1000, accumulating an SNR of several hundreds.
As a final remark, we stress that there are additional sources of errors that we did not account for here, among which a redshift dependent GW source distribution and the completeness of the detected GW catalogues.These can impact especially the high redshifts, where sources are fewer and more difficult to detect.We partially compensate for the former using wider bins at higher redshifts, although a more accurate determination of its influence might further lower the expected SNR at z ≥ 2.0.However, for lower redshifts we do not expect these effects to have a significant impact on our conclusions.

Summary and conclusions
We have derived the theoretical expression for the number count of objects living in luminosity distance space, in the broad context of scalar-tensor theories of gravity with luminal propagation of the tensors, including the full set of relativistic effects.We have shown that, contrary to the GR case, the theoretical form of the number count differs if the objects considered are gravitational waves (GW) or Supernovae (SN).This, because the propagation of GW waves is perturbed differently than that of photons in the presence of a non-minimal coupling with the scalar field.
We have investigated the impact of relativistic effects on the number count in luminosity distance, for both GW and SN, by mean of the 2-point angular correlation function associated to it.Similarly to the galaxy case, we find that the effect that dominate the signal at essentially all redshifts and scales is the density contrast -which is probed directly from the number count.Following suit, the two most important effects on top of that are the weak lensing convergence and the luminosity distance space distortions (DSD).For the latter, we have explored how their contribution to the total signal depends on the choice for the tomographic bins.We have also shown that all other perturbations remain sub-dominant with respect to lensing and DSD, with Doppler contributing at most O(1%) of the correlation signal, similarly to all potential terms, time delay, ISW and the scalar field clustering combined.
Focusing on scalar-tensor models, we have explored the deviations of the GW number count angular correlations from the GR counterpart.With the spirit of investigating the potential for detectability, we have plotted these deviations in units of cosmic variance σ GR ℓ .While departures from ΛCDM are small for all models, the theories that give the most promising signatures are GBD and f (R), for which they can reach up to 6 σ GR ℓ , especially at low redshifts.On the contrary, for the class of KGB and wCDM models departures from GR are tiny, and barely above 1 σ GR ℓ only for the most extreme values still allowed for the MG parameters.
We have also investigated the synergy of GW with other EM tracers, either living in redshift space (galaxies) or in luminosity distance space (SN).In ΛCDM and in alternative theories, we find that the number count of GW and galaxies give two very well distinguished signals for z < 3.0, differing by up to 25% and up to 5 − 6 σ GAL ℓ ; therefore, their crosscorrelation can add valuable information in cosmological analyses.On the other hand, despite the theoretical differences in the GW and SN number counts, the respective auto-correlations differ, at most, by fractions of the SN cosmic variance σ SN ℓ .However, we show that it is still possible to additionally define a new observable, in terms of the difference between the SN and GW number count.This quantity isolates the MG contributions to the number count, while being zero in ΛCDM, thus its detection could serve in principle as a smoking gun of exotic physics.We give the theoretical expression for this observable, though we strongly stress that the actual feasibility of its detection is yet to be proven.
We have also explored the potential observability of the GW number count in planned ground-based and space-based observatories, concluding that it can be detected with the next generation of interferometers only in the optimal scenarios, especially at high redshifts.Crucial for the observation are a high statistic of detected mergers, a low uncertainty on the luminosity distance determination and the access to the lowest scales, i.e. a good angular localization of the sources.We find the former factor to be the biggest limitation for a LISAlike detector, while the last two affect more heavily the ground-based ET and CE.However, for the DO family, we find that the GW number count is detectable at all redshifts, with SNR that can reach several hundreds, likely allowing for optimal cosmological analyses.We conclude that, in the long run, GW have the potentialities to become an extremely effective probe to test the full extent of the cosmological model.
As a final remark, we highlight that there are further observational uncertainties that we have not considered in this work.Likely the most impactful of them is the GW magnification bias, for which we only provided the theoretical expression in ST theories.We have shown that for GW this bias is model-dependent, as the number of sources that are detected above a certain SNR threshold depends on the background cosmology and on the value of the nonminimal coupling.We deferred the full modeling of these effects to a future investigation.However, we caution that, given the importance that the magnification bias has for galaxy surveys (see e.g.[91,92]), we expect it to be as important for GW, and to be an essential ingredient in forecast analyses.
resulting equations are presented in a more general form than done in [28], with the intent of generalizing the result to cosmological models beyond ΛCDM.

A.1 Density perturbations
We begin with the density perturbations The first term on the RHS is simply the GW clustering in Newtonian gauge.For the second term we use that ρ where n( D) is the physical comoving number density of sources, and we obtain The term ∂ ln n/∂ ln a is the evolution bias of GW, totally analogous to the same term affecting the clustering of galaxies.

A.2 Volume densities
Before proceeding to the calculation of the volume perturbations, a very useful quantity that will come in handy is the volume density.This can be defined starting from the infinitesimal volume element Then, recalling the definition of our background and perturbed frames in Sec. 2, (θ s , ϕ s ) always correspond to the physical location of the source ( θ, φ), while (θ o , ϕ o ) depend on the frame in which the volume is considered.
We can compute the background volume density v( D, θ, φ) using Eq.(A.6), and noting that, at the background level, √ −g = a 4 , θ o = θ, ϕ o = φ, u = 1/a(1, ⃗ 0) and xi = χn i .Furthermore, we can write while the partial derivatives of the comoving coordinates w.r.t. the angles are (A.9) Because u i = 0 and ∂ x3 /∂ φ = 0, in Eq. (A.6)only four terms survive, namely As expected, the physical volume density does not depend on the source sky-location, only on its distance D from the observer.
The calculation of the volume density in the perturbed space is similar to the one in physical space, but with some crucial differences.In the perturbed space we have that √ −g = a 4 (1 + Ψ − 3Φ), u = 1/a(1 − Ψ, v), θ o = θ, ϕ o = ϕ and the angles are connected, to linear order, by the relation which means that the determinant of the Jacobian of the angles transformation now reads From eq. (A.6), we can split the perturbed volume density in the sum of two terms the calculation of the first term is in all analogous to that of v( D), except now we have a number of prefactors coming from u 0 , the Jacobian of the angular matrix and √ −g.We also need to be cautions that all quantities that enter the calculation are perturbed (unbarred) quantities, except the source angles, that are connected to the observed angles through (A.12).
Retracing the same steps, we get where, in the second term, we have used that Ψ is already a first order perturbation, and we neglected higher order terms.
As for the second term, v(D, n)| µ=i , considering that ∂η/∂ θ = ∂η/∂ φ = 0 and that, once again, ∂x 3 /∂ φ = 0, there are still only 4 terms that survive the contraction of ϵ iναβ with the coordinates derivatives.These terms are Now, v is already a first order perturbation and as such, since it enters all terms, we can keep only the zeroth order terms in the prefactors and drop the rest.We obtain Joining Eq. (A.15) and (A.18), we are left with the final expression for the perturbed volume density

A.3 Volume perturbations
The definition (A.5) allows to rewrite the volume perturbations in terms of perturbations to the volume density.We can write (A.20) sin θ sin θ + cos θδθ = 1 + cot θδθ , (A.21) and that χ 2 = ( χ + δχ) 2 ≃ χ2 + 2 χδχ.We also have that Finally, we expand that we re-write as where we take the derivatives of the comoving distance perturbations δχ and of the luminosity distance perturbations ∆D along the null geodetics.Plugging Eqs.(A.21)-(A.24)into Eq.(A.19), we finally obtain Then, we need to compute the second term at the RHS of Eq. (A.20), that yields We can further explicitate this expression, by using that the perturbation on the comoving distance is given by (see [8,28]) Combining the results for both the density and volume perturbations, we ultimately obtain the general expression (2.11) for the number count in luminosity distance space.

B Numerical implementation
For the numerical implementation in EFTCAMB, we want to express the 2-point correlation function of Eq. (4.6) in a convenient form.We do so following closely the original CAMB Source paper [27].We define a source object S ∆ (k) as (B.4)Each source term can be derived by performing the expansion (4.4) on each term at the RHS of Eq. (3.13) and factoring out the integration over time.

B.1 Derivation of the source terms
• We begin by deriving the source terms corresponding to the perturbations at the source position and that do not contain spatial derivatives.These are S i δ , S i Φ , S i Φ, S i Ψ and S i SF .For a generic term of the form • We then have perturbation terms that are integrated over the GW path, such as the ISW and the Shapiro time-delay terms.These only introduce a small extra complication to the calculation.The terms are in the generic form where W i (χ) ≡ dD/dχ W i (D).Before expanding in Fourier-modes, we want to manipulate these terms changing the order of integration.Doing so, transforms ∆ ŝ as where we also have renamed the silent integration variables, switching primed and unprimed η.Now, ŝ is in the proper form to perform the usual expansion in spherical harmonics, which then allows us to read off Once integrated by parts, that gives S i ∇Φ (k) (Eq.(B.22h)).
• Finally, we look at the lensing therm.First, we perform the same trick implemented above for the ISW and time delay terms, and exchange the order of integration, obtaining For this term we will work, as it is customary, in the flat sky approximation.Then, instead of the expansion of Eq. (4.4), we take directly the 2D Fourier transform of the 2D laplacian, that simply brings down a factor −ℓ(ℓ + 1), giving

Figure 1 .
Figure1.Angular correlations of the number counts of GW (left), SN (center) and galaxies (right) in the fiducial GBD model.Each curve shows the auto-correlations of a specific relativistic effect, i.e. the number count signal obtained switching off all other perturbations.The dark blue solid curve shows instead the total signal.We use GBD to show the possible difference between GW and SN, and the individual contribution of the scalar field clustering.Comparing with the ΛCDM case, we see no difference in these curves appreciable by eye, except the absence of the scalar field curve in left column panels, thus this plot can also be used as representative of ΛCDM.

Figure 2 .
Figure 2. Contribution of the individual perturbation terms to the number count signal in ΛCDM, as a fraction of the total signal.The fractional contribution is computed as |C ℓ − Ĉℓ |/C ℓ , where Ĉ is the number count signal obtained including all perturbations effects except the selected one.Lines of different color refer to different perturbation terms (with the yellow lines indicating the contribution of all term except density contrast, lensing, DSD/RSD and Doppler), while the dashed line indicates where the lensing contribution has negative sign.The panels in the first row refer to GW, while in the bottom row are the analogous plots for the galaxies.Finally, different panels correspond to different redshifts, as indicated, in the range [0.5, 3.0].

Figure 3 .
Figure3.Relative contributions to the total signal (same as Fig.2) of the density contrast, lensing and DSD.We plot these contributions at different redshifts and for different choices of the tomographic bin width σ i z .We notice that the binning choice has an important impact in determining the relative weight of each perturbation within the number count, with density contrast and DSD dominating thinner bins, while lensing becomes more and more relevant for wider bins, as both the density contrast and the DSD average out faster.

Figure 4 .
Figure 4. Deviation from ΛCDM of the auto-correlation of GW clustering for different MG theories, in units of cosmic variance.The theories explored are, from top to bottom, wCDM, GBD, KGB and f(R).The different curves are obtained varying the theories parameters.Respectively, we varied: w 0 in the interval [−0.95, −0.825] in the first row, Ω 0 in [0.01, 0.06] in the second row, γ 2 in the interval [−0.6, −0.1] in the third row and B 0 in [10 −6 , 10 −5 ] in the bottom row.The remaining parameters for each theory were left fixed at the fiducial values of Table1.

(4. 16 )Figure 5 .
Figure 5. Top row: Relative difference between GW and EM sources number count auto-correlations (C GW ℓ − C EM ℓ )/C EM ℓ , with EM = [GAL, SN], at different redshifts in the range [0.5, 3.0].We plot the difference between the GW and the galaxy signal in ΛCDM (left panel) and in KGB (central panel), and the difference between GW and SN auto-correlations in KGB (right panel).Bottom row: same as above, but this time the difference between the auto-correlations ∆C ℓ is plotted in units of cosmic variance σ EM

Figure 6 .
Figure 6.Cumulative SNR for the GW number count correlation, at different redshifts.The observational scenarios considered in terms of accuracy of the luminosity distance measure and number of detected sources are those of Tab. 2, specifically: a network of ground based 3G detectors (first row), a LISA like space based interferometer (second row) and a DO space based observatory (bottom row).Gaussian bins widths are fixed at σ bin /D i = 10% (30% for the 3G -Conservative scenario).The dash-dotted and solid horizontal lines in each panel indicate a cumulative SNR of 1 and 5 respectively.

Table 1 .
[62]cial parameter values for the cosmological models considered in this work.The background parameters h, Ω m,0 , Ω b,0 , n s and σ 8 were kept fixed at Planck18[62]best fit values.