Squeezed relic photons beyond the horizon

Owing to the analogy with the ordinary photons in the visible range of the electromagnetic spectrum, the Glauber theory is generalized to address the quantum coherence of the gauge field fluctuations parametrically amplified during an inflationary stage of expansion. The first and second degrees of quantum coherence of relic photons are then computed beyond the effective horizon defined by the evolution of the susceptibility. In the zero-delay limit the Hanbury Brown-Twiss correlations exhibit a super-Poissonian statistics which is however different from the conventional results of the single-mode approximation customarily employed, in quantum optics, to classify the coherence properties of visible light. While in the case of large-scale curvature perturbations the degrees of quantum coherence coincide with the naive expectation of the single-mode approximation, the net degree of second-order coherence computed for the relic photons diminishes thanks to the effect of the polarizations. We suggest that the Hanbury Brown-twiss correlations are probably the only tool to assess the quantum or classical origin of the large-scale magnetic fluctuations and of the corresponding curvature perturbations.


Introduction
The squeezed states of optical photons arise in a number of diverse physical situations all related (directly or indirectly) to the quantum theory of the parametric amplification [1]. The formulation of the quantum theory of optical coherence [2,3,4,5] paved the way for the first quantum description of parametric amplification [6]. Since then various complementary descriptions of quantum amplifiers have been developed through the years [7,8,9,10] both in the context of single-mode and two-mode squeezed states (see also [11,12] for an incomplete list of review articles on the subject).
After the seminal discoveries of the COBE satellite [13] (later confirmed and extended by the WMAP experiment [14,15]) it became gradually clear that the early Universe itself could be seen, from the physical viewpoint, as an effective quantum amplifier. Consequently the applications of quantum optical techniques to the analysis of large-scale inhomogeneities has been firstly suggested by Grishchuk and collaborators in a class of problems involving the evolution of the tensor and scalar modes of the four-dimensional geometry [16,17,18]. Neither the tensor [19] nor the scalar [20,21,22] inhomogeneities of a conformally flat geometry are invariant under Weyl rescaling of the four-dimensional metric. The lack of Weyl invariance implies then the formation of squeezed states of the relic gravitons and of the relic phonons [16,17,18] (see also [23] for a review article). The key physical assumption behind these attempts rests on the quantum mechanical nature of the initial conditions of the large-scale inhomogeneities, as suggested long ago by Sakharov [24] even prior to the formulation of the conventional inflationary paradigms.
The quantum theory of parametric amplification has been later applied to the case of relic photons [25] where the quantum optical analogy is even more compelling: in this case it is precisely the time variation of the susceptibility that plays the role of the laser pump often employed for the direct experimental preparation of the squeezed states in various classes of nonlinear materials (see e.g. [1,11,12] and also [26]). The quantum theory of parametric amplification of the relic photons (but also of the relic gravitons and relic phonons) is useful for treating the problem of initial data but it becomes essential for analyzing the higher-order correlations of the large-scale fluctuations, as the quantum optical analogy clearly suggests.
There are some who argue that we have already an accurate control of the protoinflationary dynamics; along this prespective a consistent model suffices for claiming that the large-scale fluctuations have a quantum origin. In spite of this belief, it would be nice (and probably even mandatory) to develop a more objective set of sufficient criteria enabling us to infer the quantum origin of large-scale fluctuations of any spin from some sort of observational evidence. The first idea coming to mind, in this respect, it is to analyze the quantum coherence of the fluctuations in the spirit of the Glauber theory [2,3,4]. Only by looking at the higher-order correlations we shall be able, at least in principle, to establish if the large-scale curvature perturbations have a classical or a quantum origin as speculated by Sakharov [24].
A first step along this direction relies on the idea of studying (and eventually measuring) the correlation functions of the intensities of the curvature perturbations rather than the correlations of the corresponding amplitudes [27]. This concept has been originally proposed by Hanbury Brown and Twiss [28] and their 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. The applications of the HBT effect range from stellar astronomy [28] (see also [29]) to subatomic physics [31] where the interference of the intensities has been used to determine the hadron fireball dimensions [32] corresponding, in rough terms, to the linear size of the interaction region in proton-proton collisions.
In this paper the quantum theory of optical coherence is applied to the scrutiny of the statistical properties of the relic photons produced thanks to the pumping action of the susceptibility during an inflationary stage of expansion. The idea is to define the Glauber correlation functions and to focus the attention on their large-scale limit. The first and second degrees of quantum coherence correspond, in the quantum optical analogy, to the Young interferometry and to the HBT interferometry. In the zero-time delay limit the degree of second-order coherence (conventionally denoted by g (2) in quantum optics [1]) can be used to infer the statistical properties of the quantum state. In the standard lore, based on the so-called single mode approximation [1], g (2) → 1 for a coherent state (also referred to as the Poissonian limit because of the well known statistical properties of the coherent states). Conversely in the chaotic (or thermal) case we would have g (2) → 2; finally in the case of twomode squeezed states g (2) → 3 signalling a super-Poissonian but also superchaotic statistics. By comparing the the Hanbury-Brown Twiss correlations computed in the scalar case (and, more precisely, for the large-scale curvature fluctuations) with the case of relic photons we find specific physical differences which are traced back to the role of the polarizations.
The plan of the present paper is the following. In section 2 we shall discuss the squeezed states of the relic photons. In section 3 the essentials of the Glauber approach will be introduced. The large-scale limits of the correlation functions will be studied in section 4. In section 5 the physical meaning of the degrees of quantum coherence will be specifically computed and contrasted with the single-mode approximation. Section 6 contains our concluding remarks. To avoid digressions, various useful details have been relegated to the appendix.

