Glauber theory and the quantum coherence of curvature inhomogeneities

The curvature inhomogeneities are systematically scrutinized in the framework of the Glauber approach. The amplified quantum fluctuations of the scalar and tensor modes of the geometry are shown to be first-order coherent while the interference of the corresponding intensities is larger than in the case of Bose-Einstein correlations. After showing that the degree of second-order coherence does not suffice to characterize unambiguously the curvature inhomogeneities, we argue that direct analyses of the degrees of third and fourth-order coherence are necessary to discriminate between different correlated states and to infer more reliably the statistical properties of the large-scale fluctuations. We speculate that the moments of the multiplicity distributions of the relic phonons might be observationally accessible thanks to new generations of instruments able to count the single photons of the Cosmic Microwave Background in the THz region.


Quantum fluctuations and their coherence
The Cosmic Microwave Background radiation (CMB in what follows) has been repeatedly scrutinized by a number of ground based observations and various spaceborne missions starting with the COBE satellite [1]. The first data released by the WMAP collaboration [2] used the newly discovered cross-correlations between the temperature and the polarization to analyze the initial conditions of the Einstein-Boltzmann hierarchy. The subsequent releases of the WMAP experiment and the Planck explorer results [3] confirmed (and partially refined) the early determinations of the first WMAP data [2] so that today we can say, with a fair degree of confidence, that the initial conditions of the CMB temperature and polarization anisotropies are predominantly adiabatic and Gaussian. Every deviation from this paradigm leads to a number of entropic solutions 2 (see e.g. [4]). When the CMB data are combined with other cosmological data sets (such as the large-scale structure data [6] and the supernova data [7]) the typical parameters describing the large-scale curvature modes can be slightly (but not crucially) modified.
The temperature and the polarization anisotropies of the CMB are believed to originate in the early Universe from the fluctuations of the spatial curvature. Sakharov [8] was presumably the first one to raise the question of the quantum mechanical origin of density perturbations in the early Universe suggesting that the complicated patterns observed in the galaxy counts could actually have some plausible origin in the zero-point fluctuations of certain quantum fields. The adiabatic and Gaussian nature of the initial conditions of the Einstein-Boltzmann hierarchy is compatible with a quantum mechanical origin of the curvature perturbations but the observed patterns of the temperature and polarization anisotropies are not sufficient to establish their quantum origin. The quest for accuracy in the determinations of the parameters of the concordance paradigm does not tell us anything about the quantum state of large-scale perturbations and of their multiplicity distribution 3 .
The quantum state of the scalar and tensor modes of the geometry is determined by the dynamical features of the model, by the its initial conditions and, at least in the standard lore, by the overall duration of the inflationary phase. Even if there are some who suggest that we have an accurate control of the protoinflationary dynamics it would be nice to develop a set of sufficient criteria enabling us to infer the quantum origin of large-scale curvature perturbations from some observational evidence. A first step along this direction relies on 2 The current data [2,3,5] seem to suggest that a small fraction of anticorrelated entropic modes in the presence of a dominant adiabatic mode may even improve the fit of the temperature autocorrelations accounting for potential large-scale suppressions of the corresponding angular power spectra. 3 In what follows we shall often refer to the statistical properties of quantum states and of stochastic processes with discrete state space. By this we mean their degree of coherence encoding the moments of the multiplicity distribution. See, in this respect, sections 2 and 3.
the idea of studying (and eventually measuring) the correlation functions of the intensities of the curvature perturbations [9] rather than the correlations of the corresponding amplitudes. This concept has been originally proposed by Hanbury Brown and Twiss [10] and the analysis of the intensity correlations is often dubbed Hanbury Brown-Twiss (HBT) interferometry as opposed to the standard Young-type interference where only amplitudes (rather than intensities) are concerned 4 . The Glauber theory has been used through the years to reach a fuller understanding of the relations between the classical and quantum theories of light [13]. The first experimental detection of non-classical light was made nearly fifty years ago [14]. Since then the secondorder coherence effects have been widely used to infer the statistical properties of light sources. For similar purposes the Glauber theory can be applied to the analysis of the largescale curvature fluctuations. In the conventional view the questions raised by this approach are dismissed since the quantum origin of large scale curvature fluctuations is, according to the theoretical prejudice, an indisputable fact of nature. In this paper we would like to take a pragmatic attitude which is incidentally also more modest: we ought to understand in what sense the HBT interferometry [9] is sufficient to infer the statistical properties of large-scale curvature perturbations. Even if the theoretical results addressed here are encouraging, a number of experimental questions remain. While these questions are beyond the scopes of this investigation and will be left unanswered now, it is not excluded they might find satisfactory answers in the future.
The layout of this investigation is therefore the following. In section 2 we shall summarize the essentials of the Glauber approach by stressing the physical analogies with the treatment of large-scale curvature perturbations. The first and second-order coherence effects will be analyzed in section 3. The ambiguities of the second-order coherence will be addressed in section 4. In section 5 we shall then argue that the degrees of third and fourth-order coherence are sufficient to infer the statistical properties of large-scale curvature inhomogeneities. In section 6 we shall discuss some future observational perspectives. Finally, section 7 contains our concluding remarks. To avoid lengthy digressions some of the formulas instrumental in obtaining the results presented in the bulk of the paper have been rederived, in a selfcontained perspective, in the appendices A, B and C.

