Non-Gaussianity from the cross-correlation of the astrophysical Gravitational Wave Background and the Cosmic Microwave Background

Since the first LIGO/Virgo detection, Gravitational Waves (GWs) have been very promising as a new complementary probe to understand our Universe. One of the next challenges of GW search is the detection and characterization of the Stochastic Gravitational Wave Background (SGWB), that is expected to open a window on the very early Universe (cosmological background) and to provide us new information on astrophysical source populations (astrophysical background). One way to characterize the SGWB and to extract information about its origin is through the cross-correlation with other cosmological probes. To this aim, in this paper, we explore the cross-correlation between the astrophysical background anisotropies and the Cosmic Microwave Background (CMB) ones. Such a signal is sensitive to primordial non-Gaussianity (nG) through the GW bias. Thus, we study the capability of next generation space-based interferometers to detect such a cross-correlation signal and to constrain primordial nG.


Introduction
One of the future goals for Gravitational Wave (GW) astronomy is the detection of the Stochastic Background of Gravitational Waves [1][2][3][4].After three runs by the LIGO-Virgo-KAGRA (LVK) collaboration we are very close to the detection of the stochastic background produced by astrophysical objects [5,6].On the other hand a claim of a possible detection has been already done by the NANOGrav collaboration [5], which in the 12.5 years data has found the presence of a common red-noise process [7].Current operative ground-based interferometers should reach their design sensitivity in a few years and this will be probably crucial for the detection [5].
Two different contributions are expected to contribute to the background: an astrophysical component (AGWB), originated from the superposition of a large number of unresolved gravitational wave signals coming from very far or very faint astrophysical sources [8][9][10][11][12][13], and a cosmological one (CGWB), coming from the very early stages of the Universe, generated by mechanisms as inflation, first-order phase transition, or scalar-induced GWs [14][15][16][17].The former will provide us information about population properties of astrophysical sources, while the latter will give us access to the early Universe's stages.Since we are close to the detection of a SGWB, it is important to understand to which extent such signal can be used to extract information about our Universe, eventually also in combination with other cosmological probes.Hence, in this paper, we try to achieve a challenging task, aiming to use the AGWB to constrain primordial nG, typically parametrized in terms of the parameter f NL [18].More specifically we focus on primordial nG of the local type (e.g., see [19][20][21]).Typically, measurements of nG from Large-Scale Structure (LSS) have been performed by means of electromagnetic signals (e.g., photons emitted by galaxies), while in this work we consider a different signal, i.e.GWs.
Although the isotropic component of the AGWB (i.e. the monopole) already contains plenty of astrophysical information [1,4,14,22], in our study we focus on its anisotropic contribution.Such anisotropies are imprinted both at the GW production time and during the propagation of the GW signal [23][24][25][26][27].They share many analogies with CMB ones [28][29][30] and they are also generated by the same common seeds (i.e., gravitational potentials) 1 .
Since they are affected by the same large-scale (scalar and tensor) perturbations and that gravitons propagate on the same geodesics of CMB photons induces, for-free, a crosscorrelation signal.Despite the two signals are generated at different epochs, such a correlation is non-negligible, as shown in [31].This mainly arises from the interplay between the density contribution to the anisotropies of the astrophysical background and the Integrated Sachs-Wolfe (ISW) of the CMB [32].
On the astrophysical side, we consider GWs coming from unresolved Binary Black Holes (BBHs), expected to dominate the signal and very likely to be observed by Advanced LIGO [33].The quantity that links the distribution of GW sources to the underlying dark matter distribution is the bias.GW mergers, in fact, reside in galaxies, related to the distribution of halos who themselves trace the dark matter distribution.Furthermore, clusters correspond to very rare high peaks of the initial density field [34]; thus from the abundance of such objects it is possible to get information on the shape of the underlying distribution and it has been shown that the bias is affected by the presence of deviation from Gaussianity in the primordial Universe [35][36][37].In this work we try to extract information about nG from its effect on the GW bias.
Primordial nG is an extremely powerful observable to understand the level of interaction of the inflaton field, and to distinguish among different inflationary models.The Planck collaboration put strong constraints on nG from CMB observations, which happens to be very small [38].Other cross-correlation signals have been used in the past to constrain nG; for instance, [39] cross-correlates CIB anisotropies at different frequencies, showing that in the future it will be possible to constrain local nG down to |f NL | < 1; [40] put constraints on local nG by exploiting the cross-correlation between galaxy clustering and the ISW (see also [41]).
In this work we focus on the large-scale effect of local primordial nG on the GW bias.The presence of local nG, in fact, induces an additional scale-dependent term in the bias scaling as 1/k 2 and produces relevant effects on the AGWB spectrum (mainly on large scales), both on the auto-and on the cross-correlation.
To extract both the AGWB auto-correlation and the AGWB×CMB cross-spectrum we modify the CLASS code [42].We implement all the contributions to the AGWB anisotropies (density, Kaiser term, Doppler effect and gravitational potential contributions) and also the non-Gaussian contribution to the bias.We develop an external python module to model the astrophysical dependencies that are used to extract both the monopole AGWB energy density and the anisotropies, taking into account the latest LVK astrophysical constraints [5,43].In order to obtain the binary merger rate we start from the star formation rate (SFR) proposed by [44] and we convolve it with a time-delay distribution.The modelling of the waveforms is done through the phenomenological fit proposed by [45][46][47].Then we perform a signal-tonoise (SNR) analysis focusing on next generation space-based interferometers like LISA [2] and BBO [48].Finally, we forecast the capability of such detectors to constrain f NL through a Fisher matrix analysis.
The paper is organized as follows.In Section 2 we introduce the formalism of the SGWB and we briefly report the steps for obtaining the anisotropies of the AGWB.Then in Section 3 we briefly review the CMB anisotropies, fundamental part of this work.In Section 4 we discuss about the effects of primordial non-Gaussianity on the bias and finally in Section 5 we introduce the cross-correlation formalism.In Section 6 we discuss our results and we conclude in Section 7.