Squeezed states of relic photons
The conformally invariant coupling of the Abelian gauge fields is broken in different situations that can be usefully recapitulated in terms of the general action [33]: where M ρ σ (ϕ, ψ) and N ρ σ (ϕ, ψ) may depend on a number of different scalar fields and on their covariant derivatives. In a complementary perspective they can be constructed directly from fluid variables (i.e. fluid velocities, vorticities and shear). In spite of their specific form, when M ρ σ = N ρ σ the system is characterized by different electric and magnetic susceptibilities; in this situation Eq. (2.1) includes, as a special case, the derivative couplings arising in the relativistic theory of Casimir-Polder and Van der Waals interactions [34]. We shall be assuming, consistently with the observations, that the evolution of the large-scale magnetic fields takes place in a conformally flat background geometry g µν = a 2 (τ )η µν where η µν denotes the Minkowski metric, a(τ ) is the scale factor and τ denotes the conformal time coordinate. If M ρ σ = N ρ σ the comoving electric and magnetic fields obey the following set of equations: The electric and magnetic couplings are, respectively, g E = (4π/Λ E ) 1/2 and g B = (4π/Λ B ) 1/2 . Under the exchange and inversion of the susceptibilities ( where F = χ /χ, χ = √ λ is the susceptibility and the prime denotes a derivation with respect to the conformal time coordinate. The components of the Abelian field strength of Eq. (2.1) are defined as Y 0i = e i /a 2 and Y ij = − ijk b k /a 2 . The canonical electric and magnetic fields appearing in Eq. (2.5) are then given by B = a 2 √ λ b and E = a 2 √ λ e. Note that the two equations appearing in Eq. (2.5) are left invariant by the duality transformations [35] χ → 1/χ (i.e. F → −F) provided E → − B and B → E. The continuous evolution of F will define an effective horizon for the gauge modes related to E and B.
In time-dependent (conformally flat) backgrounds the Coulomb gauge (i.e. Y 0 = 0 and ∇ · Y = 0) is preserved (unlike the Lorentz gauge condition) under a conformal rescaling of the metric. For the quantum mechanical description of the problem we can therefore start with the canonical Hamiltonian (see appendix A for a derivation) (2.6) where ξ = iF/2. Equation (2.6) is reminiscent of the toy model of parametric amplifier analyzed, for the first time by Mollow and Glauber [6]. The free part of Eq. (2.6) and the two components of the interacting Hamiltonian satisfy the usual commutation relations of the SU (1, 1) Lie algebra, as we shall see in a moment. Equation (2.6) describes an interacting Bose gas at zero temperature. In this case the free Hamiltonian corresponds to the kinetic energy while the interaction terms account for the two-body collisions with small momentum transfer [43,44]. The Hamiltonian (2.6) is invariant under duality that transforms χ in its inverse, i.e. χ → 1/χ. Under this transformation we have that F → −F while the creation and annihilation operators transform as:â Recalling the notations discussed in appendix A, the Fourier representation of the field operators and of the momentâ transform asÂ k α →π k α k ,π k α → −kÂ k α (2.10) if we use Eqs. (2.7) and (2.8). In the present discussion the vacuum corresponds to the state minimizing the Hamiltonian at the onset of the dynamical evolution. 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 [43,44].
The evolution ofâ k α andâ † k α can be obtained from Eq. (2.6) and from the evolution equations in the Heisenberg description: The formal solution of Eq. (2.11) iŝ 12) where u p (τ ) and v p (τ ) satisfy |u p (τ )| 2 − |v p (τ )| 2 = 1. From Eq. (2.11) the equations obeyed by u p and v p can be written as: The solution for the evolution equations of u p (τ ) and v p (τ ) can be obtained in two complementary regions, namely for the wavelengths larger than the effective horizon (i.e. p/F 1) and for wavelengths shorter than the effective horizon (i.e. p/F 1). In the short wavelength region the solutions of Eq. (2.13) are plane waves e ±ipτ while in the long wavelength regime the solution becomes: (2.15) where A k (τ, τ ex ) and B k (τ, τ ex ) are given by: The two dimensionless integrals I B (τ ex , τ ) and I E (τ ex , τ ) are given by Thanks to Eqs. (2.14) and (2.15) the initial conditions for the evolution can be directly expressed at τ ex and can be written in terms of the values of the mode functions at the corresponding epoch (i.e. u p (τ ex ) ≡ u p and v * p (τ ex ) ≡ v * p ).
We can now remark that the two complex functions u p (τ ) and v p (τ ) (subjected to the constraint |u p (τ )| 2 − |v p (τ )| 2 = 1) can the be parametrized in terms of three real functions. The evolution of u k and v k can then be rephrased in terms of the so-called squeezing parameters [1,11,12] (see also [7,8,9,10]): 19) where ϕ p , r p and γ p are all functions of the conformal time coordinate τ even if the arguments of the functions will be dropped for the sake of conciseness. Using Eq. (2.19), Eq. (2.12) can be rewritten as:â Equation (2.20) 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 [1,45], Eq. (2.21) implieŝ where Ξ(ϕ) and Σ(ζ) (with ζ = re iγ ) are, respectively, the rotation operator and the twomode squeezing operators defined in terms of the generators of the SU (1, 1) Lie algebra: These two operators describe the evolution of the states in the Schrödinger representation; their use has been pioneered by Grishchuk and Sidorov [16] (see also [25] in the case of the relic photons). Using Eq. (2.20) into Eqs. (2.13), the evolution of the squeezing amplitude r k and of the phase ϕ p becomes: where α p = 2ϕ p − γ p and the relation between γ p and ϕ p is given by: By combining Eqs. (2.24) and (2.25) it is immediate to obtain α p = 2p + 2F sin α p tanh 2r p . (2.26) 3 Glauber description of quantum coherence

