Characterizing the Cosmological Gravitational Wave Background Anisotropies and non-Gaussianity

A future detection of the Stochastic Gravitational Wave Background (SGWB) with GW experiments is expected to open a new window on early universe cosmology and on the astrophysics of compact objects. In this paper we study SGWB anisotropies, that can offer new tools to discriminate between different sources of GWs. In particular, the cosmological SGWB inherits its anisotropies both (i) at its production and (ii) during its propagation through our perturbed universe. Concerning (i), we show that it typically leads to anisotropies with order one dependence on frequency. We then compute the effect of (ii) through a Boltzmann approach, including contributions of both large-scale scalar and tensor linearized perturbations. We also compute for the first time the three-point function of the SGWB energy density, which can allow one to extract information on GW non-Gaussianity with interferometers. Finally, we include non-linear effects associated with long wavelength scalar fluctuations, and compute the squeezed limit of the 3-point function for the SGWB density contrast. Such limit satisfies a consistency relation, conceptually similar to what found in the literature for the case of CMB perturbations.


Introduction
The current ground based interferometers are close to reach the expected sensitivity to detect the Stochastic Gravitational Wave Background (SGWB) from unresolved astrophysical sources [1]. Future space-based (such as LISA [2] and DECIGO [3]) and earth-based (like Einstein Telescope [4,5] and Cosmic Explorer [6]) interferometers have the potential to detect the SGWB of cosmological origin (see [7][8][9][10] for reviews of possible cosmological sources). It is likely that a detection of a cosmological SGWB background will require the ability to discriminate it against the astrophysical signal. Astrophysical GW background (AGWB) arises from the superposition of the signals emitted by a large population of unresolved sources that are mainly dominated by two types of events: (i) the periodic long lived sources (e.g. the early inspiraling phase of binary systems) where the frequency is expected to evolve very slowly compared to the observation time; (ii) the short-lived burst sources, e.g. core collapse to neutron stars or black holes, oscillation modes, r-mode instabilities in rotating neutron stars, magnetars and super-radiant instabilities (for example, see [11,12]). Several techniques have been developed to distinguish among the various backgrounds. The most obvious tool for this component separation is the frequency dependence [13], as several cosmological mechanisms are peaked at some given characteristic scale. However, future detectors will allow for a better angular resolution of anisotropies of the astrophysical background. Therefore, another tool could be the directionality dependence of the SGWB [14][15][16][17][18][19] and, as we explore here, its statistics.
In this work, we discuss graviton propagation through a Boltzmann approach [15] as it is typically done for the CMB. Specifically, we construct and evolve the equation for the distribution f of gravitons in a FLRW background, plus first order scalar and tensor perturbations (we also consider how non-linear effects for the specific case of squeezed non-Gaussianity, as we discuss at the end of this Introduction). At the unperturbed level, following the isotropy and homogeneity of the background, the distribution depends only on time and on the GW frequency p/2π (where p is the physical momentum of the gravitons) through the combination q ≡ p a, where a is the scale factor of the universe. Namely, the gravitons freely propagate, and their physical momentum redshifts during the propagation. This property is shared by any free massless particles, and, in particular, also by the CMB photons. On the other hand, differently from the photon distribution, the initial population of gravitons is not expected to be thermal (as we have in mind production mechanisms, such as inflation [20,21], phase transitions [22], or enhanced density perturbations leading to primordial black holes (PBH) [23,24,31], which occur at energies well below the Planck scale) which leaves in the distribution a sort of "memory" of the initial state. As we show, the fact that the spectrum is non thermal generically results in angular anisotropies that have an order one dependence on the GW frequency. This is in contrast with the CMB case, for which this dependence only arises at second order in perturbation theory.
This initial state will in general be anisotropic, as no mechanism of GW production can be perfectly homogeneous. Additional anisotropies are induced by the GW propagation in the perturbed universe. As we are interested in large scale, we work in a regime of a large hierarchy q k between the GW (comoving) momentum q and the (comoving) momentum k of the large scale perturbations. We confirm that in the angular power spectrum, the Sachs Wolfe (SW) effect is dominating on large scales also for gravitons, while the Integrated Sachs-Wolfe (ISW) contribution is subdominant.
We employ this approach to study the non-Gaussianity of the SGWB energy density. Although we are not aware of any dedicated analysis in this sense, it is reasonable to expect that the SGWB produced by incoherent astrophysical sources is Gaussian, due to the central limit theorem. In light of this fact, a measurement of non-Gaussianity would be a signal of large scale coherency, that would likely point to a cosmological origin of the signal. Previous works showed that inflation can result in a sizeable an nonvanishing 3−point function h 3 for the graviton wave function, but that this is generically non observable in interferometers [23,24], due to the decoherence of the phase the GW wave-function h induced by the GW propagation, and due to the finite duration of the measurement (see [26] for a possible exception to this conclusion, occurring for a very specific shape of the bispectrum). Since the phase does not affect the GW energy density, we argue that the energy density is a much better variable to study the statistics of the SGWB. Also in this case, the of non-Gaussianity can be induced both by the production mechanism and the propagation. As an example of the former, in ref. [31] we recently computed the 3−point function of the SGWB energy density that arises in presence of non-Gaussianity of the scalar perturbations of the local shape (in presence of this non-Gaussianity, a long-scale mode of momentum k can modulate the power of the short-scalr scalar perturbations that are responsible for the PBH formation). Here we study the 3−point function induced by the GW propagation. This is also proportional to the non-Gaussianity of the scalar perturbations. In this sense, the SGWB can be used as a novel probe (beyond the CMB and the LSS) of the non-Gaussianity of the scalar perturbations.
Although in most of this work we limit our attention to linearized fluctuations, in Section 6 we consider non-linear effects induced by long-wavelength scalar perturbations, which modulate correlation functions involving short-wavelength modes. We make use of a powerful method first introduced by Weinberg in [27], which focusses on adiabatic systems, and identifies the effects of long modes with an appropriate coordinate transformation. Applying this method to our set-up, we compute how non-linearities induce a non-vanishing squeezed limit of the 3-point function for the SGWB density contrast. We determine how such squeezed limit depends on the scale-dependence of the spectrum of primordial scalar fluctuations; on the momentum dependence of the background SGWB distribution; and on the time, scale, and direction dependence of the scalar transfer functions connecting primordial to late-time adiabatic scalar fluctuations.
The paper is organized as follows. In Section 2 we present the computation and the formal solution of the Boltzmann equation for GW propagation. In Section 3 we decompose the formal solution in spherical harmonics, paralleling a treatment that is familiar in the study of CMB perturbations. In Section 4 we compute the angular power spectrum and bispectrum of the SGWB perturbations. In Section 5 we review one physical mechanism that can result in a sizeable cosmological SGWB with some degree of anisotropy. In Section 6 we study non-linear effects on the squeezed bispectrum. These results are discussed and summarized in Section 7. The paper is concluded by three appendices. Appendix A contains the details of the computation of the anisotropies due to the large-scale tensor perturbations. Appendix B provides some intermediate steps on the computation of the GW bispectrum induced by tensor modes. Finally, Appendix C presents an immediate connection between our formal solutions and the CMB results obtained in the case of initial thermal state.
Part of the results contained in the present work were also summarized in the Rapid Communication [58].

