Glauber theory and the quantum coherence of curvature inhomogeneities

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:

$$\begin{eqnarray}&&{{\mathcal{G}}^{(n,m)}}\left({{x}_{1}},...\,{{x}_{n}},\,{{x}_{n+1}},...,{{x}_{n+m}}\right)=\text{Tr}\left[\hat{\rho}\,{{\hat{q}}^{(-)}}\left({{x}_{1}}\right)...\,{{\hat{q}}^{(-)}}\left({{x}_{n}}\right)\,{{\hat{q}}^{(+)}}\left({{x}_{n+1}}\right)...\,{{\hat{q}}^{(+)}}\left({{x}_{n+m}}\right)\right],\end{eqnarray} \tag{ 2.2 }$$

where ${{x}_{i}}\equiv \left({{\vec{x}}_{i}},\,{{\tau}_{i}}\right)$ and $\hat{\rho}$ is the density operator representing the (generally mixed) state of the field $\hat{q}$. The field $\hat{q}\left(\vec{x},\tau \right)$ can always be expressed as $\hat{q}(x)={{\hat{q}}^{(+)}}(x)+{{\hat{q}}^{(-)}}(x)$, with ${{\hat{q}}^{(+)}}(x)={{\hat{q}}^{(-)\,\dagger}}(x)$. By definition we will have that ${{\hat{q}}^{(+)}}(x)|\text{vac}\rangle =0$ and also that $\langle \text{vac}|\,{{\hat{q}}^{(-)}}(x)=0$; the state $|\text{vac}\rangle$ denotes the vacuum⁷. 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 equation (2.2) and let us observe that an operator of the type

$$\begin{eqnarray}&&\hat{O}\left({{x}_{1}},...\,{{x}_{n}}\right)={{\hat{q}}^{(-)}}\left({{x}_{1}}\right)...\,{{\hat{q}}^{(-)}}\left({{x}_{n}}\right)\,{{\hat{q}}^{(+)}}\left({{x}_{1}}\right)...\,{{\hat{q}}^{(+)}}\left({{x}_{n}}\right),\end{eqnarray} \tag{ 2.3 }$$

is needed to describe n-fold delayed coincidence measurements of the field at the space-time points $\left({{x}_{1}},...\,{{x}_{n}}\right)$. If $|\,b\rangle$ is the state before the measurement and $|\,a\rangle$ is the state after the measurement, the matrix element corresponding to the absorption of the quanta of $\hat{q}$ at each detector and at given times is $\langle a\,|{{\hat{q}}^{(+)}}\left({{x}_{1}}\right)...\,{{\hat{q}}^{(+)}}\left({{x}_{n}}\right)|\,b\rangle$. The rate at which such absorptions occur, summed over the final states, is therefore proportional to:

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \underset{a}{\sum}\,|\langle a\,|{{{\hat{q}}}^{(+)}}\left({{x}_{1}}\right)...\,{{{\hat{q}}}^{(+)}}\left({{x}_{n}}\right)|\,b\rangle {{|}^{2}}=\underset{a}{\sum}\,\langle b|{{{\hat{q}}}^{(-)}}\left({{x}_{1}}\right)...\,{{{\hat{q}}}^{(-)}}\left({{x}_{n}}\right)|a\rangle \langle a|{{{\hat{q}}}^{(+)}}\left({{x}_{1}}\right)...\,{{{\hat{q}}}^{(+)}}\left({{x}_{n}}\right)|b\rangle =\langle b|\hat{O}|b\rangle, \end{array}\end{eqnarray} \tag{ 2.4 }$$