Glauber theory of coherence
Consider the field operatorq( x, τ ) where x is the position and τ is the (conformal) time coordinate. In quantum opticsq is often identified with a single polarization of the electric field (as in the original Glauber approach [15]), with a single polarization of the vector potential (as argued by Mandel and Wolf 5 [13]) or even with a vector field itself, if suitable vector indices are included.
In the present context the field operatorq coincides either with the canonically normalized curvature perturbations of a conformally flat 6 Friedmann-Robertson-Walker space-time or with a single tensor polarization. In both cases the evolution of the quantum Hamiltonian can be written as:Ĥ whereq andπ are the two canonically conjugate and Hermitian field operators; the pump field λ depends on the conformal time and it satisfies λ * = −λ. According to the specific dictionary developed in the appendix A, Eq. (2.1) may describe the parametric amplification of either the scalar or of the tensor modes of the geometry. Equation (2.1) is equivalent to Eq. (A.6) which has been analyzed for the first time by Mollow and Glauber [16] in the framework of the quantum theory of the parametric amplification. When expressed in terms of the corresponding creation and annihilation operators the free part of the Hamiltonian and the two components of the interacting Hamiltonian satisfy the usual commutation relations of the SU(1, 1) Lie algebra. In the case of the tensor modes of the geometry the canonical Hamiltonian of Eq. (2.1) follows from the corresponding gauge-invariant action obtained long ago by Ford and Parker 7 [17]. The scalar fluctuations of the metric will be described here within the standard Bardeen formalism [19]; the corresponding gauge-invariant action has been discussed, for the first time, by Lukash [20] in a background model filled by a perfect and irrotational fluid. Various authors applied a similar analysis to the case of scalar field matter and to its fluctuations [21,22]. The Hamiltonians for the scalar and tensor modes have been discussed, in a unified perspective, in appendix A (see also [23]). Indeed, the same process leading to the curvature 5 The authors of Refs. [13] suggested indeed that the definition in terms of the vector potential is more convenient when dealing with the fotoelectric detection of light fluctuations. 6 Consistently with observational determinations [2,3], the background metric will be conformally flat, i.e. g µν = a 2 (τ )η µν where a(τ ) is the scale factor, τ denotes the conformal time coordinate and η µν is the flat metric with signature mostly minus, i.e. (+, −, −, −). 7 In the seminal papers of Grishchuk [18] the emphasis was on the breaking of Weyl invariance of the equations of motion rather than on the gauge-invariant action.
inhomogeneities also produces to a stochastic background of relic gravitons extending over three decades in frequency [24].

