Quantum coherence of relic gravitons and Hanbury Brown-Twiss interferometry

The coherence of the relic gravitons is investigated within a quantum mechanical perspective. After introducing the notion and the properties of the generalized Glauber correlators valid in the tensor case, the degrees of first- and second-order coherence are evaluated both inside and beyond the effective horizon. The inclusive approach (encompassing the polarizations of the gravitons) is contrasted with the exclusive approximation where the total intensity is calculated either from a single polarization or even from a single mode of the field. While the relic gravitons reentering the effective horizon after the end of a quasi-de Sitter stage of expansion are first-order coherent, the Hanbury Brown-Twiss correlations always exhibit a super-Poissonian statistics with different quantitative features that depend on the properties of their initial states and on the average over the tensor polarizations.


Introduction
As gravitational wave astronomy is opening a new observational window, the potential implications of the current developments for the stochastic backgrounds of relic gravitons are more accurately investigated. In a wide range of scenarios the early evolution of the spacetime curvature induces a stochastic background of primordial gravitational waves with a spectral energy density extending today from frequencies O(aHz) (i.e. 1 aHz = 10 −18 Hz) up to frequencies O(GHz) (i.e. 1 GHz = 10 9 Hz). While the specific features of different models will necessarily produce a variety of spectral amplitudes, all the current and planned experiments aiming at a direct (or indirect) detection of the relic gravitons are (or will be) sensitive to the average intensity of the gravitational field. In the language of the quantum theory of optical coherence the mean intensity of the field (or the average multiplicity of gravitons) is related to the degree of first-order coherence. With the goal of inspiring some of the future endeavours it is interesting to analyze the degrees of quantum coherence of the relic gravitons in a systematic perspective similar to the one already attempted in the case of large-scale curvature inhomogeneities.
The quantum theory of optical coherence [1,2,3] is customarily formulated in the context of vector fields but it can be generalized to the tensor and scalar cases by appropriately including (or excluding) the relevant polarizations. This is, after all, the logic already followed in quantum optical analyses where the scalar analog of the electromagnetic field is often scrutinized by focussing the attention on a single polarization (see e.g. [4]). The quantum treatment of the problem (not really mandatory prior to the celebrated series of experiments of Hanbury Brown and Twiss [5,6]) stems from the inadequacy of the classical description of the degree of second-order coherence of certain optical fields. According to the same perspective we can argue that the Young two-slit experiment (i.e. first-order correlations) is not a valid criterion to infer the quantum or classical nature of a given radiation field be it a vector field (as in the case of the photons) or a tensor field (as in the case of the gravitons). The interferometry originally developed by Hanbury Brown and Twiss can be then applied to relic phonons (i.e. quanta corresponding to large-scale curvature perturbations) and relic gravitons as firstly suggested some time ago [7,8]. While various approximations have been attempted, we intend to generalize the Glauber theory to the case of tensor fields.
The relic gravitons are potentially produced in the early Universe thanks to the pumping action of the gravitational field, as suggested in Refs. [9,10,11,12] even prior to the formulation of conventional inflationary models. The quantum theory of parametric amplification, originally developed for photons [13,14], has been generalized to the case of fields of different spins and, in particular, to the case of relic gravitons (see e.g. [7,8] and references therein). The gravitons produced from the vacuum or by stimulated emission (i.e. from a specific initial state) typically have opposite comoving momenta and lead to squeezed states [15]. Similar patterns arise in different scenarios including the conventional inflationary models where the spectral energy density has the usual quasi-flat slope [16] with a low-frequency break around 100 aHz [17]. Even if a direct detection of relic gravitons is not behind the corner both for technical and conceptual reasons, the degree of quantum coherence of the large-scale correlations could be used to disambiguate their origin, at least in principle [7,18,19]. In the light of these ambitious targets the present investigations are therefore mandatory. Indeed, the astrophysical events observed by the Ligo/Virgo collaboration (e.g. the three-detector observations of gravitational waves from black-hole coalescence [20], the evidence of gravitational waves from neutron star inspiral [21], and the observation of a 50-solar-mass binary black hole coalescence for a redshift z = 0.2 [22]) are qualitatively different from the potential signals coming from relic gravitons. Even if the spectra of relic gravitons are rightfully advertised as potential snapshot of the early Universe, their spectral energy density between few mHz (i.e. 1 mHz = 10 −3 Hz) and the kHz is rather minute 2 . The signal may sharply augment when the spectral energy density increases for frequencies larger than the mHz as it happens when the tensor modes of the geometry inherit a refractive index [24,25] or in the presence of stiff phases. In these cases it can happen that h c = O(10 −25 ) [24], while the chirp amplitudeh c corresponding to the astronomical signals detected so far by the Ligo/Virgo collaboration is O(10 −21 ) [20,21,22]. Even if the current upper limits on stochastic backgrounds of relic gravitons are still far from the final targets [27,28], the terrestrial interferometers in their advanced version [29,30] will hopefully probe chirp amplitudes O(10 −25 ) corresponding to spectral amplitudes h 2 0 Ω gw = O(10 −11 ). In the foreseeable future the Japanese Kagra (Kamioka Gravitational Wave Detector) [36,37] (effectively a prosecution of the Tama-300 experiment [35]) and the Einstein telescope [38] should be both operational in the audio band. The GEO-600 detector [39,40] is already progressing towards a further reduction of the quantum noise that will probably be essential for the third generation of terrestrial wide-band interferometers.
The layout of this investigation is the following. In section 2 the Glauber correlators are introduced in the tensor case. In section 3 the degree of first-order coherence is evaluated both inside and beyond the effective horizon. Section 4 is devoted to the estimate of the degree of second-order coherence while the role of the initial states will be analyzed in section 5. Section 6 contains the concluding discussion. To make the analysis self-contained the most relevant technical results have been relegated to the appendices A and B.