Boltzmann equation for gravitational waves
We consider first order perturbations around a Friedmann-Lemaitre-Robertson-Walker (FLRW) background in the Poisson gauge where a(η) is the scale factor as a function of the conformal time η. Φ and Ψ are scalar perturbations while χ ij represent the transverse-traceless (TT) tensor perturbations. We neglect linear vector modes since they are not produced at first order in standard mechanisms for the generation of cosmological perturbations (as scalar field inflation), and we consider tensor modes at linearised order. Given the statistical nature of the GW we can define a distribution function of gravitons as f = f (x µ , p µ ), which is function of their position x µ and momentum p µ = dx µ /dλ, where λ is an affine parameter along the GW trajectory. As we will see, observables as number density, spectral energy density, and flux (directions) can be derived from the distribution function. The graviton distribution function obeys the Boltzmann equation where L ≡ d/dλ is the Liouville term, while C and I account, respectively, for the collision of GWs along their path, and for their emissivity from cosmological and astrophysical sources [15]. The collision among GWs affects the distribution at higher orders (in an expansion series in the gravitational strength 1/M P , where M P is the Planck mass) with respect to the ones we are considering, and they can be disregarded in our analysis (see [29] and references therein for a discussion of collisional effects involving gravitons). The emissivity can be due to astrophysical processes (such as black hole merging) in the relatively late universe, as well as cosmological processes, such as inflation or phase transitions. In this work we are only interested in the stochastic GW background of cosmological origin, so we treat the emissivity term as an initial condition on the GW distribution. This leads us to study the free Boltzmann equation, df /dη = 0 in the perturbed universe wheren ≡p is the GW direction of motion, and where we have used the comoving momentum q ≡ | p|a (as opposed to the physical one, used in [15,30]). This simplifies the equations by factorizing out the universe expansion. The first two terms in (2.3) encode free streaming, that is the propagation of perturbations on all scales. At higher order this term also includes gravitational time delay effects. The third term causes the red-shifting of gravitons, including the Sachs-Wolfe (SW), integrated Sachs-Wolfe (ISW) and Rees-Sciama (RS) effects. The fourth term vanishes to first order, and describes the effect of gravitational lensing. We shall refer to these terms as the free-streaming, redshift and lensing terms, respectively, as customarily done in CMB physics.
Keeping only the terms up to first order in the perturbations, Eq. (2.3) gives where we have followed the standard procedure developed for the CMB in [30,39]. The distribution function f can be expanded as f η, x i , q, n i =f (q) + f (1) η, x i , q, n i + .... ≡f (q) − q ∂f ∂q Γ η, x i , q, n i + .... , (2.5) where the dominant, homogeneous and isotropic contributionf (q) solves the zeroth order Boltzmann equation. The function f (1) (η, x i , q, n i ) is the solution of the first order equation, and the ellipses denote the higher order solutions in a perturbative expansion. In this expression we have parameterized the first order solution in terms of the function Γ, so to simplify the first order Boltzmann equation [15]. For a thermal distribution with temperature T , one finds Γ = δT /T . This is particularly the case for the CMB, for which, due to the thermalization, the temperature anisotropies are frequency-independent up to second order in the perturbations. For gravitons, as we already mentioned, the collisional term is extremely small, and, for a generic production mechanism, Γ generically retains an order one dependence on frequency (as we show below, also for the GW case the propagation effects induce frequency-independent perturbations at linear order). The zeroth order homogeneous Boltzmann equation simply reads ∂f /∂η = 0, and it is solved by any distribution that is function only of the comoving momentum q, namely f =f (q). In our approach this solution is simply given as the homogeneous part of the initial condition. As a consequence, the physical momentum of the individual gravitons redshifts proportionally to 1/a, and the physical graviton number density n ∝ d 3 pf (q) is diluted as a −3 as the universe expands. This is also the case for CMB photons, whose distribution functionf CM B = (e p/T − 1) −1 is only controlled by the ratio p/T ∝ a p = q, where T is the temperature of the CMB bath. We see that these rescalings with a are a consequence of the free particle propagation in the expanding FLRW background, and they do not rely on the distribution being thermal.
As anticipated, from the graviton distribution function, evaluated at the present time η 0 , we can compute the SGWB energy density where we have introduced the spectral energy density Ω GW and the critical density ρ crit = 3H 2 M 2 p . Here H ≡ (1/a 2 ) da/dη is the Hubble rate. Following standard conventions, the suffix 0 denotes a quantity evaluated at the present time.
Contrary to most of the studies of the SGWB, that assume a homogeneous Ω GW , in our case the GW energy density depends on space. We denote the homogeneous component of For the full spectral energy density, we define and we introduce the SGWB density contrast In terms of the function Γ, the first order Boltzmann equation reads [15] ∂Γ ∂η n i n j ∂χ ij ∂η is the source function which includes the physical effects due to cosmological scalar and tensor inhomogeneities. We note that the source is q−independent (thus showing that the anisotropies arising at first order from propagation effects are frequency-independent, as we anticipated).
To solve this equation, it is convenient to Fourier transform with respect to spatial coordinates, and analogously for the other variables (we use the same notation for a field and for its Fourier transform, as the context always clarifies which object we are referring to). This leads to Γ + i k µ Γ = S(η, k,n) , (2.12) where from now on prime denotes a derivative with respect to conformal time, and where we denote by µ ≡k ·n , (2.13) the cosine of the angle between the Fourier variable k and the direction of motionn of the GW. In Fourier space the source term reads (2.14) With this information in mind, Eq. (2.12) is readily integrated to give Γ η, k, q,n = e ikµ(η in −η) Γ η in , k, q,n We integrate the second term in the second line by parts, and obtain with the last two terms in the first line being the boundary terms of this integration. In the following section, we decompose then-dependence of the solution (representing the arrival direction of the GW on our sky) in spherical harmonics. As we are not interested in the monopole term, we can disregard the −Φ(η, k) contribution to the solution, and write The first term, which was disregarded in [15], carries the "memory" of the initial conditions. Due to this term, the GW energy density anisotropies are generically dependent on the frequency q. We discuss an example of this fact in Section 5, where we study the SGWB produced in axion inflation. Generally, this term has also a dependence onn. This implies that the solution has a dependence on the directionn, which is more general than the one arising from the projection of k on the line of sightn. (Indeed, the remaining terms in Eq. (2.17) depend onn only through the µ ≡k ·n combination. Thanks to this fact, they result in angular correlators that are statistically isotropic (as we show in the next two sections).) On the other hand, the angular dependence present in the first term of (2.17) could result in statistically anisotropic correlators (specifically, 2-point and 3-point correlators) that have a more general dependence on the multipoles coefficients i and m i than Eqs. (4.3). This would indicate an overall anisotropy of the mechanism responsible for the GW across the entire universe, and, ultimately, a departure from an exact FLRW geometry. While we believe that this can be an interesting topic for future exploration, the present work focuses on the statistically isotropic case, and we assume an initial condition of the form Γ in = Γ(η in , k, q), which guarantees such a condition.