The Glauber correlation function
The parametric amplification described by the Hamiltonian (2.1) produces specific quantum correlations in the final state. Since the degrees of quantum coherence provide a systematic approach to the statistical properties of the final state [15], it is rather plausible to analyze the large-scale curvature perturbations in the light of the Glauber approach. To comply with this program, the first step is to introduce the Glauber correlation function in its most general form: andρ is the density operator representing the (generally mixed) state of the fieldq. The fieldq( x, τ ) can always be expressed asq( . By definition we will have thatq (+) (x)|vac = 0 and also that vac|q (−) (x) = 0; the state |vac denotes the vacuum 8 . Provided the total duration of inflation exceeds the minimal number of about 65 efolds [3,23,27], the vacuum initial data are the most plausible, at least in the conventional lore.
Let us now consider in greater detail the expression of Eq. (2.2) and let us observe that an operator of the typê is needed to describe n-fold delayed coincidence measurements of the field at the space-time points (x 1 , . . . x n ). If | b is the state before the measurement and | a is the state after the measurement, the matrix element corresponding to the absorption of the quanta ofq at each detector and at given times is a |q (+) (x 1 ) . . .q (+) (x n )| b . The rate at which such absorptions occur, summed over the final states, is therefore proportional to: In the discussion of section 3 the vacuum corresponds to the state minimizing the Hamiltonian at the onset of the dynamical evolution as discussed in appendix B. This state can be explicitly constructed by diagonalizing the Hamiltonian in terms of an appropriate canonical transformation. A similar procedure is used to derive the ground state wavefunction of an interacting Bose gas at zero temperature [25,26].
whereÔ has been given in Eq. (2.3) and the second equality of Eq. (2.4) follows from the completeness relation. It is clear from Eq. (2.4) that when b|Ô|b is averaged over the ensemble of the initial states of the system it becomes identical with Eq. (2.2) for x n+r = x r (with r = 1, 2, . . ., n and n = m). Since this is the case that will be studied hereunder we shall denote the Glauber correlation function as

The degrees of quantum coherence
Recalling the result of Eq. (2.5), the coherence properties of the quantum fieldq(x) can be discussed by introducing the normalized version of the n-point Glauber function [15]: . (2.6) While, by definition, |g (1) (x 1 , x 2 )| ≤ 1 the higher order correlators are not restricted in absolute value as it happens for g (1) (x 1 , x 2 ). A fully coherent field must therefore satisfy the following necessary condition: for n = 1, 2, 3, . . .. If only a limited number of normalized correlation functions will satisfy Eq. (2.7) we shall speak about partial coherence. So for instance if g (1) (x 1 , x 2 ) = 1 and g (2) (x 1 , x 2 , x 3 , x 4 ) = 1 (but g (3) (x 1 , x 2 , x 3 , x 4 , x 5 , x 6 ) = 1) we shall say that the radiation is second-order coherent. We shall be specifically interested in the first, second, third and fourth-order correlators. The degrees of first-and second-order coherence are: Similarly the third-and fourth-order degrees of coherence are: . (2.11) The degree of first-order coherence appears naturally in the Young two-slit experiment. Whenever the degree of first-order coherence is equal to 1 the visibility is maximized [13].
The degree of second-order coherence enters the discussion of the Hanbury Brown-Twiss effect [10] and its different applications ranging from stellar interferometry [13] to high-energy physics [11,12]. The degree of second-order coherence arises naturally when discussing the correlations of the intensities of the fieldq (see section 3). Notice that the intensity correlators relevant to the HBT interferometry can be easily obtained from Eq. (2.6) by identifying the space-time points as follows: (2.12) In this case the original Glauber correlator will effectively be a function of n points and and it will describe the correlation of n intensitiesÎ( Equations (2.10) and (2.11) describe, respectively, the degrees of third-order and fourthorder coherence. It has been recently suggested, in quantum optical applications, that the degree of second-order coherence might not always be sufficient to specify completely the statistical properties of the radiation field [28,29,30,31]. A specific discussion of these points can be found in section 4.
Since a fully coherent field must satisfy Eq. (2.7) the Poissonian case will be used as a benchmark value to analyze the higher-order correlators. From the definitions given above it can be argued that while the degree of first-order coherence cannot exceed 1 (see Eq. (2.7) and discussion thereafter), the degrees of higher-order coherence do not have an obvious upper bound. Various inequalities satisfied by degrees of coherence of various orders have been derived by Titulaer and Glauber [32] but the simpler results given above are adequate, at least for the present discussion.

Super-Poissonian degree of coherence
The degrees of first-order and second-order coherence given Eqs. (2.8) and (2.9) will now be specifically analyzed when the evolution of the field operators is governed by the Hamiltonian (2.1) and the corresponding pump fields are determined by the quasi-de Sitter evolution. The initial conditions imposed on the evolution of the field operators will imply that the Hamiltonian of Eq. (2.1) is both diagonal and minimized. Equation (2.5) in the cases n = 1 and n = 2 implies: The correlators G (1) (x 1 , x 2 ) and G (2) (x 1 , x 2 , x 3 , x 4 ) obey the following relations: The quantum theory of the parametric amplification will be discussed in the Heisenberg description 9 . When the pump field λ is determined by the quasi-de Sitter the degree of firstorder coherence goes always to 1; this happens in the limit of wavelengths larger than the Hubble radius and in spite of the correlation properties of the initial state. In the same limits the degree of second-order coherence goes to 3; the curvature phonons (and the gravitons) are therefore bunched [13,37], their statistics is super-Poissonian and their degree of correlation exceeds the typical value of a chaotic source.

Quasi-de Sitter evolution: first-order coherence
The evolution of the field operators depends on two complex functions v k (τ ) and u k (τ ) whose explicit expression can be obtained by solving Eqs. (C.1) and (C.2). In the case of a quasi-de Sitter evolution v k (τ ) and u k (τ ) are given by Eqs. (C.3) and (C.4); n the pure de Sitter case their expression is given instead by Eq. (C.5). The first-order correlation function of Eq.
(3.1) at separate space-time points: Consequently, from Eqs. (2.8) and (3.4) the normalized degree of first-order coherence is: We are interested in the value of the first-order coherence when the relevant wavelengths are larger than the Hubble radius and, in this limit, the explicit expressions of Eqs. (C.3) and where µ is given, in the scalar and tensor case, by Eq. (C.4). Equation (3.6) has been obtained by expanding the Hankel functions of Eq. (C.3) in the limit of small arguments [38]. The first subleading term appearing in Eq. (3.6) is suppressed by 3−2µ while the second correction is suppressed, in comparison with the leading term, by |kτ | 1/2 which is small in 9 The quantum theory of parametric amplification has been originally developed in [16] and subsequently analyzed in connection with the dynamics of the squeezed quantum states [33,34]. An explicit discussion of the evolution of cosmological perturbations in the Schrödinger representation has been presented long ago by Grishchuk and collaborators Ref. [35] and subsequently discussed by various authors (see for instance [9,36]); see also the discussion at the end of appendix C.
the large-scale limit 10 . In the limit µ → 3/2 we have, from Eq. (3.6), that v k (τ ) = −i/(2kτ ) which coincides with the exact result obtainable by setting µ = 3/2 in Eq. (C.3). Performing the angular integrals we have that Eq. (3.6) implies where j 0 (k 1 r) denotes the spherical Bessel function of zeroth order [38]. The numerator and the denominator of Eq. (3.7) depend on τ 1 and τ 2 . However, owing to the specific form of v k (τ ), the dependence on τ 1 and τ 2 simplifies and the final form of Eq. (3.7) becomes where the integrals are evaluated over all the modes larger than the Hubble radius and the second relation clearly holds in the limit k 1 r ≪ 1 (corresponding to large angular separations). The same result implied by Eq. (3.8) can be obtained when kτ ≫ 1 in Eq. (C.3). In this case the Hankel functions become plane waves and in the limit k 1 r ≪ 1 we have that g (1) ( r, τ 1 , τ 2 ) always tend to 1 implying the first-order coherence of the underlying fluctuations.

Quasi-de Sitter evolution: second-order coherence
The general form of the second-order correlation function involves four separated space-time points. As in the case of the first-order coherence Eq. (3.2) can be written in terms of v k (τ ) and u k (τ ): If we identify the space-time points two by two (i.e x 1 = x 3 and x 2 = x 4 ) the correlation function of Eq. (3.9) describes the interference of two beams with intensitiesÎ( x 1 , τ 1 ) and where r = x 1 − x 2 . Inserting Eqs. (3.9) and (3.10) into Eq. (2.9) the degree of second-order coherence becomes: According to Eq. (3.11) the degree of second-order coherence involves four field operators at two separated spatial points [13,37]. From the observational viewpoint the correlations of two intensities (i.e. two beams) translates into the correlation of the output of four distinct brightness perturbations (see also the discussion in section 6). Equation (3.11) can be further simplified by appreciating that the mean number of produced particles per Fourier mode is simply n k 1 (τ 1 ) = |v k 1 (τ 1 )| 2 and n k 2 (τ 2 ) = |v k 2 (τ 2 )| 2 . Furthermore in the pure de Sitter case Eq. (C.4) imply: The same relation holds also in the quasi-de Sitter case to leading order in the slow-roll expansion and for scales larger than the Hubble radius (i.e. kτ ≪ 1). Thanks to the preceding two observations we have that the degree of second-order coherence, to leading order 11 in 1/n k is given by: The normal-ordered expectation values dominate the degrees of coherence since the average multiplicity of the produced phonons and gravitons is large; this observation defines the 1/n k expansion (see, in this respect, the final part of appendix C). In the pure de Sitter case, from Eq. (C.5), n k (τ ) = (4k 2 τ 2 ) −1 ≫ 1 for kτ ≪ 1.
The large-scale limit the spherical Bessel functions go to 1 and therefore Eq. (3.13) becomes 12 : 14) The result of Eq. (3.14) holds provided the total number of efolds N tot is larger than the maximal number of efolds today accessible by observations 13 denoted by N max . The ambition of the present analysis is not to endorse a particular initial state of the curvature inhomogeneities but rather to stress the need of a systematic statistical scrutiny of the large-scale correlations. Given the larger arbitrariness of the initial data, the analysis of the degrees of higher-order coherence seems even more compelling when the initial state is not the vacuum. The methods employed in subatomic physics involve the construction of two-particle (and eventually three-particle) correlation functions from the distribution of particles radiated from a hot (and spatially localized) source [11]. While the concepts are very much the same, the notations may be slightly different so that the pivotal variable is often defined as R( r, τ ) = g (2) All in all the large-scale curvature fluctuations are only partially coherent. They are firstorder coherent but, according to the Glauber terminology, their statistics is super-Poissonian with a degree of second-order coherence which is thrice the one of a fully coherent field (see Eq. (2.7) and discussion thereafter).

Potential ambiguities
In conventional HBT interferometry the space-time dimensions of the emitters need to be determined but the statistical properties of the source are known. For large-scale curvature perturbations the reverse is true and the gross uniformity of the temperature fluctuations at last scattering implies that curvature perturbations, prior to matter-radiation equality, had typical wavelengths larger than the Hubble radius at the corresponding epoch. In this situation the degrees of first and second-order coherence are however not sufficient 12 The result of Eq. (3.14) holds also in the case τ 1 → τ 2 ; in the zero time-delay limit τ 1 − τ 2 = 0 and in this case the degree of second-order coherence will simply depend on τ and r. 13 The value of N max depends in turn on the post-inflationary history and conservative estimates suggest N max = 63 ± 15 [9,27]. If the post-inflationary history is standard (with sudden reheating) N max = 63.6 + 1 4 ln ǫ. In this case N max coincides approximately with the minimal number of efolds N min needed to solve the kinematic problems of the standard cosmological model. All in all the result of Eq. (3.14) is pretty robust except when N tot ∼ N max ∼ N min . In this case the correlations of the initial state may affect the degree of quantum coherence. There are some who think that for the consistency of the inflationary scenarios we must anyway demand that N tot exceeds N max .
to disambiguate the statistical properties of large-scale curvature inhomogeneities. This aspect can be already appreciated by noticing that the degree of second-order coherence of large-scale curvature inhomogeneities turns out to be numerically very close to the case of Bose-Einstein correlations.

Bose-Einstein correlations
Indeed the result of Eq. (3.14) can be usefully compared with the one of Bose-Einstein correlations in the same large-scale limit. Let us then suppose, for the sake of comparison, that the fieldq( x, τ ) is in a thermal state in standard Minkowski space-time. In this case the density matrix can be written, in the Fock basis, as: The multimode probability distribution appearing in Eq. (4.1) is given by: where n k = Tr[ρd † kd k ] is the average occupation number of each Fourier mode. Furthermore, following the standard notation, |{n} = |n k 1 |n k 2 |n k 3 ... where the ellipses stand for all the occupied modes of the field. Even though the density matrix describes a mixed state, the average multiplicity in each Fourier mode(i.e. n k ) does not need to coincide with the Bose-Einstein occupation number 14 . Using Eq. (2.9) and performing the appropriate averages in by means of the density matrix of Eq. (4.1) the normalized degree of second-order coherence becomes, in this case, While the calculation leading to Eq. (4.3) is a bit lengthy we just recall that the expectation value contains creation and annihilation operators with four different momenta; three typical leading terms will arise: in the first term all the momenta coincide while in the two remaining terms the momenta are paired two by two. The expectation value becomes 15 whered i andd † j denote the annihilation and creation operators related two generic momenta, i.e. for instanced q andd † p ; furthermore, following the same shorthand notation, δ i j denotes the delta functions over the three-momenta (i.e. δ q, p ). Performing each of the averages of Eq. (4.4) in terms of the density matrix (4.1) leads to Eq. (4.3) after careful use of the various delta functions over the momenta. We can finally integrate over the angular coordinates and the result will be 16 : In the large-scale limit Eq. (4.5) will then imply that g (2) ( r, τ ) → 2 [13]; conversely Eq. (3.14), in the same physical limit, implies g (2) ( r, τ ) → 3.

Why are we so close to Bose-Einstein correlations?
The curvature quanta follow a Bose-Einstein multiplicity distribution 17 . This is, in a nutshell, the reason why the degree of second order coherence of curvature inhomogeneities is very close to the Bose-Einstein case (i.e. 3 rather than 2). Consider the case of a two-mode squeezed state which corresponds to a single k-mode of the field. We are interested in computing m m |ζ where |ζ = Σ(ζ)|0 0 and Σ(ζ) is actually given in Eq. (C.9). We will have that the operator Σ(ζ) can be factorized as [41]: where K + , K − and K 0 obey the standard commutation relations of the SU(1, 1) Lie algebra 18 .
Using the previous equation we will have that the multiplicity distribution is given by: | m m |ζ | 2 = 1 n + 1 n n + 1 m , n = sinh 2 r. (4.7) If we now take into account that we have a two-mode state for each Fourier mode the form of the density matrix becomeŝ (4.8) 16 As in the case of Eq. (3.14) the result of Eq. (4.5) holds also in the in the zero time-delay limit τ 1 → τ 2 . 17 This occurrence can be traced back to a property of the SU (1, 1) group. The connection between the multiplicity distributions and the SU(1,1) group structure can be neatly expressed by computing, in explicit terms, the Wigner matrix element of the positive discrete series (see, in this respect, Ref. [40]) 18 Recalling the explicit expressions of Eq. (C.8) it is simple to show that [K − , K + ] = 2 K 0 and that whose diagonal elements define the multiplicity distribution which has precisely the Bose-Einstein form: (4.9) accounting for the way curvature quanta of each Fourier mode (i.e. n k ) are distributed as a function of their average multiplicity (i.e. n k = sinh 2 r k ). The behaviour of the off-diagonal elements of Eq. (4.8) is dictated by the phases α k = (2ϕ k − γ k ) whose evolution is specifically discussed in appendix C (see, in particular, Eq. (C.7) and discussion therein). This explains the nature of the multiplicity distribution which is of Bose-Einstein type but with a larger amount of correlations (hence the value 3 instead of 2 for the normalized degree of second-order coherence). Such a difference also persists to higher order. By appropriately reducing the density matrix, the multiplicity distribution can become exactly to the one of a thermal state 19 .

The quantum mechanical correspondence
The numerical values of the degrees of second-order coherence in the large-scale limit are solely determined by the statistical properties of the quantum state. This statement can be made more precise: the degrees of first-and second-order coherence in the large-scale limit coincide with the results obtainable if only a single mode of the field is excited in the limit of zero time delay between the quantum operators 20 .
Indeed, for a single mode of the field the degrees of first-and second-order coherence are defined as: 11) where the overline at the left hand side distinguishes Eqs. (4.10) and (4.11) from Eqs. (2.7) and (2.8)-(2.9) holding in the general case. Equations (4.10) and (4.11) define, respectively, 19 Even if this point might be related to the dynamical evolution of α k (see Eq. (C.7) and discussion therein) we just point out that by averaging over α k Eq. (4.8), the density matrix can be reduced (i.e. ρ red k = 1 2π 2π 0 dα kρ k ) to a the density matrix of a chaotic state. A similar observation has been made long ago in a related context [42]. 20 In actual interferometry the electric field is first split into two components through the beam splitter, then it is time-delayed and finally recombined at the correlator. The limit of zero time delay between the signals is commonly used, in both cases, to characterize the statistical properties of the source. We are here using a similar logic.
the degrees of first and second-order temporal coherence: in the zero time-delay limit τ 1 − τ 2 → 0 and in this case the degree of second-order coherence will be denoted by g (2) . Recalling the analysis of section 4 we have already demonstrated, in practice, that: lim kr≪1 g (2) (r, τ ) = g (2) . (4.12) Equation (4.12) is verified in the case of a single-mode coherent state (i.e.â|α = α|α ) where Eqs. (4.10) and (4.11) imply g (1) = g (2) = 1. Similarly for a chaotic state with statistical weights provided by the Bose-Einstein distribution the density matrix has the same form quoted before but for a single mode. Therefore we have that g (1) = 1 but g (2) = 2. In the case of a single Fock state g (2) = (1 − 1/n) < 1 showing that Fock states lead always to sub-Poissonian behaviour and they are anti-bunched [13,37]. Finally a single-mode of the field in the quasi-de Sitter case would correspond toâ = cosh rb − sinh rb † where the phases have been fixed to zero; taking the limit of zero time-delay and inserting these expressions in Eq. (4.11) we have that: g (2) = 3 + 1/n with n = sinh 2 r. In the quasi-de Sitter case the quantum correspondence holds provided the number of produced quanta in each Fourier mode is large so that terms of the order of 1/n (and higher) are all negligible since n. In this 1/n expansion the normal-ordered expectation values dominate the degrees of coherence (see also the last part of appendix C). Chaotic light is an example of bunched quantum state (i.e. g (2) (0) > 1 implying more degree of second-order coherence than in the case of a coherent state). Fock states are instead antibunched (i.e. g (2) < 1) implying a degree of second-order coherence smaller than in the case of a coherent state. The correspondence of Eq. (4.12) holds as well for the higher-order degrees of coherence, as we shall see. In the zero time-delay limit the degree of second-order coherence has a simple relation with the variance of the probability distribution associated with a given quantum state. More specifically Eq. (4.11) implies g (2) = (D 2 − n)/n 2 where n = N and D 2 = N 2 − N 2 . It is sometimes useful to define the so-called Mandel parameter whose expression is Q = n[g (2) − 1] = D 2 /n − 1. The parameter Q is directly related to the way second-order correlations are parametrized in subatomic physics [11]. In the case of a coherent state [15] we have that D 2 = N while Q = 0; this means that the distribution underlying this state is just the Poisson distribution (i.e. the variance coincides with the mean value).

