Anisotropies and non-Gaussianity of the cosmological gravitational wave background

The stochastic gravitational wave background (SGWB) is expected to be a key observable for gravitational wave (GW) interferometry. Its detection will open a new window to early Universe cosmology and to the astrophysics of compact objects. Using a Boltzmann approach, we study the angular anisotropies of the GW energy density, which is an important tool to disentangle the different cosmological and astrophysical contributions to the SGWB. Anisotropies in the cosmological background are imprinted both at its production and by GW propagation through the large-scale scalar and tensor perturbations of the Universe. The first contribution is not present in the cosmic microwave background radiation (as the Universe is not transparent to photons before recombination), causing an order 1 dependence of the anisotropies on frequency. Moreover, we provide a new method to characterize the cosmological SGWB through its possible deviation from Gaussian statistics. In particular, the SGWB will become a new probe of the primordial non-Gaussianity of the large-scale cosmological perturbations.


I. INTRODUCTION
Operating ground-based interferometers is not far from reaching the sensitivity to detect the stochastic gravitational wave background (SGWB) from unresolved astrophysical sources [1,2]. On the other hand, future space-based gravitational wave detectors, like LISA [3] and DECIGO [4], and earth-based gravitational wave detectors, like the Einstein Telescope [5] and Cosmic Explorer [6], may be able to detect a stochastic background of cosmological origin, generated by early Universe mechanisms of production of gravitational waves (GWs) [7][8][9][10][11][12][13]. The most immediate way to differentiate the two backgrounds is by their frequency profiles [14]. However, given that the SGWB is the sum of different contributions whose profiles are not fully known, it is important to also develop other means to characterize them. In this work, we study the statistics of the angular anisotropies in the energy density of the GWs, which are either produced primordially or imprinted in the GWs as they propagate in the perturbed Universe [15][16][17][18][19][20].
This approach has several analogies with the well established formalism developed for cosmic microwave background (CMB) anisotropies. Following [16], we study, as is commonly done for the CMB, the GW phasespace distribution function f, which can immediately be related to their energy density. We solve the collisionless Boltzmann equation for this distribution and compute the two-point and three-point correlators of the GW energy density anisotropies on our sky. We focus on one crucial difference from the CMB: while the CMB temperature anisotropies are generated only at the last scattering surface, or afterward, the Universe is instead transparent to GWs at all energies below the Planck scale. Therefore, the SGWB provides a snapshot of the Universe right after inflation, and its anisotropies retain precious information about the primordial Universe and the mechanisms for the GW formation. In particular, the primordial signal may be characterized by a significant (i.e., order 1) dependence of the anisotropies on frequency. On the contrary, this dependence is very small in the CMB case since any initial condition is erased by collisions before recombination, and any frequency dependence of the anisotropies is generated only at second order in perturbations [21][22][23]. We show, through a representative example of sourced GWs during axion inflation [24,25], that a primordial GW signal visible at interferometer scales can indeed lead to anisotropies with a large frequency dependence.
Secondly, we study another important tool to characterize the cosmological SGWB-namely, its non-Gaussianity. Recent works, starting with [26], have investigated whether the GW three-point function hh 3 i can be tested at interferometers. The measurement of this signal requires the measurement of phase correlations of the GW wave functions. As shown in [27,28], two effects make such a measurement unfeasible: (i) the GW propagation in the perturbed Universe destroys any hh 3 i correlation possibly present in the primordial signal, and (ii) modes of nearby frequencies get confused with one another due to the finite duration of the experiment, also resulting in a large phase decorrelation. There is, however, another type of non-Gaussianity that can be observed, and it is the one present in the spatial distribution of the GW energy density. This does not involve initial phase correlation of the GW field itself: here we present the first steps for the computation and study of the three-point function (the bispectrum) of the GW energy density. Notice that Planck set the tightest limits from CMB data on deviations from Gaussian statistics for cosmological fluctuations [29]. Still, this does not rule out the possibility of primordial non-Gaussian signatures: the observable we discuss (the angular bispectrum of GW energy density anisotropies) relies on future measurements that might be sensitive enough to probe primordial non-Gaussian signals. It might be the case for models of inflation (like those where the inflaton interacts with additional fields), producing a GW signal with a peak at interferometric scales (see [9]) which would favor the measurement of the two-and three-point angular correlators.
For brevity reasons, this paper contains results under the simplest conditions only. In a companion paper [30], we shall present the details of these computations, extend them to include the GW propagation to second order in perturbations, and develop a more extended analysis of the GW bispectrum.