General form of the Glauber correlation function
The statistical properties of any quantum state and its degrees of quantum coherence can be used to reconstruct, at least in partially, the physical nature of the source [1,2,3,4].
In quantum optics the Glauber theory is often used in an exclusive manner: the statistical properties of visible light are reduced to the study of a single mode of the field. This is what goes under the name of single-mode approximation. Conversely, in the analysis of the large-scale cosmological fluctuations of different spin, a more inclusive approach is needed since the correlation functions contain all the modes of the field. In its most general form the Glauber correlation function can be written as [2,4]: . By definition we will have thatÂ (+) i (x)|vac = 0 and also that vac|Â (−) i (x) = 0; the state |vac denotes the vacuum. The vacuum corresponds to the state minimizing the Hamiltonian at the onset of the dynamical evolution. 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 [43,44]. Provided the total duration of inflation exceeds the minimal number of about 65 efolds, the vacuum initial data are the most plausible, at least in the conventional lore (see, however, Ref. [27] for different choices in a related context). The correlation function defined in Eq. (3.1) depends on the polarizations as the free indices clearly show. It is useful, for future convenience, to introduce the Glauber correlation function for a scalar degree of freedom. In this case Eq. (3.1) simply becomes: The quantum fieldq(x) defines, for instance, the normalized curvature perturbations on comoving orthogonal hypersurfaces.
It is relevant to remark that Eq. (3.1) (and, similarly, Eq. (3.2)) contain an operator that can be written as: 3) The operator of Eq. (3.3) 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 ofÂ i at each detector and at given times is a |Â The rate at which such absorptions occur, summed over the final states, is therefore proportional to [1,2,4]: where the second equality of Eq. (3.4) follows from the completeness relation.