Correlations from compound Poisson processes
In the perspective of the present investigation the degrees of first and second-order coherence are not sufficient to disambiguate the statistical properties of the quantum state and need to be complemented by the analysis of some higher-order degrees of coherence. We shall now show, in a specific class of examples, that very different quantum states lead to the same degrees of first-order and second-order coherence.
Let us suppose that the quantum state of cosmological perturbations is unknown but characterized by a given density matrix whose statistical weights will be denoted by P n , i.e. ρ = ∞ n=0 P n |n n|, ∞ n=0 P n = 1. (4.13) We work in the zero-delay limit and exploit the results of Eqs. (4.10), (4.11) and (4.12). The degrees of correlations can then be expressed as 21 : 14) where n = Tr[â †âρ ]; the function A(s) appearing in Eqs. (4.14)-(4.17) is the probability generating function defined as: One of the simplest situations coming to mind is the one where P n is a Poisson distribution.
In this case we have that P n = e −n n! n n , A(s) = e n(s−1) (4.19) and thanks to Eqs. (4.14)- (4.18) it is immediate to verify that the degree of nth-order coherence in the case of a Poisson distribution is always 1. Note that Eq. (4.19) does not define a coherent state but rather a mixed state with Poisson distribution. Another possible case which can be easily handled with Eqs. (4.14)-(4.18) is when P n is a Bose-Einstein distribution. In this case g (n) = n! and, as already discussed, g (2) = 2.
To reproduce the degree of second-order coherence of cosmological perturbations the simplest idea is to deviate slightly from the Bose-Einstein distribution and to consider the sum A N of a number N of mutually independent random variables, i.e.
A N = X 1 + X 2 + .... + X N . (4.20) If N is fixed the probability generating function P(s) of the sum A N is simply given by the product of the generating functions. If N is itself a random number the generating function is given by the compound distribution. Suppose that X i are identically and independently distributed random variables with probability distribution P k and generating function P(s). If the distribution of N is G k its corresponding generating function will be, as usual G(s) = k s k G k . The generating function of describing the compound process of Eq. (4.20) will then be given by A(s) = G[P(s)] [43].
If N is distributed as a Poissonian variable with multiplicity N c the probability generating function will be given by exp {N c [P(s) − 1]}. If we take P(s) to be the generating function of a logarithmic distribution, i.e.
The generating function of the random sum of random variables will then be given by: (4.22) In the case k = 1 Eq. (4.22) reproduces the probability generating function of a Bose-Einstein distribution while, in the limit k → ∞ A(s) becomes exactly the Poisson generating function. Consider now, for the sake of concreteness, the case k = 1/2: the normalized degrees of first and second-order coherence reproduce the ones of a squeezed state to leading order in 1/n, namely from Eq. (4.22) Even assuming the the initial state of cosmological perturbations is exactly the state that minimizes the Hamiltonian of the fluctuations, the degree of second-order coherence can be easily reproduced by the density matrix of an appropriately mixed state whose statistical weights correspond to a compound Poisson process where each of the independent random variables of the sum follow a logarithmic distribution. The probability generating function of Eq. (4.22) appears in the so-called pure death stochastic process when the evolution equation for the probability distributions are solved for a specific set of initial data [43]. In this case, however, the interpretation of Eq. (4.22) would be slightly different and the value of k bound to be integer.