Stochastic Gravitational Wave Background
The SGWB is usually characterized in terms of the GW energy density per logarithmic observed frequency f o , defined as where Ω o is the solid angle along the line of sight n.Here, ρ c = 3H 2 0 /(8πG) represents the critical density of the Universe, with H 0 and G the Hubble parameter and the Newton's constant respectively.The total energy density includes both a background contribution, by definition homogeneous and isotropic, and a directional dependent one, which is the main focus of this work.Typically, one can also decompose the total energy density into a sum over the contributions to the GW energy density coming from different GW sources where the index "i" runs over different sources that contribute to the SGWB.As anticipated, we will focus on the contribution coming from BBHs.

Anisotropies
Defining with Ω GW the isotropic (background) energy density, the total relative fluctuation can be written as Analogously, one can define the single contribution coming from each source, so that the total one can be obtained as a weighted sum of the single ones, with weights w The anisotropic contribution was firstly derived in [27,49,50] and in a general gauge by [25] applying the Cosmic Rulers formalism [51] (see also the follow up paper [52]).In the following, we will report such a contribution in the Poisson gauge, following [25], since it is easier to distinguish the different terms and it is the way they are implemented in CLASS.The line element can be written as where we indicated with η the conformal time, x i the comoving coordinates, a the scale factor and Ψ and Φ the Bardeen potentials.It results [25,52,53] Again the index "i" runs over different sources that contribute to the AGWB.We indicated with χ(z) the comoving distance at redshift z, with H the comoving Hubble rate.The prime indicates the conformal time derivative, while the symbol ∂ ∥ a derivative along the line of sight; v is the peculiar velocity.It is important to specify that the quantity δ GW in equation (2.5) is expressed in the comoving gauge and not in the Poisson gauge, where it would not be gauge invariant.Thus, before linking it to the underlying matter density contrast as explained in the introduction, we map δ GW from the Poisson to the comoving gauge (where it is gauge invariant) as δ [i](P) evo is the evolution bias and accounts for the formation of new sources in time.It is defined as [54,55] b with N [i] the comoving number density of sources.The function W(z) is the astrophysical kernel that will be described in the Appendix B. It contains all the astrophysical information regarding the sources under consideration; it acts as a weight function for each redshift considered.On the first line of Eq. (2.5) we report the "density" term.Actually it is an instrinsic anisotropy and it is due to the inhomogeneous and anisotropic distribution of GW sources in the sky.This typically gives a dominant contribution in the auto-and cross-correlation.On the second line we can distinguish the Kaiser and the Doppler terms [25], on the third line the gravitational potential terms and, finally, on the fourth line the terms that are evaluated at the observer (indicated with the subscript o).Actually these latter give contribution only to the monopole and to the dipole, so in this analysis they have been disregarded.Before concluding this part we introduce the bias in our formalism.Black Hole mergers, in fact, are tracers of the underlying dark matter distribution.This translates in having a bias parameters linking their overdensity to the underlying dark matter distribution one as (2.8) In the following we will rename b BBH GW = b.We conclude adding that even if in general the bias is allowed a time dependence, in this work we consider an effective bias for the BBHs fixing b = 1.5, following [52,56] (see also e.g.[57,58]).

