Many cosmological and astrophysical scenarios predict the existence of stochastic gravitational wave backgrounds (SGWBs) with a sufficiently large amplitude to be detectable at interferometer scales, see e.g. [1][2][3] for reviews. A detection of a SGWB signal would provide important information about its origin, hence it is essential to characterize its properties in detail, and to devise observables that will enable us to distinguish among different possible sources of gravitational waves (GWs).
The statistical features of a SGWB can offer various observables that can be measured with interferometers. While so far the attention has been mainly focussed on the two-point function (power spectrum) of the GW signal (see e.g. the recent review [4]), in this work we study for the first time how LISA [5] (and for this matter any other interferometer) can probe the three-point function (bispectrum) of the SGWB. A non-vanishing GW three-point function is associated with the non-Gaussian features of the SGWB. Some of the key ingredients of early universe models (in particular inflation and topological defects) predict that the primordial SGWB can have large non-Gaussianity, whose features can be calculated from first principles, and can be used to discriminate among different models. On the other hand, a SGWB due to a combination of a large number of uncorrelated astrophysical sources, or due to sub-horizon cosmological processes such as phase transitions or preheating, is Gaussian to a high degree, due to the central limit theorem (it remains to be studied to which degree the central limit theorem applies to the distribution of astrophysical sources that contribute to the stochastic GW background measured by LISA). Thus a detection of a non-Gaussian SGWB would be a distinctive indication of a signal of cosmological origin, and the study of its properties would provide crucial information on the physics of the very early universe.
The Fourier transform of a GW three-point function defines the non-Gaussian tensor bispectrum, which depends on the frequencies of the three GWs involved, as well as on their polarization. Its scale and chirality dependence can be very rich, in contrast to the GW two-point function, whose features are very constrained by the symmetries of the underlying background geometry. This implies that a study of the GW bispectrum can lead to a large number of new observables that can be used to differentiate among different models of cosmology. So far, the physics of primordial tensor fluctuations has been mainly investigated at Cosmic Microwave Background (CMB) scales, and indeed the non-Gaussianity of primordial GWs produced in certain models of inflation have been observationally constrained by the Planck satellite [6].
In this work, we investigate for the first time the capabilities of interferometers, in particular LISA, to test the non-Gaussian features of a SGWB via its three-point functions 1 . We do so by developing a formalism to compute the interferometer cubic response function to a non-Gaussian SGWB, connecting the GW bispectrum to the statistics of the actual signal outputs obtained from LISA (see Sections 2 and 3). We show that cubic correlators of the interferometer signal are sensitive to the properties of the GW bispectrum -its dependence on the GW wavenumber and chirality indexes -and also depend on specific features of the interferometer, like its arm configurations. Measurements of the cubic interferometer signal can distinguish among different cosmological models for a SGWB, and we show how different examples of primordial bispectra, with distinctive dependence on the momenta of the incoming GWs, lead to qualitatively distinguishable features in the interferometer cubic correlators. We also demonstrate that any measurement of the GW bispectrum at LISA is invariant under parity. Nonetheless, we show that in the presence of non-Gaussianity the analysis of LISA data in different frequency intervals can allow to extract nontrivial dependence on the correlators of the left-and right-handed gravitational waves, such as the relative amplitudes of the ( RRR + LLL ) and ( RRL + LLR ) correlation functions.
These results allow us to build the optimal estimator for the GW bispectrum observed at LISA, generalizing results previously developed for studying the GW two-point function [7] to the three-point function, see Section 4. Our final expression for the Signal-to-Noise-Ratio (SNR) is physically transparent, and we apply it to a specific scenario of a primordial non-Gaussian GW signal whose bispectrum is amplified at a specific scale to quantitatively demonstrate the ability of LISA to measure tensor non-Gaussianity, depending on the amplitude of the GW bispectrum.
Our results demonstrate that the statistics of the GW signal measured at interferometer is sensitive to various distinctive properties of the bispectrum of primordial GWs from the early universe: Section 5 surveys existing cosmological models capable of producing large non-Gaussianity of the SGWB, analyzing the features of the corresponding bispectra, and briefly discussing prospects of detection with the LISA interferometer and with CMB experiments, in light of our findings.
Section 6 contains our conclusions, with a summary of our results and suggestions for future studies, while six appendices contain technical results used in the main text.
In this section we build a formalism for describing tensor non-Gaussianity, which can be used to conveniently analyze how to probe this observable with interferometers. We consider a stochastic background of gravitational waves associated with the transverse-traceless metric perturbation h ab (t, x) of the background metric (where a, b correspond to spatial indexes, while λ denotes tensor polarization). We decompose the tensor modes as
where k ≡ | k|, k denotes a unit vector in the direction of the vector k. In this expression, we sum over the two transverse polarizations λ of a GW, with e ab,λ ( k) being the polarization operator for the polarization λ. For a more detailed discussion on the definition and construction of different polarization basis, and the properties they satisfy, we refer the reader to Appendices A, B.
As we require that h ab (t, x) is real, it must follow that
holds. This condition, however, does not fix completely the choice of the polarization tensors e ab,λ ( k). Two basis for the GW polarizations are commonly employed in the literature, the {+, ×} basis, and the {R, L} chiral basis. In the explicit computations performed in this paper, the chiral basis is used. Furthermore, we choose chiral polarization operators that specifically satisfy e * ab,λ ( k) = e ab,λ (-k) [we discuss the freedom of choosing different polarization operators in Appendices A, B]. This implies that the same property is also satisfied by the momentum space variable h λ (t, k), so that the reality of h ab (t, x) is ensured. Given our (2.3)
In particular, the second equality in the first expression states that e ab,R ↔ e ab,L under a parity transformation.
The statistical properties of the GW background are controlled by its correlation functions in Fourier space. In this work we consider the 2-point correlator h 2 , and, more in detail, the 3-point correlator h 3 , which is non-vanishing for a non-Gaussian SGWB. Assuming statistical isotropy, the equal-time momentum-space correlator is given by
where P λ (k) is the power spectrum of the helicity λ, and the numerical factor at the righthand side has been fixed imposing that the combination of eqs. (2.1) and (2.4) leads to the real space correlator
For studying the 3-point function, we use an ansatz analogous to the one used for describing the statistics of primordial scalar perturbations [8][9][10][11][12][13], see for example [14][15][16]. Specifically, we assume a small departure from Gaussianity, so that a tensor mode is the sum of a dominant Gaussian component, and its quadratic convolution
This ansatz is characterized by a kernel K λλ λ , which depends on the GW momenta and polarizations. We shall see that the kernel defines the properties of the non-Gaussian tensor bispectrum. We assume that, in general, the two different polarizations can be coupled in the convolution. The kernel is symmetric in the last two arguments (λ ↔ λ together with
The transformation of the kernel under a rotation is discussed in Appendix A. The dependence of the kernel on the magnitude of the three momenta controls the so-called shape of the non-Gaussianity [11], namely how the bispectrum changes according to different triangular configurations in Fourier space. The simplest form of non-Gaussianity is the so called local shape, enhanced in the squeezed limit of the bispectrum, for which (assuming also that the different helicities are not mixed in the convolution) K λλ λ k, p, q ∝ δ λλ δ λλ . More in general, we normalize K λλ λ (k p v1 , k p v2 , k p v3 ) = 1 for a reference unit triangle formed by v1 , v2 and v3 = -v 1 -v2 (see Eq. (3.25)) and some given pivot scale k p . In this way, the size of non-Gaussianity is controlled by the nonlinear parameter f λ,λ ,λ NL . In the concrete example of non scale-invariant non-Gaussianity that we study in Section 4.3, the pivot scale is chosen to be the scale at which the bispectrum is maximal, see eq. (4.22).
It is important to note that the mode functions appearing in the relation (2.6) are evaluated today. If the GW has a cosmological origin, we need to account for its evolution. It is conventional to encode this in a cosmological transfer function
where t 0 indicates the present time, and h pr is the primordial value of the GW mode. For an adiabatic tensor mode produced during inflation, h pr λ is constant (time independent) at super-horizon scales. For GW produced inside the horizon after inflation, we take the value of the mode at the end of GW production. Correspondingly, eq. (2.6) is changed into
× δ (3) k -p -q K λλ λ k, -p, -q .
(2.9)
This leads to the relation
between the parametrization of non-Gaussianity in terms of the primordial vs. the present day GW mode functions. 2 While the product of the last two factors on the rhs is more directly related to what we measure, the lhs is more immediately connected to the theory that provides the origin of the non-Gaussianity.
From the ansatz (2.6) we obtain the equal-time three-point function, to linear order in the nonlinear parameters, as
with
(2.12)
This expression for the bispectrum can describe the different shapes of non-Gaussianity, and the dependence on chirality. The expression for the tensor three-point function (2.11) evaluated at non-equal times, which will be used in Section 4, is discussed in Appendix D.
Our aim in this Section is to build the tools to study the primordial tensor three-point functions with LISA: we connect the theoretical results of the previous section with the actual LISA design. We analyze how the signal cubic correlator can be used to probe the tensor bispectrum described in Section 2, deriving for the first time the interferometer cubic response function to tensor non-Gaussianity. As we shall see, the interferometer response function is sensitive to the shape and polarization dependence of the primordial tensor bispectrum.
We first analyze how a SGWB influences the time a photon takes for traveling along a single arm of an interferometer. We then study the signal measured at LISA. We follow the derivation done in [4,17] based on the frequency basis (E1) that we introduce in the appendix E, but reformulated in the basis (2.1).
We consider an interferometer arm at rest, with mass 1 and mass 2 at the two ends. Mass 1 and mass 2 are located, respectively, at x 1 and at x 2 = x 1 + L l12 , where L is the length of the arm, and l12 the unit vector in the direction from mass 1 to mass 2. The change in the light-travel time for a photon emitted at mass 1 at the time t 1 and arriving at mass 2, due to a passing gravitational wave h ab is given by
Inserting the decomposition (2.1) and performing the integral, we find
with
where n is the direction of propagation of the wave. Let us consider now the signal s 12 (t, x 1 ), measured at time t by a detector of mass 1 located at x 1 , as the total change in the light-travel time for a photon emitted at mass 1 (at the time t -2L), arriving to mass 2 (at the time t -L), and then coming back to mass 1, i.e.,
where n 1 (t) = n 1 (t, x 1 ) is the noise measured by the detector at the vertex x 1 . Using eq. (3.2) we obtain
where the detector transfer functionfoot_2 is
as well as the 3-point correlation function
where we have assumed that signal and noise are uncorrelated (i.e., s n = s n 2 = s 2 n = 0).
Equation (3.5) describes the signal generated in a single arm of an interferometer. From this result we can construct the response functions of the full instrument, where the phase measurements in the individual arms are combined to minimize the instrumental noise (see e.g. [18]). Combining two arms of the equilateral triangular LISA configuration with a common mass at x 1 yields a Michelson interferometer:
Cyclic permutation of the endpoints x 1 , x 2 , x 3 results in a total of three Michelson interferometers, which we label X, Y and Z. Specifically
Defining the three arm directions as
Labeling of the LISA arms used in this work, cf. eq. (3.12).
(see Figure 1), the three measurements (3.11) acquire the form
with
From these, the standard LISA output channels A, E and T are constructed as [18]
We write these through the compact expression
where O = {A, E, T } labels the detector channel and the matrix c is given by
(3.17)
We now have the tools to study the two-point and three-point correlation functions for the signal, and analyze how they depend on the statistical properties of the SGWB, namely its power spectrum and bispectrum as described in Section 2. We derive the quadratic and cubic interferometer response functions, and discuss how the interferometer measurements allow to probe properties of the primordial tensor non-Gaussianity. The wavenumber k of the GW has been chosen to be 0.1/L. The response functions for the left-handed polarizations are obtained by a parity transformation about the plane of the detector (the xz plane, in our choice), which for our case corresponds to φ → -φ (as discussed in the text).
The quadratic auto-correlation of the interferometer signal reads
where
Here we use standard spherical coordinates, where θ ∈ [0, π] is the polar angle, and φ ∈ [0, 2π] is the azimuthal angle to express the direction of k. To obtain these expressions we have inserted eq. (2.4) and used the fact that Q i ab e ab,λ -k = Q i * ab e * ab,λ k [recall that this property is specific to our choice of polarization operators eq. (2.3)].
The quantity R OO λ (k) gives the (scale-dependent) detector response function to the GW of polarization λ for the 2-point correlation function of the channel O. This is obtained after
over all the directions k of the incoming GW.
The angle-dependent integrand is shown in Fig. 2, for the specific choice k = 0.1/L of the GW wavenumber and for the two channels A (left panel) and E (right panel), respectively. We disregard the T channel, as its sensitivity to the GW background is well below that of the A and E channels [18] (we have verified that this is the case also for the 3-point functions that we study below). For definiteness, we set the detector in the xz-plane, with the three masses, respectively at the {x, z} locations given by {0, 0} , {0, L}, and
(we verified that rotating the plane of the interferometer does not change the final value of the response function obtained after integrating over the angular directions of k). As expected GWs arriving from the perpendicular direction to the detector (the y axis, characterized by θ = π 2 and φ = π 2 , 3π 2 ) are those that provide the largest contribution to the detector twopoint function. For equal length interferometer arms the response function of the A -E cross-correlation R AE vanishes [18]. In Fig. 3 we show the detector response function (after integrating over the GW directions). The same result is obtained for both polarizations (we discuss this below) and in the two channels A and E. The response function is nearly constant up to k ∼ 0.2/L, and then it strongly decreases with higher k's. By Taylor expanding the relations (3.14) in the limit of small k, we obtain the analytic result
This agrees with the first of eqs. (18) in [19].
The response function for short wavelengths kL 0.1, contrary to the long wavelength limit, becomes more supressed the shorter the wavelength. Simple inspection of eq. (3.14) leads to predict that its envelope amplitude scales as ∝ (k * /k) 2 . This power law suppression at large k can be clearly appreciated in Fig. 3.
Starting from eq. (3.16), we obtain the three point correlation function of the measurement as
With respect to the two-point function, the computation of the three-point function has the additional complication that the bispectrum changes under a rotation of the triangle formed by the three vectors k i . We discuss this technical point in Appendix A, in particular see eq. (A19). To perform the angular integration, we fix the orientation of a reference triangle, formed by vectors k * i satisfying the momentum conservation condition k * 1 + k * 2 + k * 3 = 0, and compute the reference bispectrum B λ 1 ,λ 2 ,λ 3 k * i . For any generic orientation we have then
where R is a rotation matrix, and where Φ λ i k * i , R is the combination of phases given in eq. (A19). We note that this phase is eliminated by the corresponding opposite phase in the transformation of the helicity operators, cf. eqs. (A17) and (A18). Namely,
The disappearance of Φ from the product is due to the fact that the measurement is proportional to h ab h cd h ef , which is a product of three tensors.
Among the possible choices for the orientation of the reference vectors k * i , in what follows we adopt this one:
where
Starting from the configuration (3.24), a generic orientation of the bispectrum is obtained through the rotation matrix
The angles in this rotation matrix can be used to parametrize the independent angular integrals in eq. (3.21). Namely, we note that the angle between k 1 and k 2 is fixed by the fact that
As a rotation preserves angles between vectors, this angle has been imposed in eq. (3.25). The right matrix in eq. (3.26) then represents a rotation (by a generic angle φ 2 ) of k 2 around the axis of k 1 , at a fixed value of θ k 1 , k 2 . (We note that this matrix leaves v1 unchanged.)
The left matrix in eq. (3.26) then represents a common rotation of k 1 and k 2 . Therefore, the matrix R [θ 1 , φ 1 , φ 2 ] parametrizes the most generic orientation for k 1 and k 2 , compatible with the fact that the three wavevectors close into a triangle. No additional independent angle is required to parametrize k 3 , as this vector is simply -
With this in mind, eq. (3.21) can be rewritten as
where sin θ 1 is the Jacobian of the transformation from the integrals in eq. (3.21) to those in (3.28) (after the d 3 k 3 integration has been eliminated through the δ (3)
We then note that
where
Therefore, eq. (3.28) can be written as
with the three-point response function
We use this expression to numerically evaluate the response function. To verify that the response function is real, it is convenient to rewrite it as where the quantities γ are defined in (A15). One can verify by direct inspection that all terms in this integrand go to their complex conjugate under parity (namely, under
. Therefore, the response function is real. As we will discuss in Section 3.2.3, measurements of correlation functions by LISA are invariant under parity. This implies that R RRR = R LLL and R RRL = R LLR , for any choice of channels and momenta. For this reason, we show the response function in the RRR and RRL cases only. In the same Section, we will also demonstrate that the only nonvanishing 3-point response functions are the one among three E channels, and the one among two A channels and one E channel. We will show that they are opposite to each other, namely EEE = -AAE (and permutations).
We now use the expression (3.32) to evaluate the response function numerically. We do so for two cases:
-An equilateral configuration. This is the case of equal wavenumbers,
In Figure 4 we show the nonvanishing response functions for the RRR (left panel) and RRL (right panel) cases. The figure confirms that EEE = -AAE . We also see that the response functions vanish both in the kL 1 and kL 1 limit. In the large k limit, the same behavior is obtained for the two-point function, cf. Figure 3, and it is due to the inability of the interferometer to resolve scales much smaller than its size. The small k limit is instead studied in Appendix C.
-A squeezed isosceles configuration. For the isosceles squeezed case,
Starting from eq. (3.28), we find that, in this limit, the combinations R aA R bB e AB,λ 1 (v 1 ) and R cC R dD e CD,λ 2 (-v 1 ) can be expressed as e ∓2iλ 1 φ 2 times a factor that is independent of φ 2 , respectively. We further note that in this limit, the quantities Q i ab x i ; Ûp ; k 1 Rv 1 and Q j cd x j ; Ûq ; -k 1 Rv 1 are independent of φ 2 , and that
is also φ 2 -independent (the matrix c O k in the formula above is given in eq. (3.17). The φ 2 integration then gives (all other quantities are φ 2 independent, and factorize out of it)
namely we find that the result is nonvanishing only for λ 1 = λ 2 , in which case the λ 1,2 dependence drops, and it is also independent of λ 3 . Specifically, we find
The remaining quantities collect into
Apart from the F ef (θ 1 , φ 1 ) factor, the last two lines in this expression are the two-point response function R OO λ times 4π. This is a consequence of the fact that the squeezed bispectrum is nothing but the power spectrum of the short wavelength modes modulated by the long wavelength mode. Numerical evaluation of this expression results in the response functions shown in Figure 5. As for the equilateral case, we find the only nonvanishing response functions shown in the figure, and we find that the response function for two A channels and one E channel is of equal magnitude but opposite sign compared to the one for three E channels. The reason for this is explained in the next Subsection.
In this Subsection we describe how the geometrical features of the LISA configuration are sufficient to characterise some of the properties of the interferometer three point function. . Squeezed isosceles bispectrum, for λ 1 = λ 2 ≡ λ, and λ 3 ≡ λ (the same result is obtained for all helicities). The third channel is taken in the squeezed limit,
Besides the dependence on the shape and orientation of the triangle formed by the three momenta, an additional feature that can characterize the tensor bispectrum is its behaviour under a parity transformation. Certain models of early universe cosmology (see Section 5), predict tensor bispectra which violate parity, and values of the bispectrum components which are not invariant under a parity transformation. We now show that LISA -and any planar interferometer -cannot distinguish among components of the bispectrum that differ only by a parity transformation, and consequently cannot detect parity breaking effects. Reference [7] noted that the LISA response function for the two-point function is the same for the left and right polarizations of the GW, namely R OO L = R OO R , due to a mirror symmetry with respect to the plane of the interferometer. As we show now, this result can be readily extended also to the three-point function, implying the equalities R RRR = R LLL and R RRL = R LLR . To understand this, consider a triangular configuration of test masses in the xz plane. The mirror symmetry then corresponds to changing the y component of a vector. We would like to understand what this implies for the response functions. Denoting by V and by V ⊥ the component of a vector V , respectively, within and perpendicular to xz plane, we note that the GW momentum, and the basis vectors expressed by eqs. (A4) change as
As the unit vectors in eq. (3.12) (expressing the displacements between the test masses) lie in the plane of the instrument, only the parallel components contribute to the Q i ab e ab,λ contraction. Due to this, and thanks to the properties of (3.39), we see that Ûa Ûb e ab,λ goes to Ûa Ûb e ab,-λ under this mirror transformation (and analogously for the contractions with V and Ŵ ). On the other hand, the transfer function T and the phase in the first term in front of the square brackets in eqs. (3.14) are invariant under this transformation. This implies that e ab,λ Q j ab k , k ⊥ = e ab,-λ Q j ab k , -k ⊥ in eq. (3.13), and analogously for the linear combinations in eq. (3.15). These combinations form the response functions, which, therefore, satisfy
(3.40) After integrating over the angles, this implies the equalities R RR = R LL for the two-point function found in [7], as well as the equalities R RRR = R LLL and R RRL = R RLL for the three-point function. This implies that parity violation in the tensor bispectrum can not be detected using a planar interferometer, like LISA.
In eqs. (3.11) we described how the three Michelson interferometer signals X, Y , Z associated with the interferometer vertices are constructed. Since LISA forms an equilateral triangle, we should expect that physical results are independent on the vertex labelling, and should be invariant exchanging vertices. For example, exchanging the second and third vertices, one finds
and physical results are invariant under this transformation. Analogously, exchanging the first and second vertex, one has
and lastly, exchanging first and third vertex, one finds invariance of results under the simultaneous transformations
We can derive some consequences of these relations for the bispectra. First, the previous formulas imply that
Moreover, only one component of the three-point function is independent, say σ X σ 2 Y , and all the remaining ones coincide with this correlator up to a sign. We can write these relations as (we use the notation:
This also implies that three-point functions computed in terms of E, A channels are particularly simple to express. One finds
while all the other three-point functions vanish. We note that the identity
A is confirmed by the numerical computations shown in Figures 4 and5. Among the vanishing correlators, we note that Σ 3 A = 0. We also note that all the bispectra involving the T channel vanish. For example
On the one hand, this prediction for the null-channel, together with the relation (3.46) provides a non-trivial consistency check of the underlying assumptions of our analysis, such as statistical isotropy and non-gaussian instrument noise. On the other hand, these relations predict that for any given frequency interval, (cross-)correlating different interferometer channels only allows us to extract a single measurement of the bispectrum, a serious obstacle in measuring its different helicity contributions. However, with mild assumptions on the frequency dependence of the bispectrum, the different frequency dependences of the R RRR = R LLL and R RRL = R LLR response functions (see Fig. 4) can be used to measure the chirality structure by comparing measurements in different frequency intervals.
4
The optimal signal-to-noise ratio
In this section we construct a frequency-dependent estimator of the stochastic gravitational wave bispectrum. Our procedure extends the arguments developed in [7] to estimate the optimal signal-to-noise ratio in the gravitational wave power spectrum to the case of a signal three-point function. In most of this section we assume for generality that all the correlators AAA , AAE , AEE and EEE can be nonvanishing and independent from each other; only at the end of Section 4.2 we specialize to the case of an equilateral LISA configuration, where only the correlators EEE = -AAE are nonvanishing.
Since our goal is to generate a frequency-dependent estimator, our first task is to compute the three-point function of the signal in frequency space, defined as (for a return flight along the first arm)
where ñ12 (f ) denotes the corresponding noise in frequency space. Defining
where now O denotes a LISA channel, we obtain the following three-point function
the calculation of which requires the knowledge of the time correlator for the gravitational waves at unequal times. The expression for the latter quantity is given in Appendix D.
Using eq. (D5) we obtain
so that we can express the result in terms of quantities evaluated for a reference closed triangle with sides k * i , as discussed in Section 3.2.2 (see in particular eq. (3.31))
where
is the signal three-point response function computed in Section 3.2.2.
Following [7] we define a frequency-dependent estimator for the three-point function as
where the filter function W ijk (f 1 , f 2 , f 3 ) is totally symmetric under simultaneous permutations of both the indices and the arguments,
and satisfies the reality condition
with expectation value
Under the assumptions that the noise is Gaussianfoot_3 (so that its three-point function vanishes) and uncorrelated with the signal, one has si (f 1 ) sj (f 2 ) sk (f 3 ) = Σi (f 1 ) Σj (f 2 ) Σk (f 3 ) , where Σi (f ) is defined in eq. (4.2). This implies that the expectation value of the estimator is
where we have defined the real quantity
We next compute the variance of F assuming that the signal is noise dominated, with
(for brevity, from now on we omit the absolute value in the frequency dependence of the noise) so that
where the factors of 6 and 9 originate from the symmetry properties of W ijk . The signal-tonoise ratio (SNR) is given by F / F2 .
We note that for any given pair of objects A ijk (f 1 , f 2 , f 3 ) and B ijk (f 1 , f 2 , f 3 ) with the properties of our filter function W ijk we can define a scalar product
so that we can write the SNR as
In writing this formula we make the hypothesis that W ijk (0, f 2 , f 3 ) = 0, whose validity will be checked in short, which implies that ijk df
, the second term in the scalar product (4.13) does not contribute to the numerator in eq. (4.14). The SNR is thus maximized for
up to a multiplicative constant that cancels out from the expression of SNR. This last equation shows that our hypothesis W ijk (0, f 2 , f 3 ) = 0 is valid, since the noise diverges at low frequencies.
To sum up, for this optimal estimator the SNR is given by
where T is the duration of the experiment, the indices i, j and k can take only the values A and E and the sum will contain 8 terms. This relation simplifies in the case in which the channels involved in the correlation have the same noise. The noises of the measurement (3.10) at the vertex X and at the two other vertices Y and Z of the instrument satisfy
where, due to symmetry 5 ,
The linear combinations (3.15) have diagonal noise (4.11), with spectral dependence
which shows that the A and E channels have indeed identical noise. The factor (2L) 2 has been inserted in the definition of P n (f ) so to convert from time displacement to strain, by dividing ∆T by the round-trip light travel distance 2L [20]. In this way, P n (f ) has the dimension of time.
Using these two channels, the relation (4.16) can be explicitly written as
5
The cross-correlation arises because the different interferometers share one common arm.
where we used the results of Section 3.
As an explicit example of application of these results, in this subsection we evaluate the SNR (4.21) for a specific shape of non-Gaussian bispectrum. We make the following choice of kernel entering in eq. (2.6), depending on a reference closed triangle k * i and a reference scale k:
where we assume that σ k, and normalize the shape function to one at its maximum. This localized, chiral bump in non-Gaussianity well approximates the predictions of certain models of early universe cosmology, for instance the signal obtained in the model of [21] that we will review in the subsection 5.2. For later convenience it is also useful to introduce the energy density of the universe per logarithmic wavenumber interval. Using the expression
ḣij ḣij for the energy density in GW, one finds
where H 0 3.24 × 10 -18 h Hz is the present value of the Hubble rate.
At this point we can use the kernel function of Eq. (4.22) and the expressions for the SNR we developed in Section 4. We relegate technical steps of the calculations to Appendix F, and write here the final expression for the SNR for an equilateral configuration:
where k * i denote the unit vector along the direction of k * i and equilateral transfer functions are plotted in the left panel of Fig. 4.
For reference we consider the optimal scale kL = 0.028 at which the ratio
is maximum (for the expressions of the noise power spectrum P n see [20]). For this value we find |R EEE LLL | 0.025 and thus substituting into the previous expression, we find
where in the second line we normalize the various quantities to useful reference values, to be able to more easily appreciate the relevance of the result. For definiteness, we normalize k to 3.4 × 10 -3 Hz, which is the frequency corresponding to k; we then normalize the noise functions P n (k) and P n (2k) to 7 × 10 -41 Hz -1 and 3.9 × 10 -41 Hz -1 respectively, which are (parametrically) the values obtained at k (see for instance Fig. 3 of [20], where the square root of P n is shown; we note that the noise power spectrum P n is related to the sensitivity S n by P n = R 2 S n , where R 2 is the two-point response function shown in Fig. 3.). We normalize the time to the nominal LISA mission time of 4 years times a 75% duty factor. As discussed in Section 2, the value of the non-linear parameter entering in this relation is the present one, related to the primordial one by eq. (2.10). A detailed discussion of the transfer functions entering in this relation can be found for instance in [22]. The precise behavior is not needed for the present estimate. When they are well outside the horizon, the GW have constant amplitude; instead the amplitude decreases as the inverse power of the scale factor while it is well inside the horizon. For the present estimate can simply take
We recall that, in the case of GWs produced during inflation, a k * is the value of the scale factor when the mode of frequency k * re-entered the horizon in the Friedmann-Lemaître-Robertson-Walker (FLRW) stage after inflation (we normalize the scale factor to one today). The same relation applies for GW modes that are produced when their size is comparable to the horizon side during radiation domination. Assuming that the universe is radiation dominated at that moment, and that radiation domination continues until the recent stage matter-radiation equality (at z eq ∼ 3, 400) we find (see Appendix G)
Notice for instance that for f LLL,primordial NL = 10 3 this equation implies that a detection of non-Gaussianities can only take place for Ω GW k h 2 10 -11 . If however, Ω GW k h 2 10 -13 this implies that a detection of non-Gaussianities requires at least an f LLL,primordial NL 10 7 . As we discuss in Appendix G, modes that have presently the frequency k = 3.4 × 10 -3 Hz reentered the horizon at the temperature of ∼ 50 TeV. A different value for a * can be obtained if we make instead the unconventional hypothesis that the universe had a different equation of state than radiation, before, but close, to Big-Bang Nucleosynthesis.
A different value for a * is instead obtained if we assume that GW are produced inside the horizon by a sudden episode that took place during radiation domination at the temperature T * . In this case we find (see Appendix G)
The previous result is recovered if we take T * = 50 TeV.
In both these relations, Ω GW k h 2 has been normalized to the threshold level that can be detected in the 2-point correlation function at LISA. Finally, we recall that these results assume kL = 0.028 Hz, as this scale maximizes the SNR for a flat spectrum Ω GW (k). In order to evaluate the SNR at arbitrary scale, all the k dependent quantities (in particular |R EEE LLL |, and Ω GW (k) if it is scale dependent) must be evaluated consistently.
In the previous sections, we developed techniques for investigating the non-Gaussianity of a SGWB by studying the three-point correlation function of signals detectable with LISA. We have shown that LISA is in principle able to measure specific properties of the SGWB bispectra, as the dependence on the amplitude of the momentum and polarization of the GW signal.
In this section we review and discuss the current theoretical understanding of non-Gaussian features of a primordial SGWB sourced by Early Universe physics, in particular inflation and cosmological defects. Our aim is to provide theoretical motivations for the results derived in the previous Sections, showing that they can be used to distinguish different sources for primordial SGWBs. Besides analyzing models, we also briefly review current and perspective constraints on tensor non-Gaussianity from the physics of CMB, that is able to probe primordial SGWBs at frequency scales much smaller than interferometers.
Cosmological inflation predicts the existence of a stochastic background of tensor modes, produced by quantum fluctuations of the metric spin-2 degrees freedom during the phase of inflationary expansion. CMB experiments constrain the amplitude of the primordial SGWB power spectrum at large CMB scales in terms of the tensor-to-scalar ratio r: the current upper bound from BICEP2/Keck and Planck is r < 0.07 at 95% confidence level [23,24] (assuming the consistency relation r = -8n T ), and future CMB polarization experiments [25][26][27] can lower this bound down to around 10 -3 in absence of a detection.
Tensor self-interactions are expected to make the non-Gaussianity of the primordial SGWB relatively large, already within Einstein gravity. Moreover, as we shall review next, there are models of inflation which exploit specific couplings between fields present during the inflationary era, in order to enhance the amplitude of the tensor power spectrum at scales that can be probed by the future generation of gravitational interferometers. This gives the opportunity to probe inflation at scales much smaller than CMB scales: in these scenarios, the large couplings among the fields involved can enhance tensor non-Gaussianity, making it a very useful observable for distinguishing among different scenarios. We briefly survey the topic of inflationary tensor non-Gaussianity, considering in succession models of increasing complexity.
The simplest situation to investigate -as a toy model for inflation -is pure General Relativity (GR) in de Sitter space. The metric for de Sitter space can be expressed as ds 2 = -dt 2 + e 2Ht d x 2 , with H the constant Hubble parameter. As we have seen in Section 2, the transversetraceless spin-2 tensor fluctuations can be decomposed in two helicity modes, λ = L (left) and λ = R (right). The Einstein-Hilbert action
can be straightforwardly expanded around the de Sitter background up to third order in fluctuations [12,28]. From the second order action one can compute the two-point function for the Fourier transform of the tensor fluctuations, which defines the primordial tensor power spectrum
From the third order action one obtains the three-point function for the fluctuations [12,28], which defines the primordial tensor bispectrum
and can be obtained using for example the in-in formalism. The tensor bispectrum is the lowest order statistics providing information on non-Gaussianity of tensor fluctuations. Its dependence on the three wave-vectors characterizes the shape of tensor non-Gaussianity, one of the properties that allows one to distinguish among different models. We list here some properties of the tensor bispectrum and the corresponding tensor non-Gaussianities for GR in de Sitter space:
1. The amplitude of tensor non-Gaussianity -parameterized by the ratio between the tensor bispectrum and the square of the tensor power spectrum -are of order one in pure GR (but not larger). For example, one finds for equilateral configurations [29] |B RRR | P 2 R 3.6 .
(5.5) This is a difference with respect to scalar curvature fluctuations in single field vanilla models of inflation, where scalar non-Gaussianities are suppressed by slow-roll parameters, and are then at most of order a few percent [8,9,12,13]. This is due to the fact that Einstein gravity is a non-linear theory, and cubic interactions are not suppressed with respect to quadratic ones by small (e.g. slow-roll) parameters. In Einstein gravity the shape of tensor bispectra is peaked in squeezed configurations.
2. The amplitude of tensor bispectra depends on chirality, that is on the polarization indexes L, R, and in general bispectra characterized by different indexes have different amplitudes. For example, in the case of GR in de Sitter space one finds [29,30]
This is different with respect to the power spectrum, where L and R modes have the same amplitude (see eq. ( 5.3)) and the cross correlation vanishes, h L h R = 0. This has interesting phenomenological consequences since we have many independent bispectrum components we can use to build observables for distinguishing among different models of inflation. As we have learned in Section 3, LISA can in principle probe different chiral components of primordial bispectra.
3. Pure Einstein gravity around de Sitter space preserves parity, in the sense that the amplitudes of power spectrum and bispectrum components obtained interchanging all the L, R indexes is the same. For example,
etc. This is potentially not true in more complex inflationary scenarios, as we are going to discuss in the following subsection.
We can go beyond a pure Einstein-Hilbert action, and investigate non-Gaussian tensor fluctuations around de Sitter space in theories of gravity including higher order curvature invariants [28,31,32], schematically denoted as W 3 and W W 2 (W µνρσ being the Weyl tensor). Such perturbative contributions modify the third-order tensor action, and lead to new, parity preserving shapes for the tensor bispectra. However, their amplitudes can not be (much) larger than the ones one finds within General Relativity, due to unitarity constraints.
An enhanced amplitude of tensor bispectra can be obtained in single field inflation, with the time variable controlled by a scalar field, the inflaton. The power spectrum of tensor modes in single-field slow-roll models of inflation is straightforward to compute, and one finds the same result as in de Sitter space, but with an additional scale dependence. In the simplest models of inflation the tensor spectrum is red, with a spectral tilt given by
with r the tensor-to-scalar ratio 6 . The analysis of the tensor bispectrum is instead less straightforward: the scalar can have non-minimal derivative couplings with the metric -see e.g. the model [34] built in terms of the Horndeski theory -which can complicate the analysis. The complete classification of tensor non-Gaussianity in single field models of inflation based on Horndeski scalar-tensor theories of gravity has been done in [30], finding that besides the Einstein-Hilbert part [12], the third order action for tensor modes acquires a parity-preserving contribution proportional to the cube of time derivatives of tensor fluctuations:
where G is a function of the scale factor and of the homogeneous profile of the scalar field and its first time derivatives (it is related with the function G 5 controlling the quintic Horndeski action). Interestingly, there are no obvious unitarity constraints on the size of this correction to the Einstein-Hilbert action, hence -depending on the inflationary model and the profile for the scalar field -this contribution can lead to a tensor bispectrum whose amplitude is parametrically larger than the GR result (although it might be hard to find an explicit model satisfying all the CMB constraints in the scalar sector). The shape of the tensor bispectrum has been investigated in [30], finding that it is maximized for equilateral configurations, hence it is different from the GR one, maximized in the squeezed limit. In the equilateral limit, the bispectrum components associated with action (5.8) satisfy the relation
(5.9) This is different from the GR result of eq (5.6): we learn that the chiral structure of the tensor bispectra -i.e. a distinctive hierarchy for the amplitudes of the tensor bispectra depending on the polarization indexes L, R -might be used to build observables to distinguish among different models.
Tensor non-Gaussianity in single field models of inflation based on generalization of Horndeski actions [35][36][37][38][39][40] is a topic still under investigation, see e.g. [41], and also [42,43] for an analysis carried on using techniques based on effective field theory for inflation. Also, tensor non-Gaussianity in parity-violating scalar-tensor theories [44,45] that spontaneously break Lorentz invariance [46] is particularly interesting to motivate a search of parity violating effects in the tensor bispectrum: an analysis of the tensor bispectrum in parity violating gravitational theories can be found in [31,32,47], although their generalization to the extended set-up introduced in [46] is still an open question.
Since the simplest single field models of inflation predict a red tensor tilt, the amplitude of the primordial SGWB is too small to be detected with interferometers. On the other hand, there are more complex inflationary scenarios that might allow one to enhance the tensor spectrum at interferometer scales [33,48], and that may lead to large tensor non-Gaussianity with distinctive features. We briefly discuss here two frameworks, whose predictions for the tensor bispectrum have been studied so far:
If other sources of GWs are present during inflation, the primordial SGWB can have richer features. A possibility to source primordial GWs is to couple fields driving inflation with additional scalars [49][50][51][52][53][54], U (1) gauge vectors [55][56][57][58], non-Abelian vector fields [59][60][61][62][63][64][65][66][67][68][69], or Standard Model fields [70]. These mechanisms usually exploit instabilities for the additional source fields during inflation. Such instabilities feed in the evolution of tensor modes through higher order contributions to the anisotropic stress, and affect both power spectra and bispectra of fluctuations: see e.g. [21,[71][72][73]. Given the nature of the couplings, the resulting tensor power spectrum and bispectra can be parity violating, and the bispectra are usually enhanced in equilateral configurations. The tensor power spectrum profile can acquire a feature (a 'bump' profile) and can also get sufficiently enhanced at small scales to be detected with interferometers [48]. Models involving couplings between (pseudo)scalars and gauge fields are theoretically well understood, and their observational prospects in CMB polarization and interferometer experiments are being developed in great detail, in particular for what respect the parity violating features of power spectra and bispectra, see e.g. the recent papers [66,74] and references therein.
Another possibility for enhancing the bispectrum of tensor modes is to break space-time symmetries during inflation. One way to do so are scenarios of (super)solid inflation: see e.g. [75][76][77][78][79][80][81][82][83][84], where the vacuum expectation value of additional scalar fields spontaneously breaks space diffeomorphisms during inflation, or more in general in models where spatial diffeomorphisms are broken spontaneously during inflation. Tensor bispectra are maximized in squeezed configurations, and their amplitude can be parametrically larger than in Einstein gravity. In the scenarios explored so far in the literature, power spectrum and bispectrum components are all parity preserving.
Alternatively, space-time symmetries can be broken explicitly as for example a violation of Lorentz symmetry in Hořava-Lifshitz gravity: in this case inflationary tensor fluctuations are automatically chiral [85] and the tensor bispectra can be chiral and enhanced with respect to GR results [86,87].
For these systems, some additional mechanisms (as growth of perturbations due to instabilities) might be necessary to amplify the tensor spectrum at interferometer scales.
After the general survey of the previous subsection, we describe in detail the predictions of a concrete model able to produce a signal detectable at interferometers with large tensor non-Gaussianity. We assume that during inflation, which is driven by a scalar φ, an auxiliary axion field χ interacts with a gauge field A µ : the Lagrangian density for the system is
where 1/f is a coupling constant with the dimension of a length, and F µν ( F µν ) is the field (dual) strength tensor. GW non-Gaussianities in different helicity channels become larger than those of the scalar perturbations by different orders of magnitudes. This can be understood easily noticing that both scalar and tensors are sourced with similar efficiency (and source) but the scalars are suppressed by the helicity conservation, while the tensors can be large and also chiral, meaning that the contribution of the different correlator are different for different combinations of the helicity.
In such a model the gauge field A µ acts as a source for both scalar and tensor metric perturbations, and, for the latter, this translates into a source term in the equation of motion of
where
∆ and where T EM lm represents the spatial part of the stress-energy tensor of the gauge field. The source term arises not from the F F term but from the F F term in eq. (5.10). Parityviolating information of the gauge field is transmitted through this source term to the tensors that acquire a chiral component. The GW solution for such equation can be found using the Green function method, see [88,89].
Using the same procedure as described in [55,88] we can compute both the two-point function and the three-point function at super-horizon scales (kτ → 0) for the tensor modes. The GW two-point function receives contributions from the parity conserving amplification of vacuum fluctuations and from the excited electromagnetic modes which source the parity violating parts. These two contributions are uncorrelated and the overall right-and lefthanded power spectra are
P l e 4πξ ξ 6 , (5.12)
where the parameter ξ encodes the velocity of the field χ, ξ = χ/(2f H). For the three-point function, at late times (kτ → 0), we have
where the function F lm λ ( k, q) is the same as in [88]. The largest contribution to the three-point function of the GW, given by the h R h R h R correlator, was already computed in [88,90], where was also shown that the shape of such a bispectrum peaks in the equilateral limit
Here we extend such computations to estimate the amplitude of the other mixed contributions to the GW bispectrum, which are
(5.14)
Note that the signs and reality condition are specific to our choice of chiral polarization operators obeying the properties in eq. ( 2.3). We find that the mixed helicity components are non-zero and in Section 4 we have seen how LISA is sensitive to such components. The shape of the GW three-point functions is close to equilateral and this can be explained by the fact that the mechanism of generation of perturbations happens at sub-horizon scales.
To put our investigation in a wider context, we discuss in this subsection the current and perspective constraints on tensor non-Gaussianity from CMB physics. CMB observations at large scales impose stringent constraints on the non-Gaussianity of scalar curvature (density) perturbations, and have the potential to constrain tensor non-Gaussianity as well, but at frequencies much smaller than interferometer scales. The tightest constraints come from the Planck collaboration [6], setting strong limits on scalar non-Gaussianity for a variety of wellmotivated models. For example, for the "standard" local, equilateral and orthogonal shapes the analysis of Planck temperature (T) and E-mode polarization data constrain the level of these types of non-Gaussianity to f local N L = 0.8 ± 5.0, f equil N L = -4 ± 43 and f ortho N L = -26 ± 21 (68% CL, statistical) from the measurements of the angular CMB bispectrum. These values point into the direction of consistency with the predictions of the standard cosmological model based on single-field models of slow-roll inflation, meaning that the structure that we observe today have been sourced by (almost) Gaussian seed perturbations. However there are also much less studied observables that could open up a new window into the physics of the early universe in the near future. Among these there are non-Gaussian signatures coming from tensor (gravitational waves) 3-point correlations, and also non-Gaussianity coming from mixed (tensor-scalar) 3-point correlators. Interesting bounds on these kind of signals already exist. A first observational limit on the tensor non-linearity parameter f tens N L has been obtained by Planck (from temperature data): f tensor N L /10 2 = 4 ± 15 (68%CL), which is consistent with a previous WMAP data analysis [91]. The tensor bispectra analyzed derive from some models of inflation predicting a specific parity violation in the tensor sector. However, future observations of CMB will try to improve constraints on tensor fluctuations focusing on B-mode polarization. The reason is that in this case primordial tensor modes are not hidden by scalar perturbations. The main contamination to the B mode polarization signal is polarized dust [92] and lensing of the E-modes to B-modes [93]. However, precise B-mode measurements have just started and the prospects of improvement are large. For example the constraints on the BT T CMB bispectrum arising from mixed correlators tensor-scalarscalar can improve by an order of magnitude with a CMB-Stage IV mission w.r.t the ones achievable from temperatura data alone (see, e.g., [94]). Observational constraints on such a primordial tensor-scalar-scalar bispectrum already exist and have been obtained recently in [95] using WMAP temperature data finding a constraint g tss = -48 ± 28 (68%CL) on the amplitude of such interactions. The inclusion of polarization information will certainly improve such a constraint. On the other hand it might be crucial to have the possibility to probe such signals on different scales, like the ones which are probed by direct measurements from interferometers. This might serve not only as a possible consistency check, but also to test models of inflation that eventually produce a relevant bispectrum signature only at interferometric scales.
So far in Section 5, we have considered GW backgrounds from inflation which can lead to a significant deviation from Gaussian statistics. Inflation, however, is not the only source of GWs from the early Universe. Various post-inflationary mechanisms can be also responsible for the generation of GW backgrounds having today a large amplitude. Non-linear field dynamical processes after inflation can generate non-trivial quadrupolar field distributions, resulting in an active and efficient source of GWs. The most representative processes are non-perturbative particle production [96][97][98][99][100][101][102][103], strong first order phase transitions [104][105][106][107][108][109][110] and cosmic defect networks [111][112][113][114][115][116][117]. For a recent review on early Universe GW sources see [3].
The stochastic background of GWs from strong first order phase transitions is expected to exhibit a single "bump" spectrum, with a rather large amplitude, depending on the ef-ficiency of the dynamical processes involved in the emission of GWs. In particular, if the electroweak phase transition is sufficiently strong due to the presence of beyond the Standard Model particle physics, the peak frequency of the associated GW background lies naturally within the LISA frequency window. For a discussion on the ability of LISA to measure a GW background from first order phase transitions (and in particular from the electroweak phase transition) see [3,109]. The stochastic background of GWs from non-perturbative phenomena like periodic [102,118], tachyonic [100,101] or other [103,119,120] excitation of fields (typically expected in preheating), is also characterized by a single "bump" spectrum, with an amplitude that can be typically large but peaked at very high frequencies (unless extreme fine tuned couplings are considered). In both cases, phase transitions and preheating, the production of GWs is due to a causal process. This implies that today's background is formed by the superposition of billions of independent signals (corresponding to the GWs emitted from the many uncorrelated regions at the time of the background generation), and hence it must have necessarily a Gaussian distribution, due to the central limit theorem. The angular resolution needed to probe any of the individual signals, beyond their effective stochastic Gaussian nature, is far beyond the reach of any reasonable GW detector. GW backgrounds from preheating and phase transitions are therefore expected to be highly Gaussian, and hence their 3-point function is expected to vanish at the data of any interferometer. For a discussion on this aspect, see Section 3.1 of Ref. [3].
The case of GW emission from cosmic defect networks is however very different. The GW signal expected from cosmic defects can be the sum of two components. The first component is the emission of GWs produced around the horizon at each time t, by the anisotropic stress of the network [116,[121][122][123][124]. This first component is expected to be emitted by any network of cosmic defects in scaling, independently of the topology and origin of the defects [116]. The second contribution is only expected for a network of cosmic strings -local gauge one-dimensional topological defects -and is given by the superposition of GWs emitted from sub-horizon strings chopped off from the string network all along cosmic history, as well as from super-horizon strings with superimposed small-scale structure [112] due to string interconnections leading to kink formation. For a cosmic string network, this second contribution is expected to be the dominant one. However, the argument of causality discussed above, remains valid, even if the GWs are continuously emitted during several Hubble times. Thus, this background cannot be resolved beyond its stochastic nature, and is expected to be Gaussian 7 . We therefore expect that any non-Gaussianity in the continuous stochastic background sourced by a defect network, can only be due to the first contribution (even if this is sub-dominant in the case of cosmic strings, in terms of the power spectrum).
Let us thus focus on the first contribution mentioned above, produced by the anisotropic stress of a defect network. This contribution represents an irreducible emission of GWs sourced by any type of viable defect network that reached scaling, meaning that the energy density 7 In reality, on top of the continuous stochastic Gaussian background from cosmic strings, there can be individual bursts emitted by nearby strings or a "popcorn" discontinuous noise. The former results from the fact that at small redshifts there are not many sources and events do not overlap, thus the time interval between events is longer than the duration of a single event. The latter is a noise generated by rare bursts at high redshifts leading to unresolved GW sources. Considering the number of sources at a given frequency as a Poisson process, one may thus expect either to get a signal or to get a superposition of signals, leading to what is called a popcorn-like noise [125]. These signals due to bursts represent therefore, in a sense, a temporal deviation from Gaussianity, that can be measured from the two-point function. However they do not correspond to the type of non-Gaussianity that we discuss in this paper, as they do not form a continuous stochastic background.
in the defect network is constant fraction of the overall energy density of the universe. Such defects can be either one-dimensional topological defects -global, Abelian, non-Abelian, semi-local strings -or (viable) global defects -domain walls, monopoles or non-topological textures [116]. Indeed, in the case of local strings, scaling is reached by the emission of gravitational waves. Global defects reach scaling mainly through long-range interactions and emission of Goldstone bosons [126]. A defect network in scaling produces an irreducible background of GWs emitted continuously around the horizon scale at every moment of the cosmic history. Its spectrum is predicted to be exactly scale-invariant for the modes that are emitted during radiation domination [116]. At the level of the power spectrum, this background mimics therefore the shape of the inflationary background due to quantum fluctuations. The irreducible emission of GWs from a defect network is however expected to be highly non-Gaussian. This is simply due to the fact that the source of the GWs is the (transverse-traceless part of the) energy-momentum of the network, and this is a function bilinear in the amplitude (modulo derivatives) of the fields that the cosmic defects are made of. Since the GW source is bilinear in field amplitudes, this implies automatically that any correlator of an odd number of tensor perturbations will be characterized by the correlation of an even product of fields. This is non-vanishing even if the fields were Gaussian. Thus, odd tensor correlation functions are non-vanishing, i.e. the GW background is not Gaussian. As the GWs are continuously emitted around the horizon scale at every moment, we cannot apply now the previous causal argument based on the superposition of many domains from the past.
Even though the shape of the power spectrum of the irreducible GW background from defects is well understood theoretically, its ultimate amplitude depends on the fine details of the so called unequal-time-correlator of the network's energy-momentum tensor. Unfortunately, this correlator that can only be obtained accurately from sufficiently fine lattice simulations of defect networks. It is therefore difficult to assess at this point whether this background can be detectable with LISA, much less whether its 3-point function can be measured. The GW signal can be however estimated analytically in a simplified case, known as the large N limit of a global phase transition due to the spontaneous symmetry breaking of O(N ) into O(N -1). This kind of phase transition leads to the formation of so called "selfordering scalar field" configurations, which correspond to non-topological defects (textures). The scale-invariant GW power spectrum due to the dynamics of such global defects has been estimated, in the limit N 1, by Refs. [122,123], whereas the 3-point function (in the equilateral configuration) has been presented in Ref. [124]. Order of magnitude calculations in the large N limit have found a GW bispectrum peaked in the equilateral configuration [124]
where P h is the total power spectrum (summing over the two polarizations) and N 1 is the number of components of the symmetry-breaking field. Even though this is a rough estimate, it clearly indicates that we should expect a large departure from Gaussianity for the irreducible GW background from any defect network. A proper assessment of the ability of LISA to detect the power spectrum and bispectrum of this stochastic background, requires however further work; namely lattice simulations of defect networks with a large dynamical range.
Let us finally comment on string gas cosmology [127,128] a string theory motivated early universe scenario, which does not involve a period of cosmological inflation. The model is based upon T-duality, making use of (fundamental) string oscillatory modes and winding modes. In this scenario, the universe may have started in a a quasi-static Hagedorn phase, during which thermal fluctuations of closed (fundamental) strings generate the density perturbations. The obtained spectrum of cosmological perturbations [129], and of GWs [130] are nearly scale invariant. The power spectrum of scalar metric fluctuations has a slight red tilt, while the spectrum of gravitational waves has a slight blue tilt. The string gas model has one free parameter and one free function, the former being the ratio of the string to the Planck length, and the latter the wavenumber dependence of the temperature. In this scenario, the spectrum of cosmological fluctuations may have large non-Gaussianities due to the thermal origin of the initial perturbations leading to strong three-point correlations [131]. Such thermal effects have been shown to be important for some inflationary models, such as chain inflation [132] or warm inflation [133]. In string gas cosmology, the non-Gaussianities of the spectrum of cosmological fluctuations depend linearly on the wavenumber, while the amplitude depends sensitively on the string scale. In slow-roll inflation, to leading order in perturbation theory, matter fluctuations do not couple to tensor ones, whereas in string gas cosmology matter fluctuations induce both scalar and tensor fluctuations. Hence, one may expect also non-Gaussianities in the tensor modes, for which the formalism discussed in the previous sections may be applicable.
In this work we studied the three-point correlation function of the SGWB expected to be measured by LISA, and we developed the formalism required for this analysis. The threepoint correlation function is a key observable to gain information on the statistics of this background and crucial in order to study its departure from Gaussianity. This can be an important discriminant for a signal of cosmological origin, since the stochastic background due to a large number of uncorrelated astrophysical sources is Gaussian to a very high degree, due to the central limit theorem. The three-point correlator (also known as bispectrum), if present in the data, is a richer observable than the standard two-point correlator (the power spectrum), due to several reasons. Firstly, this is due to the frequency of the three modes present in the correlator. The wave-vectors of the modes entering in the correlator must add up to zero, namely k + k = 0 for the power spectrum, and k + k + k = 0 for the bispectrum, implying that the momentum vectors form a closed triangle. For a statistically isotropic SGWB, the power spectrum only depends on the signal wavenumber (i.e. | k| = | k |), while the bispectrum depends on the momentum orientation (see eq. (A19)), on an overall scaling, and on two relative sizes. This last dependence is known as the "shape" of the bispectrum. Secondly, additional combinations of LISA channels can be measured using the three-point signal, some of which may be useful as null tests for non-Gaussianity from unresolved astrophysical point sources. Thirdly, the planar nature of the instrument (together with the assumption of statistical isotropy) results in the insensitivity of the twopoint function to the polarization of the two GW modes. This situation is ameliorated for the three-point function, since more combinations of GW polarizations can be considered (see below).
In eq. (2.6) we provided a very general parametrization for the case of a mildly non-Gaussian SGWB. This ansatz is borrowed from the wide literature on the non-Gaussianity of the primordial density perturbations probed through the fluctuations of the CMB radiation. This gives rise to a bispectrum of rather general shape (namely, the dependence on the wavenumbers of the three GW involved in the correlation, and, ultimately, on the frequencies of the three measured signals), see eq. (2.12). Different cosmological mechanisms produce non-Gaussianity of different shapes, so that this study can potentially be a powerful discriminant between them. In Section 5 we provided a survey of several possibilities present in the literature.
In Section 3 we computed all the non-vanishing three-point LISA response functions R OO O λλ;λ k, k , k which, together with the GW correlation function h λ ( k)h λ ( k )h λ ( k ) , provide the signal s O s O s O measured by the different LISA channels (we study the response for all possible combinations of the A and E channels). We provided the formulae to compute the response functions for generic wavenumbers, and we then provided the explicit results for two cases of particular relevance (as they probe distinct physical mechanisms for the generation of non-Gaussianity): the equilateral configuration k = k = k and the isosceles squeezed configuration k k = k . As we already mentioned, the planarity of the instrument and the fact that it has equal arms results in identities between different response functions, in the case of an isotropic SGWB. It is most immediate to discuss this in terms of circular polarizations, as they change into each other under a parity transformation. In the 2-point functions we have R OO LL = R OO RR (where L and R stand, respectively, for the left handed and the right handed circular polarization), resulting effectively in an inability to probe a chiral GW signal. On the other hand, for the three-point function we have
RRR , and
RRL , but these two objects have different frequency dependences. Hence by comparing measurements of the signal three-point function in different frequency ranges, one can hope to discriminate between a chiral and a non-chiral SGWB.
Moreover, in Section 4 we constructed the frequency-dependent estimator for the SGWB bispectrum that maximizes the signal-to-noise ratio in the measurement. As an example, we evaluated the general formula for the estimator for the specific case of a non-Gaussian signal narrowly localized at one given scale, k = k = k = k * .
To summarize, we have provided the complete formalism to compute the three-point correlation function of the signal measured at LISA for any given theoretical SGWB threepoint correlator. This formalism can be readily extended to any other combination of GW interferometers. Moreover, if the primordial bispectrum is amplified on a squeezed shapehence coupling modes of different frequencies -we can use the three-point function to correlate signals from different experiments probing distinct frequency ranges, as LISA, LIGO, or PTA. This can allow us to 'break' the planarity condition of a single interferometer, and measure effects of parity violation also in small frequency ranges that cannot be probed by other means (see also the discussion in Ref. [7]). If non-vanishing, measurements of three-point correlation functions can provide a wealth of information on different cosmological mechanisms, due to its dependence on the scale, the shape, and, possibly, the polarization of the SGWB. Even a null measurement can set limits on specific models that would be otherwise be unconstrained.
Let us first of all reproduce for clarity the tensor decomposition given in Eq. (2.1),
expressed in terms of an arbitrary polarization basis e ab,λ ( k). A basis for the polarization operators often considered in the literature is the {+, ×} basis. To define it one can introduce an auxiliary unit vector, e.g. in the z-direction of our frame of coordinates, i.e. êz , and introduce the orthogonal unit vectors
The set { k, û, v} forms an orthonormal basis. Note that k and û are odd under parity, while v is even. In standard spherical coordinates k = (sin θ cos φ, sin θ sin φ, cos θ) , (A3) the unit vectors û and v acquire the form û = (sin φ, -cos φ, 0) , v = k × û = (cos θ cos φ, cos θ sin φ, -sin θ) .
Choosing an auxiliary unit vector different from êz does not change any physical results, but induces a rotation of û, v around k, hence transforming the polarization operators and the GW mode functions [1]. We discuss this in appendix B. In this appendix (like in the bulk of the text) we stick to identifying the auxiliary vector with êz , and then study how the polarization operators, and the mode functions, change under a rotation of the wave vector k. These two different transformations, changing êz vs. changing k, produce a similar effect on the polarizations operators and on the mode function, as can be understood from the structure of (A2).
Out of these vectors, it is customary to define and they satisfy the properties (2.3) given in the main text. We are interested in studying how the polarization operators change under a rotation. As we will see they acquire an additional factor with respect to the transformation of a rank 2 tensor, which is due to the fact that the fixed vector êz used in (A2) does not change under the rotation. To see this, it is convenient to rewrite e (R)
To find how the polarization operators transform under a rotation, we will now study the transformation properties of a,L/R (this is facilitated by the fact that they have one spacial index, rather than two). Specifically, we want to obtain the quantities R -1 û R k and R -1 v R k , where R represents a spatial rotation. To do so, we compute the coordinates of these two quantities in the k, û k , v k basis:
We note that
and analogously for v. The various coordinates c uu , . . . , c vk can be then computed via scalar products, such as
It is convenient to define the metric
on the plane transverse to k, and the product with respect to this metric
In terms of this product, the expression (A10) acquires the compact form
Computing the other coordinates in an analogous way, we find that c uk = c vk = 0 (namely, the rotated û and v remain orthogonal to the rotated k, as one would expect), and that
where
Combining this with (A7), we can write
where λ = +1 denotes the R helicity, while λ = -1 denotes the L one. This leads to
As it is clear from eq. (A1), the product between the polarization operator and the mode function is a rank two tensor. We then learn that, under the rotation, the GW mode function transforms as
We note from (A15) that γ -k, R = -γ k, R . As a consequence, the power spectrum P λ (k) is invariant under a rotation, since the two modes involved in the correlation have momenta k = -k, cf. eq. (2.4). This is not the case for the bispectrum. Denoting by
An identical relation is satisfied by the Kernel function
As mentioned in the previous appendix, the decomposition in eq. (A1) [c.f. eq. (2.1)], depends on the choice of orientation of {û, v} within the plane transverse to k. For a given GW mode propagating in the direction k, our choice of the orthonormal basis {û, v} in eq. (A4), based on choosing êz as an auxiliar fixed unit vector, was actually arbitrary (even though very convenient, as we will see). Our canonical choice of {û, v} can be actually rotated by an arbitrary angle α around k, into a new triad {û(α), v(α), k}, where û(0) ≡ û, v(0) ≡ v are given by eq. (A4). This can be implemented explicitly by a rotation matrix around k as
from where it is easy to verify that |û(α)| 2 = |v(α)| 2 = 1 and û(α) × v(α) = k, as it should.
The rotated triad {û(α), v(α), k} form therefore an equally valid orthonormal basis, from which to build up new +, × polarizations tensors as
and from here a set of chiral polarization tensors in the rotated basis,
Basic algebra leads to relate the new chiral polarization tensors with respect to the canonical 'unrotated' basis {û, v} eq. (A4), as
where the chiral polarization tensors on the right hand side are those from eq. (A6), and on the left hand side those from eq. (B3). This implies that if we were to express the GW in eq. (A1) [c.f. eq. (2.1)] in terms of the new e (B5)
These transformation rules simply reflect, actually, the spin-2 nature of the tensor field h ij representing the GWs. Equations (B4), (B5) are relevant as they determine the complex phase of the response function we have calculated in Section 3. Let recall that the response functions presented along the main text were always computed for a tensor basis obeying the property expressed in eq. (2.3), namely that under a "parity" tranformation k ←→ -k, the chirality polarization tensors transform as e ab,λ (-k) = e * ab,λ ( k). That property is certainly true when the polarization tensors are computed in the canonical basis {û, v} of eq. (A4). However Eq. (B4) implies that such relations are not universal: on the contrary, they are only verified by a very small subset of orthonormal basis {û(α), v(α)}. In particular, from eq. (B4), and using eq. ( 2.3), we deduce immediately that, in general, under a flip of momentum, the rotated basis of polarization tensors transform as
As expected, for α = 0 we recover the result expressed in eq. ( 2.3) that that used in the main text, i.e. e ab,λ (-k) = e * ab,λ ( k). We see that this result also holds for α = π, which is nothing else but the flip of the canonical vectors as û -→ -û, v -→ -v. For any angle α for which α mod(π) = 0, the property e ab,λ (-k) = e * ab,λ ( k) will not hold. For instance for α = π/2, this is equivalent to û -→ v and v -→ -û, and hence it holds that e ab,λ (-k) = -e * ab,λ ( k). In general, the result of flipping the momentum will introduce a phase, as indicated in eq. (B6).
It is clear then that the choice of a 'rotated' vector basis other than the canonical one(s), will determine the complex phase of the polarization tensor via eq. (B4), and hence this will have necessarily an impact on the final expression of the 3-point response functions that we have derived in Section 3. The reason is that R λλ λ is always proportional to three polarization tensors, schematically ∼ e λ * * e λ * * e λ * * . Denoting as R λλ λ [α 1 , α 2 , α 3 ] the response functions calculated in a orthonormal vector basis rotated an angle α i around each vector k i , with respect the canonical basis {û, v} eq. (A4), we deduce the following relation(s)
where the response functions R λλ λ [0] on the right hand side should be identified with the response functions we presented in the main text. Due to eq. (B7), we see that we can have pure real, imaginary or simply complex response functions, simply depending on the choice of the phases α i . The fact that our 3-point response function is real applies because of the choice of polarization operators built from the canonical basis given in eq. (A4). For example, let us take all equal angles, α i = α ∀i = 1, 2, 3, and consider α = 0 or α = π. Then we find that all 3-point response functions are real (like in the main text), while for α = π/6 the response functions R RRR , R LLL are still real (but different in a phase), while R LRR , R RLL become complex (and differing in a phase). The condition to be purely real, purely imaginary, or simply complex, depends on the choice of α i 's according to eq. (B7). Of course, the final 3-point signal Eq. (3.31) [c.f. Eqs. (3.21,3.28)] in the GW detector does not depend on the choice of the û, v basis. The left hand side of eq. (3.31) is the same no matter the choice of α. However, depending on α, one must modify our expressions for the response functions by the phases indicated in eq. (B7), and at the same time correct the tensor bispectra B λλ λ by exactly the opposite phases, so that at the end the physical signal eq. (3.31) remains invariant.
In order to use our formalism, special care must be put into what vector basis is being used in the computation of the tensor bispectrum, and hence what properties of the polarization operators and response functions hold.
C Vanishing of the 3-point response function in the kL → 0 limit
In this Appendix we show that the LISA three-point response function vanishes in the long wavelength limit Lk → 0. In this limit, the three quantities in eq. (3.14) simplify to
We can then write the expressions for the three measurements in a compact way as
and
We recall that the index O = 1 refers to the A channel, while O = 2 to the E channel. We instead find Σ T = 0 (so α 3 = β 3 = γ 3 = 0).
This gives the response function
The quantity contracted with D O ab D O cd D O ef in the rhs of the above expression, carries 6 spatial indices, none of which can coincide with the spatial indices of any of the ki directions (otherwise the response function could be different for different orientations of the LISA arms). In fact, the direct evaluations confirms that
where only the number N λ 1 ,λ 2 ,λ 3 depends on the orientation of the reference plane used to evaluate the bispectrum. For the choice given by eqs. (3.24) and (3.25) we obtain N RRR = -9π 2 35 and N RRL = -π 2 35 . The structure in the curly parenthesis is compatible with the GW origin of the signal, and it vanishes when we trace over the ab, the cd, or the ef pair. Taking into account the fact that D O aa = 0, this expression simplifies to
For brevity, we denote the sum of the two traces as T OO O . Using the values in (C3), this expression gives
Finally, recalling the definitions (3.12) of the arm directions, we have Û1
As we anticipated, this shows (through a fully analytical computation) that the threepoint response function vanishes in the large wavelength limit, in agreement with what obtained numerically in Figures 4 and5.
To provide the expressions for the 2-point and 3-point correlators, we introduce a standard Fourier space decomposition
assuming that h λ k is a stochastic variable that satisfies
where statistical isotropy has been assumed, and where we allow for the possibility of a chiral GW background, P R (k) = P L (k). Then the two-point function at non-equal times reads
To compute the three-point function at non-equal times, recall our ansatz (2.6) for small departures from a Gaussian stochastic background, h λ t, k = h λ,g t, k + λ ,λ f λ,λ ,λ NL d 3 p d 3 q h λ ,g (t, p) h λ ,g (t, q) × δ (3) k -p -q K λλ λ k, -p, -q , (D4) with the properties of kernel function given in Eq. (2.7). From the ansatz (D4), and the 2-point correlator (D3), one then obtains the unequal time 3-point correlation function
(D5)
It is instructive to compare the GW decomposition (2.1) used in this paper with the decomposition in terms of positive and negative frequencies that is often encountered in the literature of the stochastic GW background [134] h ab (t, x) = where we are considering left and right polarizations, and where the requirement of a real GW background is used to define the negative frequency field, h λ (-f, n) = h * -λ (f, n) (in the {+, ×} basis, one has instead h σ (-f, n) = h * σ (f, n)). Using this property, the decomposition (E1) can be rewritten as where we note that we also sent {n, λ} → {-n, -λ} in the second term of the second line. Comparing this with the first line of (2.1) allows to relate the operators in the two decompositions as
Using the relation between the two lines of (2.1) we can also relate the operator in (E1) to the full GW Fourier transform as
where dot denotes time differentiation. We can use this relation to also relate correlators between modes in the two basis. For instance, for the 2-point correlator we have
Eq. (2.4) in the main text provides the equal time correlator in the Fourier basis. This is obtained from
These relations lead to the unequal time correlator
(which immediately reduces to (2.4) at equal times). Inserting this relation in (E5) leads to
As a check on our algebra, we verified that combining this result with the decomposition (E1) leads again to the equal-time real-space correlator (2.5).
In this Appendix we provide some technical steps that were omitted in the computation of Section 4.3. Inserting the kernel of eq (4.22) relation into eqs. (4.4) and (4.10), we obtain
The curly bracket has three terms. In each term, two dk i integrals are immediately performed thanks to the δ-functions. The third integral is performed in the limit of narrow width of the bump, in which the functions multiplying the exponential term are evaluated at k. The third integral is then also immediately done by extending it from -∞ to ∞ (which is also appropriate in the narrow width approximation). We obtain
We need to square this quantity, and insert it in (4.16). The mixed products between the three different terms in the curly bracket provide a negligible contribution to this result. Let us for instance discuss the product between the first and the second term. In the narrow peak approximation, the first term is mostly supported at |f 2 | = |f 3 | = k, while the second term at |f 1 | = |f 3 | = k. So the mixed product is supported only around the points where the magnitude of all the frequencies is equal to k * . But this is incompatible with the δ-function δ (f 1 + f 2 + f 3 ) present in (4.16).
Therefore, for the purpose of computing the signal to noise ratio, only the sum of the square of the three terms in the curly bracket is relevant. In this sum we also approximate all functions multiplying the exponential terms by their value at the center of the bump. This leads to We perform the df 2 and the df 3 integrals in the narrow width limit. The second line of this expression then becomes
The contribution proportional to δ (f 1 ) can be disregarded, since the noise diverges in that limit. This leads to eq. (4.24) written in the main text.
The parameterization (2.6) of non-Gaussianity in terms of a nonlinear parameter is subject to the cosmological evolution, giving rise to the relation (2.10). When the momenta in the convolution have equal magnitude, we have a single rescaling, that can be roughly approximated as
where k * is the magnitude of the momentum and a k * denotes the value of the scale factor when the modes re-entered the horizon. For GW generated after inflation inside the horizon, a k * denotes instead the value of the scale factor at the moment of the GW production. In this relation, the scale factor is normalized to one at the present time. In this appendix we estimate the value of a k * for both cases. We start from GW generated during inflation. In this case, the modes re-enter the horizon during radiation domination. We can write
a eq a k * = 1 a eq t eq t * = 1 a eq a * H * a eq H eq , (G2)
where we have used the fact that a ∝ t 1/2 and H ∝ 1 t during radiation domination. In this relation, the suffix 'eq' refers to the moment of matter-radiation equality, while the star denotes the moment at which the mode k * re-enters the horizon.
At horizon re-entry, a * H * = 2πk * , where the 2π is a consequence of the Fourier transform convention adopted in this work (see for instance eq. (2.1)). The Hubble rate scales as the square root of the energy density. At equality, the energy density was twice that of radiation, ρ eq = 2 × ρ rad,eq = 2 × a -4 eq ρ rad,0 = 2 × a -4 eq Ω rad,0 ρ 0 , where ρ 0 is the present value of the energy density. Therefore,
eq Ω rad,0 .
In scaling the energy density of radiation as a -4 we have disregarded the fact that massive neutrinos become non-relativistic, so we consistently evaluate the present fractional density in radiation as if neutrinos were massless, leading to Ω rad,0 4.18 × 10 -5 h -2 . This leads to Let us instead assume that the GW are produced inside the horizon, when the temperature of the universe was T * . For convenience, we normalize the temperature at 100 GeV. We then have (disregarding the variation of relativistic degrees of freedom for the purpose of our estimate)
, production inside horizon at temperature T * .
(G5)
Here we are interested in probing the intrinsic non-Gaussian statistics of the SGWB, and not in describing 'pop-corn-like' non-Gaussian bursts sometimes present in the unresolved stochastic background due to combining many individual sources. See[4] for a detailed review on searches for this kind of non-Gaussian signatures using the SGWB two-point function.
To be precise, both sides in this relation are multiplied by δ(3) ( k + p + q). Eq. (2.10) indicates in an unambiguous way how the nonlinear parameter fNL and the kernel function K should be independently rescaled, once we demand that both the present day and the primordial kernel function are normalized to K λλ λ (kpv1, kpv2, kpv3) = K λλ λ ,pr (kpv1, kpv2, kpv3) = 1 (see text below eq. (2.7)).
The detector transfer function T encodes the response of the detector to a gravitational wave, and should not be confused with the cosmological transfer function introduced in eq. (2.8).
This turned out to be the case, for a broad range of frequencies around 1 mHz, in LISA Pathfinder, and is being considered as a working assumption in the LISA Data Challenge (Carlos Fernández Sopuerta, private communication).
Equation (5.7) holds at leading order in slow-roll for inflationary models with unit tensor sound speed. See[33] for a systematic discussion of more general scenarios that can change this relation.