where $\hat{O}$ has been given in equation (2.3) and the second equality of equation (2.4) follows from the completeness relation. It is clear from equation (2.4) that when $\langle b|\hat{O}|b\rangle$ is averaged over the ensemble of the initial states of the system it becomes identical with equation (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

$$\begin{eqnarray}&&{{\mathcal{G}}^{(n)}}\left({{x}_{1}},...\,{{x}_{n}},\,{{x}_{n+1}},...,{{x}_{2n}}\right)=\text{Tr}\left[\hat{\rho}\,{{\hat{q}}^{(-)}}\left({{x}_{1}}\right)...\,{{\hat{q}}^{(-)}}\left({{x}_{n}}\right)\,{{\hat{q}}^{(+)}}\left({{x}_{n+1}}\right)...\,{{\hat{q}}^{(+)}}\left({{x}_{2n}}\right)\right].\end{eqnarray} \tag{ 2.5 }$$

Recalling the result of equation (2.5), the coherence properties of the quantum field $\hat{q}(x)$ can be discussed by introducing the normalized version of the n-point Glauber function [15]:

$$\begin{eqnarray}&&{{g}^{(n)}}\left({{x}_{1}},...\,{{x}_{n}},\,{{x}_{n+1}},...,{{x}_{2n}}\right)=\frac{{{\mathcal{G}}^{(n)}}\left({{x}_{1}},...\,{{x}_{n}},\,{{x}_{n+1}},...,{{x}_{2n}}\right)}{\sqrt{ \Pi _{j=1}^{2n}\,{{\mathcal{G}}^{(1)}}\left({{x}_{j}},{{x}_{j}}\right)}}.\end{eqnarray} \tag{ 2.6 }$$

While, by definition, $|{{g}^{(1)}}\left({{x}_{1}},\,{{x}_{2}}\right)|\leqslant 1$ the higher order correlators are not restricted in absolute value as it happens for ${{g}^{(1)}}\left({{x}_{1}},{{x}_{2}}\right)$. A fully coherent field must therefore satisfy the following necessary condition:

$$\begin{eqnarray}&&{{g}^{(n)}}\left({{x}_{1}},...\,{{x}_{n}},\,{{x}_{n+1}},...,{{x}_{2n}}\right)=1,\end{eqnarray} \tag{ 2.7 }$$

for $n=1,\,2,\,3,...$. If only a limited number of normalized correlation functions will satisfy equation (2.7) we shall speak about partial coherence. So for instance if ${{g}^{(1)}}\left({{x}_{1}},\,{{x}_{2}}\right)=1$ and ${{g}^{(2)}}\left({{x}_{1}},\,{{x}_{2}},\,{{x}_{3}},\,{{x}_{4}}\right)=1$ (but ${{g}^{(3)}}\left({{x}_{1}},\,{{x}_{2}},\,{{x}_{3}},\,{{x}_{4}},\,{{x}_{5}},\,{{x}_{6}}\right)\ne 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:

$$\begin{eqnarray}\begin{array}{*{35}{l}} {} & {} & {{g}^{(1)}}\left({{x}_{1}},{{x}_{2}}\right)=\frac{{{\mathcal{G}}^{(1)}}\left({{x}_{1}},{{x}_{2}}\right)}{\sqrt{{{\mathcal{G}}^{(1)}}\left({{x}_{1}},{{x}_{1}}\right)\,{{\mathcal{G}}^{(1)}}\left({{x}_{2}},{{x}_{2}}\right)}}, \end{array}\end{eqnarray} \tag{ 2.8 }$$

$$\begin{eqnarray}\begin{array}{*{35}{l}} {} & {} & {{g}^{(2)}}\left({{x}_{1}},{{x}_{2}},{{x}_{3}},{{x}_{4}}\right)=\frac{{{\mathcal{G}}^{(2)}}\left({{x}_{1}},{{x}_{2}},{{x}_{3}},{{x}_{4}}\right)}{\sqrt{{{\mathcal{G}}^{(1)}}\left({{x}_{1}},{{x}_{1}}\right)\,{{\mathcal{G}}^{(1)}}\left({{x}_{2}},{{x}_{2}}\right)\,{{\mathcal{G}}^{(1)}}\left({{x}_{3}},{{x}_{3}}\right)\,{{\mathcal{G}}^{(1)}}\left({{x}_{4}},{{x}_{4}}\right)}}. \end{array}\end{eqnarray} \tag{ 2.9 }$$

Similarly the third- and fourth-order degrees of coherence are:

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{g}^{(3)}}\left({{x}_{1}},\,{{x}_{2}},\,{{x}_{3}},\,{{x}_{4}},\,{{x}_{5}},\,{{x}_{6}}\right) & = & \frac{{{\mathcal{G}}^{(3)}}\left({{x}_{1}},\,{{x}_{2}},\,{{x}_{3}},\,{{x}_{4}},\,{{x}_{5}},\,{{x}_{6}}\right)}{\sqrt{{\prod}_{i=1}^{6}{{\mathcal{G}}^{(1)}}\left({{x}_{i}},\,{{x}_{i}}\right)}}, \end{array}\end{eqnarray} \tag{ 2.10 }$$

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{g}^{(4)}}\left({{x}_{1}},\,{{x}_{2}},\,{{x}_{3}},\,{{x}_{4}},\,{{x}_{5}},{{x}_{6}},\,{{x}_{7}},\,{{x}_{8}}\right) & = & \frac{{{\mathcal{G}}^{(4)}}\left({{x}_{1}},\,{{x}_{2}},\,{{x}_{3}},\,{{x}_{4}},\,{{x}_{5}},\,{{x}_{6}},\,{{x}_{7}},\,{{x}_{8}}\right)}{\sqrt{{\prod}_{i\,=\,1}^{8}{{\mathcal{G}}^{(1)}}\left({{x}_{i}},\,{{x}_{i}}\right)}}. \end{array}\end{eqnarray} \tag{ 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 field $\hat{q}$ (see section 3). Notice that the intensity correlators relevant to the HBT interferometry can be easily obtained from equation (2.6) by identifying the space-time points as follows:

$$\begin{eqnarray}&&{{x}_{1}}\equiv {{x}_{n+1}},\quad \quad {{x}_{2}}\equiv {{x}_{n+2}},\quad...\,\quad \quad {{x}_{n}}\equiv {{x}_{2n}}.\end{eqnarray} \tag{ 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 $\hat{\mathcal{I}}\left({{x}_{i}}\right)=\hat{q}\left({{x}_{i}}\right)\hat{q}\left({{x}_{i}}\right)$ where ${{x}_{i}}={{x}_{1}},\,{{x}_{2}},...{{x}_{n}}$.

Equations (2.10) and (2.11) describe, respectively, the degrees of third-order and fourth-order 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–31]. A specific discussion of these points can be found in section 4.

Since a fully coherent field must satisfy equation (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 equation (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.

The evolution of the field operators depends on two complex functions ${{v}_{k}}(\tau )$ and ${{u}_{k}}(\tau )$ whose explicit expression can be obtained by solving equations (C.1) and (C.2). In the case of a quasi-de Sitter evolution ${{v}_{k}}(\tau )$ and ${{u}_{k}}(\tau )$ are given by equations (C.3) and (C.4); in the pure de Sitter case their expression is given instead by equation (C.5). The first-order correlation function of equation (3.1) at separate space-time points is:

$$\begin{eqnarray}&&{{\mathcal{G}}^{(1)}}\left({{x}_{1}},\,{{x}_{2}}\right)=\frac{1}{2{{(2\pi )}^{3}}}{\int}^{}\frac{{{\text{d}}^{3}}{{k}_{1}}}{{{k}_{1}}}\,v_{{{k}_{1}}}^{\ast}\left({{\tau}_{1}}\right)\,{{v}_{{{k}_{1}}}}\left({{\tau}_{2}}\right)\,{{\text{e}}^{-\text{i}{{{\vec{k}}}_{1}}\centerdot \left({{{\vec{x}}}_{1}}-{{{\vec{x}}}_{2}}\right)}}.\end{eqnarray} \tag{ 3.4 }$$

Consequently, from equations (2.8) and (3.4) the normalized degree of first-order coherence is:

$$\begin{eqnarray}&&{{g}^{(1)}}\left({{\vec{x}}_{1}},\,{{\vec{x}}_{2}};\,{{\tau}_{1}},\,{{\tau}_{2}}\right)=\frac{{\int}^{}{{\text{d}}^{3}}{{k}_{1}}\,v_{{{k}_{1}}}^{\ast}\left({{\tau}_{1}}\right)\,{{v}_{{{k}_{1}}}}\left({{\tau}_{2}}\right)/{{k}_{1}}\,{{\text{e}}^{-\text{i}{{{\vec{k}}}_{1}}\centerdot \left({{{\vec{x}}}_{1}}-{{{\vec{x}}}_{2}}\right)}}}{\sqrt{{\int}^{}{{\text{d}}^{3}}{{k}_{1}}\,v_{{{k}_{1}}}^{\ast}\left({{\tau}_{1}}\right)\,{{v}_{{{k}_{1}}}}\left({{\tau}_{2}}\right)/{{k}_{1}}}\sqrt{{\int}^{}{{\text{d}}^{3}}{{k}_{2}}\,v_{{{k}_{2}}}^{\ast}\left({{\tau}_{2}}\right)\,{{v}_{{{k}_{2}}}}\left({{\tau}_{2}}\right)/{{k}_{2}}}}.\end{eqnarray} \tag{ 3.5 }$$

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 equations (C.3) and (C.4) imply:

$$\begin{eqnarray}&&{{v}_{k}}(\tau )=-\text{i}\,\frac{ \Gamma (\mu )}{\sqrt{2}\pi}\,{{\text{e}}^{-\text{i}\pi (\mu +1/2)/2}}\,{{\left(-\frac{k\tau}{2}\right)}^{1/2-\mu}}\left[1+\mathcal{O}(3-2\mu )\right]+\mathcal{O}\left(|k\tau {{|}^{1-\mu}}\right),\end{eqnarray} \tag{ 3.6 }$$

where μ is given, in the scalar and tensor case, by equation (C.4). Equation (3.6) has been obtained by expanding the Hankel functions of equation (C.3) in the limit of small arguments [38]. The first subleading term appearing in equation (3.6) is suppressed by $3-2\mu$ while the second correction is suppressed, in comparison with the leading term, by $|k\tau {{|}^{1/2}}$ which is small in the large-scale limit⁹. In the limit $\mu \to 3/2$ we have, from equation (3.6), that ${{v}_{k}}(\tau )=-\text{i}/(2k\tau )$ which coincides with the exact result obtainable by setting $\mu =3/2$ in equation (C.3). Performing the angular integrals we have that equation (3.6) implies

$$\begin{eqnarray}&&{{g}^{(1)}}\left(\vec{r};\,{{\tau}_{1}},\,{{\tau}_{2}}\right)=\frac{{\int}^{}{{k}_{1}}\text{d}{{k}_{1}}\,v_{{{k}_{1}}}^{\ast}\left({{\tau}_{1}}\right)\,{{v}_{{{k}_{1}}}}\left({{\tau}_{2}}\right)\,{{j}_{0}}\left({{k}_{1}}r\right)}{\sqrt{{\int}^{}{{k}_{1}}\text{d}{{k}_{1}}\,v_{{{k}_{1}}}^{\ast}\left({{\tau}_{1}}\right)\,{{v}_{{{k}_{1}}}}\left({{\tau}_{2}}\right)}\sqrt{{\int}^{}{{k}_{1}}\text{d}{{k}_{1}}\,v_{{{k}_{1}}}^{\ast}\left({{\tau}_{1}}\right)\,{{v}_{{{k}_{1}}}}\left({{\tau}_{2}}\right)}},\end{eqnarray} \tag{ 3.7 }$$

where ${{j}_{0}}\left({{k}_{1}}r\right)$ denotes the spherical Bessel function of zeroth order [38]. The numerator and the denominator of equation (3.7) depend on ${{\tau}_{1}}$ and ${{\tau}_{2}}$. However, owing to the specific form of ${{v}_{k}}(\tau )$, the dependence on ${{\tau}_{1}}$ and ${{\tau}_{2}}$ simplifies and the final form of equation (3.7) becomes

$$\begin{eqnarray}&&{{g}^{(1)}}(r)=\frac{{\int}^{}\text{d}{{k}_{1}}\,{{k}^{2-2\mu}}\,{{j}_{0}}\left(k{{r}_{1}}\right)}{{\int}^{}\text{d}{{k}_{1}}\,{{k}^{2-2\mu}}}\to 1,\end{eqnarray} \tag{ 3.8 }$$

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\ll 1$ (corresponding to large angular separations). The same result implied by equation (3.8) can be obtained when $k\tau \gg 1$ in equation (C.3). In this case the Hankel functions become plane waves and in the limit ${{k}_{1}}r\ll 1$ we have that ${{g}^{(1)}}\left(\vec{r},{{\tau}_{1}},{{\tau}_{2}}\right)$ always tend to 1 implying the first-order coherence of the underlying fluctuations.

The general form of the second-order correlation function involves four separated space-time points. As in the case of the first-order coherence equation (3.2) can be written in terms of ${{v}_{k}}(\tau )$ and ${{u}_{k}}(\tau )$:

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{\mathcal{G}}^{(2)}}({{x}_{1}},\,{{x}_{2}},\,{{x}_{3}},\,{{x}_{4}}) & = & {\int}^{}\frac{{{\text{d}}^{3}}{{k}_{1}}}{2{{k}_{1}}{{(2\pi )}^{3}}}{\int}^{}\frac{{{\text{d}}^{3}}{{k}_{2}}}{2{{k}_{2}}{{(2\pi )}^{3}}} \\ {} & \, & \times \,[v_{{{k}_{1}}}^{\ast}({{\tau}_{1}})\,v_{{{k}_{2}}}^{\ast}({{\tau}_{2}})\,{{v}_{{{k}_{1}}}}({{\tau}_{3}})\,{{v}_{{{k}_{2}}}}({{\tau}_{4}}){{\text{e}}^{-\text{i}{{{\vec{k}}}_{1}}\centerdot ({{{\vec{x}}}_{1}}-{{{\vec{x}}}_{3}})}}{{\text{e}}^{-\text{i}{{{\vec{k}}}_{2}}\centerdot ({{{\vec{x}}}_{2}}-{{{\vec{x}}}_{4}})}} \\ {} & \, & +\,v_{{{k}_{1}}}^{\ast}({{\tau}_{1}})\,v_{{{k}_{2}}}^{\ast}({{\tau}_{2}})\,{{v}_{{{k}_{1}}}}({{\tau}_{4}})\,{{v}_{{{k}_{2}}}}({{\tau}_{3}}){{\text{e}}^{-\text{i}{{{\vec{k}}}_{1}}\centerdot ({{{\vec{x}}}_{1}}-{{{\vec{x}}}_{4}})}}{{\text{e}}^{-\text{i}{{{\vec{k}}}_{2}}\centerdot ({{{\vec{x}}}_{2}}-{{{\vec{x}}}_{3}})}} \\ {} & \, & +\,v_{{{k}_{1}}}^{\ast}({{\tau}_{1}})\,u_{{{k}_{1}}}^{\ast}({{\tau}_{2}})\,{{u}_{{{k}_{2}}}}({{\tau}_{3}})\,{{v}_{{{k}_{2}}}}({{\tau}_{4}}){{\text{e}}^{-\text{i}{{{\vec{k}}}_{1}}\centerdot ({{{\vec{x}}}_{1}}-{{{\vec{x}}}_{2}})}}{{\text{e}}^{-\text{i}{{{\vec{k}}}_{2}}\centerdot ({{{\vec{x}}}_{3}}-{{{\vec{x}}}_{4}})}}]. \end{array}\end{eqnarray} \tag{ 3.9 }$$

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 equation (3.9) describes the interference of two beams with intensities $\hat{\mathcal{I}}\left({{\vec{x}}_{1}},{{\tau}_{1}}\right)$ and $\hat{\mathcal{I}}\left({{\vec{x}}_{2}},{{\tau}_{2}}\right)$, i.e.

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{\mathcal{G}}^{(2)}}\left({{x}_{1}},\,{{x}_{2}}\right) & = & \langle \hat{\mathcal{I}}\left({{{\vec{x}}}_{1}},{{\tau}_{1}}\right)\,\hat{\mathcal{I}}\left({{{\vec{x}}}_{2}},{{\tau}_{2}}\right)\rangle ={\int}^{}\frac{{{\text{d}}^{3}}{{k}_{1}}}{2{{k}_{1}}{{(2\pi )}^{3}}}{\int}^{}\frac{{{\text{d}}^{3}}{{k}_{2}}}{2{{k}_{2}}{{(2\pi )}^{3}}} \\ {} & {} & \times \left\{|{{v}_{{{k}_{1}}}}\left({{\tau}_{1}}\right){{|}^{2}}\,|{{v}_{{{k}_{2}}}}\left({{\tau}_{2}}\right){{|}^{2}}\left[1+{{\text{e}}^{-\text{i}\left({{{\vec{k}}}_{1}}-{{{\vec{k}}}_{2}}\right)\centerdot \vec{r}}}\right] \right. \\ {} & {} & \left. +\,v_{{{k}_{1}}}^{\ast}\left({{\tau}_{1}}\right)\,u_{{{k}_{1}}}^{\ast}\left({{\tau}_{2}}\right)\,{{u}_{{{k}_{2}}}}\left({{\tau}_{1}}\right)\,{{v}_{{{k}_{2}}}}\left({{\tau}_{2}}\right)\,{{\text{e}}^{-\text{i}\left({{{\vec{k}}}_{1}}+{{{\vec{k}}}_{2}}\right)\centerdot \vec{r}}}\right\}, \end{array}\end{eqnarray} \tag{ 3.10 }$$

where $\vec{r}={{\vec{x}}_{1}}-{{\vec{x}}_{2}}$. Inserting equations (3.9) and (3.10) into equation (2.9) the degree of second-order coherence becomes:

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{g}^{(2)}}\left(\vec{r},{{\tau}_{1}},{{\tau}_{2}}\right) & = & \frac{\langle \hat{\mathcal{I}}\left({{{\vec{x}}}_{1}},{{\tau}_{1}}\right)\,\hat{\mathcal{I}}\left({{{\vec{x}}}_{2}},{{\tau}_{2}}\right)\rangle}{\langle \hat{\mathcal{I}}\left({{{\vec{x}}}_{1}},{{\tau}_{1}}\right)\rangle \langle \hat{\mathcal{I}}\left({{{\vec{x}}}_{2}},{{\tau}_{2}}\right)\rangle} \\ {} & = & 1+\frac{{\int}^{}{{k}_{1}}\text{d}{{k}_{1}}|{{v}_{{{k}_{1}}}}\left({{\tau}_{1}}\right){{|}^{2}}\,{{j}_{0}}\left({{k}_{1}}r\right)\,\,{\int}^{}{{k}_{2}}\text{d}{{k}_{2}}|{{v}_{{{k}_{2}}}}\left({{\tau}_{2}}\right){{|}^{2}}\,{{j}_{0}}\left({{k}_{2}}r\right)}{{\int}^{}{{k}_{1}}\,\text{d}{{k}_{1}}|{{v}_{{{k}_{1}}}}\left({{\tau}_{1}}\right){{|}^{2}}\,{\int}^{}{{k}_{2}}\,\text{d}{{k}_{2}}|{{v}_{{{k}_{2}}}}\left({{\tau}_{2}}\right){{|}^{2}}} \\ {} & {} & +\,\frac{{\int}^{}{{k}_{1}}\text{d}{{k}_{1}}\,u_{{{k}_{1}}}^{\ast}\left({{\tau}_{2}}\right)v_{{{k}_{1}}}^{\ast}\left({{\tau}_{1}}\right)\,{{j}_{0}}\left({{k}_{1}}r\right)\,\,{\int}^{}{{k}_{2}}\text{d}{{k}_{2}}\,{{u}_{{{k}_{2}}}}\left({{\tau}_{1}}\right){{v}_{{{k}_{2}}}}\left({{\tau}_{2}}\right)\,{{j}_{0}}\left({{k}_{2}}r\right)}{{\int}^{}{{k}_{1}}\text{d}{{k}_{1}}|{{v}_{{{k}_{1}}}}\left({{\tau}_{1}}\right){{|}^{2}}\,\,{\int}^{}{{k}_{2}}\text{d}{{k}_{2}}|{{v}_{{{k}_{2}}}}\left({{\tau}_{2}}\right){{|}^{2}}}. \end{array}\end{eqnarray} \tag{ 3.11 }$$