Spherical harmonics decomposition
We separate the solution (2.17) in three terms Γ η, k, q,n = Γ I η, k, q,n + Γ S η, k,n + Γ T η, k,n , where I, S, and T stand for Initial, Scalar and Tensor sourced terms respectively and they are given by Similarly to what is usually done for the CMB, in order to compute the angular power spectrum, in an all-sky analysis we decompose the fluctuations using spin-0 or spin-2 spherical harmonics. Since Γ is a scalar, we can express it as where we recall thatn is the direction of motion of the GWs. More specifically:

Initial condition term and q−dependent anisotropies
Let us first evaluate the initial condition term Following the standard treatment for CMB anisotropies [30], we make use of the identity (where j and P are, respectively, spherical Bessel functions and Legendre polynomial) so to obtain Here x 0 denotes our location (that can be set to the origin), η 0 denotes the present time, and η in the initial time. Once again we stress the peculiar property of the initial condition, namely its dependence on the frequency q. In Section 4 we discuss how this imprints the SGWB angular spectrum.

Scalar sourced term
A second source of anisotropy is due to the GW propagation in the large-scale scalar perturbations of the universe (the wavenumber of these perturbations k is many order of magnitudes smaller than the GW frequency q, and the GW acts as a probe of this large-scale background). As long as the scalar perturbation is in the linear regime (which is the case for the large-scale modes that leave an impact on the large-scale anisotropies of our interest), we can express it [30] as a transfer function (a deterministic function that encodes the time-dependence of the perturbations) times a stochastic variable ζ. This assumes the absence of isocurvature modes, and, in particular, of anisotropic stresses, as for example those due to the relic neutrinos. This also assumes that the statistical properties of ζ have been set well before the propagation stage that we are considering (for instance during inflation, or during some early phase transition). Therefore, the scalar perturbations are Under the above assumptions, T Φ (η, k) = T Ψ (η, k). However, we keep these two terms as distinct in our intermediate computations, so that the present analysis can be most easily generalized, if one wishes to introduce more general sources. With this in mind, the scalar sourced term becomes and we note that we are assuming a single adiabatic mode (i.e. ζ( k) is the operator associated with the conserved curvature perturbation at super-horizon scales). Proceeding as above, (3.10) As we can see, also the SGWB, feels, similarly to the CMB, a Sachs-Wolfe and integrated Sachs-Wolfe effect, which are represented by the first and the second term in (3.10), respectively.