Angular decomposition
In order to characterize the AGWB anisotropies we decompose them on the two dimensional sky, i.e. on the sphere.So, we can work in harmonic space, decomposing the energy density in spherical harmonics as All the physics is contained in the coefficients of the expansion a GW ℓm that are given by (2.10) The index α runs over the various anisotropies defined in subsection 2.1, while the functions S [i]α ( k) are called source functions; they contain the information about the anisotropies and will be an important piece in the construction of the cross-correlation.Their full expressions can be found in [25].We report as an example the one related to the density anisotropies, with j ℓ (x) the spherical Bessel function and T δm is the transfer function linking δ m to the primordial curvature perturbation ζ and η is the conformal time.Here T δm is proportional to D(η)T (k), where D(η) is the growth function and T (k) is the usual matter transfer function obtained with CLASS [28,42].The angular power spectrum is defined as It is straightforward to show that where P ζ (k) is the primordial power spectrum (i.e. the two-point correlation function in Fourier space).The dimensionless power spectrum can be built from P ζ (k) as [28,59] where A s and n s are the scalar perturbations amplitude and tilt, respectively, and k p is the pivot scale.

Cosmic Microwave Background
The Cosmic Microwave Background has been found in 1964 by Penzias and Wilson, that were able to see just the monopole [60].Up to now we have been able to characterize both the monopole and the anisotropies of the signal with very high accuracy (up to ℓ ∼ 2500) [61].
Since the description and characterization of the CMB is out of the purposes of this work, but at the same time is an important part of the work, we will discuss just the main points, reporting the reader to specific works for a more detailed analysis.We start from defining the CMB temperature as [62] T Here the anisotropies are contained in the function Θ.The three main contributions to the CMB anisotropies on large scales are the Sachs-Wolfe, the Doppler and the Integrated Sachs-Wolfe (ISW).Actually the first one is imprinted at the moment of the decoupling and is a gravitational redshift effect.The Doppler depends on the peculiar velocity of the Earth with respect to the cosmic fluid, while the third one is still a gravitational redshift effect, but imprinted during the propagation of the signal (we will discuss this latter effect in more detail later).Also for the CMB the typical procedure is to expand the perturbation in spherical harmonics as with x o position of the observer.At the end one gets Here we indicated with T ℓ ( k) the source functions and also in this case there will be one for each term of the anisotropies.We report as an example the source function for the ISW, expected to dominate the cross-correlation with the density anisotropies described before, i.e.
In this equation T Φ and T Ψ are transfer functions that relate the gravitational potential to the primordial curvature perturbation ζ (for the explicit expression see e.g.[28]).The j ℓ 's are the spherical Bessel functions, η the conformal time and τ the optical depth.