Symmetric form of the correlation function
According to Eq. (3.4), when b|Ô|b is averaged over the ensemble of the initial states of the system it becomes identical with Eq. (3.1) 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 symmetric form of the Glauber correlation function as: Thanks to Eq. (3.5), the coherence properties of the quantum fieldÂ i (x) can be discussed by introducing the normalized version of the n-point Glauber function [2,4]: . (3.6) While, by definition, |g 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 [1,2,4]: for n = 1, 2, 3, . . .. If only a limited number of normalized correlation functions will satisfy Eq. (3.7) we shall speak about partial coherence. The degrees of first-and second-order coherence can be written as: , (3.9) where, in agreement with the general definitions of Eq. (3.1), the correlation functions appearing in Eqs. (3.8) and (3.9) are given by: In a similar manner it is possible to define, for instance the third-and fourth-order degrees of coherence , (3.12) where, for the sake of simplicity, we just suppressed the polarization indices. If g (1) the quantum field is second-order coherent. We shall be interested in the first and second degrees of coherence even if 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 [46,47,48,49].

Electric and magnetic correlation functions
The Glauber correlation function of Eq. (3.5) has been originally defined not in terms of the vector potentials but rather using the electric fields: From Eq. (3.13) the corresponding degrees of second-order coherence can also be defined. Equation (3.5) has been instead proposed as basic correlator in the approach of Mandel and Wolf [1]. Both approaches are somewhat convenient for applications to questions relating to photoelectric detection of light fluctuations. In the present context exactly the same discussion can be carried on in the case of the magnetic correlator defined as: in (x 2n ) . (3.14) From Eqs. (3.13) and (3.14) the normalized degrees of quantum coherence can be easily defined from the expressions already derived 3 using Eq. (3.5).
The degree of first-order coherence of Eq. (3.8) appears naturally in the Young two-slit experiment and whenever the degree of first-order coherence is equal to 1 the visibility is maximized [1]. The degree of second-order coherence of Eq. (3.9) enters the discussion of the Hanbury Brown-Twiss effect [28] and its different applications ranging from stellar interferometry [1] to high-energy physics [31,32]. The degree of second-order coherence arises naturally when discussing the correlations of the intensities of the fieldsÂ i ,Ê i and B i . Notice that the intensity correlators relevant to the HBT interferometry can be easily obtained from Eqs. (3.6) and (3.9) by identifying the space-time points as follows: (3.15) 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. The same observation can be made in the case of Eqs. (3.13) and (3.14). The explicit expressions of the HBT correlators can then be written from Eqs. (3.1), (3.13) and (3.14) with the help of Eq. (3.15): where the sum over repeated indices is pleonastic since the usual convention of sum over repeated indices has been adopted throughout. Nonetheless the explicit form of Eqs. (3.16), (3.16) and (3.18) can be revealing when compared with the explicit form of Eq. (3.2) in the case of HBT correlations: The difference between Eqs. (3.16)-(3.18) and Eq. (3.19) will have a direct repercussion on the large-scale limits of the degree of quantum coherence, as we shall see in the following section.

Quantum correlators beyond the effective horizon
The correlation functions introduced in section 3 will now be computed in the case of the squeezed quantum states associated with the Hamiltonian of Eq. (2.6). To avoid digressions some of the relevant details have been relegated in the appendices B and C.

Explicit form of the correlators
In the case n = 1, Eqs. (3.5) and (3.13)- (3.14) give the the explicit expressions of the first-order correlators: ii (x 1 , x 2 ). (4.3) Within the notations Eq. (4.3), the corresponding degrees of first-order electric and magnetic coherence are, respectively, .

(4.5)
As a consequence of Eq. (4.2) we also have that g B (x 1 , x 2 ). Equations (3.5), (3.13) and (3.14) give the degree of second-order coherence when written in the case n = 2. More specifically, when n = 2 Eq. (3.14) is given by Eq. (B.4) of the appendix; then, after making explicit the expectation values (see Eq. (B.5)) the final result is: For the present ends and as a preparation for the discussion of the last part of section 5, it is relevant to contrast Eq. (4.6) with the degree of second-order coherence obtainable in the case of a scalar field [27]. The Hamiltonian coincides, in this case, with Eq. (2.6) but the sum over the polarizations and the polarization dependence of the creation and annihilation operators are absent. When m = n = 2 Eq. (3.2) implies: where the results of Eqs. (B.6) and (B.7) have been taken into account. Equations (4.6) and (4.7) are similar but the polarizations introduce a quantitive difference which is even more apparent when Eq. (4.6) is written in explicit terms: While the electric and the magnetic correlators of Eqs. (3.13) and (3.14) lead to the same results (i.e. Eq. (4.8)), if we use the vector potential as pivotal variable (as suggested, for instance, in [1]) we get, formally, a different correlator. However, the expressions of Eqs. (4.6) and (4.9) are equivalent and only differ in the contribution of the phase space. Furthermore these differences are immaterial when estimating the degree of second-order coherence in the large-scale limit (see section 5).