According to equation (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 ${{\bar{n}}_{{{k}_{1}}}}\left({{\tau}_{1}}\right)=|{{v}_{{{k}_{1}}}}\left({{\tau}_{1}}\right){{|}^{2}}$ and ${{\bar{n}}_{{{k}_{2}}}}\left({{\tau}_{2}}\right)=|{{v}_{{{k}_{2}}}}\left({{\tau}_{2}}\right){{|}^{2}}$. Furthermore in the pure de Sitter case equation (C.4) imply:

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \frac{u_{{{k}_{1}}}^{\ast}\left({{\tau}_{1}}\right)\,v_{{{k}_{1}}}^{\ast}\left({{\tau}_{1}}\right){{u}_{{{k}_{2}}}}\left({{\tau}_{2}}\right)\,{{v}_{{{k}_{2}}}}\left({{\tau}_{2}}\right)}{|{{v}_{{{k}_{1}}}}\left({{\tau}_{1}}\right){{|}^{2}}\,|{{v}_{{{k}_{2}}}}\left({{\tau}_{2}}\right){{|}^{2}}}=1+4{{k}_{1}}{{\tau}_{1}}{{k}_{2}}{{\tau}_{2}}+2\text{i}\,\tau \left({{k}_{2}}{{\tau}_{2}}-{{k}_{1}}{{\tau}_{1}}\right). \end{array}\end{eqnarray} \tag{ 3.12 }$$

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\tau \ll 1$). Thanks to the preceding two observations we have that the degree of second-order coherence, to leading order¹⁰ in $1/{{\bar{n}}_{k}}$ is given by:

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{g}^{(2)}}\left(\vec{r},{{\tau}_{1}},{{\tau}_{2}}\right) & = & 1+2\frac{{\int}^{}{{k}_{1}}\text{d}{{k}_{1}}{{j}_{0}}\left({{k}_{1}}\,r\right)\,{{{\bar{n}}}_{{{k}_{1}}}}\left({{\tau}_{1}}\right)\,{\int}^{}{{k}_{2}}\text{d}{{k}_{2}}{{j}_{0}}\left({{k}_{2}}\,r\right)\,{{{\bar{n}}}_{{{k}_{2}}}}\left({{\tau}_{2}}\right)}{{\int}^{}{{k}_{1}}\text{d}{{k}_{1}}{{{\bar{n}}}_{{{k}_{1}}}}\left({{\tau}_{1}}\right)\,{\int}^{}{{k}_{2}}\text{d}{{k}_{2}}\,{{{\bar{n}}}_{{{k}_{2}}}}\left({{\tau}_{2}}\right)} \\ {} & {} & +\,\frac{{\int}^{}{{k}_{1}}\text{d}{{k}_{1}}{{j}_{0}}\left({{k}_{1}}\,r\right)/\sqrt{{{{\bar{n}}}_{{{k}_{1}}}}\left({{\tau}_{1}}\right)}\,\,{\int}^{}{{k}_{2}}\text{d}{{k}_{2}}{{j}_{0}}\left({{k}_{2}}\,r\right)/\sqrt{{{{\bar{n}}}_{{{k}_{2}}}}\left({{\tau}_{2}}\right)}\,}{{\int}^{}{{k}_{1}}\text{d}{{k}_{1}}{{{\bar{n}}}_{{{k}_{1}}}}\left({{\tau}_{1}}\right){\int}^{}{{k}_{2}}\text{d}{{k}_{2}}\,{{{\bar{n}}}_{{{k}_{2}}}}\left({{\tau}_{2}}\right)}. \end{array}\end{eqnarray} \tag{ 3.13 }$$