Non-Gaussian Bias
The quantity that links the distribution of GW sources to the underlying dark matter distribution is the bias.In this work we focus on the effects of primordial local nG on it.In this case the primordial gravitational potential can be modeled as [20,21]  As expected, the density contribution (blue line) dominates the spectrum, even if the velocity terms (i.e., Kaiser and Doppler) provide a non-negligible contribution, especially at the large scales.The green line represents the gravitational potential contribution.We recall that the plot has been obtained modifying the public code CLASS [42] considering only the BBH contribution.
where Φ L is the Gaussian potential, while f NL is a parameter that quantifies the amount of primordial nG.Actually, starting from (4.1) one can show that the effects of nG can be recasted in an additional scale-dependent contribution to the bias ∆b such that [35][36][37] Here D(z) is the growth function, T (k) the matter transfer function and δ c = 1.686 is the spherical collapse linear over-density [34,63,64].Such a correction is proportional to f NL , meaning that in principle we could use this signal to constrain this parameter.Moreover we observe that the scale dependence goes like k −2 meaning that the effect of the correction will be important on large scales (small k).This is an important remark since actually we expect our cross-correlation signal to be sizeable mainly on large-scales and so sensitive to this correction.This strengthens the idea that we could use this signal to constrain nG.Actually, before proceeding, it is interesting to underline that the authors of [37] specify how equation (4.2) is not always valid.Summarizing, one can re-write the bias as  with p depending on the halo properties, e.g. if it is the result of a recent merger, if it is an old halo etc.In our case we are accounting for astrophysical sources of GW.Being formed typically at the end of a star life, they are located in galaxies with a very intense and recent stellar activity.As argued by [37] (see also e.g., [65,66]), finding the right value for p would make a difference in constraining f NL .In order to understand this better, we evaluate σ f NL /f NL for different values of p, i.e we consider the two extreme cases p = 1 and p = 1.6 in [37], since they are obtained through analytical estimation (rather than simulations) and are independent of the halo mass, to check how the constraints on f NL would vary.In general the proper value of p depends on the link between the properties of the tracer considered and the properties of the host halo.We further specify that in our analysis we fix the value of p instead of choosing a prior and then marginalize over it, due to the degeneracy between the p and f NL .Moreover, marginalizing over p is quite challenging, because its (but more generally b ϕ ) degeneracy with f NL makes it difficult to constrain both of them simultaneously.On the other hand, also the marginalization is not trivial.As shown in [67], for example, the choice of the prior even if quite wide, could bias the final result for f NL .Following [52] we work with time-and host-independent PDFs for the astrophysical parameters, since interested in giving an order of magnitude estimate of the signal.We perform the analysis only with BBO, since the SNR resulting from LISA is too low to perform a Fisher analysis, as we show in Appendix A. We conclude underlying that we perform our analysis in a full GR context, in agreement with [54].

Cross-Correlation formalism
We evaluate the cross-correlation of the AGWB with the CMB, including the effect of nG in the bias.Actually the cross-correlation arises mainly from the anisotropy related to the inhomogeneous distribution of GW sources in the sky and the late-ISW effect for the CMB.Let us better understand the physics behind this: BBH are relatively recent sources in the Universe evolution and as remarked before they trace the underlying dark matter distribution; so these sources reside in the large-scale gravitational potential.On the other hand, the late-ISW effect is first of all very recent (it is non-vanishing since the cosmological constant started to dominate the energy density) and also dominant on large scales.It is due to the time-variation of the gravitational potentials, that are constant during matter domination, but start varying again during Dark Energy domination (and this explains why it is recent).Moreover, as the perturbation re-enters the horizon it is damped by the expansion of the Universe and this justifies why the larger the scale, the higher the amplitude of the perturbation is expected.
Let us now think for a moment on the effects we expect when adding the nG correction to the bias (equation (4.2)) in the analysis.Assuming to be in the case in which (b − p) > 0, the effects of nG are expected just to increase (positive f NL ) or suppress (negative f NL ) the cross-correlation spectrum, while in the opposite case (b − p) < 0 the effects are just inverted (in the range of f NL allowed by Planck).Physically, this is due to a variation of the clustering properties of the tracers considered induced by primordial nG itself and it is expected to have an impact on the cross-correlation, but also on the constraints on f NL .

Prescription
In order to extract the cross-correlation signal we start by considering the single a ℓm .For the AGWB we have where s ℓm is the contribution coming from the actual signal, while n ℓm refers to the noise of the detector.In the case of the CMB we have neglected the instrumental noise of the CMB experiment, since at the multipoles of interest (ℓ ≤ 100), available CMB data are completely cosmic variance dominated, i.e. a CMB ℓm = s CMB ℓm .The angular cross-correlation spectrum results Note that the term ⟨n GW ℓm s CMB ℓm ⟩, vanishes since the noise of the instrument and the CMB signal are uncorrelated.In general the C ℓ coming from (5.2) are given by (notice the presence of the source functions defined before) Of course one should write the C ℓ as a sum of all the possible combinations among the various contributions seen before.Here we report the explicit expression for the dominant one, i.e.2 (5.4) On the second line, the integral runs from 0 to today.Actually the astrophysical kernel acts as a window function that cuts the integral to recent times (star formation has not been active since the beginning of the Universe).Thus the lower bound can be simply replaced with some η min , depending on the model chosen for the star formation.T δm is the transfer function for the density contrast, T Ψ and T Φ are the transfer functions for the gravitational potentials defined before.