II. BOLTZMANN EQUATION FOR GWs
We consider a distribution f ¼ fðη; x i ; q;n i Þ of GWs as a function of their position x μ and momentum p μ ¼ dx μ =dλ, where λ is an affine parameter along the GW trajectory. This distribution obeys the Boltzmann equation L½f ¼ C½fðλÞ þ I½fðλÞ, where the Liouville term is L ≡ d=dλ, while C and I account, respectively, for the collision of GWs along their patch and for their emissivity from cosmological and astrophysical sources [16]. The collision among GWs affects the distribution at higher orders (in an expansion series in the gravitational strength 1=M planck ) with respect to the ones we are considering, and it can be disregarded. The emissivity may be due to astrophysical processes (such as black holes merging) in the relatively late Universe, as well as cosmological processes, such as inflation or phase transitions. In this work, we are interested only in the stochastic GW background of cosmological origin, so we treat the emissivity term as an initial condition on the GW distribution (see [31] and the references therein for a discussion of collisional effects involving gravitons). This leads us to study the free Boltzmann equation, df=dη ¼ 0, in the perturbed Universe. Specifically, we consider scalar (Φ and Ψ) and tensor (h ij ; taken to be transverse and traceless) perturbations in the so-called Poisson gauge, around a homogeneous and isotropic background, giving the line element where aðηÞ is the scale factor, and η is conformal time.
Dividing the free Boltzmann equation by p 0 leads to wheren ≡p is the direction of motion of the GWs, while q ≡ j ⃗pja is the comoving momentum that we use (as opposed to the physical one that was used in [16], following the standard computation done for the CMB photon propagation [32]), as it simplifies Eq. (3). The first two terms in Eq. (2) 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 redshifting of gravitons, including the Sachs-Wolfe (SW), integrated Sachs-Wolfe (ISW), and Rees-Sciama 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, in a way similar to CMB physics. Keeping only the terms up to first order in the perturbations, Eq. (2) gives In analogy to the split in Eq. (1), we also assume that the GW distribution has a dominant, homogeneous, and isotropic contribution, with distribution functionf, plus a subdominant contribution δf. The two functions are obtained by solving Eq. (3) at zeroth and first order in perturbations. Doing so, one immediately finds that any functionfðqÞ of the comoving momentum solves Eq. (3) at zeroth order. As a consequence, the associated number density n ∝ R d 3 pfðqÞ is diluted as a −3 as the Universe expands. This is also the case for CMB photons, whose distribution functionf CMB ¼ ðe p=T − 1Þ −1 is controlled only by the ratio p=T ∝ pa ¼ q, where T is the temperature of the CMB bath. This is a consequence of the free particle propagation in an expanding background, and it does not rely on the distribution being thermal.
The subdominant anisotropic component δf can be present as an initial condition. However, even if it is initially absent, Eq. (3) shows that this anisotropy is produced by the propagation of the isotropic component f in the perturbed background. Assuming that ∂f=∂q ≠ 0 (otherwise the solution of δf also becomes trivial), it is convenient to rescale the perturbed part of the distribution function as In this variable and in Fourier space, Eq. (3) gives where from now on prime will denote a derivative with respect to conformal time, μ is the cosine of the angle between ⃗ k andn, while the source function is S ¼ Ψ 0 − ikμΦ − 1 2 n i n j h 0 ij . As we now show, the quantity Γ can be immediately related to the anisotropic component of the GW energy density, ρ GW ≡ R d 3 ppf. It is customary to parametrize the GW energy density measured at the time η at the location ⃗ x in terms of its fractional contribution Ω GW through where ρ crit ¼ 3H 2 M 2 p is the critical energy density of the Universe, and H is the Hubble rate. Nearly all studies assume Ω GW to be homogeneous. Since we are interested in its inhomogeneous and anisotropic component, we have allowed Ω GW to depend on space. We account for the anisotropic dependence by defining ω GW through Ω GW ¼ R d 2n ω GW ðη; ⃗ x; q;nÞ=4π, and by introducing the density contrast δ GW ≡ δω GW ðη; ⃗ x; q;nÞ=ω GW ðη; qÞ. Using Eq. (4), one then finds that withΩ GW being the homogeneous, isotropic component of Ω GW .
In the CMB case, by inserting definition (4) into the Planck distribution and expanding to first order, one finds that Γ CMB ¼ δT=T. The main difference between the CMB and the GW case is that, before recombination, the collision term between photons and baryons suppresses any existing temperature anisotropy, thus removing any memory of the initial state. The observed temperature anisotropies δT=T arise since recombination, following an equation analogous to Eq. (5), with a source that, to first order, is independent from the energy of the CMB photons. While in the CMB this dependence arises only to second order in perturbations, a significantly greater dependence can be present in the GW distribution as an initial condition. In the following, we first compute and discuss the cosmological correlators of the GW anisotropies, then show through a concrete example that they can indeed have a significant dependence on frequency.