Tensor sourced term
Finally, the third contribution Γ m,T is due to the GW propagation in the large-scale tensor modes To evaluate such term we decompose the tensor modes in right and left-handed (respectively λ = ±2) circular polarizations (see e.g. [43]), The three factors involved in each term are, respectively, the tensor circular polarization operator, the tensor mode function (equal for the two polarizations), and the stochastic variable for that tensor polarization (that is the analog of ζ we discussed in the scalar case). Inserting this decomposition in Eq. (3.11), a lengthy algebra, that we report in Appendix A, leads to which is formally analogous to Eq. (3.10), with the productζ Y * m replaced by the combination λ=±2ξ λ ( k) −λ Y * m (Ω k ), involving the spin-2 spherical harmonics, and with the scalar transfer function replaced by the tensor one.

Summary of the three contributions
The results derived in the three previous subsections can be written in the (slightly) more compact form where we have introduced the linear transfer function T X(z) , with X = S, T which represents the time evolution of the graviton fluctuations originated from the primordial perturbation

Correlators of GW anisotropies and SGWB non-Gaussianity
We now compute the 2-point Γ m Γ * m and the 3-point Γ 1 m 1 Γ 2 m 2 Γ 3 m 3 angular correlators of the solutions (3.14). The statistical operators entering in these solutions are the four momentum-dependent quantities Γ(η in , k, q), ζ( k), ξ R ( k), and ξ L ( k), while the other terms encode deterministic effects such has the time evolution of the large-scale modes (in the linearized theory of the cosmological perturbations) and the projection of the GW anisotropies in the harmonic space. In this study, we assume that the stochastic variables are nearly Gaussian, with the 2-point functions and a subdominant 3-point component The assumption of nearly Gaussian modes is experimentally verified for the large-scale perturbations of ζ and of ξ λ , as obtained from the CMB data [44]. We assume that this is the case also for the initial condition term. The expressions (4.1) and (4.2) can be readily used to compute the angular correlators of the solutions in (3.14). Moreover, for simplicity of exposition, we have here assumed that the various terms are not cross-correlated. These results in separate sets of correlators for the three terms in (3.14). This assumption can be easily relaxed, and in fact, we did so in [31] where we studied the anisotropic distribution of the GW originated in models with primordial black holes, as we review in Section 5.
The computations performed so far assume statistical isotropy (recall the discussion at the end of Section 2). Correspondingly, when we combine (4.1) and (4.2) with (3.14) we obtain angular correlators with well specific dependence on the multipole indices. Specifically, the two point correlators have the dependence while, under the above assumption, the angular power spectrum and the reduced bispectrum consists of the three separate contributions We recall that the form of the bispectrum factorizes the Wigner-3j symbols [59], which are nonvanishing only provided that i m i = 0 and that the three i satisfy the triangular inequalities.
In the following we provide the explicit expression for the various contributions to the power spectrum and the reduced bispectrum.

Angular power spectrum of GW energy density
We start with the computation of the two-point function of the initial condition term. From the first of (3.14) we can write The correlator of the initial condition term is then given by the first of (4.1). Using this, and the orthonormality condition of the spherical harmonics, d 2n which indeed is of the form dictated by statistical isotropy. The other two terms are obtained analogously. Altogether, we find We know from the CMB that the large-scale tensor modes have a power smaller than the scalar ones. At large scale, the scalar contribution is dominated by the term proportional to the initial value of Φ in T (0) , which is the analog of the SW contribution for the CMB. The large-scale modes that we are considering re-entered the horizon during matter domination. For these modes, ignoring the late time dark energy domination, T Φ = T Ψ = 3/5 [30]. So, for scale invariant power spectra, The second term can be compared to the SW contribution to the CMB anisotropies. In that case, the final temperature anisotropy is 1/3 times the scalar perturbation at the last scattering surface, while Φ at that moment decreased by a factor 9/10 in the transition from radiation to matter domination [30]. With this in mind, the second term in (4.8) leads to C SW = (3/10) 2 C ,S , in agreement with the CMB literature. On the other hand, if the two contributions are correlated, as it would be the case for adiabatic initial condition for Γ I , then both terms in (4.8) contribute to the SW effect for the SGWB.