The quantum coherence of relic gravitons
The canonical theory of optical coherence is formulated in terms of vectors [1,2,3,4] but it is not uncommon to consider the scalar analog of the electromagnetic field by discussing a single polarization. The aim of this section is to extend the Glauber approach to the case of the divergenceless and traceless tensor fields describing the evolution of the relic gravitons.

Glauber correlation functions
The transverse and traceless fluctuations of the metric are conventionally denoted by The background geometry will be taken to be a conformally flat g µν = a 2 (τ )η µν where η µν is the Minkowski metric with signature (+, −, −, −) and a(τ ) is the scale factor in the conformal time parametrization. For practical purposes the correlation functions will be defined via the rescaled tensor amplitudê µ ij (x) =ĥ ij (x)a(τ ) and the hats will denote throughout the quantum field operators as opposed to their classical analog. The operatorsμ ij (x) consist of a positive and of a negative frequency part, i.e.μ ij (x) =μ If |vac is the state that minimizes the tensor Hamiltonian when all the modes are inside the effective horizon (for instance at the onset of inflation) the operatorμ ij (x) = 0). In the tensor case the Glauber correlation function is given by: whereρ is the density operator representing the (generally mixed) state of the fieldμ ij . Equation (2.1) generalizes the Glauber correlator (normally written in the case of photons) to the case of gravitons 3 . In quantum optics an exclusive perspective is often invoked by purposely neglecting one of the two polarizations of the photon [4]. This choice, motivated by specific empirical requirements 4 , amounts to expunging from the Glauber correlators the 3 Instead of (n + m) vector indices [2], in Eq. (2.1) we have (n + m) pairs of tensor indices (i.e. (i 1 j 1 ) ... (i n j n ) ... (i n+m j n+m )). 4 In quantum optics the single-polarization approximation is motivated by various experiments dealing with a single polarization (for instance in a cavity); this approach is exclusive since the experiments are typically conceived by only considering a single polarization. Since in the case of relic gravitons the initial conditions are not observationally accessible, it is useful consider also an inclusive approach where the sum over the polarizations is not neglected in the definition of the intensity and the quantum state of relic gravitons is unpolarized.
vector indices. The same logic and notations will be used in the case of the gravitons so that we define the scalar analog of Eq. (2.1): with the proviso that now Eq. (2.2) holds in the case of a single tensor polarization. It is finally not uncommon to treat the Mach-Zehnder and Hanbury Brown-Twiss interferometry in terms of a single mode of the field [41]. The single-mode experiments use plane parallel light beams whose transverse intensity profiles are not important for the measured quantities. In these situations is often sufficient to consider the light beams as exciting a single mode of the field. This viewpoint is even more exclusive than the one described by Eq. (2.2) where all the modes of the field are taken into account. In the explicit estimates of the following sections we shall consider, in this order, the inclusive description of Eq. (2.1), the exclusive approach of the single-polarization approximation (i.e. Eq. (2.2)) and finally the single-mode approximation.

Physical interpretations of the Glauber correlators
Equations (2.1) and (2.2) arise when considering the n-fold delayed coincidence measurement of the tensor field at the space-time points (x 1 , . . . x n , x n ). Let us focus, in particular, on the operator appearing inside the trace of Eq. (2.1), namelŷ and let us define |{ a} as the state of the field after the measurement and |{ b} the state of the field before the measurement. The matrix element corresponding to the absorption of gravitons at different times and at different locations of the hypothetical detectors can then be expressed as {a} |μ To obtain the rate at which the absorptions occur we must sum over the final states, i.e.  to Eq. (2.1) for x n+r = x r and (i n+r j n+r ) = (i r j r ) with r = 1, 2, . . ., n and n = m. The Glauber correlation function in the tensor case will therefore correspond to: The analog of Eq. (2.5) in the single-polarization approximation follows instead from Eq. (2.2) and it is given by: Having extended the Glauber correlators to the tensor case, it is now useful to introduce the corresponding degrees of quantum coherence.

Degrees of coherence for tensor fields
The first-order Glauber correlation function follows from Eq. (2.5) for n = 1 and it is: From Eq. (2.6) we can similarly obtain the analog of Eq. (2.7) in the single-polarization approximation: When (i 1 j 1 ) = (i 2 j 2 ) Eq. (2.7) describes the intensity averaged over the tensor polarizations and shall be denoted by Equation (2.5) in the case n = 2 determines the second-order correlation function which is relevant when discussing the Hanbury Brown-Twiss interferometry in the tensor case: Similarly from Eq. (2.6) we can obtain the analog of Eq. (2.10) in the single-polarization approximation: Equations (2.10) and (2.11) can describe the correlations of the intensities at two separate space-time points. This choice corresponds to the interferometric strategy pioneered by Hanbury Brown and Twiss (HBT) [5,6] as opposed to the standard Young-type experiments where only amplitudes (rather than intensities) are allowed to interfere. The applications of the HBT ideas range from stellar astronomy [5,6] to subatomic physics [42]. The interference of the intensities has been used to determine the hadron fireball dimensions [43,44] corresponding to the linear size of the interaction region in proton-proton collisions. To disambiguate the possible origin of large-scale curvature perturbations [7,8] and of relic gravitons, probably the only hope is the analysis of the degree of second-order coherence, as we shall argue. Since the intensity must be Hermitian [41] the standard HBT correlators follow from Eqs. (2.10) and (2.11) by requiring Thus with the identifications (2.12) the intensity correlators will be given by: (2.14) Note that Eq. (2.13) does depend on the polarizations through the sum over the tensor indices while the sum does not appear in Eq. (2.14). From the results of Eqs. (2.7)-(2.9) and of Eqs. (2.12)-(2.14) the corresponding degrees of quantum coherence can be easily obtained. More specifically the degrees of first-order coherence are: .