III. CORRELATORS OF GW ANISOTROPIES AND NON-GAUSSIANITY
As is standard [32], we express each of the sources appearing in Eq. (5) as a mode function times an initial variable that is constant at large scales, assuming for simplicity adiabatic scalar perturbations, and whose statistical properties have been set well before the propagation stage that we are considering (for instance during inflation, or during an early phase transition). Therefore, the scalar modes are (disregarding anisotropic stresses as, for example, those due to the relic neutrinos) Ψ ¼ Φ ≡ T Φ ðη; kÞ ζð ⃗ kÞ; we then decompose the tensor modes as h ij ≡ P λ¼AE2 e ij;λ ðkÞhðη; kÞξ λ ðk i Þ, where the sum is over rightand left-handed (respectively, λ ¼ AE2) circular polarizations, and the polarization operators are constructed as in [26]. We insert these expressions in the source function into Eq. (5) and solve for Γ. We then follow the treatment done for CMB perturbations, and we expand the solution in spherical harmonics, ΓðnÞ ¼ P l P l m¼−l Γ lm Y lm ðnÞ, where we recall thatn is the direction of motion of the GWs, and thus the direction at which the GWs arrive at our sky. The multipoles Γ lm are the sum of three contributions. The first contribution arises from the initial conditions, where η 0 denotes the present time, and we set our location to ⃗ x 0 ¼ 0. We also remark that this term in general depends on q. The second contribution is due to the scalar sources in Eq. (5), where the scalar transfer function T ð0Þ l is the sum of a term analogous to the SW effect for CMB photons, T Φ ðη in ; kÞ j l ½kðη 0 − η in Þ, and the analog of the ISW term, Finally, the third contribution Γ lm;T is due to the tensor modes in Eq. (5), and it is formally analogous to Eq. (9), with the productζY Ã lm replaced by the combination involving the spin-2 spherical harmonics, and with the scalar transfer function replaced by the tensor one T ðAE2Þ l ðk; η 0 ; η in Þ, given by We are interested in statistical correlators of the anisotropies. Under the assumption of statistical homogeneity and isotropy, the two-point and three-point correlators of ζ are expressed in terms of, respectively, the scalar power spectrum and bispectrum through hζð ⃗ kÞζ Ã ð ⃗ k 0 Þi 0 ¼ð2π 2 =k 3 Þ P ð0Þ ðkÞ and hζ 3 ð ⃗ k i Þi 0 ¼ B ð0Þ ðk i Þ [we use the standard notation of the prime to eliminate the momentum conservation Dirac delta and the ð2πÞ 3 coefficient]. Analogously, correlators P ðλÞ and B ðλÞ can also be defined for the two tensor polarizations. Moreover, we impose correlators of the same structure for the initial conditionsnamely, hΓðη in ; ⃗ k; qÞΓ Ã ðη in ; ⃗ k 0 ; qÞi 0 ¼ ð2π 2 =k 3 ÞP ðIÞ ðkÞand for the bispectrum B ðIÞ . In this work, we assume that the different contributions are uncorrelated. Under these assumptions, one obtains hΓ lm Γ Ã l 0 m 0 i ≡ δ ll 0 δ mm 0C l ¼ δ ll 0 δ mm 0 ½C l;I ðqÞ þC l;S þC l;T , where we denote the correlators with a tilde to distinguish them from the CMB case. The contribution from the initial condition readsC where again we stress the possible frequency dependence. The other two terms arẽ At large scales, this contribution is dominated by the term proportional to the initial value of Φ in T ð0Þ l (the analog of the SW contribution for the CMB). For modes that reenter the horizon during matter domination (as is the case for those that give the large-scale anisotropies that we are considering), T Φ ¼ 3=5 at early times [32]. 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 [32]. Therefore,C l;S ¼ ð10=3Þ 2 C SW l . If instead the initial condition term Γ I and the scalar propagation term Γ S are correlated (as for instance under the assumption that Γ I is controlled by the adiabatic scalar perturbation) then the sum of both terms should be compared to the SW for the CMB.
The structure of the bispectrum is forced by statistical isotropy to be the product of an l i -dependent term and Gaunt integrals [33], h Q 3 i¼1 Γ l i m i i ≡b l 1 l 2 l 3 G m 1 m 2 m 3 l 1 l 2 l 3 . The initial condition term leads tõ b l 1 l 2 l 3 ;I ¼ The scalar termb l 1 l 2 l 3 ;S is analogous, with the first spherical Bessel function replaced by the transfer function T ð0Þ l i , and with the scalar bispectrum B ð0Þ as the last term. In particular, for a primordial bispectrum of the local form [34], B ð0Þ ðk 1 ; k 2 ; k 3 Þ ¼ ð6f NL =5Þ½ð2π 2 Þ 2 = ðk 3 1 k 3 2 ÞP ð0Þ ðk 1 ÞP ð0Þ ðk 2 Þ þ 2 perm:, when applied to the CMB result [34,35], the same rescaling done after Eq. (13) gives the dominant SW contribution at large scales:b l 1 l 2 l 3 ;S ≃ 2f NL ½C l 1 ;SCl 2 ;S þ 2 perm:: Finally, the tensor term reads with the Wigner 3-j symbols being employed in defining For the case of purely adiabatic fluctuations, the formalism developed here allows us to determine consistency relations for the squeezed limit of the bispectrum B δ ð ⃗ Þi 0 . Such a squeezed limit is determined by nonlinear effects coupling long and short modes, and it can be computed using well-known techniques developed in the context of cosmic inflation [36] and the CMB [37][38][39]. Focusing on a matter dominated Universe, neglecting second-order tensor fluctuations, and considering for simplicity only isotropic contributions to the density contrast, we find that Hence the squeezed limit of the bispectrum, in this specific situation, depends on the tilt of the scalar fluctuations, and also on derivatives of the background distribution function fðqÞ, which is modulated by long modes. We plan to further explore this subject in [30].