Higher-order correlations
The degree of second-order coherence should be complemented by some informations on the higher-order correlators. In this way the statistical properties of the quantum state can be disambiguated.
The importance of the higher-order correlation functions of the detector's readings has been emphasized in [28] with the purpose of determining the coherence of the states of the field in an optical cavity. Similarly in [29] the third-and fourth-order autocorrelation functions have been used to detect the nonclassical character of the light transmitted through a photonic-crystal nanocavity containing a strongly coupled quantum dot probed with a train of coherent light pulses. The value of g (3) has been contrasted with the conventional diagnostic based on g (2) . It has been demonstrated that, in addition to being necessary for detecting two-photon states emitted by a low-intensity source, g (3) provides a more clear indication of the nonclassical character of a light source. In [29] preliminary data have been presented to demonstrate bunching in the fourth-order autocorrelation function g (4) as the first step toward detecting three-photon states. Higher-order autocorrelations are necessary to characterize the multiphoton nature of nonclassical light also in the context of the so called quantum dots [30,31]. Because of their strong interaction with light and ease of integration into optoelectronic devices, self-assembled quantum dots are promising candidates for quantum light sources, i.e. light sources leading, in our language, to g (2) ≪ 1. All in all we can say that the specific analysis of g (2) is often not sufficient to determine the statistical properties of light so that higher-order correlations have been recently analyzed [30,31].
The logic of the HBT interferometry will now be generalized by considering the correlations of three and four intensities in quasi-de Sitter space. Recalling Eq. (2.5) in the case n = 3 the third order correlator corresponding to the three-point function of the intensity of the field evaluated at the same conformal time is given by: where the six points appearing in Eq. (2.5) in the case n = 3 have been identified using the same procedure already explained when passing from Eq. (3.9) to Eq. (3.10). The various correlators have been made explicit by eliminating three out the six integrals over the momenta thanks to the three delta functions arising from each of the four expectation values appearing at the right hand side of Eq. (5.1). The result of this manipulation is given by: where, by definition, r ij = x i − x j . The degree of third order coherence can be immediately computed and it is: By evaluating Eq. (5.3) in the large-scale limit, after some algebra the above result can be expressed as It can be observed that the integral d 3 k 1 |v k 1 (τ )| 2 /k 1 appearing in the denominator of Eq. (5.4) can be factorized in the numerator so that the expression for the degree of third-order coherence can be finally expressed as: In the large-scale limit, the degree of third-order coherence becomes therefore: The reason why this is a numerical value and does not depend on the momenta shows that, in this limit, the degree of coherence only depends on the statistical properties of the quantum state, as explained in section 4. Indeed, Eq. (5.5) coincides with the result obtainable in the case of a single degree of freedom (recall the discussion after Eqs. (4.10) and (4.11)) neglecting the higher orders in the 1/n expansion: Process k g (1) g (2) g (3) g (4) Poisson Compound 1/2 1(1) 3 (3) 15 (11) 105 (93) Bose-Einstein 1 1 2 6 24 Poisson ∞ 1 1 1 1 Table 1: The degrees of coherence in various processes. In the first row we report inside the brackets the degrees of coherence computed for the quantum state of curvature perturbations.
Finally going back to Eq. (2.5) in the n = 4 case, the fourth-order correlator is given by: Equation (5.8) must be evaluated in the large-scale limit and the degree of fourth-order coherence becomes: . (5.9) After making explicit all the expectation values and after rearranging the integrations over the momenta one of the integrals simplifies between the numerator and the denominator.
Inserting the explicit expressions for u k (τ ) and v k (τ ) in the large-scale limit we have that (5.10) where N (4) (τ ) and D (4) (τ ) are defined as Again Eq. (5.10) agrees with the quantum mechanical correspondence in the 1/n expansion, i.e. g (4) → g (4) where g (4) = â †â †â †â †ââââ n 4 = 93 + 76 n + 7 n 2 . (5.12) The third and fourth-order degrees of coherence complement the informations coming from the standard HBT interferometry. This analysis shows that the statistical properties of the quantum states can then be disambiguated by examining the higher degrees of coherence as exemplified in Tab. 1. In the case of the compound Poisson process of section 4 the degrees of first-order and second-order coherence coincide with the ones of large-scale curvature modes (see the case k = 1/2 in Tab. 1). By going to higher order the results differ substantially (see the values in brackets in the last two column). In the second and third lines we illustrated, for comparison, the Bose-Einsteins case (i.e. k → 1 in Eq. (4.22)) and the Poisson case (i.e. k → ∞ in Eq. (4.22)).