Continuity of the effective horizon
For a reliable implementation of the large-scale limit of the degrees of quantum coherence, a continuous evolution of the extrinsic curvature, of the susceptibility and of the effective horizon is mandatory. For this purpose we shall consider the following expressions for the scale factors across the inflationary transition 4 : The continuous evolution of χ can then be parametrized in two complementary ways. In the first case the susceptibility approaches exponentially the constant asymptote and the evolution of χ(τ ) across the boundary τ = −τ i will then be parametrized as 5 : From the explicit expressions of Eqs. (4.11) and (4.12) we have that χ i (−τ i ) = χ r (−τ i ) and, similarly, χ i (−τ i ) = χ r (−τ i ) implying that both the functions and their first derivatives 4 Note that the γ appearing in Eq. (4.10) has nothing to do with the γ p appearing in Eqs. (2.19)-(2.25).
This remark avoids potential confusions. 5 If the solution (4.11) is simply matched to a constant value of χ for τ > −τ i the first derivative will be discontinuous while the second derivative of χ at the transition will be singular. All the parametrizations must then contain a transition regime (as in Eqs. (4.11)-(4.12) and (4.13)-(4.14)) which can be studied, though, in the sudden limit (i.e., respectively, for β 1 and α 1).
are continuous. The continuity of the susceptibility and of its first derivative implies the continuity of F = χ /χ. In the cosmic time parametrization we shall have that F = aF where F =χ/χ and the overdot denotes a derivation with respect to the cosmic time coordinate t. In Eq. (4.12) the rate with which the constant value χ 1 is approached is controlled by β. The interesting physical limit will be the one where β 1: in this limit the transition is continuous but it occurs suddenly.
The same sudden limit can be studied using a power-law parametrization for the transition regime, like, for instance: In Eq. (4.14) the parameter α ≥ 1 plays the same role of β in Eqs. (4.11) and (4.12): it controls the rate of the transition in the intermediate regime and as α increases the transition gets more sudden. The expressions of Eqs. (4.13) and (4.14) are continuous and differentiable, as it can b explicitly checked i.e. χ inf (−τ i ) = χ rad (−τ i ) and χ inf (−τ i ) = χ rad (−τ i ). In spite of the different analytical details, the parametrizations of Eqs. (4.11)-(4.12) and (4.13)-(4.14) lead to the same results in the sudden limit. In numerical studies of the problem (see e.g. third paper of [37]) the continuous evolution of the susceptibility and of the effective horizon have been always enforced even if there are some who confuse the sudden approximation (i.e. the regime β 1 or α 1) with a discontinuity of the effective horizon.

Evolution of the squeezing parameters
According to Eqs. (2.19) and (2.24)-(2.26) the evolution r p , γ p and α p follows directly from u p and v p : Eqs. (2.24)-(2.26) have been derived from Eq. (2.13) by means of Eq. (2.19). However, instead of solving Eqs. (2.24)-(2.26) it is more practical to derive u p (τ ) and v p (τ ), rephrase the result in terms of the squeezing parameters and take, when needed, the largescale limit. In this procedure Eq. (2.13) and Eqs. (2.24)-(2.26) can be used interchangeably in order to simply some of the asymptotic expressions.
The same strategy leading to Eq. (4.15) could also be employed in the regime τ > −τ i ; the idea would be to insert Eqs. (4.12) (or (4.14)) inside Eq. (2.13) and then deduce the corresponding solutions. However, if χ scales with (τ /τ i ) (i.e. χ = χ(z) with z = τ /τ i ) the equation for (u k + v * k ) obeys, in spite of the functional form of χ(z) where z = τ /τ i is the scaling variable. Provided the transition occurs through a scaling period where χ = χ(τ /τ i ), the first term inside the square bracket of Eq. (4.16) is always negligible: kτ i is at most of order 1 since the largest amplified wavenumber is O(1/τ i ). In similar terms we also have The solution of Eqs. (4.16) and (4.17) to lowest order in k 2 τ 2 i can be written as 7 : For an analytically tractable solution it is practical to use an explicit profile such as the one of Eq. (4.12). The full solution for τ > −τ i is therefore given by 8 : ν (x i ) + 7 Equations (4.18) and (4.19) hold under the condition kτ i ≤ 1 which is is verified for all the amplified modes of the spectrum; this condition is less stringent than the usual requirement that the modes are larger than the effective horizon (i.e. kτ < 1). 8 While this solution holds in the case of the profile (4.11)-(4.12) a similar result can be obtained in the case of Eqs. (4.13)-(4.14) but, for the sake of conciseness, the details will be skipped. (4.21) where, for simplicity, we defined C = 1 − (1 − 2ν)/(2β) and D = (1 − 2ν)/(2β); for τ ≥ −τ i the solutions u (rad) k (τ ) and v (rad) k (τ ) have been denoted, respectively, by u k (τ ) and v k (τ ). It follows from Eqs. (4.20) and (4.21) that |u k (x i , τ )| 2 −|v k (x i , τ )| 2 = 1. Note that the obtained solution, as required, is continuous and differentiable everywhere and, in particular, at the transition point τ = −τ i (recall, for this purpose, that C + D = 1).