Angular bispectrum of GW energy density
The characterization of the non-Gaussian properties of the SGWB is a potential tool to discriminate whether a SGWB has a primordial or astrophysical origin. The primoridal 3point function of the GW field, h 3 , is unobservable, due to the decoherence of the associated phase (because of the propagation, and the finite duration of the measurement [23,24]), with, possibly, the exception of very specific shapes [26,33]. It is more convenient to consider the non-gaussianity associated to the GW energy density angular distribution, which is not affected by this problem [31]. This gives rise to the bispectra in (4.4), which we evaluate now.
As we did for the power spectrum, also in this case we start from the initial condition term. Combining the first of (3.14) and the first of (4.2) leads to We then use the representation of the Dirac δ−function, 10) and the orthonormality of the spherical harmonics, to arrive to where we have introduced the Gaunt integrals We remark that also the bispectrum from the initial condition also generally as an O (1) dependence on the GW frequency. An analogous computation leads to the contribution from the scalar modes For the tensor sourced contribution we have (4.14) Following [60], in Appendix B we show that also this contribution can be cast in a similar form to the previous two terms: (4.16) We remark once again that we have neglected for simplicity all the mixed scalar-tensor correlators.

Reduced Bispectrum and estimation
The three contributions to the bispectrum found above have the correct form (4.3) as dictated by statistical isotropy. For convenience, we collect here the explicit form of the reduced bispectra contributing to (4.4) To estimate the SGWB bispectrum, we consider only the scalar source contribution B 1 2 3 ,S and we assume the simplest small non-linear coupling local ansatz for the curvature perturbation where ζ g ( x) denotes the linear Gaussian part of the perturbation. With the local ansatz, the bispectrum of the scalar perturbations assumes the form [46,47] We insert this in the second of (4.17) and we assume a matter transfer function T Φ (η, k) = T Ψ (η, k) = 3/5 g (η) with the growth factor g(η) = 1 and a scale invariant spectrum for the primordial curvature fluctuations. We can then integrate over one of the internal momenta k i , 2 π dk k 2 j (k η 0 ) j (k r) (4.20) The relation (4.20) is exact if k ranges up to infinity, which is not the case for the innermost momentum (as the integral (4.20) is performed first, this will necessarily be the the momentum that we order to be the innermost one), due to the triangular inequalities associated to the bispectrum. The condition 1 ensures that the support of the integration occurs at sufficiently small k, so that the relation (4.20) becomes exact at large . The result then allows to immediately perform the integral over r. We then find that the reduced bispectrum from the scalar contribution, assuming that the SW is the dominant contribution, is This result can also be written in terms of the two-point functions found in Eq. (4.7): which resembles the one for the CMB angular bispectrum in the Sachs-Wolfe regime [46]. So, the SGWB bispectrum is specified by the f NL parameter and the angular spectrum. Also in this estimate we neglected a possible correlation between the initial and scalar source contributions that should be taken into account when, for instance, Γ I is controlled by the adiabatic scalar perturbation (see [31] for an example).

An example: the axion-inflation case
The goal of this section is to understand under which conditions the initial term Γ I (q) has a nontrivial q−dependence, that distinguishes it from the other contributions to the anisotropy. There are several mechanisms for the generation of a cosmological GW signal visible at interferometer scales (see [8][9][10] for recent review). In this section we focus on a specific mechanism: we consider the case where an axion inflaton φ sources gauge fields, which in turn generates a large GW background. In particular we consider the specific evolution shown in Figure 4 of [42], where the inflaton potential is chosen so to lead to a peak in the GW signal at LISA frequencies, without overproducing scalar perturbations and primordial black holes. The amount of GWs sourced in this mechanism is controlled by the parameter ξ ≡ (φ/2f φ H), where f φ is the decay constant of the axion inflaton. The present fractional energy in GW, Ω GW (η 0 , q), is related to the primordial GW power spectrum P λ (η in , q) by This relation, taken from [10], interpolates between large and small scales. Since we are interested in the modes with q q eq , that entered the horizon during radiation domination, we consider only the second term in the square bracket, and we find and, as we will see, the constant term is not relevant for our computation. We are interested in the contribution from the initial condition Γ in . So we can set the long modes ζ( k) = h λ (k) = 0 in this discussion. We therefore assume that the value of the energy density that arrives to the location x from the directionn is controlled by the parameter ξ =ξ + δξ ( x + dn) .
In this relation, ξ is the value that this parameter had during inflation at the location x + dn, where d is the distance covered by the gravitons between the initial and the present time (equal for all directions, since we are disregarding the effect of the long scale modes ζ; we note that these modes will contribute to the term Γ S , that we are not discussing in this section).
In writing this relation, we have assumed that the parameter ξ is in turn controlled by a dynamical field (the rolling axion, in the example of [42]), which results in the background valueξ, and in the perturbation δξ.
We then generalize the relation (5.1) to which has the background valueΩ GW (η 0 , q) = constant × λ P λ q,ξ . The constant factor drops in the ratio as well as in δ GW (η 0 , x, q,n) = P λ (q, ξ (η 0 , x,n)) − P λ q,ξ P λ q,ξ = ∂ ln λ lnP λ q,ξ ∂ξ δξ ( x + dn) , (5.6) where we have expanded the GW primordial power spectrum to linear order in δξ. In this wat, the relation (2.9) can be recast in the form where we have also made use of the standard definition of the tensor spectral tilt n T . The question of whether we have or have not spectral distortion depends on whether the quantity F(q,ξ) is or is not q−dependent. This provides an immediate criterion for evaluating whether and how much the GW anisotropies depend on frequency (as, in principle, one could imagine a GW power spectrum for which the dependence on q of F vanishes, or is extremely suppressed). This conclusion only assumes that the primordial GW signal is function of some additional parameter ξ which has small spatial inhomogeneities, and therefore it likely applies to several other mechanisms.
We show in Figure 1 the evolution of the function F corresponding to the GW production shown in Figure 4 of [42]. We see that indeed this quantity presents a nontrivial scale dependence, and therefore the correlators of the anisotropies will be different at different frequencies. f Hz ℱ(f) Figure 1. Quantity F as a function of the frequency f = q/2π of the GW signal for the model of axion inflation described in the text.