Counting detectors for single CMB photons
The considerations developed here suggest the interesting possibility that counting detectors for single CMB photons could provide the observational counterpart for the considerations developed in the previous sections. While it would naively seem totally hopeless to conceive and build this kind of detectors, a variety of cryogenic technologies are being developed [44,45,46] for CMB observation at THz (1 THz = 10 3 GHz) frequencies. Furthermore the idea suggested in Ref. [9] of applying the HBT interferometry to CMB photons has been taken seriously in [45,46]. According to these references it is not unconceivable that, in the future, bunched CMB photons could be identified and collected in the THz range. It is not our purpose to review here the ideas related to photon counting THz interferometry (PCTI in the language of Ref. [46]) which is, in oversimplified terms, an advanced version of the HBT interferometry [46]. In what follows we shall instead focus the discussion on the definition of the various observables (i.e. the degrees of coherence). The sole purpose of these considerations will be to relate the degrees of quantum coherence to actual CMB observables. If we consider the curvature inhomogeneities as quantum fields, the brightness perturbations will also become operators. More specifically the curvature perturbation R is related to the field operatorq (see appendix A). The intensity fluctuations of the radiation field (i.e. the warmer and the cooler regions in the CMB sky) are given, in turn, by ∆ I ( x, τ dec ) ≃ −R( n, τ dec )/5 under the approximation of sudden decoupling (i.e. assuming that the visibility function is a narrow Gaussian centered at τ dec ) [47,48]. This correspondence will now be swiftly outlined.