Crossing of the effective horizon
The condition defining the time when a given mode reenters the effective horizon is obtained by requiring χ rad /χ rad = k 2 ; the latter condition implies: where Eq. (4.12) has been explicitly used. Equation (4.22) defines the crossing of the effective horizon as a function of x i = kτ i . Since kτ i ≤ 1 (kτ i 1 for the typical scale of the gravitational collapse) we will have that (4.23) where kτ i = k/(a i H i ). To get an idea of the accuracy of this expansion we can compute k/(a i H i ) in terms of the fiducial parameters of the concordance scenario: 4.24) where A R is the amplitude of the power spectrum of scalar fluctuations at the pivot scale k p = 0.002 Mpc −1 .
To compute the degrees of quantum coherence we must fix a reference time and we shall take this reference time to coincide with τ re . Alternatively one can keep the time-scale generic and expand the relevant correlation functions in the limit x i 1. Inserting then Eq. (4.22) into Eq. (4.20) and (4.21) we obtain 9 These equations are still exact but they can be expanded around the effective horizon. From Eqs. (4.25) and (4.26) the squeezing parameters can be obtained, as we shall now show. Since, by definition n k (x i , ν, β) = |v k (x i , ν, β)| 2 the average multiplicity can be computed by expanding, at once, the whole expression: where we extensively used that C + D = 1 (and hence that C 2 − D 2 = C − D). The same result of Eq. (4.27) can be obtained if we expand around the effective horizon but keep the Hankel functions in their exact form. The result of this procedure is: . (4.28) Recalling that e −iϕ k (τ ) |u k (τ )| = u k (τ ) we can express ϕ k (τ re ) in a closed form: Equation (4.29) can be obtained by writing u k (x i , ν, β) as 30) where Q(x i , ν, β) and P (x i , ν, β) are both real and given by: 4.31) Exactly with the same strategy we can compute γ k which is given by . (4.32) By combining Eqs. (4.29) and (4.32) we also have that With the results obtained so far we shall be able to discuss in detail the degrees of first-order and second-order coherence.

Degrees of coherence in the large-scale limit
The degrees of first-order and second-order coherence will now be computed. We shall then contrast the results with the benchmark values obtained in the context of the single-mode approximation.