(2.15)
When the degree of first-order coherence has an overline it means that it is evaluated in the single-polarization approximation. The single-mode approximation for the degrees of coherence will be distinguished by a subscript (i.e. g (1) s ). For the sake of conciseness in Eq. (2.15) the following notations have been used: 16) and the same notations spelled out in Eq. (2.16) will be used throughout. Finally the degrees of second-order coherence for the relic gravitons follow from Eqs. (2.13) and (2.14) and they are given by: . (2.17) As in the case of the degree of first-order coherence the overline refers to the case of a single polarization; as usual g (2) s will denote the degree of second-order coherence in the single-mode approximation.

First-order coherence of relic gravitons
The field operators describing the positive and negative frequency parts can be expressed as:μ where P = √ 8πG and e where λ = iH = ia /a and, as already mentioned, the prime will denote throughout the discussion a derivation with respect to the conformal time coordinate. Equation (3.3) includes the sum over the two polarizations of the gravitons and is the continuous-mode generalization of the quantum mechanical Hamiltonian of Mollow and Glauber [13] (see also [46]). In the analysis of the evolution of the scalar and tensor modes of the geometry the quantum optical analogy has been firstly pointed out in Ref. [15,47].

General expressions for the degree of first-order coherence
The evolution ofâ k, α (τ ) andâ † − k, α (τ ) is determined by the Hamiltonian of Eq. (3.3) via the corresponding equations in the Heisenberg representation (see Eq. (A.9)). Their solution is then expressed in terms of the values of the creation and annihilation operators at the reference time τ i :â where the subscript p denotes the comoving three-momentum while the subscript α refers to the polarization. All the wavelengths that are today of the order of (or smaller than) the Hubble radius were presumably inside the effective horizon at τ i (i.e. kτ i 1) as it happens, for instance, in the case of conventional inflationary models. The two complex functions u p, α (τ, τ i ) and v p, α (τ, τ i ) appearing in Eqs. (3.4) and (3.5) are then solely determined by the specific dynamical evolution of the pump field λ appearing in Eq. (3.3) but they are also subjected to the condition |v p, α (τ, τ i )| 2 − |u p, α (τ, τ i )| 2 = 1 since the commutation relations between the two sets of creation and annihilation operators must be preserved. For each of the two tensor polarizations, u p, α (τ, τ i ) and v p, α (τ, τ i ) depend upon one amplitude and two phases and in spite of their specific unitary evolution they can always be parametrized as: where δ k, α , θ k, α and r k α are all real and depend on τ and τ i . The canonical transformation of Eqs. (3.4) and (3.5) is generated by the squeezing operator S(z) and by the rotation operator R(δ): where σ(z) and n(δ) involve both an integral over the modes and a sum over the two tensor polarizations 6 : Note that δ k,λ is real while z k, λ = r k, λ e iθ k, λ ; the parametrization of Eq. (3.6) follows from Eq. (3.4) and (3.5) by appreciating that 7 : Recalling the specific form of Eqs.  6 These two operators are typically expressed for a discrete set of modes but, in the present context, their continuous-mode generalization must be considered as it happens in the derivation of the ground state wavefunction of an interacting Bose gas at zero temperature [48,49]; the same approach has been used to describe the superfluid ground state [50,51] 7 Note, incidentally, that δ k, α denotes the phase δ with modulus of three-momentum k and polarization α; this quantity has noting to do with the Kroeneker δ αβ where α and β are instead two generic tensor polarizations. With this specification no confusion is possible.
the explicit form of the first-order Glauber correlation function given in Eq. (2.7): If the operatorb k, α (τ i ) annihilates the initial state at τ i the expectation value appearing in the last line of Eq. (3.12) corresponds to δ (3) ( k + p) δ αβ . However, to account for the possible presence of a finite number of gravitons at τ i , the expectation value shall be modified as denotes the average multiplicity of the initial state. With these specifications, Eq. (3.12) becomes: In the standard situation the relic graviton background is not polarized so that the u k,α and v k,α are the same for each of the two polarizations: Inserting the condition (3.14) into Eq. (3.13) we have that the first-order correlator becomes: The degree of first-order coherence is determined by Eq. (3.16) with k = i and = j: (3.17) where r = | x 1 − x 2 | and j 0 (kr) = sin kr/(kr) is the zeroth-order spherical Bessel function [52,53]. In the single-polarization approximation the analog of Eq. (3.17) reads In Eq. (3.18) the factor 4 comes from the sum over the polarizations that is counted in T (1) ( r; τ 1 , τ 2 ) but not in S (1) ( r; τ 1 , τ 2 ). As a consequence the degree of first-order coherence is The degree of first-order coherence computed in the single-polarization approximation, i.e. g (1) (r; τ 1 , τ 2 ) coincides with g (1) (r; τ 1 , τ 2 ) because of Eq. (3.18). When τ 1 → τ 2 and r → 0 we also have that Thus, in the zero-delay limit and for spatially coincident points the relic gravitons are always first-order coherent. The result implied by Eq. (3.20) is actually more general and it holds in all relevant physical regimes, as we shall discuss in the following three subsections.  22) where N = π/2e iπ(ν+1/2)/2 and ν = (β + 1/2). In the exact de Sitter case case β → 1 and ν → 3/2 and the results of Eqs. (3.21) and (3.22) then imply