Squeezed limit and consistency relations of the SGWB
Non-linear effects associated with the propagation of interacting GWs in a non-linear universe lead to non-vanishing connected n-point functions even in absence of intrinsic, primordial non-Gaussianity. In particular, the squeezed limit of bispectra associated with GW observables should acquire a non-vanishing value, and satisfy consistency relations that resemble Maldacena's consistency relations [54]. This is analogous to what happens for CMB [55][56][57].
In this Section we compute the squeezed limit of the bispectrum for the graviton distribution function in the case of adiabatic fluctuations. As in Section 2, we write in momentum space ω GW η, k i , q, n j =ω GW (η, q) 1 + δ GW η, k i , q, n i , (6.1) whereω GW (η, q) is associated with the energy density of the isotropic SGWB. This quantity depends on time η and on the GW momentum q. Small anisotropies of the SGWB are controlled by the quantity δ GW given in Eq. (2.9). We re-write it here, expressing it in terms of the functionf (q): where recall that Γ S controls the fluctuations in the distribution function (see the definitions in Section 2). In this Section we focus on the contribution due to scalar fluctuations. We assume there is no anisotropic stress, and that scalar perturbations in Newtonian gauge satisfy the adiabaticity condition: where g(η) is a function mapping the superhorizon seed (controlled by ζ( k)) to the scalar fluctuations at small scales (see e.g. [61,62]). It is generally time dependent although it is equal to unity in pure matter domination. Then the contribution Γ S reads (see eq (3.9)) where µ =n ·k and T S is the definition for the scalar transfer function we adopt here. In matter domination this becomes Notice that Γ S does not depend on q. Assembling the definitions above, we can then write Indicating with P Γ the power spectrum, we can write the 2-point correlators in momentum space as where a prime corresponds to correlators understanding the (2π) 3 δ( k i ) factor. Then: In matter domination, as we learned above, |T S | 2 = 9/25, but in general |T S | 2 can depend on η, k,n.
In what follows, we study how the two-point correlation functions of SGWB anisotropies, when evaluated at small scales k, are modulated by the presence of a long-scale mode ζ L ≡ ζ( k L ), with | k L | | k|. Such modulation induces a non-vanishing squeezed limit for the three-point function of δ GW . The anisotropies δ GW depend on various quantities, (η, k, q, µ), which can be sensitive in a different way to the long mode. We use the systematic approach pionereed by Weinberg [27] that unambiguously associates the effects of a long mode with an appropriate coordinate transformation. We shall closely follow the treatment of [56], which develops the arguments of [27] for the case of CMB, applying it to the SGWB (for similar approaches see also [55,57]).