The large-scale limit the spherical Bessel functions go to 1 and therefore equation (3.13) becomes¹¹:

$$\begin{eqnarray}&&{{g}^{(2)}}\left(r,{{\tau}_{1}},{{\tau}_{2}}\right)\to 3,\quad \quad \underset{{{\tau}_{1}}\to {{\tau}_{2}}}{{\lim}}\,{{g}^{(2)}}\left(r,{{\tau}_{1}},{{\tau}_{2}}\right)={{g}^{(2)}}(r,\tau ).\end{eqnarray} \tag{ 3.14 }$$

The result of equation (3.14) holds provided the total number of efolds N_tot is larger than the maximal number of efolds today accessible by observations¹² 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\left(\vec{r},\tau \right)={{g}^{(2)}}\left(\vec{r},\tau \right)-1$.

All in all the large-scale curvature fluctuations are only partially coherent. They are first-order 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 equation (2.7) and discussion thereafter).

Indeed the result of equation (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 field $\hat{q}\left(\vec{x},\tau \right)$ is in a thermal state in standard Minkowski space-time. In this case the density matrix can be written, in the Fock basis, as:

$$\begin{eqnarray}&&\hat{\rho}=\underset{\{n\}}{\sum}\,\,{{P}_{\{n\}}}\,|\{n\}\rangle \langle \{n\}|,\quad \quad \underset{\{n\}}{\sum}\,\,{{P}_{\{n\}}}=1.\end{eqnarray} \tag{ 4.1 }$$

The multimode probability distribution appearing in equation (4.1) is given by:

$$\begin{eqnarray}&&{{P}_{\{n\}}}=\underset{{\vec{k}}}{\prod}\,\frac{\bar{n}_{k}^{{{n}_{{\vec{k}}}}}}{{{\left(1+{{{\bar{n}}}_{k}}\right)}^{{{n}_{{\vec{k}}}}+1}}},\end{eqnarray} \tag{ 4.2 }$$

where ${{\bar{n}}_{k}}=\text{Tr}\left[\hat{\rho}\,\hat{d}_{{\vec{k}}}^{\dagger}\,{{\hat{d}}_{{\vec{k}}}}\right]$ is the average occupation number of each Fourier mode. Furthermore, following the standard notation, $|\{n\}\rangle =|{{n}_{{{{\vec{k}}}_{1}}}}\rangle \,|{{n}_{{{{\vec{k}}}_{2}}}}\rangle \,|{{n}_{{{{\vec{k}}}_{3}}}}\rangle...$ 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. ${{\bar{n}}_{k}}$) does not need to coincide with the Bose–Einstein occupation number¹³. Using equation (2.9) and performing the appropriate averages by means of the density matrix of equation (4.1) the normalized degree of second-order coherence becomes, in this case,

$$\begin{eqnarray}&&{{g}^{(2)}}\left(\vec{r},{{\tau}_{1}},{{\tau}_{2}}\right)=\frac{{\int}^{}{{\text{d}}^{3}}{{k}_{1}}{{{\bar{n}}}_{{{k}_{1}}}}\left({{\tau}_{1}}\right)/{{k}_{1}}\,{\int}^{}{{\text{d}}^{3}}{{k}_{2}}\,{{{\bar{n}}}_{{{k}_{2}}}}/{{k}_{2}}\,\left[1+{{\text{e}}^{-\text{i}\left({{{\vec{k}}}_{1}}+{{{\vec{k}}}_{2}}\right)\centerdot \vec{r}}}\right]}{{\int}^{}{{\text{d}}^{3}}{{k}_{1}}{{{\bar{n}}}_{{{k}_{1}}}}\left({{\tau}_{1}}\right)/{{k}_{1}}\,{\int}^{}{{\text{d}}^{3}}{{k}_{2}}\,{{{\bar{n}}}_{{{k}_{2}}}}\left({{\tau}_{2}}\right)/{{k}_{2}}}.\end{eqnarray} \tag{ 4.3 }$$