Degree of first-order coherence beyond the effective horizon
When the given wavelength is larger than the effective horizon (i.e. kτ The squeezing parameters and the squeezing phases can be determined directly from Eqs. (A.15) and (A.16) but the same result also follows from the explicit expressions of u k (τ ) and v k (τ ) with the help of Eq. (3.6). For instance in the case of an exact de Sitter background of Eq. (3.23) the squeezing parameter and the two corresponding phases are: where y = H/(2k) = 1/(−2kτ ). The conformal time coordinate τ is negative during an exact de Sitter phase so that y → −∞ for typical wavelengths larger than the effective horizon. In the same limit it also follows from Eqs. (3.25) and (3.26) that the combination (δ k + θ k /2) is practically vanishing to an excellent approximation: While the result of Eq. (3.27) holds in the exact de Sitter case, in the quasi-de Sitter case we In the quasi-de Sitter case the results of Eq. (3.22) can then be expanded by recalling that for the range of parameters characteristic of the slow-roll dynamics (i.e. < 1) we have where the plus and the minus apply to the Hankel functions of second and first kind respectively [52,53]. Since to leading order in |kτ | there is is a cancellation in v k (τ ) and u k (τ ) the correct asymptotic result follows from Eq. (3.28) by keeping the next-to-leading correction in |kτ |: which coincides with Eq. (3.23) for ν = 3/2 and in the limit |kτ | 1. According to Eq. (3.29) the asymptotic values of the phases appearing in Eqs. (3.25) and (3.26) are given by δ k = π(1/2 − ν)/2 + O(kτ ) and by (θ k + δ k ) = π(5/2 − ν)/2 + O(kτ ). Conversely in the quasi-de Sitter case the combination ( From the results obtained so far it then follows that the degrees of first-order coherence of Eq. (2.15) can be explicitly computed from Eq. (3.29) when the relevant wavelengths are larger than the effective horizon: Equation (3.31) shows that the degrees of first-order coherence defined in Eq. (2.15) coincide and the single-polarization approximation gives the same result of the exact Glauber correlator. Since the dependence on τ 1 and τ 2 disappears from Eq. (3.31) the degree of first-order coherence goes always to 1 beyond the effective horizon. Equation (3.31) explains and justifies the analog result already mentioned in the zero-delay limit (see Eq. (3.20)).

Degree of first-order coherence inside the effective horizon
When the expansion rate exceeds the wavenumber a mode is said to be beyond the effective horizon: this does not necessarily have anything to do with causality [45]. The qualitative description of the evolution of the tensor modes of the geometry stipulates that a given wavelength exits the effective horizon (also dubbed sometimes Hubble radius) at a typical conformal time τ ex (for instance during an inflationary stage of expansion) and approximately reenters at τ re , when the Universe still expands but at a decelerated rate. Inside the effective horizon, i.e. in the limit |kτ 1, Eqs. (3.24)-(3.26) the squeezing parameters become By taking the large argument limits of the corresponding Hankel functions in Eqs. (3.21) and (3.22), u k (τ ) and v k (τ ) can be determined when |kτ | 1: implying δ k kτ and (θ k + δ k ) kτ . These results also imply that δ k + θ k /2 kτ 1 and apply when the modes of the fields are inside the effective inflationary horizon. Equation (3.33) is however not applicable during the radiation of matter phases but only describe the modes inside the effective horizon during the inflationary phase.
To compute the degrees of coherence inside the effective horizon after the end of inflation it is simpler to avoid specific exact solution such as the ones discussed in Eqs. (3.21)-(3.22) and directly work within an appropriate WKB approximation where the evolution of u k and v k will be approximate but more generally applicable. The strategy will be to ensure the correct normalization of u k (τ ) and v k (τ ) in the limit kτ 1 and then compute their form when the modes reenter either during the radiation-dominated phase or during the matter epoch. This analysis has been relegated to appendix B and it can be found in the current literature. As a consequence the functions u k (τ ) and v k (τ ) can be expressed as: where c ± (k) have been determined in Eq. (B.4) from the exact matching across the turning points τ ex and τ re (see also [15,24]). As a function of kτ the explicit expressions of Eqs.
(3.34)-(3.35) hold in the limit kτ 1. The values of c ± (k) depend on the reentry and on the exit of the given mode and they can be usefully approximated from Eq. (B.4) as 10 where ν ± (k) = k(τ ex ∓ τ re ) and the relation |c + (k)| 2 − |c − (k)| 2 = 1 holds explicitly. Equation (3.36) follows from the results of appendix B and it implies that, for wavelengths shorter than the effective horizon, the degree of first-order coherence becomes The result of Eq. (3.37) goes to 1 in the zero-delay limit; conversely if τ 1 = τ 2 the degree of first-order coherence computed from Eq. (3.37) is always smaller than 1, i.e. |g (1) (0, τ 1 , τ 2 )| ≤ 1 where the sign of equality holds for τ 1 = τ 2 .

Single-mode approximation
If the integrals over the modes and the sum over the polarizations are replaced by a single quantum oscillator the two terms appearing in the exponent of the operator S(z) become: This is, in a nutshell, the idea of the single-mode approximation which is not fully realistic since, in Eq. (3.38) the squeezing parameter and the phases must anyway follow the dynamics coming from the full Hamiltonian. With these caveats the single-mode approximation is crude but interesting, at least for comparison. We shall denote with the calligraphic style the continuous-mode operators while the single-mode operators will be denoted by the standard capital letter in roman style; so for instance where z = r e iθ , α = |α|e iϕ and so on and so forth. Using the properties of Eq. (3.39) we will have, for instance, The degree of first-order coherence in the single-mode case will then be given by where the subscript s reminds that we are here discussing the single-mode approximation.
In spite of the statistical properties of the state g (1) s → 1 in the zero-delay limit 11 (i.e. for τ 1 → τ 2 ). The single-mode limit can also be implemented in a slightly different way by introducing two different oscillators [ĉ,ĉ † ] = 1 and [d,d † ] = 1 (with [ĉ,d] = 0). In this case the two-mode squeezing and rotation operators [54,55] (analogs of R(δ) and S(z) but valid in the two-mode case) imply  The Hanbury Brown-Twiss correlations preliminary presented in Eq. (2.13) will now be estimated. For the sake of convenience the discussion mirrors exactly the same steps of the previous section. After determining the correlation functions, the degrees of second-order coherence will be explicitly discussed in various limits. From the expressions ofμ (−) (x) and µ (+) (x) Eq. (2.13) becomes: The expectation value appearing in the last line of Eq. (4.1) must be first referred to the initial time τ i when, by definition, all the relevant modes are inside the effective horizon (i.e. kτ i 1). For this purpose, using Eqs. (3.4) and (3.5) we can write Since the relic graviton background is not polarized, as in Eq. (3.13) the standard quantummechanical initial conditions imply u k,⊕ (τ, If the initial state is not the vacuum 12 at τ i the expectation value can be expressed as is the average multiplicity of the initial state and the standard case vacuum result will be recovered from the final HBT correlations in the limit n k (τ i ) → 0. After inserting into Eq. (4.1) the results of Eq. (4.2), the explicit form of the HBT correlations is: Equation (4.3) differs from the single-polarization approximation described by Eqs. (2.11) and (2.14). In the single-polarization approximation Eq. (4.3) must then be replaced by the result that follows, in the same physical situation, from Eqs. (2.11) and (2.14): (2.17) only depend on r = | x 1 − x 2 |: ,

Degree of second-order coherence beyond the effective horizon
When the wavelengths of the gravitons exceed the effective horizon, the functions u k (τ, τ i ) and v k (τ, τ i ) are given by the result of Eqs. (3.29). Using Eq. (4.3), to leading order in kτ 1 1 and pτ 2 1 the HBT correlations become: where the function E(k, p, r) appearing in Eq. (4.7) can be expressed as: where this time the analog of Eq. (4.8) is given, with the same notations, by: The degrees of second-order coherence defined in Eqs. (4.5) and (4.6) follow then from Eqs. (4.8) and (4.9) by recalling the explicit form of T (1) (τ ) and S (1) (τ ) given in Eqs. (3.17) and (3.18). Bearing in mind the results for the first-order correlations, it turns out that the intensity correlations are factorized as follows: (4.12) when the wavelengths exceed the effective horizon. The approximate equalities remind that Eqs. (4.11)-(4.12) hold in the limits |kτ 1 | 1 and |pτ 2 | 1. Thanks to Eqs. (4.5) and (4.6) the results of Eqs. (4.11) and (4.12) imply that 13 According to Eq. (4.13) the degree of second-order coherence is always super-Poissonian when the relevant wavelengths exceed the effective horizon since both g (2) and g (2) are larger than 1. In the single-polarization approximation the degree of second-order coherence goes to 3 while the presence of the polarization reduces the degree of second-order coherence.

Degree of second-order coherence inside the effective horizon
The initial conditions of the Einstein-Boltzmann hierarchy (both for the scalar and for the tensor modes of the geometry) are set before matter-radiation equality when the relevant wavelengths are larger than the Hubble radius at the corresponding epoch [7,8]. For an experiment probing the degree of second-order coherence in the CMB, the results of Eq. (4.14) The last to terms of Eq. (4.14) can be rewritten by factoring |c − (k)| 2 |c − (p)| 2 and by using Eqs. (3.36) in the obtained expression; the result of this step is given by: (4.15) Equation (4.15) can be explicitly estimated in two complementary limits. In the first case kτ re = pτ re 1; this choice corresponds to the situation where, in the vicinity of the turning points, |a /a| = 0. If the modes reenter when the condition |a /a| → 0 then kτ re < 1 and pτ re < 1 (see also appendix B and discussion therein). In both situations the results are similar and Eq. (4.15) can be approximated as: The result of the angular integration appearing in Eq. (4.16) is a complicated function of kr and pr that goes to a constant for kr < 1 and pr < 1. Therefore the degrees of second-order coherence will receive the dominant contribution for kr ∼ pr ∼ O(1): where ∆τ = τ 1 − τ 2 . Equations (4.17) and (4.18) coincide with Eq. (4.13) in the zero-delay limit. When τ 1 = τ 2 it can be demonstrated that |g (2) (τ )| < g (2) (0) and |g (2) (τ )| < g (2) (0) which implies, in a quantum optical language, that the degree of second-order coherence is not only super-Poissonian but also bunched [1,41].

Single-mode approximation and its physical interpretation
The degree of first-order coherence has been analyzed in the single-mode approximation at the end of section 3 and in full analogy with the definition of Eq. (3.41) the degree of second-order (temporal) coherence will be defined as: In the zero-delay limit Eq. (4.19) becomes: where, for simplicity, we will employ the notation g (2) s = g (2) s (0). After using the commutation relations Eq. (4.20) can be expressed in terms ofN =â †â and of the dispersion σ 2 : Equation (4.21) is often presented in terms of the so-called Mandel Q parameter [1] defined, within our notations, as : In the case of a single-mode squeezed state we have that the previous quantities can be all expressed in terms of a single parameter which is the average multiplicity of the state denoted hereunder by N = n sq = sinh 2 r: In the case of a coherent state the Mandel parameter vanishes so that the result of Eq. (4.22) can be dubbed by saying that the degree of second-order coherence is super-Poissonian. The results obtained in the present section suggest the following hierarchy: The first equality follows from the comparison between the single-mode approximation of Eqs. (4.21)-(4.22) and single-polarization approximation discussed in Eqs. (4.12), (4.13) and (4.18). Equation (4.23) suggests that the interference of the intensities of a singlepolarization can be approximated with the interference of the intensity of a single mode of the field. The effect of the polarizations is a progressive reduction of the degree of secondorder coherence. This reduction preserves the super-Poissonian character of the quantum state so that the Poissonian limit (typical of the coherent state) is never reached.

Stimulated versus spontaneous emission
The stimulated emission of relic gravitons does not reduce the degrees of coherence below the Poissonian limit provided the average multiplicity of the initial state does not dominate against the average multiplicity of the produced gravitons. This conclusion partly follows from Eqs. (3.37), (4.17) and (4.18) where the average multiplicity of gravitons at τ i has been already considered. When the average multiplicity of the initial state does not vanish the previous results describe the interference of the intensities of the relic gravitons produced by stimulated emission. Conversely in the limit n k (τ i ) → 0 the same expressions hold in the case of spontaneous emission whereb k (τ i ) annihilates the vacuum at τ i . While the previous analyses show that the super-Poissonian behaviour is not altered by more general parametrizations of the initial state, it is also true that the parametrization of the initial state suggested above is not the most general one. This complementary aspect of the present analysis will now be clarified by considering the initial conditions provided by a coherent state in the continuous mode representation. Since a given multiparticle density matrix can be projected on the coherent state basis via the so-called Klauder-Sudarshan P -representation [4] this analysis seems appropriate and sufficiently conclusive, at least for the present purposes.

Squeezed coherent states
For a continuum of modes the Glauber displacement operator is defined as [56] with the same notations already employed for the squeezing and rotation operators S(z) and R(δ) of Eqs. (3.7) and (3.8). The squeezed coherent states of relic gravitons can be introduced in two complementary perspectives mirroring their quantum optical analogs originally discussed by Caves [57] (sometimes also referred to as ideal squeezed state) an by Yuen [58] (the so-called two-photon coherent states). In the Caves representation the initial state is rotated, squeezed and finally dispalced 14 i.e. |{α z δ} = D(α)|{z δ} . In the Yuen representation [58] the squeezed-coherent states of relic gravitons are instead defined as |{z δ β} = S(z)R(δ)|{β} . According to the strategy of Ref. [58] (appropriately extended to the continuous-mode description of Eq. (5.1) the creation operators transform as: If the action of the displacement operator precedes the squeezing and the rotation, as suggested in Ref. [57], the creation operators transform instead as: +e −i δ q λ cosh r q λâ q λ − e i(θ q λ +δ q λ ) sinh r q λâ † − q λ . α q,λ = e −i δ q,λ cosh r q,λ β q,λ − e i(θ q,λ +δ q,λ ) sinh r q,λ β * − q,λ .

Degrees of coherence
The first-order and second-order Glauber correlators for a squeezed-coherent are The action of D(α) overμ (−) (x) andμ (+) (x) displaces the field operators by their classical value [56]: where we stress that µ c (x) is not an operator but a classical field. Inserting Eq. (5.7) into Eq. (5.5) the first-order Glauber correlator is: The first term in Eq. (5.8), analog to the condensate arising in the theory of superfluidity [50,51], depends on x 1 and x 2 . Conversely the second term W 0 (r, τ 1 , τ 2 ) is the quantum contribution which is a function of the distance. The explicit form of the HBT correlations is given by Besides the contribution of the condensate and of the quantum fluctuations (first line of Eq. (5.9)), we can identify a mixed contribution (second line of Eq. (5.9)) and two interference terms (third and fourth lines of Eq. (5.9)) that depend on the mutual phases characterizing the displacement, squeezing and rotation operators. The degree of second-order coherence can then be expressed in the following manner , (5.10) where, from Eq. (5.8), we defined W 0 (τ ) = W 0 (0, τ, τ ) and also According to Eq. (5.10) and thanks to the explicit form of Eqs. (5.11)-(5.13) the sign of g (2) (x 1 , x 2 ) − 1 may become negative and this possibility is well known from quantum optical studies where a squeezed state with a strong coherent component may have a sub-Poissonian statistics [1,41,57,58] provided the average multiplicity of the coherent component dominates against the the squeezing contribution, as already anticipated at the beginning of this section. A similar conclusion will be reached hereunder and to investigate the analog phenomenon in the present case we consider the regime x 1 → x 2 where Eq. (5.10) becomes: (5.14) The explicit expressions of W 0 (τ ), W 1 (τ ), W 2 (τ ) and W 3 (τ ) appearing in Eq. (5.14) is: While the term appearing in the first line at the right hand side of Eq. (5.14) is always positive semidefinite, the remaining two terms (in the second line of the same equation) do not have a definite sign. We can therefore conclude that g (2) (x) > 1 as long as the average multiplicity of the produced gravitons exceeds the coherent component of the initial state.
In particular when µ c (x) → 0 the case treated in the previous section is recovered and g (2) (x) → 3.

Wavelengths inside and beyond the effective horizon
If the spectrum of the initial fluctuations is characterized by a given wavelength (for instance a thermal wavelength) at τ i , the present value of this characteristic scale will be much larger than the Hubble radius at the present time unless the total number of efolds N t is very close 15 to the critical number of efolds N crit = O(66) [59,60]. Even assuming (or tuning) N t ∼ N crit it seems difficult to conceive an initial state that could make the statistics sub-Poissonian. For this purpose we can investigate more carefully the sign of g (2) (x) − 1 when the coherent component dominates over the average multiplicity of the produced gravitons. It is useful to introduce the quantities 0 (x) = W 0 (τ )/|µ c (x)| 2 < 1 and 3 (x) = W 3 (τ )/|µ c (x)| 2 < 1 that are both smaller than one when the coherent component exceeds the squeezed contribution; Eq. (5.14) can then be written as: where µ c (x) = e iϕ(x) |µ c (x)| has been separated in its modulus and phase. The first term at the right hand side of Eq. (5.18) contains contributions O( 2 0 ) and O( 2 3 ) that are subleading in comparison with the remaining two terms (of order 0 and 3 respectively). The sign of [g (2) (x) − 1] will then be determined by the balance of the dominant contributions at the right hand side. If we now recall Eq. (3.6) and assume, for simplicity, that ϕ is constant we can rephrase the dominant contributions of Eq. (5.18) as: where ζ k = (ϕ − θ k /2). This result implies that g (2) (x) − 1 < 0 when cos 2ζ k > |v k |/|u k | and provided the coherent component exceeds the squeezing contribution. In this limit, however, the large-scale fluctuations will simply correspond to the coherent contribution. For this possibility it will then be essential to tune the total number of efolds around their critical value. If we now recall Eq. (3.6) the condition (|v k (τ )| 2 − cos 2ζ k |u k (τ )| |v k (τ )|) < 0 becomes sinh r k [sinh r k − cosh r k ] cos 2 ζ k + sinh r k [sinh r k + cosh r k ] sin 2 ζ k < 0. (5.20) After simple algebra, Eq. (5.20) can also be expressed as: showing that if ζ k = 0 the inequality is always verified even in the limit r k 1. The results of the Caves approach [57] discussed so far can be translated into the Yuen description [58] by making use of Eq. (5.4) which can also be written, in the single-polarization approximation, . From Eq. (5.4) the relation between ζ k and γ k is given by: which also demands e 4r k (1 − e −2r k ) sin 2 γ k < e −2r k (1 − e −2r k ) cos 2 γ k . In analogy with Eq. (5.21), if γ k → 0 the previous inequality is always verified even in the limit r k 1. But unfortunately the limit γ k → 0 can only be realized in the exact de Sitter case (see Eq. (3.23)) or if we consistently tune χ k → (δ k + θ k /2) for all modes that exceed the Hubble radius. Conversely, without fine-tuning, Eq. (3.30) implies γ k −π /2 in the quasi-de Sitter case and Eq. (5.23) demands In the second term we can always neglect the 2 correction which is small with respect to 1; from the remaining terms we have 1 < r k < (1/6) ln [4/(π 2 2 ) − 1/2] i.e. r k < 2.15 for = 0.001. This condition is therefore unphysical since the average multiplicity of the gravitons produced during inflation would be negligibly small and this would imply that the tensor modes are not amplified in comparison with the scalar modes. The result of Eq. (5.24) is, in some sense, pleonastic since the relic gravitons potentially observable today (for instance in the audio band) are all inside the effective horizon (i.e. kτ 1) and, in this limit, γ k = 0 for independent reasons. Again the statistics can never become sub-Poissonian unless the average multiplicity of gravitons produced during inflation is negligible 16 . Inside the Hubble radius the phases are determined by the pair of conditions where z = a re /a ex , q = k/H re and ν ± (k) = k(τ ex ∓ τ re ). In Eqs. (5.25) and (5.26) the corresponding expressions have been simplified in the limit z 1 and k/H re < z 2 so that the squeezing phases are given by But this means δ k + 2θ k [ν + + ν − ]/2 = kτ ex = O(1) implying that, for squeezed coherent states, the statistics of the relic gravitons is never sub-Poissonian. Following the ideas conveyed in this section, different initial states such as squeezed-number states or even mixed states (e.g. squeezed thermal states) can be analyzed with similar qualitative results. In some cases the degree of second-order coherence can be reduced by the presence of an appropriate initial state. While we leave the explicit analysis of this interesting point to a more topical discussion, we can safely conclude that the properties of the initial states may very well interfere with the squeezing contribution but do not affect the super-Poissonian character of the degree of second-order coherence and of the HBT correlations especially inside the effective horizon. 16 The present results are at odd with the claim of Ref. [19] (obtained in the single-mode approximation) where the authors suggest that the analog of γ k is generically vanishing. This is only true provide the phase of the coherent state is tuned in a way that the squeezed contribution is exactly cancelled.

Concluding discussion
To assess the classical or quantum origin of the relic gravitons the only sound strategy is a careful scrutiny of the higher degrees of quantum coherence. The first step along this direction is a proper extension of the Glauber theory of quantum coherence to the case of the tensor modes of the geometry and this has been the aim of the present investigation. The degree of first-order coherence of relic gravitons always tends to 1 when the corresponding wavelengths are either larger or smaller than the effective horizon. In standard Young interferometry the degree of first-order coherence goes to 1 when the interference fringes are maximized while it goes to 0 in the opposite case when the field is incoherent. Classical configurations and quantum states of a given optical field lead to comparable degrees of first-order coherence and this conclusion remains practically unchanged in the case of relic gravitons. The analysis of the Hanbury Brown-Twiss correlations in their canonical form shows instead that the degree of second-order coherence is always larger than 1. In the quantum optical jargon the relic gravitons are therefore bunched and their statistics is super-Poissonian. The results are physically similar if more exclusive approaches are adopted when, for instance, a single tensor polarization or even a single mode contributes to the total intensity of the field. The obtained conclusions do not change if we consider stimulated (rather than spontaneous) emission of relic gravitons. While the super-Poissonian nature of the degree of second-order coherence is a necessary condition if we want to infer the quantum origin of the relic gravitons, such a requirement is not sufficient since other states (for instance mixed) may lead to a super-Poissonian degree of second-order coherence. Even if, according to some, the quantum origin of the relic gravitons and of large-scale curvature perturbations is indisputable, the spirit of the present analysis is more modest and it aims at formulating a set of criteria which could independently confirm (or infirm) the conventional viewpoint. The obtained results suggest the Hanbury Brown-Twiss interferometry and the scrutiny of the higher degrees of quantum coherence could be a valid (and probably unique) tool in these matters.

A Quantum theory of parametric amplification
The action describing the parametric amplification of the relic gravitons can be compactly expressed as [10] S (t) = 1 8 2 where γ ± p α is related to the presence of the anisotropic stress and my lead to a coherent complex functions u p α (τ, τ i ) and v p α (τ, τ i ) must therefore depend, for each polarization and for each mode, upon three real functions (i.e. two phases and one amplitude) and have been parametrized in terms of the a squeezing amplitude r p,α supplemented by two squeezing phases (i.e. θ p,α and δ p,α ), as already mentioned in Eqs. In the limit |kτ | 1 we have, approximately, that δ k + θ k /2 0 which is indeed the correct result. However, by summing up Eqs. (A.15) and (A.16) and by taking the limit of the obtained expression we would be led to conclude that the same result also holds for |kτ | 1 (i.e. inside the effective horizon). This however does not happen as it can be verified from the exact solution of Eqs. (3.24)-(3.26) implying that the wanted combination is given by: If we now take the limit kτ 1 we will have from Eq. (A.17) that δ k + θ k /2 4k 3 τ 3 /3 1 that coincides with Eq. (3.27). In the opposite limit (i.e. |kτ | 1) Eq. (A.17) does not imply δ k + θ k /2 1 but rather δ k + θ k /2 kτ 1. This result indeed agrees with the exact solution (3.24)-(3.26) in the limit |kτ | 1 (see also Eq. (3.32)).

B Evolution inside the effective horizon
The relic gravitons potentially observable today are all inside the effective horizon and to estimate their degrees of quantum coherence the full expressions of u k (τ ) and v k (τ ) must be evaluated for τ > τ re (see e.g. Eq. (3.32) and discussion thereafter). For this purpose the idea is to estimate the amplification of the mode functions and then relate the obtained result to the asymptotes of u k (τ ) and v k (τ ) in the limit kτ 1 and τ > τ re . As already mentioned in the bulk of the paper, Eqs. (A.12) and (A.13) can be decoupled. In particular the combination (u k − v * k )/ √ 2k = f k obeys the standard equations f k + [k 2 − a /a]f k = 0 whose solution is well known in the different asymptotic regimes (see e.g. [24]). A given wavelength exits the effective horizon (also dubbed sometimes Hubble radius [45]) at some typical conformal time τ ex during an inflationary stage of expansion and approximately reenters at τ re , when the Universe still expands but in a decelerated manner. The two typical times τ ex and τ re are the turning points of the WKB approximation [15] and are both determined from the condition k 2 = |a /a|. If |a /a| = 0 in the vicinity of the turning point, then kτ 1 and this is what normally happens at τ ex . However, when |a /a| → 0 in the region where the turning point is located, the situation is different. Since |a /a| → 0 during radiation, if the given mode reenters during the radiation epoch (or anyway in a regime where the pump field vanishes either approximately or exactly) then kτ re 1 (even if, for τ > τ re , |kτ | > 1). In the WKB approximation we therefore have Note that g k = f k − Hf k and g k (τ ) = a ex a g k (τ ex ) − k 2 a(τ ) τ τex a(τ 1 )f k (τ 1 ) dτ 1 .