Brightness perturbations
Indeed the degrees of coherence can be directly assessed by studying higher order temperature and polarization correlations in the limit of large angular scales. In the concordance scenario with no tensors the temperature anisotropies are customarily described, in real space by ∆ I (n, τ ) (the brightness perturbation of the intensity). The polarization anisotropies are instead discussed in terms of the E and B-modes, denoted, respectively, by ∆ E (n, τ ) and ∆ B (n, τ ). Both ∆ E and ∆ B are spin 0 quantities, exactly as ∆ I [49]. Recalling that µ = cos ϑ we have, in the concordance framework and in the absence of gravitons, The brightness perturbations ∆ I (n, τ ) and ∆ P (n, τ ) obey, in Fourier space, the following pair of equations: where S P (k, τ ) can be expressed as the sum of the quadrupole of the intensity, of the monopole of the polarization and of the quadrupole of the polarization, i.e., respectively, S P (k, τ ) = (∆ I2 + ∆ P0 + ∆ P2 ); note that, in Eqs. (6.2) and (6.3), ε ′ and ε(τ, τ 0 ) denote, respectively, the differential optical depth and the optical depth itself [47,48] ε ′ = x eñe a σ γe , ε(τ, τ 0 ) = τ τ 0 x eñe a σ γe dτ. (6.4) The line of sight solution of Eqs. (6.2) and (6.3) can be written, respectively, as: where K(τ ) = ε ′ e −ε(τ,τ 0 ) is the visibility function and x(τ ) = k(τ 0 − τ ). The visibility function can be approximated as a double Gaussian with two peaks roughly corresponding to the redshifts of recombination and reionization [47,48]. In the present analysis the focus will be on the large-angular scales corresponding to typical multipoles ℓ ≤ √ z rec where the finite width of the visibility function is immaterial and the opacity suddenly drops at recombination. This implies that the visibility function presents a sharp (i.e. infinitely thin peak at the recombination time). Thus, since K(τ ) is proportional to a Dirac delta function and e −ε(τ,τ 0 ) is proportional to an Heaviside theta function. Under the latter approximations, Eq. (6.2) leads to the well known pair of separated contributions, i.e. the Sachs-Wolfe (SW) and the integrated Sachs-Wolfe (ISW) contributions: (6.8) where, by definition, x(τ rec ) = y rec . The SW and the ISW contributions of Eqs. (6.7) and (6.8) can be separately evaluated. In particular, as anticipated, the ordinary SW contribution becomes: 10) while the ISW contribution is: Both in Eqs. (6.9) and (6.11), R( k, τ ) denotes the constant value of curvature perturbations at τ i < τ eq . By further approximating the integrand in Eq. (6.11) the whole large-scale contribution can be written, for the present purposes, as where q rec = a rec /a eq ≃ 3.04 for h 2 0 Ω M0 ≃ 0.134 (as usual Ω M0 denotes the total matter density in critical units at the present time). Note that if we neglect the ISW term and assume that q ≫ 1 we obtain the conventional simplified form of the SW contribution.
If the curvature perturbations are described by a quantum field operator, also the brightness perturbations must be treated quantum mechanically. The degrees of coherence defined in Eqs. (2.9), (2.10) and (2.11) can therefore be written as: 14) , (6.15) and similarly for g (4) (m,n,r,ŝ). In Eqs. (6.14) and (6.15) we decomposed the positive and negative frequency parts of the operators as∆(n) =∆ (+) (n) +∆ (−) (n); note finally that ∆(n) can denote, indifferently,∆ I of∆ E (recall that we are here discussing the case where the B-mode polarization is absent). This means that we can define separate degrees of coherence for the temperature, for the E-mode polarization and for their cross-correlation.
To give an example we can compute, for instance, the degree of second-order coherence of the temperature. By substituting∆ with∆ I in Eq. (6.14) we shall have that , A ℓ (q,m ·n) = (2ℓ + 1)j 2 ℓ (kτ 0 )P ℓ (m ·n), (6.16) where j ℓ (kτ 0 ) are the spherical Bessel functions and P ℓ (m ·n) the Legendre polynomials. Equation (6.16) directly relates the curvature perturbations to the brightness perturbations. All the inequalities established for the degree of second-order coherence and all the considerations presented before are also applicable to Eq. (6.16): in the limitm ·n → 1, using well known identities, the degree of coherence coincides with the results obtained before. We shall skip these discussions which can be found elsewhere [9].