IV. ANISOTROPIES IN THE PRIMORDIAL SGWB AND THEIR FREQUENCY DEPENDENCE
Several mechanisms for the generation of a cosmological GW signal visible at interferometer scales have been studied in the literature [8][9][10]. Here we comment on a specific mechanism where an axion inflaton ϕ sources gauge fields, which in turn generates a large GW background. 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 inflaton background value then results in a backgroundξ, and thus in a homogeneous and isotropic GW component. The inflaton fluctuations result in a perturbation δξ, and thus in the inhomogeneities of the primordial GW background. The anisotropy in the GW energy density arriving today at our location from a directionn is controlled by the value assumed by the parameter ξ during inflation at the position ⃗ x 0 þnd, where d is the distance traveled by the GWs from their production during inflation to today. GW modes observable at interferometers reentered the horizon during the radiation-dominated era. The present fractional energy density Ω GW of these modes is equal to their primordial power spectrum P GW times a q-independent factor. Then, by linearizing the primordial power spectrum in δξ, relation (7) can be recast in the form Γ I ðη 0 ; ⃗ x 0 ; q;nÞ ¼ F ðq;ξÞδξð⃗ x 0 þ dnÞ, with We have then provided 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 assumes only that the primordial GW signal is a function of some additional parameter ξ which has small spacial inhomogeneities, and therefore it likely applies to several other mechanisms. For axion inflation, we consider the specific evolution shown in Fig. 4 of [25], where the inflaton potential is chosen so as to lead to a peak in the GW signal at LISA frequencies, without overproducing scalar perturbations and primordial black holes. We show in Fig. 1 the corresponding evolution of the parameter F. We see that this quantity indeed presents a nontrivial scale dependence, and therefore the correlators of the anisotropies will be different at different frequencies.

V. FUTURE WORK
We plan to extend the results presented here to analyze several additional physical effects, including the effects of neutrinos on the GW amplitude [40], the possible direct dependence of Γ I onn, tests of nonstandard expansion in the early Universe, possible mixed bispectra among the three contributions to Γ that we have discussed, and the feasibility of measuring the frequency dependence of the two-point function and the bispectra at GW interferometers.