Numerical Implementation
To obtain the cross-correlation angular spectra one should solve the Boltzmann equations for the perturbations from inflation up to today and then evaluate the various integrals as shown in (5.4).For this purpose we modify the publicly available code CLASS [42,68,69] where the CMB anisotropies are already implemented.We add the AGWB ones including the presence of the primordial nG correction to the bias.We also implement the possibility to use a time-dependent bias as an input and the possibility to vary the p parameter in the nG additional contribution.We fix the bias to 1.5 and we allow for two values of p (p = 1 and p = 1.6).
Moreover we develop a python module to model and evaluate the astrophysical dependencies.Actually such kernel is built with different moduli, in such a way it can be easily adapted to include other GW sources and astrophysical dependencies.More details on the astrophysical part are given in Appendix B.

Results
We describe in the following subsections the plots and the results we obtained in the evaluation of the cross-correlation spectra and in the constraining of nG.

Auto-and Cross-spectra
We plot in Fig. 1 and Fig. 2 the AGWB auto-correlation and the AGWBxCMB crosscorrelation angular power spectra.We plot the spectra up to ℓ ∼ 100, since it is unlikely that next generation GW detectors will be able to reach smaller scales.In the auto-correlation we observe how the density term dominates the spectrum, even if the velocity terms (Kaiser and Doppler) are expected to give a valuable correction [25].We specify that this behaviour is of course model dependent; more specifically the information about the source population comes mainly in the bias parameter, but also in the evolution bias described before.Furthermore, the behaviour of the auto-correlation and in particular of the density term is useful for a qualitative prediction of the cross-correlation spectrum, as will be explained below.Let us now focus on the cross-correlation spectrum.As anticipated before, we expect that the main contribution to such a spectra comes from the ISW-density term.The ISW, in particular its late contribution, results to be higher at the largest scales, to then decrease at smaller ones.On the other hand, the density anisotropy grows as the scales decrease.So it is clear that a cross-correlation among the two signals is expected to grow at the lowest ℓs (following the density behaviour) and then, after reaching a maximum, decrease following the dominating ISW behaviour.Fig. 2 shows the result of the numerical computation performed with the modified version of CLASS that confirms the expectations.The cumulative effects of the density cross-correlations, i.e. the correlation of all the CMB terms with the density anisotropies described before (blue line), clearly dominates the signal.The total contribution (black line) and the density almost coincide on all the scales, even if we observe that on the largest scales the correlation with the velocity anisotropies (red line) is non-negligible, giving a visible contribution (few percent) on the total signal.Among the density anisotropies, indicated with dash-dotted lines in the Figure , we observe that actually the density-ISW one prevails (cyan line), almost coinciding with the total density term (blue line), especially on large scales.
We specify that in principle the contribution coming from the cross-correlation of the density anisotropies with the Sachs-Wolfe could be important on scales k ≪ H/c and should be accounted [70].In this analysis we integrate starting at k ∼ 0.1 H/c (see also [55]) and for low values of f NL such a contribution can be neglected, at first approximation.We leave a further detailed analysis of this effect to a future work.
In Figure 3 we fix the bias at a reference value b = 1.5 and p = 1, and we plot the percent difference of the cross-correlation spectrum with respect to the Gaussian case due to different values of f NL .For positive values the spectrum increases, up to 40% in the case in which f NL = 10, while for the negative ones we observe a suppression of the spectrum on the largest scales of the same amount 3 .This symmetry in the behaviour of the spectrum when varying from positive to negative values of f NL is due to the peculiar linear dependence of the cross-correlation on such a parameter.The auto-correlation, for example, will present both a quadratic term in f NL ∝ (f NL (b − p)) 2 that will always give a positive contribution and a linear one that goes like f NL (b − p).This latter on the other hand will behave as in the cross-correlation case.It is clear that the combination of the two leads to a break of the symmetry of the nG effects on the auto-correlation.
Finally we specify that the effects on the spectrum depend also on the values chosen for b and p.At fixed p, if b > p and f NL is positive we get the behaviour described above, but for a bias b < p, actually the effects on the power spectrum are reversed.We stress that what really matters is the whole f NL (b − p) contribution and not only f NL : in the case in which b ≃ p, it would not be possible to appreciate the effects of nG.This further underlines the importance of a precise modelling of the bias.Moreover, accounting for a time-dependence in the bias could make the difference, since in principle it could be greater than p at a certain time and less than p in another, leading to interesting effects in the behaviour of the power spectra.This of course would need to take into account the evolution of the GW sources but also of the halos in which they reside.We leave such a detailed analysis to a future work.