Multiplicity distributions and all-order observables
We conclude with some speculation concerning the multiplicity distributions of the total number of phonons. If observed these distributions might usefully complement the interferometric approach discussed in the previous subsection.
In the hypothesis that curvature quanta are the sole source of temperature inhomogeneities the distribution P (n) of the total number of phonons n = k n k must reflect the distribution of the warmer and cooler regions. The multiplicity distribution P (n) accounts for the way the total number of phonons n is distributed as a function of its mean value; P (n) can be very different from the P {n k } which has been discussed in Eq. (4.9). Denoting with p({n}) the joint probability distribution of the set of phonon occupation numbers {n} of the field, we shall have that [13] p({n}) = k 1 (1 + n k )(1 + 1/n k ) n k . (6.17) For any mode for which n k = 0, the corresponding factor must be interpreted as δ n k 0 . In the following we shall suppose, quite generally, that only a subset consisting of ǫ modes of the field is actually occupied and we shall restrict the attention to this subset of modes.
If n is the total number of phonons and P (n) is the multiplicity distribution of n, then P (m) = {n} p({n})δ m n . In quantum optics an analog of the multiplicity distribution P (m) describes the statistical properties of (unpolarized) chaotic light beams [13]. The evaluation of P (m) can be in general difficult but it becomes easy in the physical case when the average occupation number n k of all the ǫ occupied modes become equal 22 ; in this case the joint probability distribution of the occupied modes becomes: Every non-vanishing term in the summation P (m) = {n} p({n})δ m n has the same value and the required probability P (m) is simply p({n}) given by Eq. (6.18) multiplied by a well known combinatorial factor accounting for the way the n phonons are distributed among the ǫ modes: P n (n, ǫ) = Γ(n + ǫ) Γ(ǫ)Γ(n + 1) n n + ǫ n ǫ n + ǫ ǫ . (6.19) How could we estimate the probability distribution from actual data? Is it somehow related to the statistics of cold and hot spots? Leaving aside these important questions it is clear that the correct limit where Eq. (6.19) should be analyzed is given by n ≫ 1 and n ≫ 1. More specifically, the realizations of the distribution (6.19) with different n. In the limit n ≫ ǫ the negative binomial distribution becomes a Gamma distribution with a characteristic scaling law: nP n (n, ǫ) ≃ ψ ǫ (w), ψ ǫ (w) = ǫ ǫ Γ(ǫ) w ǫ−1 e −ǫw , w = n/n. (6.20) discussed in quantum optics and also relevant to in high-energy physics [50]. The scaling limit is defined by the asymptotic behaviour of P n (n) as n → ∞, n → ∞ while n/n is kept fixed. For fixed ǫ the distribution of Eq. (6.19) scales as ψ ǫ (w). 22 In the case of a thermal light beam which is either fully polarized or fully unpolarized the use of a rectangular spectral density is an excellent approximation in the derivation of the photocounting statistics which is also experimentally accessible [13].

Concluding remarks
There are some who think that the quantum fluctuations are the only plausible initial data for the evolution of the metric perturbations. The complementary approach pursued here has been to scrutinize the coherence properties of large-scale curvature inhomogeneities when the relevant wavelengths are larger than the Hubble radius. This step is both plausible and mandatory if we ought to distinguish the quantum nature of large-scale curvature inhomogeneities from more mundane statistical correlations.
Since the pioneering attempts of Hanbury Brown and Twiss it has been realized that the analysis of first order interference between the amplitudes cannot be used to distinguish the nature of different quantum states of the radiation field. Young interferometry is not able, by itself, to provide information on the statistical properties of the quantum correlations since various states with diverse physical properties (such as laser light and chaotic light) may lead to comparable degrees of first-order coherence. The Hanbury Brown-Twiss interferometry demand the analysis of the second-order correlations of corresponding intensities. The degree of second-order coherence, however, does not describe unambiguously the correlation properties of large-scale curvature perturbations. The applications of the Glauber theory allow for a systematic scrutiny of the higher-order correlations. More specifically the analysis of the degrees of third-and fourth-order coherence is necessary to assess the statistical properties of curvature inhomogeneities and their plausible quantum origin. Similar conclusions have been recently drawn in a quantum optical context where higher-order autocorrelations are mandatory for characterizing the multiphoton nature of nonclassical light.
Since we showed that the statistical properties of the quantum states can be disambiguated by examining the higher degrees of coherence, it would be interesting to tailor specific observational strategies along this direction. The new generations of CMB detectors and the Hanbury Brown-Twiss interferometry in the THz region could be a reasonable hope in this respect. While this possibility seems rather interesting, there are some who might object that an improved statistical accuracy in the determination of the parameters of the concordance paradigm should be the primary target of future experimental endeavours. A full account of this potentially interesting debate is beyond the scopes of the present investigation.
− λ * (τ )u k (τ )v k (τ ) (b † pb p +b − pb of the conformal time coordinate τ . Using Eq. (C.6) into Eqs. (C.1)-(C.2), the evolution of the squeezing amplitude r k and of the two phases becomes: where α k = 2ϕ k − γ k ; combining the previous equations it is also simple to obtain α ′ k = 2k + 2F sin α k / tanh 2r k (see also Refs. [9,23]). Equation (C.6) can be swiftly obtained by considering a single p-mode and by noticing that the operators K ± and K 0 obey the commutation relations of the SU(1, 1) Lie algebra: Using the the standard Campbell-Baker-Hausdorff theorem [13], Eq. (C.8) implieŝ a = Σ † (ζ) Ξ † (ϕ)b 1 Ξ(ϕ)Σ(ζ) = e −iϕ cosh rb 1 − e iγ sinh rb † 2 , (C.9) where Ξ(ϕ) and Σ(ζ) (with ζ = re iγ ) are, respectively, the rotation operator and the twomode squeezing operators [33,34] defined in terms of the generators of the SU(1, 1) Lie algebra: These two operators play an important role in the evolution of the states in the Schrödinger representation; their use has been pioneered by Grishchuk and Sidorov [35] (see also [36]). We conclude by proving that the normal-ordered expectation values dominate the degrees of coherence provided the average multiplicity of the produced phonons and gravitons is large; this has been dubbed in the bulk of the paper as the 1/n expansion. Consider, for this purpose, the following normalized expectation values: where the superscripts denote, respectively, the normal, antinormal and symmetric ordering. Using the commutation relations the first two relations of Eq. (C.11) can be written in terms of the normal-ordered correlator. For the antinormal correlator we have g (A) = g (N ) + 3/n + 2/n 2 while for the symmetric correlator we get g (S) = g (N ) + 2/n where, as usual, n = N is the average multiplicity. Since n ≫ 1 all the terms O(1/n) can be neglected. In quantum optics it is natural to impose the normal ordering in the higher-order correlators since the detection of light quanta (i.e. in the optical range of frequencies) occurs by measuring a photo-current, i.e. a current induced by the absorption of a photon. Notice finally that using the commutation relations we can easily obtain g (N ) = 1 + (D 2 − n)/n 2 where D 2 = N 2 − N 2 denotes in the bulk of the paper the variance of the underlying multiplicity distribution.