While the calculation leading to equation (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¹⁴

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \langle \hat{d}_{i}^{\dagger}\,\hat{d}_{j}^{\dagger}\,{{{\hat{d}}}_{k}}\,{{{\hat{d}}}_{\ell}}\rangle =\langle \hat{d}_{i}^{\dagger}\,\hat{d}_{i}^{\dagger}\,{{{\hat{d}}}_{i}}\,{{{\hat{d}}}_{i}}\rangle {{\delta}_{i\,j}}\,{{\delta}_{j\,k}}\,{{\delta}_{\ell \,k}}+\langle \hat{d}_{i}^{\dagger}\,\hat{d}_{j}^{\dagger}\,{{{\hat{d}}}_{i}}\,{{{\hat{d}}}_{j}}\rangle \,{{\delta}_{i\,k}}\,{{\delta}_{j\,\ell}}\left[1-{{\delta}_{ij}}\right]+\langle \hat{d}_{i}^{\dagger}\,\hat{d}_{j}^{\dagger}\,{{{\hat{d}}}_{j}}\,{{{\hat{d}}}_{i}}\rangle \,{{\delta}_{i\,\ell}}\,{{\delta}_{j\,k}}\left[1-{{\delta}_{ij}}\right], \end{array}\end{eqnarray} \tag{ 4.4 }$$

where ${{\hat{d}}_{i}}$ and $\hat{d}_{j}^{\dagger}$ denote the annihilation and creation operators related two generic momenta, i.e. for instance ${{\hat{d}}_{{\vec{q}}}}$ and $\hat{d}_{{\vec{p}}}^{\dagger}$; furthermore, following the same shorthand notation, ${{\delta}_{i\,j}}$ denotes the delta functions over the three-momenta (i.e. ${{\delta}_{\vec{q},\,\vec{p}}}$). Performing each of the averages of equation (4.4) in terms of the density matrix (4.1) leads to equation (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¹⁵:

$$\begin{eqnarray}&&{{g}^{(2)}}\left(\vec{r},\tau \right)=1+\frac{{\int}^{}{{k}_{1}}\text{d}{{k}_{1}}{{{\bar{n}}}_{{{k}_{1}}}}(\tau )\,{{j}_{0}}\left({{k}_{1}}r\right){\int}^{}{{k}_{2}}\text{d}{{k}_{2}}\,{{{\bar{n}}}_{{{k}_{2}}}}(\tau ){{j}_{0}}\left({{k}_{2}}r\right)}{{\int}^{}{{k}_{1}}\text{d}{{k}_{1}}{{{\bar{n}}}_{{{k}_{1}}}}(\tau )\,{\int}^{}{{k}_{2}}\text{d}{{k}_{2}}\,{{{\bar{n}}}_{{{k}_{2}}}}(\tau )}.\end{eqnarray} \tag{ 4.5 }$$

In the large-scale limit equation (4.5) will then imply that ${{g}^{(2)}}\left(\vec{r},\tau \right)\to 2$ [13]; conversely equation (3.14), in the same physical limit, implies ${{g}^{(2)}}\left(\vec{r},\tau \right)\to 3$.

The curvature quanta follow a Bose–Einstein multiplicity distribution¹⁶. 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 $\vec{k}$-mode of the field. We are interested in computing $\langle m\,m\,|\zeta \rangle$ where $|\zeta \rangle = \Sigma (\zeta )|0\,0\rangle$ and $\Sigma (\zeta )$ is actually given in equation (C.9). We will have that the operator $\Sigma (\zeta )$ can be factorized as [41]:

$$\begin{eqnarray} &&\Sigma (\zeta )=\exp \left[-\frac{\zeta}{|\zeta |}\tanh |\zeta |\,{{K}_{+}}\right]\times \exp \left[-2\ln \cosh |\zeta |\,{{K}_{0}}\right]\times \exp \left[\frac{{{\zeta}^{\ast}}}{|\zeta |}\,\tanh |\zeta |{{K}_{-}}\right],\end{eqnarray} \tag{ 4.6 }$$

where K₊ , K₋ and K₀ obey the standard commutation relations of the SU(1, 1) Lie algebra¹⁷. Using the previous equation we will have that the multiplicity distribution is given by:

$$\begin{eqnarray}&&|\langle m\,m\,|\zeta \rangle {{|}^{2}}=\frac{1}{\bar{n}+1}{{\left(\frac{{\bar{n}}}{\bar{n}+1}\right)}^{m}},\quad \quad \bar{n}={{\sinh}^{2}}r.\end{eqnarray} \tag{ 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 becomes

$$\begin{eqnarray}&&{{\hat{\rho}}_{{\vec{k}}}}=\frac{1}{{{\cosh}^{2}}{{r}_{k}}}\underset{{{n}_{{\vec{k}}}}=0}{\overset{\infty}{\sum}}\,\underset{{{m}_{{\vec{k}}}}=0}{\overset{\infty}{\sum}}\,{{\text{e}}^{-\text{i}{{\alpha}_{k}}\left({{n}_{{\vec{k}}}}-{{m}_{{\vec{k}}}}\right)}}{{\left(\tanh {{r}_{k}}\right)}^{{{n}_{{\vec{k}}}}+{{m}_{{\vec{k}}}}}}|{{n}_{{\vec{k}}}}\,\,{{n}_{-\vec{k}}}\rangle \langle {{m}_{-\vec{k}}}\,\,{{m}_{{\vec{k}}}}|,\end{eqnarray} \tag{ 4.8 }$$

whose diagonal elements define the multiplicity distribution which has precisely the Bose–Einstein form:

$$\begin{eqnarray}&&{{P}_{\left\{{{n}_{{\vec{k}}}}\right\}}}=\underset{{\vec{k}}}{\prod}\,{{P}_{{{n}_{{\vec{k}}}}}},\quad \quad {{P}_{{{n}_{{\vec{k}}}}}}\left({{\bar{n}}_{k}}\right)=\frac{\bar{n}_{k}^{{{n}_{{\vec{k}}}}}}{{{\left(1+{{{\bar{n}}}_{k}}\right)}^{{{n}_{{\vec{k}}}}+1}}},\end{eqnarray} \tag{ 4.9 }$$

accounting for the way curvature quanta of each Fourier mode (i.e. ${{n}_{{\vec{k}}}}$) are distributed as a function of their average multiplicity (i.e. ${{\bar{n}}_{k}}={{\sinh}^{2}}{{r}_{k}}$).

The behaviour of the off-diagonal elements of equation (4.8) is dictated by the phases ${{\alpha}_{k}}=\left(2{{\varphi}_{k}}-{{\gamma}_{k}}\right)$ whose evolution is specifically discussed in appendix C (see, in particular, equation (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¹⁸.

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¹⁹.

Indeed, for a single mode of the field the degrees of first- and second-order coherence are defined as:

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{{\bar{g}}}^{(1)}}\left({{\tau}_{1}},{{\tau}_{2}}\right) & = & \frac{\langle {{{\hat{a}}}^{\dagger}}\left({{\tau}_{1}}\right)\,\hat{a}\left({{\tau}_{2}}\right)\rangle}{\sqrt{\langle {{{\hat{a}}}^{\dagger}}\left({{\tau}_{1}}\right)\,\hat{a}\left({{\tau}_{1}}\right)\rangle}\,\sqrt{\langle {{{\hat{a}}}^{\dagger}}\left({{\tau}_{2}}\right)\,\hat{a}\left({{\tau}_{2}}\right)\rangle}}, \end{array}\end{eqnarray} \tag{ 4.10 }$$

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{{\bar{g}}}^{(2)}}\left({{\tau}_{1}},{{\tau}_{2}}\right) & = & \frac{\langle {{{\hat{a}}}^{\dagger}}\left({{\tau}_{1}}\right){{{\hat{a}}}^{\dagger}}\left({{\tau}_{2}}\right)\,\hat{a}\left({{\tau}_{2}}\right)\,\hat{a}\left({{\tau}_{1}}\right)\rangle}{\langle {{{\hat{a}}}^{\dagger}}\left({{\tau}_{1}}\right)\,\hat{a}\left({{\tau}_{1}}\right)\rangle \langle {{a}^{\dagger}}\left({{\tau}_{2}}\right)\,a\left({{\tau}_{2}}\right)\rangle}, \end{array}\end{eqnarray} \tag{ 4.11 }$$

where the overline at the left hand side distinguishes equations (4.10) and (4.11) from equations (2.7) and (2.8)–(2.9) holding in the general case. Equations (4.10) and (4.11) define, respectively, the degrees of first- and second-order temporal coherence: in the zero time-delay limit ${{\tau}_{1}}-{{\tau}_{2}}\to 0$ and in this case the degree of second-order coherence will be denoted by ${{\bar{g}}^{(2)}}$. Recalling the analysis of section 4 we have already demonstrated, in practice, that:

$$\begin{eqnarray}&&\underset{kr\ll 1}{{\lim}}\,{{g}^{(2)}}(r,\tau )={{\bar{g}}^{(2)}}.\end{eqnarray} \tag{ 4.12 }$$

Equation (4.12) is verified in the case of a single-mode coherent state (i.e. $\hat{a}|\alpha \rangle =\alpha |\alpha \rangle$) where equations (4.10) and (4.11) imply ${{\bar{g}}^{(1)}}={{\bar{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 ${{\bar{g}}^{(1)}}=1$ but ${{\bar{g}}^{(2)}}=2$. In the case of a single Fock state ${{\bar{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 $\hat{a}=\cosh r\hat{b}-\sinh r{{\hat{b}}^{\dagger}}$ where the phases have been fixed to zero; taking the limit of zero time-delay and inserting these expressions in equation (4.11) we have that: ${{\bar{g}}^{(2)}}=3+1/\bar{n}$ with $\bar{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/\bar{n}$ (and higher) are all negligible since $\bar{n}$. In this $1/\bar{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. ${{\bar{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. ${{\bar{g}}^{(2)}}<1$) implying a degree of second-order coherence smaller than in the case of a coherent state. The correspondence of equation (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 equation (4.11) implies ${{\bar{g}}^{(2)}}=\left({{D}^{2}}-\bar{n}\right)/{{\bar{n}}^{2}}$ where $\bar{n}=\langle \hat{N}\rangle$ and ${{D}^{2}}=\langle {{\hat{N}}^{2}}\rangle -{{\langle \hat{N}\rangle}^{2}}$. It is sometimes useful to define the so-called Mandel parameter whose expression is $\mathcal{Q}=\bar{n}\left[{{\bar{g}}^{(2)}}-1\right]={{D}^{2}}/\bar{n}-1$. The parameter $\mathcal{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}}=\langle \hat{N}\rangle$ while $\mathcal{Q}=0$; this means that the distribution underlying this state is just the Poisson distribution (i.e. the variance coincides with the mean value).

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.

$$\begin{eqnarray}&&\hat{\rho}=\underset{n=0}{\overset{\infty}{\sum}}\,{{P}_{n}}\,|n\,\rangle \langle n|,\quad \quad \underset{n=0}{\overset{\infty}{\sum}}\,\,{{P}_{n}}=1.\end{eqnarray} \tag{ 4.13 }$$

We work in the zero-delay limit and exploit the results of equations (4.10)–(4.12). The degrees of correlations can then be expressed as²⁰:

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{{\bar{g}}}^{(1)}} & = & \frac{\langle {{{\hat{a}}}^{\dagger}}\,\hat{a}\rangle}{{\bar{n}}}=\frac{1}{{\bar{n}}}\frac{\text{d}\mathcal{A}}{\text{d}s}{{|}_{s=1}}, \end{array}\end{eqnarray} \tag{ 4.14 }$$

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{{\bar{g}}}^{(2)}} & = & \frac{\langle {{{\hat{a}}}^{\dagger}}\,{{{\hat{a}}}^{\dagger}}\hat{a}\,\hat{a}\rangle}{{{{\bar{n}}}^{2}}}=\frac{1}{{{{\bar{n}}}^{2}}}\frac{{{\text{d}}^{2}}\mathcal{A}}{\text{d}{{s}^{2}}}{{|}_{s=1}}, \end{array}\end{eqnarray} \tag{ 4.15 }$$

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{{\bar{g}}}^{(3)}} & = & \frac{\langle {{{\hat{a}}}^{\dagger}}\,{{{\hat{a}}}^{\dagger}}\,{{{\hat{a}}}^{\dagger}}\hat{a}\,\hat{a}\,\hat{a}\rangle}{{{{\bar{n}}}^{3}}}=\frac{1}{{{{\bar{n}}}^{3}}}\frac{{{\text{d}}^{3}}\mathcal{A}}{\text{d}{{s}^{3}}}{{|}_{s=1}}, \end{array}\end{eqnarray} \tag{ 4.16 }$$

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{{\bar{g}}}^{(4)}} & = & \frac{\langle {{{\hat{a}}}^{\dagger}}\,{{{\hat{a}}}^{\dagger}}\,{{{\hat{a}}}^{\dagger}}\,{{{\hat{a}}}^{\dagger}}\hat{a}\,\hat{a}\,\hat{a}\,\hat{a}\rangle}{{{{\bar{n}}}^{4}}}=\frac{1}{{{{\bar{n}}}^{4}}}\frac{{{\text{d}}^{4}}\mathcal{A}}{\text{d}{{s}^{4}}}{{|}_{s=1}}, \end{array}\end{eqnarray} \tag{ 4.17 }$$

where $\bar{n}=\text{Tr}\left[{{\hat{a}}^{\dagger}}\,\hat{a}\,\hat{\rho}\right]$; the function $\mathcal{A}(s)$ appearing in equations (4.14)–(4.17) is the probability generating function defined as:

$$\begin{eqnarray}&&\mathcal{A}(s)=\underset{n=0}{\overset{\infty}{\sum}}\,\,{{s}^{n}}\,{{P}_{n}},\quad \quad \mathcal{A}(1)=1.\end{eqnarray} \tag{ 4.18 }$$

One of the simplest situations coming to mind is the one where P_n is a Poisson distribution. In this case we have that

$$\begin{eqnarray}&&{{P}_{n}}=\frac{{{\text{e}}^{-\bar{n}}}}{n!}\,{{\bar{n}}^{n}},\quad \quad \mathcal{A}(s)={{\text{e}}^{\bar{n}(s-1)}}\end{eqnarray} \tag{ 4.19 }$$

and thanks to equations (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 equation (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 equations (4.14)–(4.18) is when P_n is a Bose–Einstein distribution. In this case ${{\bar{g}}^{(n)}}=n!$ and, as already discussed, ${{\bar{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}_{\mathcal{N}}}$ of a number $\mathcal{N}$ of mutually independent random variables, i.e.

$$\begin{eqnarray}&&{{A}_{\mathcal{N}}}={{X}_{1}}+{{X}_{2}}+....+\,{{X}_{\mathcal{N}}}.\end{eqnarray} \tag{ 4.20 }$$

If $\mathcal{N}$ is fixed the probability generating function $\mathcal{P}(s)$ of the sum ${{A}_{\mathcal{N}}}$ is simply given by the product of the generating functions. If $\mathcal{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 $\mathcal{P}(s)$. If the distribution of $\mathcal{N}$ is G_k its corresponding generating function will be, as usual $\mathcal{G}(s)={{{\sum}^{}}_{k}}{{s}^{k}}{{G}_{k}}$. The generating function of describing the compound process of equation (4.20) will then be given by $\mathcal{A}(s)=\mathcal{G}\left[\mathcal{P}(s)\right]$ [43].

If $\mathcal{N}$ is distributed as a Poissonian variable with multiplicity ${{\bar{N}}_{\text{c}}}$ the probability generating function will be given by $\exp \left\{{{\bar{N}}_{c}}\left[\mathcal{P}(s)-1\right]\right\}$. If we take $\mathcal{P}(s)$ to be the generating function of a logarithmic distribution, i.e.

$$\begin{eqnarray}&&\mathcal{P}(s)=1-\frac{k}{{{{\bar{N}}}_{c}}}\ln \left[\frac{{\bar{n}}}{k}(1-s)+1\right],\end{eqnarray} \tag{ 4.21 }$$

The generating function of the random sum of random variables will then be given by:

$$\begin{eqnarray}&&\mathcal{A}(s)={{\text{e}}^{{{{\bar{N}}}_{c}}\left[\mathcal{P}(s)-1\right]}}=\frac{{{k}^{k}}}{{{\left[(1-s)\bar{n}+k\right]}^{k}}},\quad \quad {{\bar{g}}^{(n)}}=\frac{k(k+1)...\,(k+n-1)}{{{k}^{n}}}.\end{eqnarray} \tag{ 4.22 }$$

In the case k  =  1 equation (4.22) reproduces the probability generating function of a Bose–Einstein distribution while, in the limit $k\to \infty$ $\mathcal{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/\bar{n}$, namely from equation (4.22)

$$\begin{eqnarray}&&{{\bar{g}}^{(1)}}=1,\quad \quad {{\bar{g}}^{(2)}}=1+\frac{1}{k}\to 3,\quad \quad k=1/2.\end{eqnarray} \tag{ 4.23 }$$

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 equation (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 equation (4.22) would be slightly different and the value of k bound to be integer.

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–46] for CMB observation at THz ($1~\text{THz}={{10}^{3}}~\text{GHz}$) frequencies. Furthermore the idea suggested in [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 [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 $\mathcal{R}$ is related to the field operator $\hat{q}$ (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 ${{ \Delta }_{\text{I}}}\left(\vec{x},{{\tau}_{\text{dec}}}\right)\simeq -\mathcal{R}\left(\vec{n},{{\tau}_{\text{dec}}}\right)/5$ under the approximation of sudden decoupling (i.e. assuming that the visibility function is a narrow Gaussian centered at ${{\tau}_{\text{dec}}}$) [47, 48]. This correspondence will now be swiftly outlined.

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 ${{ \Delta }_{I}}\left(\hat{n},\tau \right)$ (the brightness perturbation of the intensity). The polarization anisotropies are instead discussed in terms of the E and B-modes, denoted, respectively, by ${{ \Delta }_{E}}\left(\hat{n},\tau \right)$ and ${{ \Delta }_{B}}\left(\hat{n},\tau \right)$. Both ${{ \Delta }_{E}}$ and ${{ \Delta }_{B}}$ are spin 0 quantities, exactly as ${{ \Delta }_{I}}$ [49]. Recalling that $\mu =\cos \vartheta$ we have, in the concordance framework and in the absence of gravitons,

$$\begin{eqnarray}&&{{ \Delta }_{\text{E}}}\left(\hat{n},\tau \right)=-\partial _{\mu}^{2}\left\{\left(1-{{\mu}^{2}}\right){{ \Delta }_{\text{P}}}\left(\hat{n},\tau \right)\right\},\quad \quad {{ \Delta }_{\text{B}}}\left(\hat{n},\tau \right)=0.\end{eqnarray} \tag{ 6.1 }$$

The brightness perturbations ${{ \Delta }_{I}}\left(\hat{n},\tau \right)$ and ${{ \Delta }_{P}}\left(\hat{n},\tau \right)$ obey, in Fourier space, the following pair of equations:

$$\begin{eqnarray}\begin{array}{*{35}{l}} {} & {} & \Delta _{I}^{\prime}+\left(\text{i}k\mu +{{\varepsilon}^{\prime}}\right){{ \Delta }_{I}}={{\psi}^{\prime}}-\text{i}k\mu \phi +{{\varepsilon}^{\prime}}\left[{{ \Delta }_{\text{I}0}}+\mu {{v}_{\text{b}}}+\frac{\left(1-3{{\mu}^{2}}\right)}{4}{{S}_{\text{P}}}(k,\tau )\right], \end{array}\end{eqnarray} \tag{ 6.2 }$$

$$\begin{eqnarray}\begin{array}{*{35}{l}} {} & {} & \Delta _{P}^{\prime}+\left(\text{i}k\mu +{{\varepsilon}^{\prime}}\right){{ \Delta }_{P}}=\frac{3{{\varepsilon}^{\prime}}}{4}\left(1-{{\mu}^{2}}\right){{S}_{\text{P}}}(k,\tau ), \end{array}\end{eqnarray} \tag{ 6.3 }$$

where ${{S}_{\text{P}}}(k,\tau )$ 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}_{\text{P}}}(k,\tau )=\left({{ \Delta }_{\text{I}2}}+{{ \Delta }_{\text{P}0}}+{{ \Delta }_{\text{P}2}}\right)$; note that, in equations (6.2) and (6.3), ${{\varepsilon}^{\prime}}$ and $\varepsilon \left(\tau,{{\tau}_{0}}\right)$ denote, respectively, the differential optical depth and the optical depth itself [47, 48]

$$\begin{eqnarray}&&{{\varepsilon}^{\prime}}={{x}_{\text{e}}}{{\tilde{n}}_{\text{e}}}\,a\,{{\sigma}_{\gamma \text{e}}},\quad \quad \varepsilon \left(\tau,{{\tau}_{0}}\right)={\int}_{{{\tau}_{0}}}^{\tau}{{x}_{\text{e}}}{{\tilde{n}}_{\text{e}}}\,a\,{{\sigma}_{\gamma \text{e}}}\text{d}\tau.\end{eqnarray} \tag{ 6.4 }$$

The line of sight solution of equations (6.2) and (6.3) can be written, respectively, as:

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{ \Delta }_{\text{I}}}\left(k,\mu,{{\tau}_{0}}\right) & = & {\int}_{0}^{{{\tau}_{0}}}\mathcal{K}(\tau )\left[{{ \Delta }_{\text{I}0}}+\phi +\mu {{v}_{\text{b}}}+\frac{\left(1-3{{\mu}^{2}}\right)}{4}{{S}_{\text{P}}}\right]{{\text{e}}^{-\text{i}\mu x(\tau )}} \\ {} & \, & +\,{\int}_{0}^{{{\tau}_{0}}}\text{d}\tau {{\text{e}}^{-\varepsilon \left(\tau,{{\tau}_{0}}\right)}}\left({{\phi}^{\prime}}+{{\psi}^{\prime}}\right){{\text{e}}^{-\text{i}\mu x(\tau )}}\text{d}\tau, \end{array}\end{eqnarray} \tag{ 6.5 }$$

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{ \Delta }_{\text{P}}}\left(k,\mu,{{\tau}_{0}}\right) & = & \frac{3}{4}\left(1-{{\mu}^{2}}\right){\int}_{0}^{{{\tau}_{0}}}\mathcal{K}(\tau ){{S}_{\text{P}}}(k,\tau ){{\text{e}}^{-\text{i}k\mu \left(\tau -{{\tau}_{0}}\right)}}\text{d}\tau, \end{array}\end{eqnarray} \tag{ 6.6 }$$

where $\mathcal{K}(\tau )={{\varepsilon}^{\prime}}{{\text{e}}^{-\varepsilon \left(\tau,{{\tau}_{0}}\right)}}$ is the visibility function and $x(\tau )=k\left({{\tau}_{0}}-\tau \right)$. 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 $\ell \leqslant \sqrt{{{z}_{\text{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 $\mathcal{K}(\tau )$ is proportional to a Dirac delta function and ${{\text{e}}^{-\varepsilon \left(\tau,{{\tau}_{0}}\right)}}$ is proportional to an Heaviside theta function. Under the latter approximations, equation (6.2) leads to the well known pair of separated contributions, i.e. the Sachs-Wolfe (SW) and the integrated Sachs-Wolfe (ISW) contributions:

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{ \Delta }_{\text{I}}}\left(k,\mu,{{\tau}_{0}}\right) & = & \Delta _{\text{I}}^{\left(\text{SW}\right)}\left(k,\mu,{{\tau}_{0}}\right)+ \Delta _{\text{I}}^{\left(\text{ISW}\right)}\left(k,\mu,{{\tau}_{0}}\right), \end{array}\end{eqnarray} \tag{ 6.7 }$$

$$\begin{eqnarray}\begin{array}{*{35}{l}} \Delta _{\text{I}}^{\left(\text{SW}\right)}\left(k,\mu,{{\tau}_{0}}\right) & = & {{\left(-\frac{\mathcal{R}(k,\tau )}{5}\right)}_{{{\tau}_{\text{rec}}}}}{{\text{e}}^{-\text{i}\mu {{y}_{\text{rec}}}}}, \\ \Delta _{\text{I}}^{\left(\text{ISW}\right)}\left(k,\mu,{{\tau}_{0}}\right) & = & {\int}_{{{\tau}_{\text{rec}}}}^{{{\tau}_{0}}}\left({{\phi}^{\prime}}+{{\psi}^{\prime}}\right){{\text{e}}^{-\text{i}\mu x(\tau )}}\,\text{d}\tau, \end{array}\end{eqnarray} \tag{ 6.8 }$$

where, by definition, $x\left({{\tau}_{\text{rec}}}\right)={{y}_{\text{rec}}}$. The SW and the ISW contributions of equations (6.7) and (6.8) can be separately evaluated. In particular, as anticipated, the ordinary SW contribution becomes:

$$\begin{eqnarray}\begin{array}{*{35}{l}} \Delta _{\text{I}}^{\left(\text{SW}\right)}\left(k,\mu,{{\tau}_{0}}\right) & = & -\frac{\mathcal{R}\left(\vec{k},{{\tau}_{\text{i}}}\right)}{5}\mathcal{S}(q){{\text{e}}^{-\text{i}\mu {{y}_{\text{rec}}}}}, \end{array}\end{eqnarray} \tag{ 6.9 }$$

$$\begin{eqnarray}\begin{array}{*{35}{l}} \mathcal{S}(q) & = & 1+\frac{4}{3q}-\frac{16}{3{{q}^{2}}}+\frac{16\left(\sqrt{q+1}-1\right)}{3{{q}^{3}}}, \end{array}\end{eqnarray} \tag{ 6.10 }$$

while the ISW contribution is:

$$\begin{eqnarray}\begin{array}{*{35}{l}} {} & {} & \Delta _{\text{I}}^{\left(\text{ISW}\right)}\left(k,\mu,{{\tau}_{0}}\right)=-2\mathcal{R}\left(\vec{k},{{\tau}_{\text{i}}}\right){\int}_{{{\tau}_{\text{rec}}}}^{{{\tau}_{0}}}\mathcal{T}_{\mathcal{R}}^{\,\prime}(\tau ){{\text{e}}^{-\text{i}\mu x(\tau )}}\,\text{d}\tau, \end{array}\end{eqnarray} \tag{ 6.11 }$$

$$\begin{eqnarray}\begin{array}{*{35}{l}} {} & {} & {{\mathcal{T}}_{\mathcal{R}}}(\tau )=1-\frac{\mathcal{H}(\tau )}{{{a}^{2}}(\tau )}{\int}_{0}^{\tau}{{a}^{2}}\left({{\tau}^{\prime}}\right)\,\text{d}{{\tau}^{\prime}}. \end{array}\end{eqnarray} \tag{ 6.12 }$$

Both in equations (6.9) and (6.11), $\mathcal{R}\left(\vec{k},\tau \right)$ denotes the constant value of curvature perturbations at ${{\tau}_{\text{i}}}<{{\tau}_{\text{eq}}}$. By further approximating the integrand in equation (6.11) the whole large-scale contribution can be written, for the present purposes, as

$$\begin{eqnarray}&&{{ \Delta }_{\text{I}}}\left(k,{{\tau}_{0}}\right)=-\mathcal{R}\left(\vec{k},{{\tau}_{\text{i}}}\right){{\left\{\frac{\mathcal{S}(q)}{5}+\left[\mathcal{T}_{\mathcal{R}}^{\,\prime}\right]\right\}}_{{{q}_{\text{rec}}}}}\,{{\text{e}}^{-\text{i}\mu {{y}_{\text{rec}}}}},\end{eqnarray} \tag{ 6.13 }$$

where ${{q}_{\text{rec}}}={{a}_{\text{rec}}}/{{a}_{\text{eq}}}\simeq 3.04$ for $h_{0}^{2}{{ \Omega }_{\text{M}0}}\simeq 0.134$ (as usual ${{ \Omega }_{\text{M}0}}$ denotes the total matter density in critical units at the present time). Note that if we neglect the ISW term and assume that $q\gg 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 equations (2.9)–(2.11) can therefore be written as:

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{g}^{(2)}}\left(\hat{m},\hat{n}\right) & = & \frac{\langle {{{\rm{\hat \Delta }}}^{(-)}}\left(\hat{m}\right)\,{{{\rm{\hat \Delta }}}^{(-)}}\left(\hat{n}\right)\,{{ \Delta }^{(+)}}\left(\hat{m}\right)\,{{ \Delta }^{(+)}}\left(\hat{n}\right)\rangle}{\langle {{{\rm{\hat \Delta }}}^{(-)}}\left(\hat{m}\right)\,{{{\rm{\hat \Delta }}}^{(+)}}\left(\hat{m}\right)\rangle \langle {{{\rm{\hat \Delta }}}^{(-)}}\left(\hat{n}\right)\,{{{\rm{\hat \Delta }}}^{(+)}}\left(\hat{n}\right)\rangle}, \end{array}\end{eqnarray} \tag{ 6.14 }$$

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{g}^{(3)}}\left(\hat{m},\hat{n},\hat{r}\right) & = & \frac{\langle {{{\rm{\hat \Delta }}}^{(-)}}\left(\hat{m}\right)\,{{{\rm{\hat \Delta }}}^{(-)}}\left(\hat{n}\right)\,{{{\rm{\hat \Delta }}}^{(-)}}\left(\hat{r}\right){{ \Delta }^{(+)}}\left(\hat{m}\right)\,{{ \Delta }^{(+)}}\left(\hat{n}\right)\rangle}{\langle {{{\rm{\hat \Delta }}}^{(-)}}\left(\hat{m}\right)\,{{{\rm{\hat \Delta }}}^{(+)}}\left(\hat{m}\right)\rangle \langle {{{\rm{\hat \Delta }}}^{(-)}}\left(\hat{n}\right)\,{{{\rm{\hat \Delta }}}^{(+)}}\left(\hat{n}\right)\rangle \langle {{{\rm{\hat \Delta }}}^{(-)}}\left(\hat{n}\right){{{\rm{\hat \Delta }}}^{(+)}}\left(\hat{n}\right)\rangle}, \end{array}\end{eqnarray} \tag{ 6.15 }$$

and similarly for ${{g}^{(4)}}\left(\hat{m},\hat{n},\hat{r},\hat{s}\right)$. In equations (6.14) and (6.15) we decomposed the positive and negative frequency parts of the operators as $\rm{\hat \Delta }\left(\hat{n}\right)={{\rm{\hat \Delta }}^{(+)}}\left(\hat{n}\right)+{{\rm{\hat \Delta }}^{(-)}}\left(\hat{n}\right)$; note finally that $\rm{\hat \Delta }\left(\hat{n}\right)$ can denote, indifferently, ${{\rm{\hat \Delta }}_{I}}$ of ${{\rm{\hat \Delta }}_{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 $\rm{\hat \Delta }$ with ${{\rm{\hat \Delta }}_{I}}$ in equation (6.14) we shall have that

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{g}^{(2)}}\left(\hat{m},\hat{n}\right) & = & 1+\frac{{{{\sum}^{}}_{\ell,\,{{\ell}^{\prime}}}}{\int}^{}q|{{v}_{q}}{{|}^{2}}\,\text{d}q{\int}^{}p|{{v}_{p}}{{|}^{2}}\,\text{d}p{{A}_{\ell}}\left(q,\hat{m}\centerdot \hat{n}\right){{A}_{{{\ell}^{\prime}}}}\left(\,p,\hat{m}\centerdot \hat{n}\right)}{{{{\sum}^{}}_{\ell}}{\int}^{}q\text{d}q|{{v}_{q}}{{|}^{2}}{{A}_{\ell}}\left(q,\hat{m}\centerdot \hat{n}\right)\,\,{{{\sum}^{}}_{{{\ell}^{\prime}}}}{\int}^{}p\text{d}p|{{v}_{p}}{{|}^{2}}{{A}_{{{\ell}^{\prime}}}}\left(q,\hat{m}\centerdot \hat{n}\right)} \\ {} & \, & +\,\frac{{{{\sum}^{}}_{\ell,\,{{\ell}^{\prime}}}}{\int}^{}qv_{q}^{\ast}\,u_{q}^{\ast}\text{d}q{\int}^{}p{{v}_{p}}\,{{u}_{p}}\,\text{d}p{{A}_{\ell}}\left(q,\hat{m}\centerdot \hat{n}\right){{A}_{{{\ell}^{\prime}}}}\left(\,p,\hat{m}\centerdot \hat{n}\right)}{{{{\sum}^{}}_{\ell}}{\int}^{}q\text{d}q|{{v}_{q}}{{|}^{2}}{{A}_{\ell}}\left(q,\hat{m}\centerdot \hat{n}\right)\,\,{{{\sum}^{}}_{{{\ell}^{\prime}}}}{\int}^{}p\text{d}p|{{v}_{p}}{{|}^{2}}{{A}_{\ell}}\left(q,\hat{m}\centerdot \hat{n}\right)}, \\ {{A}_{\ell}}\left(q,\hat{m}\centerdot \hat{n}\right) & = & (2\ell +1)j_{\ell}^{2}\left(k{{\tau}_{0}}\right){{P}_{\ell}}\left(\hat{m}\centerdot \hat{n}\right), \end{array}\end{eqnarray} \tag{ 6.16 }$$

where ${{j}_{\ell}}\left(k{{\tau}_{0}}\right)$ are the spherical Bessel functions and ${{P}_{\ell}}\left(\hat{m}\centerdot \hat{n}\right)$ 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 equation (6.16): in the limit $\hat{m}\centerdot \hat{n}\to 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].

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 section.

In the hypothesis that curvature quanta are the sole source of temperature inhomogeneities the distribution P(n) of the total number of phonons $n={{{\sum}^{}}_{{\vec{k}}}}{{n}_{{\vec{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}_{\left\{{{n}_{k}}\right\}}}$ which has been discussed in equation (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]

$$\begin{eqnarray}&&p\left(\{n\}\right)=\underset{{\vec{k}}}{\prod}\,\frac{1}{\left(1+{{{\bar{n}}}_{{\vec{k}}}}\right){{\left(1+1/{{{\bar{n}}}_{{\vec{k}}}}\right)}^{{{n}_{{\vec{k}}}}}}}.\end{eqnarray} \tag{ 6.17 }$$

For any mode for which ${{\bar{n}}_{k}}=0$, the corresponding factor must be interpreted as ${{\delta}_{{{n}_{{\vec{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)={{{\sum}^{}}_{\{n\}}}p\left(\{n\}\right){{\delta}_{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 ${{\bar{n}}_{{\vec{k}}}}$ of all the ε occupied modes become equal²¹; in this case the joint probability distribution of the occupied modes becomes:

$$\begin{eqnarray}&&p\left(\{n\}\right)=\frac{1}{{{\left(1+\bar{n}/\epsilon \right)}^{\epsilon}}{{\left(1+\epsilon /\bar{n}\right)}^{n}}},\quad \quad \bar{n}=\underset{{\vec{k}}}{\sum}\,{{\bar{n}}_{k}}=\epsilon \,{{\bar{n}}_{k}}.\end{eqnarray} \tag{ 6.18 }$$

Every non-vanishing term in the summation $P(m)={{{\sum}^{}}_{\{n\}}}p\left(\{n\}\right){{\delta}_{m\,n}}$ has the same value and the required probability P(m) is simply p({n}) given by equation (6.18) multiplied by a well known combinatorial factor accounting for the way the n phonons are distributed among the modes:

$$\begin{eqnarray}&&{{P}_{n}}\left(\bar{n},\epsilon \right)=\frac{ \Gamma (n+\epsilon )}{ \Gamma (\epsilon ) \Gamma (n+1)}{{\left(\frac{{\bar{n}}}{\bar{n}+\epsilon}\right)}^{n}}{{\left(\frac{\epsilon}{\bar{n}+\epsilon}\right)}^{\epsilon}}.\end{eqnarray} \tag{ 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 equation (6.19) should be analyzed is given by $n\gg 1$ and $\bar{n}\gg 1$. More specifically, the realizations of the distribution (6.19) with different $\bar{n}$. In the limit $\bar{n}\gg \epsilon$ the negative binomial distribution becomes a Gamma distribution with a characteristic scaling law:

$$\begin{eqnarray}&&\bar{n}{{P}_{n}}\left(\bar{n},\epsilon \right)\simeq {{\psi}_{\epsilon}}(w),\quad \quad {{\psi}_{\epsilon}}(w)=\frac{{{\epsilon}^{\epsilon}}}{ \Gamma (\epsilon )}{{w}^{\epsilon -1}}{{\text{e}}^{-\epsilon w}},\quad \quad w=n/\bar{n}.\end{eqnarray} \tag{ 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}}\left(\bar{n}\right)$ as $n\to \infty$, $\bar{n}\to \infty$ while $n/\bar{n}$ is kept fixed. For fixed ε the distribution of equation (6.19) scales as ${{\psi}_{\epsilon}}(w)$.