Fisher Forecast
We perform a Fisher forecast to provide estimates of the capability of future GW detectors to constrain f NL through the AGWBxCMB cross-correlation.The Fisher matrix, in fact, provides a powerful tool for describing the information content coming from observables in terms of underlying parameters; it can be written as [71] where in the last step we have defined ), C ℓ being the covariance matrix and θ α , θ β the parameters.In our case we have just one single parameter, i.e. f NL .Since the C CMB ℓ do not depend on the parameter we are willing to estimate, equation (6.1) reduces to We have included the angular spectrum of the instrumental noise N ℓ of the interferometers in the covariance matrix.Finally one obtains that the error is given by , while solid ones show the case in which also the C XC ℓ are added.In this latter case it is possible to appreciate the non-negligible improvement in the value of σ fNL .Since the exact value of p is strongly dependent on the evolution properties of the halos and of the GW sources tracing them, we show how the variation of this parameter impacts on the constraining power of the signal.
In Fig. 4 we plot the σ f NL /f NL as function of f NL for two different values of p, i.e. p = 1 (blue lines) and p = 1.6 (red lines).We keep fixed the bias to the reference value of 1.5 and we consider the BBO detector in the evaluation of the covariance for each multipole.Furthermore we show the case in which only the C AGWB ℓ are considered in the analysis (dashed lines) and the one in which also the C XC ℓ are included, to show the impact of these latter on the value of σ f NL .Before discussing the effects of different values of p, we highlight that we get a higher constraining power for higher values of f NL .The low angular sensitivity of LISA leads to a low SNR, as shown in Appendix A, suggesting, in this case, a poor constraining power.On the other hand, as we can see in Fig. 4, BBO will be able to put tight constraints on the f NL parameter, being able to reach percent error levels.
Let us now discuss the impact of different values of p on σ f NL .We can see that in the case p = 1 we are able to put tighter constraints with respect to p = 1.6 case.This behaviour is of course dependent on the value of the bias and on the parameters chosen for modelling the clustering properties of the tracers.In this case, fixing the bias b to 1.5, it is clear that for p = 1.6 the effect of nG are suppressed while for p = 1 they are enhanced and so in this latter case we achieve a better constraining power.Thus, the plot shows how a good modelling of the sources is crucial in constraining f NL .Furthermore, the difference between the dashed and the solid lines in the plot shows the importance of the cross-correlation.The inclusion of C XC ℓ leads to a significative improvement in the value of σ f NL .We can conclude, then, that not only the cross-correlation among the AGWB and the CMB is non-negligible, but it is also relevant for the f NL constraints.
Then we compute the SNR for both LISA and BBO (see Appendix A for more details on the SNR analysis).In Figure 5 we plot the SNR as a function of ℓ max for LISA and BBO.We can see that we are able to gain information mainly from the very first multipoles (up to ℓ ∼ 4 − 5), meaning that, in principle, we could achieve a greater constraining power on f NL , if a better angular resolution could be reached.We have included only the terms coming from the auto-and cross-correlation parts (i.e.C AGWB ℓ and C XC ℓ ).We can see that LISA will not be able to reach SNR larger than one, even for large (and positive) values of f NL .On the other hand, in the case of BBO, we can see that we are able to reach SNR larger than one for 5 years of observation.This reflects on the ability of constraining f NL .Furthermore, we stress that this behaviour is strongly dependent on the choices made for both the cosmological and the astrophysical parameters, underlying that a better modeling is necessary.