Long wavelength modes as coordinate transformations
We discuss how to identify the effects of a long mode with an appropriate coordinate transformation. We limit our attention to effects due to scalar fluctuations. The metric including long-wavelength scalars in Poisson gauge is (6.10) We assume that the long-scale mode depends on a momentum k L , with magnitude much smaller than that of the momentum of the short-scale modes introduced in eq. (6.17), but with a certain direction, and we discuss how the quantities (η, k, q, µ), transform under a coordinate redefinition adsorbing the long modes. We start by noticing that the following coordinate transformation preserves the Poisson gauge structure (ζ L indicates the long mode of curvature fluctuations at large scales, responsible for the modulation effects): 11) with λ constant. After performing such gauge transformation, we can 'gauge away' the long wavelength scalar modes making the gauge choice so that in the hat coordinates the metric is purely FRW with no long-wavelength perturbations. As explained in [27,56], in order to be consistent with the small k limit of Einstein equations, we need to impose the conditions (in absence of anisotropic stress) λ = 1 , where η * is some initial reference time. Eq (6.15) immediately leads to the equality = −2H + 1 . (6.16) After performing the coordinate redefinition (6.11), (6.12), we can write a metric containing short-wavelength scalar fluctuations in Poisson gauge 'on top' of long fluctuations: In fact, such metric contains the long-scale modes within the definition of the hat coordinates. We can then express the perturbations in terms of the original coordinates (η, x i ) using again relations (6.11), (6.12). Such operation teaches us how the short wavelength modes are modulated by the long wavelength ones: Importantly, the short modes acquire a second order correction due to long modes. As we shall discuss in what comes next, these non-linear, higher-order corrections modulate the 2-point function for short modes, and lead to a non-vanishing squeezed limit for the 3-point function.
As a concrete example, that we shall use in what follows, we can consider the case of constant proportionality between pressure and energy density, p = wρ. Being in this case a(η) ∝ η 2/(1+3w) , H = 2/[η(1 + 3w)] we get and which, for matter domination, gives H = 2/5. We also need to evaluate how the Fourier transform of a function f (x i ) changes under a rescaling of spatial coordinates, as in eq (6.12). We find that if we apply a constant rescaling of spatial coordinates to a function f , then its Fourier transform, given by transforms as (at first order in a ζ L expansion) This implies that up to first order in ζ L , under the coordinate transformation we are interested in, we have: As a last step, we now investigate how to transform the coordinates (q,n i ) that control the GW four-momentum. At first order, neglecting tensors, the GW four-momentum components are given by We wish to express the previous quantities in terms of hat coordinates, including the effects of the long modes. In particular, we are interested in determining the quantitiesq and n i that are contained into the GW four-momentum, when it is expressed in terms of hat coordinates. We use the fact that P µ is a vector, transforming in the usual way under coordinate transformations (in particular transformations (6.11), (6.12)). Using this fact, we findq Condition (6.26) gives, at first order in the long-scale modes, On the other hand, condition (6.27) giveŝ These are the results that we need. It is convenient to write more compact expressions aŝ with β q,n functions of time In matter domination we find β q = 2/5 and β n = 0.

Coordinate transformations and the GW distribution function
We now apply the previous results to the problem at hand. We start by re-writing the GW energy density whereω and δ GW (η, k, q, n i ) = − ∂ lnf (q) ∂ ln q Γ S η, k, q, n i . (6.34) We now transform each contribution in the previous formulas under the coordinate transformation discussed in Section 6.1. The background quantitiesω GW andf (q) transform as The quantity Γ S is mapped to that, expanded at linear order in ζ L , becomes We now assemble the results obtained. The SGWB energy density, including anisotropies, is modulated by the long mode ζ L as Eq. (6.39) is the basic expression that we need: all quantities at the RHS are evaluated in terms of the original coordinates without the hat. Notice that even in absence of intrinsic small-scale anisotropies, the GW energy density is modulated by the long mode: a dependence on ζ L is indeed still present by setting δ GW = 0 in eq (6.39). This is the effect studied by Alba and Maldacena [14]. For example, in pure matted domination, we have β q = H = 2/5. Setting δ GW = 0, eq (6.39) simply becomeŝ In this case, the modulation of ω GW is then controlled by the momentum dependence of the functionf (q), associated with the isotropic distribution function of the SGWB energy density [14].