First-order coherence
From the discussion of section 4, the degrees of coherence can be computed at any time τ i < τ ≤ τ re but the most relevant reference time is τ = O(τ re ); in this case, v k (τ ) and u k (τ ) are given by Eqs. (4.20) and (4.21). From Eqs. (4.2) and (4.3), after angular integration, the first-order correlation function at separate space-time points are where j 0 (k 1 r) denotes the spherical Bessel function of zeroth order [50]. From Eqs. (5.1)-(5.2) the normalized degree of first-order coherence defined in Eqs. (4.4) and (4.5) becomes: G ( x 1 , x 2 ; τ 1 , τ 2 ) = . 3) depend on τ 1 and τ 2 but, as a consequence of Eqs. (4.20) and (4.21), this dependence simplifies when computing the degrees of quantum coherence in the large-scale limit. Therefore the final form of Eqs. (5.3)-(5.4) can be written as: E (r) = dp p 5−2ν j 0 (pr) dp p 5−2ν → 1, G (r) = dp p 3−2ν j 0 (pr) dp p 3−2ν → 1, (5.6) where the integrals are evaluated over all the modes larger than the effective Hubble radius and the second relation clearly holds in the limit k 1 r 1 (corresponding to large angular separations). Equations (5.5) and (5.6) remain clearly valid in the zero time-delay limit (i.e. τ 1 → τ 2 ).
In quantum optics the numerical values of the degrees of first-and second-order coherence are customarily classified by considering a single mode of the field and a single polarization 10 . For a single mode of the field the degrees of first-and second-order coherence are defined as: 17) where the overline at the left hand side distinguishes Eqs. (5.16) and (5.17) from Eqs. (5.8) and (5.9) holding in the general case. Equations (5.16) and (5.17) define, respectively, 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) . For a singlemode coherent state (i.e.â|α = α|α ), Eqs. (5.16) and (5.17) imply (5.18) so that a coherent state is both first-order and second-order coherent in the single mode approximation. For a chaotic state in the single approximation the statistical weights of the the density matrix are provided by the Bose-Einstein distribution [1,30] and the results for the degrees of coherence imply: (5.19) so that the degree of second-order coherence is twice the result of a coherent state. In the case of a Fock state g (2) = (1 − 1/n) < 1 showing that Fock states lead always to sub-Poissonian behaviour and they are anti-bunched [1,30]. Let us now come to the most interesting case for the present discussion, namely the case of a squeezed state [26], corresponding 11 to 10 This approximation is often referred to as single mode quantum optics (see, e.g. chapter 5 of Ref. [30] and also [1]). The rationale for this approximation is that many experiments use plane parallel light beams whose transverse intensity profiles are not important for the measured quantities. As a consequence it is often sufficient in interpreting the data to consider the light beams as exciting a single mode of the field. 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. 11 For simplicity, the phases have been fixed to zero since they do not affect the degree of second-order coherence in the single-mode approximation. â = cosh rb − sinh rb † . Taking the limit of zero time-delay and inserting these expressions in Eq. (5.17) we have that: (5.20) Equation (5.20) also implies that in the limit n 1 the degree of second-order coherence goes to 3.
In the single-mode approximation, chaotic light is an example of bunched quantum state (i.e. g (2) > 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. Finally squeezed light is bunched and also superchaotic, meaning that the degree of second-order coherence is larger than in the case of thermal state.
Based on the single-mode approximation, we have that the degree of second-order coherence of our problem should have implied that g (2) X → 3, for X = B, E, G. We instead obtained g (2) X → 5/3 (and g (2) X → 1). The reason for this apparent disagreement stems from the contribution of the polarizations to the degree of second-order coherence.
To prove this statement let us consider the case of a scalar field. For this analysis we shall adapt the results of Ref. [27] valid in the case of the scalar modes of the geometry. Recalling the results of Eqs. (3.2), (4.7) the correlation function of Eq. (B.6) ( when x 1 = x 3 and x 2 = x 4 ) describes the interference of two beams with intensitiesÎ( x 1 , τ 1 ) andÎ( x 2 , τ 2 ), i.e. coherence in the scalar case becomes g (2) ( r, τ 1 , τ 2 ) = 1 + 2 k 1 dk 1 j 0 (k 1 r) n k 1 (τ 1 ) k 2 dk 2 j 0 (k 2 r) n k 2 (τ 2 ) k 1 dk 1 n k 1 (τ 1 ) k 2 dk 2 n k 2 (τ 2 ) The large-scale limit the spherical Bessel functions go to 1 and therefore Eq. (5.23) becomes: The result of Eq. (5.15) holds also in the zero time-delay limit τ 1 − τ 2 = 0. This analysis demonstrates that the degree of second-order coherence for the squeezed relic photons does not go to 3 in the large-scale limit but rather to 5/3. It is interesting to stress, as we close, that the single-mode approximation is perfectly sound when the fluctuations beyond the horizon are described by a scalar field as it happens for the curvature perturbations [27]. In this case we could even go to higher order and compute the degrees of third-or fourth-order coherence (see Eqs. (3.11) and (3.12)) and confirm the same result. While the lengthy details will be omitted we can say that g (3) = 11 + O(1/n) and g (4) = 93 + O(1/n): this result holds also in the scalar case when all the modes of the field are taken into account. In the case of the squeezed relic photons, however, the role of the polarizations is essential, as the comparison of Eqs. (4.6), (4.7) and (4.8) clearly shows.

Concluding remarks
Among the six fiducial parameters characterizing the concordance scenario with massless neutrinos, a single number (i.e. the scalar spectral index) accounts for the presence of largescale inhomogeneities. A further source of inhomogeneity is represented by the tensor modes of the geometry even if their amplitude is, at least, one order of magnitude smaller than the one of the scalar modes. Furthermore since we do observe magnetic fields over large distance scales we may even admit the presence of large-scale gauge inhomogeneities. In the standard lore provided by conventional inflationary models all the potential sources of large-scale perturbations could stem from the zero-point fluctuations of quantum fields of different spins. At the moment the only argument in favour for this appealing possibility is merely theoretical: since a long stage of inflation is supposed to iron efficiently all preexisting inhomogeneities, it is logically plausible that large-scale fluctuations originated quantum mechanically. Because of the various assumptions behind this suggestion, it would be highly desirable to a have a more operational way of deciding about the statistical properties of large-scale fluctuations.
As we showed a possible answer to these questions involves the application of the tenets of Glauber theory, originally developed to address the coherence properties of optical fields. This analysis can be applied to the large-scale curvature perturbations but also to the largescale fluctuations of the gauge fields. Since the pioneering attempts of Hanbury Brown and Twiss, it has been realized that the study of first order interference between the amplitudes cannot be used to distinguish the nature of different quantum states of the radiation field. Young interferometry (indirectly based on the concept of power spectrum) 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. It is only by correlating intensities that the possible quantum origin of large-scale inhomogeneities can be independently assessed. In quantum optics the Glauber approach is often used in an exclusive manner by reducing the statistical properties of light to the analysis of a single (polarized) mode of the field: this is commonly referred to as the single-mode approximation. When dealing with large-scale fluctuations of different spins in cosmology the approach can only be inclusive since the correlation functions are typically unpolarized and contain all the modes of the field.
While the overall attempt of this paper is rather pragmatic, the obtained results are potentially inspiring. The modest viewpoint conveyed in this analysis is that precision cosmology, by itself, cannot validate its own premises. If new generations of astrophysical detectors will be able to resolve single photons the analysis of second-order interference effects may become feasible, at least in the case of the Cosmic Microwave Background.

A Basic conventions and notations
In time-dependent conformally flat backgrounds and in the Coulomb gauge (i.e. Y 0 = 0 and ∇ · Y = 0) the action (2.1) can be written as: where 12 A = Λ E /(4π) Y ; we have assumed that χ E and χ B are only dependent on the conformal time coordinate τ . In terms of the canonical momentum π conjugate to A the canonical Hamiltonian is simply given by: The discussion can be carried on in the case of different susceptibilities and different gauge couplings (see e.g. [33]); however we shall now focus on the case χ E = χ B = χ so that Eq. (A.2) becomes: The vector potential and the canonical momenta are explicitly related to the canonical electric and magnetic fields as B = ∇× A and as E = − π. In Fourier space the corresponding field operators are: iπ p α (τ ) e −i p· x ,π p, α = −i p 2 (â p α −â † − p α ), (A.5) where V is a fiducial (normalization) volume. In the discussion it is practical to switch from discrete to continuous modes where the creation and annihilation operators obey [â k α ,â † p β ] = δ αβ δ (3) ( k − p) and the sums are replaced by integrals according to k → V d 3 k/(2π) 3 . This observation should be borne in mind when discussing the explicit results; in terms of Eqs.

B Four-point functions
We report here some of the explicit expressions involved in the derivations of the fourpoint functions appearing in sections 4 and 5. Let us recall that, according to the present 12 The 1/ √ 4π is purely conventional and its presence comes from the factor 16π included in the initial gauge action. conventions: The two-point functions define the degree of first-order coherence and they are: Using Eq. (3.5) in the case n = 2 we have: The degrees of quantum coherence can also be defined in terms of the electric fields themselves, as originally suggested by Glauber. Equation (3.13) in the case n = 2 becomes: Finally, if we write Eq. (3.14) in the case n = 2 the result is: If the evolution of the susceptibility is not continuous (or not differentiable) we can still write a generic form of the u k (τ ) and v k (τ ), namely: where δ k , in this context, is just an arbitrary phase possibly picked up at the transition and, as usual, x i = kτ i . While in principle c ± (x i ) cannot be determined since the evolution is not continuous we can try to fix them by imposing, artificially, the continuity of the solutions for τ < −τ i and τ ≥ −τ i . The result of this procedure will be where, for simplicity, we denoted c ∓ (x i ) = c ∓ (x i )e ∓iδ k . The magnetic and the electric power spectra are, respectively, P B (k, τ ) = k 4 4π 2 |u k (τ ) − v * k (τ )| 2 , P E (k, τ ) = k 4 4π 2 |u k (τ ) + v * k (τ )| 2 (C.5) Using Eqs. (C.3) and (C.4), Eq. (C.5) becomes: (C. 6) In the sudden approximation (i.e. β → ∞ and C → 1) Eqs. (C.1) and (C.6) give the same result for x i 1 and kτ < 1. The reverse is not always true since the technique leading to Eq. (C.6) is based on the continuity of the susceptibility which is not verified in practice. The correct junction conditions for the susceptibility and for the extrinsic curvature are therefore essential for a correct derivation of the power spectra and of the degrees of quantum coherence.