Conclusions
In this paper we analysed the effects of local primordial nG on the cross-correlation between CMB and AGWB anisotropies.We explored to which extent it is possible to use astrophysical GW sources to constrain early-universe cosmology.We introduced the theoretical formalism for describing the CMB and AGWB anisotropies and then we discussed the cross-correlation signal.This is expected to be dominated by the correlation between the density anisotropies of the AGWB and the CMB (late)-ISW effect, both effects being very recent and important on large scales.We confirmed such a behaviour by plotting the expected angular power spectrum obtained modifying the publicly available code CLASS.We added the effects of local primordial nG on the bias and the AGWB anisotropies and we modelled the population of BBHs accounting for the latest LVK astrophysical constraints.We found that the inclusion of primordial nG would lead to an enhancement (or a suppression) of the power spectrum, mainly on large scales.Then we performed a Fisher forecast on the f NL parameter to understand to which extent future space-based interferometers (e.g. and BBO) will be able to put constraints on nG.We find that BBO will be able to put very tight constraints on f NL .This behaviour is confirmed by the SNR analysis, which clearly shows that BBO can reach a SNR ∼ 4 for 5 years of observation.On the other hand the SNR for LISA results quite low.
Our analysis shows the importance of the clustering properties of GW sources and their impact on the signal.We considered the simplified case of a constant bias and b ϕ ∝ (b − p), obtaining that positive values of f NL enhance the spectrum, while negative values suppress it.We stress that this is not a general behaviour, meaning that the enhancement/suppression depends not only on the value of f NL , but also on the clustering properties of the sources.An accurate modelling of such quantities is crucial for constraining f NL .A time-dependent bias and a better modelling of the p parameter would impact our analysis, avoiding underor over-estimation of the error bars. 4e conclude stressing that, even if GW detectors are limited in angular sensitivity to the first few multipoles, they can play an important role in constraining early-universe cosmology.

A SNR Evaluation
In the main text we reported the results for the computation of the expected AGWB×CMB signal at different values of f NL .We focus in this Appendix on its detectability by means of a SNR analysis.We evaluate the SNR accounting for the contributions coming from both the auto-and the cross-correlations, i.e.
where C −1 ℓ is the inverse covariance matrix that here we explicit Here the σ 2 ℓ s are the covariances and contain both the cosmic variance and the noise contribution.We report as an example We consider the angular power spectra for the noise of LISA and BBO (normalized with respect to the monopole) and we evaluate the resulting SNR.In Fig. 5 we report the cumulative SNR.We consider just the contribution coming from the auto-and crosscorrelation, neglecting the CMB contribution (otherwise the SNR coming from the CMB, would be so high to cover the effects of the SGWB we are interested in).For both the detectors we see that only the multipoles up to ℓ ∼ 4 − 5 give a relevant contribution.In particular we see that LISA will not be able to detect such a signal, resulting in a SNR lower than 10 −1 , suggesting that the Fisher forecast will be unable to put tight constraints on f NL .On the other hand we see that BBO is able not only to reach a SNR of order one, but is so sensitive that the SNR is limited for the first multipoles just by the cosmic variance; it then saturates at a value of almost 4.This clarifies the very strong constraining power for f NL showed in the main text but at the same time underlines how promising is the GW search and the need to improve angular resolution.The plots also show the impact of nG on the SNR.We see that effect of nG are more relevant for LISA, causing and increase (positive f NL ) or a decrease (negative f NL ) of the SNR.We clarify that the suppression just explained happens only for low values of f NL .In the case of a negative f NL but with a large absolute values, assuming (b − p) positive, the second term of (4.3) would give a positive dominating contribution in the auto-correlation spectrum, leading to an increase of the SNR.For the models adopted in this work, this effect would be important at values of f NL around 10, so outside the constraints given by Planck.