The squeezed limit of 3-point correlation functions
We now apply the general result of (6.39) to study how correlation functions of small-scale GW anisotropies δ GW are influenced by the long mode. We start by studying how two-point correlation functions are modulated by ζ L ; we continue investigating the squeezed limit of the three-point correlation functions.
Using eq (6.39), we find the result 1 where the modulating factor M reads Notice that the contributions in the first line of eq (6.39) that depend only on the long mode (without being weighted by δ GW ) do not contribute to M. Therefore they do not modulate the short-mode two point function.
We now apply to the results derived above the definitions of δ GW and Γ power spectra, eqs (6.8), (6.9). We find the following expression for the modulation of the power spectrum due to a long mode: (6.43) All quantities inside the square parenthesis in the RHS are again evaluated at the same values of η,n, k; hence we understand this dependence. We find that the power spectrum of δ GW is modulated by the long mode ζ( k L ) through three (physically distinct) effects: 1. A modulation due to the scale dependence of the primordial curvature spectrum, as in Maldacena's consistency relation. This is contained in the first line of eq (6.43), second term in the RHS. (Notice that the contributions coming from derivatives of the 1/k 3 factor cancel out, as expected.) 2. A contribution due to the momentum-dependence of the background distributionf (q). This is contained in the first line of eq (6.43), third term in the RHS. This is a close relative of the effect pointed out by Alba and Maldacena [14], although it is not exactly the same result because we find contributions depending on second derivatives of the functionf (q).
3. A contribution due to the time, scale, and direction dependence of the transfer function of scalar modes. This is contained in the second line of Eq. (6.43).
In the previous discussion we learned how the long mode modulates the 2-point function. This effect is expected to lead to a non-vanishing squeezed limit for the 3-point function involving the anisotropies δ GW . Indeed, expressing a large scale limit of δ GW in terms of ζ aŝ for a small | k 3 |, we can write the schematic relation (all δ GW 's are evaluated at the same values of η, n i , q so we understand their dependence) where in the second line we used eq (6.41). This non-vanishing result gives the squeezed limit of the three-point function for δ GW . We adopt the following definition 2 for the non-linear parameter f δ GW In our case, using the previous results, we find and we can apply to this result the very same considerations made after eq (6.43).
The formula simplifies considerably in the case of pure matter domination. In this case, T S = 3/5, β q = H = 2/5. Then, Recalling thatf (q) is related with the GW isotropic energy density Ω GW by the relation  As an illustrative toy model which demonstrates this effect, we can consider a GW spectral density with the shape of a broken power law. The following parameterisation for the spectral energy changes slope at a scale q = q :

Conclusion
The amount of information extracted from the detection of GW signals by the LIGO-Virgo collaboration has shown the power of GW to study astrophysical compact object and to give relevant cosmological information on the late time universe. At the same level, the improving angular resolution of future GW detectors will allow one to extract precious information from the detection of the stochastic background of GWs generated both from the superposition of unresolved astrophysical sources and from cosmological sources, like inflation, phase transition or topological defects. However, high sensitivity alone will be not sufficient for discriminating among different contributions. So it becomes necessary to characterize such backgrounds using observables that can give a clear hint about the origin of the signals. As recently studied, a parity-violating SGWB, which represents a smoking gun for some cosmological signals, can be probed using ground and space-based interferometers [63]. Another important tool is the directionality dependence of the SGWB. As shown for astrophysical GW, the distribution of sources implies that the energy density is characterized by an anisotropic contribution beyond the isotropic one. In the same way we expect that, analogously to CMB photons, also primordial GW are charaterized by anisotropies that can be generated both at the moment of production and during their propagation. In this paper we focused on the stochastic background of cosmological origin and we studied the anisotropies due to the production mechanism (that we encode in an initial condition term) plus those generated from the propagation of GW on the perturbed universe, using a Boltzmann approach. We solved the Boltzmann equation for the graviton distribution function considering a FLRW metric with both scalar and tensor inhomogeneities. We showed that, contrary to CMB photons, at the moment of production, GWs, which are characterized by a non-thermal spectrum, generically result in angular anisotropies that have an order one dependence on the GW frequency. We provide a criterion to evaluate whether and how much the GW anisotropies depend on frequency. As an example, we evaluate this criterion in the case where an axion inflaton φ sources gauge fields, which in turn generates a large GW background. Additional anisotropies are induced by the GW propagation in the the large-scale scalar and tensor perturbations of the universe. We compute the angular power spectrum of the SGWB energy density, and, analogously to CMB photons, also the gravitons distribution function gets mainly affected by the Sachs-Wolfe effect on large scales, while the Integrated Sachs-Wolfe is subdominant.
We then focus on a second observable that can be a crucial tool in discriminating an astrophysical from a cosmological background, namely its departure from a Gaussian statistics. While we expect that the astrophysical background is Gaussian, due to central limit theorem, (some) cosmological backgrounds should shown a non-Gaussian statistics. We computed the three-point function (bispectrum) of the SGWB energy density, which is not affected by de-correlation issues, both considering the effects at generation and due to propagation. We have shown that also the SGWB bispectrum carries a memory of the initial condition and that it is proportional to the non-Gaussianity of the scalar perturbations. In this sense, the SGWB can be used as a novel probe (beyond the CMB and the LSS) of the non-Gaussianity of the scalar perturbations.
Finally we consider non-linear effects induced by long-wavelength scalar perturbations, which generate a modulation effect on the correlation functions of the short-wavelength modes. We identified the effects of long modes with an appropriate coordinate transformation and we computed the effect of non-linearities in inducing a non-vanishing squeezed limit of the SGWB three-point correlation function. We quantified the dependence of the squeezed bispectrum on the scale-dependence of the spectrum of primordial scalar fluctuations similar to Maldacena consistency relation, on the momentum dependence of the background SGWB distribution function, and on the time, scale, and direction dependence of the scalar transfer function.
In summary, in this paper we have approached the possibility to use CMB techniques to describe the cosmological SGWB trying to characterize it using peculiar features that we do not expect to have in the astrophysical background. Of course the detectability with interferometers of such effects is one crucial step to address and we plan to work on it on a future paper. At the same time we also plan to analyze several additional physical effects that we have neglected in this first paper, like the effects of neutrinos on the GW amplitude or a possible direct dependence of Γ I onn, which would give distinctive signatures useful for the characterization.
with the knowledge that, when k is decomposed according to (A.3) (namely, withk directeed along thee z−axis), We need to evaluate the integral (A.5) for a generic orientation of k. On the other hand, the explicit expression of the integrand (A.6), holds only when k is oriented along the z−axis. We cope with this by rotating the integrand of the d 2 Ω n integration into a basis in which the directionn is decomposed according to Eq. (A.3).
To achieve this, we introduce the rotation matrix S (Ω k ) ≡   cos θ k cos φ k − sin φ k sin θ k cos φ k cos θ k sin φ k cos φ k sin θ k sin φ k − sin θ k 0 cos θ k where the relation (4.12) has also been used. We collect some of the factors in this expression into the combinatioñ