B Astrophysical Kernel
As reported in Subsection 2.2, all the astrophysical information is included in the function W [i] [η(z)], i.e. the astrophysical kernel.It acts as a weight function in redshift and is defined as Here f o is the observed frequency, f e is the frequency at the emission and ⃗ θ encloses all the astrophysical parameter dependencies, like the masses of the two BHs, m 1 and m 2 , and the merger rate today, R BBH 0 , expected to be between 17.9 and 44 Gpc −3 yr −1 [43].As reported by [5], the upper bound on the monopole coming from BBH is Ω BBH GW (25Hz) = 5.0 +1.7 −1.4 × 10 −10 , while for the total contribution it results 3.4×10 −9 at 25 Hz.N [i] (z, f e , ⃗ θ) is the total comoving density of GWs, written as where R [i] is the intrinsic comoving merger rate of BBHs and dE GW /df e the energy of the GWs emitted with a frequency f e .We take into account all the stages of evolution of a binary system, i.e. the inspiral, the merger and the ringdown, following [45][46][47].We computed also the merger rate of BBHs as follows [52]: • We start considering the star formation rate following the phenomenological expression from [44] R ⋆ (z) ∝ (1 + z) with λ 1 = 2.7 and λ 2 = 2.9 and that holds for 0 ≤ z ≤ 8.
• We then assume the binary formation rate to be proportional to the star formation rate as with C 1 in the range 0 -1 since of course not all the stars end up forming a binary.
• Finally we account for the time delay between the formation and the merger of the binary, so that the merger rate can be obtained after integrating over the time delay, where p(t d ) ∝ t −1 d [5], assuming t d,min = 50 Myr for BH binaries.In order to account for the multiplicative constant, we set that the final merger rate has to be consistent with R BBH today provided by the latest constraint from LIGO-Virgo-Kagra collaboration [43].
To model the BH mass distribution we consider the Powerlaw + Peak mass model [43].The w(z, ⃗ θ) function has the role to discriminate resolved from unresolved sources in the computation, so actually it represents the efficiency of the detectors in the measurement of the AGWB.We obtain it by integrating the pdf of the SNR from 0 to the threshold that determines a detection of a resolved source, that we fix to 8 [72] so that w(z, ⃗ θ) = We conclude that a better modeling of the window function would include the generation of a catalog of sources, also to better understand their clustering properties and better discriminating when a source can be considered resolved or not.We leave this for a more dedicated future analysis.

Figure 1 .
Figure 1.The plot shows the auto-correlation angular power spectrum of the AGWB anisotropies.As expected, the density contribution (blue line) dominates the spectrum, even if the velocity terms (i.e., Kaiser and Doppler) provide a non-negligible contribution, especially at the large scales.The green line represents the gravitational potential contribution.We recall that the plot has been obtained modifying the public code CLASS[42] considering only the BBH contribution.
) with α ∝ T (k)D(z) and b ϕ = 2δ c (b − 1).This latter is called universality relation, but it is not always true.One can show that a more general relation would be b ϕ = 2δ c (b − p),

Figure 2 .
Figure 2. The plot shows the cross-correlation angular power spectrum of the AGWB with the CMB.The dominant contribution, as expected, comes from the CMB ISW with the density anisotropies of the AGWB (light blue dashed line).The black line represents the total contribution.We note that on large scales the velocity anisotropies provide a non-negligible contribution in the signal, while on smaller ones dominates the correlation with the density ones.

Figure 3 .
Figure 3.The plot shows the percent difference of the angular power spectrum for the crosscorrelation between the AGWB and the CMB in the nG cases with respect to the Gaussian one.Actually we consider different values of f NL , showing that the nG contribution leads to a suppression (or enhancement) up to 40% on large scales.The exact behaviour of the curves is actually dependent on the clustering properties of the GW sources, i.e. on the parameters b and p under consideration.In this case we consider b = 1.5 and p = 1.

Figure 4 .
Figure 4.The plot shows the variation of σ fNL /f NL as function of f NL for the BBO detector for different values of p (blue lines for p = 1.6 and red lines for p = 1).Dashed lines refer to the inclusion in the analysis of only the C AGWBℓ

Figure 5 .
Figure 5.The plot shows the cumulative SNR obtained for different values of ℓ max in the cosmic variance-limited case (orange line), accounting for the BBO noise (blue lines) and for the LISA one (red lines).We also account for the effect of different values of f NL on the SNR, showing that they potentially enhance or suppress it.