The Influence of Zero-Point Fluctuations on the Photon Wave Packet Motion in Vacuum

1. Introduction

Current investigations of photon propagation in vacuum probe fundamental physics at the intersection of quantum field theory, cosmology, and spacetime symmetry tests [1]. Two phenomena remain pivotal: (i) energy-dependent velocity dispersion [2,3], potentially detectable as time delay for photons of the different energies from distant astrophysical sources, and (ii) intrinsic decoherence mechanisms linked to quantum spacetime fluctuations [4,5,6], which may depolarize radiation over cosmological scales [7,8]. While velocity dispersion tests Lorentz invariance [1], decoherence bridges quantum gravity phenomenology and observational astrophysics [5,9].

Here, we develop a framework to quantify these effects via photon dynamics in fluctuating spacetime. Formulating an optical Dirac equation for photons in curved geometry allows us to use the density matrix formalism in the first-order differential equation. The next step is to derive an equation for the photon density matrix in a random space-time, which is assumed to be Minkowski on average. In this way, one needs to calculate the correlators of a metric arising due to zero-point fluctuations of the quantum fields. Then, we introduce some quantity to describe an electromagnetic field’s degree of coherence (i.e., purity of a system). Finally, an explicit formula for the decoherence time of a photon wave packet is deduced.

2. Optical Dirac Equation for the Photon Under Gravitational Background

The wave equation inherently involves a second-order time derivative, distinguishing it from the Schrödinger equation, which governs the dynamics of massive particles and is first-order in time. Decoherence effects for such particles have been analyzed within the Schrödinger framework [4]. Extending this approach to electromagnetic fields requires redefining Maxwell’s equations in a matrix form. This reformulation allows for constructing an optical analogue of the Dirac equation, as discussed in [10], providing a suitable foundation for investigating decoherence in the photon sector.

Because it is not possible to define a vacuum state which is invariant relative to general coordinate transformations, it is reasonable to define a unique class of metrics to describe the gravitational field arising due to zero-point fluctuations [4,11,12,13,14]:

$$d{s}^2\equiv g_{\mu \nu }d{x}^{\mu }d{x}^{\nu }=a^2{\left(1-{\partial }_m{P}^m\right)}^{2}d{\eta }^2-{\gamma }_{ij}(d{x}^i+N^{i}d\eta )(d{x}^j+N^{j}d\eta ),$$

(2.1)

where $x^{\mu }=\{\eta ,\mathit{x}\}$, $\eta$ is a conformal time, ${\gamma }_{ij}$ is a spatial metric, $a={\gamma }^{1/6}$ is a locally defined scale factor, and $\gamma =det{\gamma }_{ij}$. The spatial part of the interval (2.1) reads as

$$d{l}^2\equiv {\gamma }_{ij}d{x}^{i}d{x}^j=a^2(\eta ,\mathit{x}){\tilde{\gamma }}_{ij}d{x}^{i}d{x}^j,$$

(2.2)

where ${\tilde{\gamma }}_{ij}={\gamma }_{ij}/a^2$ is a matrix with the unit determinant. Let us set $\mathit{N}=0$ and $\mathit{P}$ and write Maxwell equations in this metric. According to [15], the Maxwell equations in the three-dimensional form can be written as

$$\begin{array}{c}\hfill {\displaystyle \frac{1}{\sqrt{\gamma }}}{\displaystyle \frac{\partial }{\partial x^i}}\left(\sqrt{\gamma }{D}^i\right)=0,\end{array}$$

(2.3)

$$\begin{array}{c}\hfill {\displaystyle \frac{1}{2\sqrt{\gamma }}}{e}^{ijk}\left({\displaystyle \frac{\partial H_k}{\partial x^j}}-{\displaystyle \frac{\partial H_j}{\partial x^k}}\right)={\displaystyle \frac{1}{\sqrt{\gamma }}}{\displaystyle \frac{\partial }{\partial \eta }}\left(\sqrt{\gamma }{D}^i\right),\end{array}$$

(2.4)

$$\begin{array}{c}\hfill {\displaystyle \frac{1}{\sqrt{\gamma }}}{\displaystyle \frac{\partial }{\partial x^i}}\left(\sqrt{\gamma }{B}^i\right)=0,\end{array}$$

(2.5)

$$\begin{array}{c}\hfill {\displaystyle \frac{1}{2\sqrt{\gamma }}}{e}^{ijk}\left({\displaystyle \frac{\partial E_k}{\partial x^j}}-{\displaystyle \frac{\partial E_j}{\partial x^k}}\right)=-{\displaystyle \frac{1}{\sqrt{\gamma }}}{\displaystyle \frac{\partial }{\partial \eta }}\left(\sqrt{\gamma }{B}^i\right),\end{array}$$

(2.6)

where $B^i={\gamma }^{ij}{H}_j/\sqrt{g_{00}}$, $D^i={\gamma }^{ij}{E}_j/\sqrt{g_{00}}$. In the conformally-unimodular metric (2.1) $\sqrt{g_{00}}=a$, $\sqrt{\gamma }=a^3$, ${\gamma }^{ij}=a^{-2}{\tilde{\gamma }}^{ij}$. By introducing new quantities $\mathcal{D}=a^3\mathit{D}$, and $\mathcal{B}=a^3\mathit{B}$, the system of Eqs. (2.3),(2.4),(2.5),(2.6) acquires the form:

$$\begin{array}{c}\hfill div\mathcal{D}=0,\end{array}$$

(2.7)

$$\begin{array}{c}\hfill rot\left(\tilde{\mathit{\gamma }}\mathcal{B}\right)={\displaystyle \frac{\partial }{\partial \eta }}\mathcal{D},\end{array}$$

(2.8)

$$\begin{array}{c}\hfill div\mathcal{B}=0,\end{array}$$

(2.9)

$$\begin{array}{c}\hfill rot\left(\tilde{\mathit{\gamma }}\mathcal{D}\right)=-{\displaystyle \frac{\partial }{\partial \eta }}\mathcal{B},\end{array}$$

(2.10)

where a matrix $\tilde{\mathit{\gamma }}$ denotes the matrix ${\tilde{\gamma }}_{ij}$ with the unit determinant. These equations could be put into the form of an optical Dirac equation [10]

$${\displaystyle \frac{\partial }{\partial \eta }}\left(\begin{array}{c}\mathcal{D}\\ i\mathcal{B}\end{array}\right)=\left(\begin{array}{cc}0& \mathit{S}\widehat{\mathit{p}}\\ \mathit{S}\widehat{\mathit{p}}& 0\end{array}\right)\left(\begin{array}{c}\tilde{\mathit{\gamma }}\mathcal{D}\\ i\tilde{\mathit{\gamma }}\mathcal{B}\end{array}\right),$$

(2.11)

where the rotor is expressed through the matrix of spin 1

$$S_x=\left(\begin{array}{ccc}0& 0& 0\\ 0& 0& -i\\ 0& i& 0\end{array}\right),S_y=\left(\begin{array}{ccc}0& 0& i\\ 0& 0& 0\\ -i& 0& 0\end{array}\right),S_z=\left(\begin{array}{ccc}0& -i& 0\\ i& 0& 0\\ 0& 0& 0\end{array}\right)$$

(2.12)

by the relation $\left(\mathit{S}\widehat{\mathit{p}}\right)\mathit{F}=i\widehat{\mathit{p}}\times \mathit{F}$ for the operator $\widehat{\mathit{p}}=-i\nabla$ and an arbitrary vector $\mathit{F}$. Denoting

$$\mathit{\alpha }=\left(\begin{array}{cc}0& \mathit{S}\\ \mathit{S}& 0\end{array}\right),\mathcal{G}=\left(\begin{array}{cc}\tilde{\mathit{\gamma }}& 0\\ 0& \tilde{\mathit{\gamma }}\end{array}\right)$$

(2.13)

allows us to write the optical Dirac equation for a six-component wave function

$$i{\displaystyle \frac{\partial \mathsf{\Psi }}{\partial \eta }}=\mathit{\alpha }\widehat{\mathit{p}}\mathcal{G}\mathsf{\Psi }.$$

(2.14)

Additionally, the conditions (2.7), (2.9) on the upper and lower components of the photon wave function have to be satisfied. For a static gravitational field, the scalar product

$$<{\mathsf{\Psi }}_1|{\mathsf{\Psi }}_2>\int {\mathsf{\Psi }}_1^{+}\mathcal{G}{\mathsf{\Psi }}_2{d}^3\mathit{r}$$

(2.15)

is conserved during evolution, but for the nonstationary case, does not. We are not interested in the photon creation by nonstationary gravitational field, so let us ad hoc modify Eq. (2.14)

$$i{\displaystyle \frac{\partial \mathsf{\Psi }}{\partial \eta }}=\mathit{\alpha }\widehat{\mathit{p}}\mathcal{G}\mathsf{\Psi }-i{\displaystyle \frac{{\mathcal{G}}^{-1}}{2}}{\displaystyle \frac{\partial \mathcal{G}}{\partial \eta }}\mathsf{\Psi }.$$

(2.16)

For the equation (2.16) the scalar product (2.15) is conserved even for a time-dependent gravitational field. Let us do a non-unitary transformation

$$\psi (\mathit{r},\eta )={\mathcal{G}}^{1/2}(\mathit{r},\eta )\mathsf{\Psi }(\mathit{r},\eta ),$$

(2.17)

which leads to the photon quantum mechanics with the "flat" scalar product $<{\psi }_1|{\psi }_2>=\int {\psi }_1^{+}{\psi }_2{d}^3\mathit{r}$. As a result, the optical Dirac equation acquires the form

$$i{\displaystyle \frac{\partial \psi }{\partial \eta }}={\mathcal{G}}^{1/2}\mathit{\alpha }\mathit{p}{\mathcal{G}}^{1/2}\psi ,$$

(2.18)

for an arbitrary time-dependent gravitational field $\mathcal{G}$. Considering the gravitational field as a perturbation reduces (2.18) to the approximate equation

$$i{\displaystyle \frac{\partial \psi }{\partial \eta }}=\mathit{\alpha }\widehat{\mathit{p}}\psi +{\displaystyle \frac{1}{2}}\left((\mathcal{G}-I)\mathit{\alpha }\widehat{\mathit{p}}+\alpha \widehat{\mathit{p}}(\mathcal{G}-I)\right)\psi ,$$

(2.19)

where the first term on the right-hand side of (2.19) contains the Hamiltonian $H_0\equiv \mathit{\alpha }\widehat{\mathit{p}}$ and the remaining terms represent the interaction

$$V={\displaystyle \frac{1}{2}}\left((\mathcal{G}-I)\mathit{\alpha }\widehat{\mathit{p}}+\mathit{\alpha }\widehat{\mathit{p}}(\mathcal{G}-I)\right).$$

(2.20)

Further, we will consider a fluctuating gravitational field originating due to the zero-point fluctuation of quantum fields.

3. Correlators of the Fluctuating Gravitational Field

Let us consider empty space-time filled only by vacuum, but taking into account its quantum properties, i.e., the gravitational field created by zero-point fluctuations of quantum fields. For simplicity, only the scalar field will be considered. Fluctuations of the gravitational field in an arbitrary metric could be significant. However, we assume that fluctuations of the gravitational field in a class of conformally-unimodular metrics (2.1) are relatively small, allowing us to consider them as perturbations. Scalar perturbations of the conformally-unimodular metric [14,16] is written as

$$d{s}^2=a{(\eta ,\mathit{x})}^2\left(d{\eta }^2-\left(\left(1+{\displaystyle \frac{1}{3}}\sum _{m=1}^3{\partial }_m^{2}F(\eta ,\mathit{x})\right){\delta }_{ij}-{\partial }_i{\partial }_{j}F(\eta ,\mathit{x})\right)d{x}^{i}d{x}^j\right),$$

(3.1)

where the perturbations of the locally defined scale factor is

$$a(\eta ,\mathit{x})=e^{\alpha \left(\eta \right)}(1+\mathsf{\Phi }(\eta ,\mathit{x})).$$

(3.2)

We will refer to $\mathsf{\Phi }$ as a gravitational potential. Stress-energy momentum tensor in the hydrodynamic approximation [17]

$$T_{\mu \nu }=(p+\rho )u_{\mu }{u}_{\nu }-p{g}_{\mu \nu }$$

(3.3)

includes the perturbations of the energy density $\rho (\eta ,\mathit{x})={\rho }_v+\delta \rho (\eta ,\mathit{x})$ and pressure $p(\eta ,\mathit{x})=p_v+\delta p(\eta ,\mathit{x})$ around the vacuum mean values, where the index v will denote a uniform component of the vacuum energy density and pressure.

The zero-order equations for a flat universe take the form [18,19,20]

$$\begin{array}{c}\hfill M_p^{-2}{e}^{4\alpha }{\rho }_v-{\displaystyle \frac{1}{2}}{e}^{2\alpha }{\alpha }^{\prime 2}=const,\end{array}$$

(3.4)

$$\begin{array}{c}\hfill {\alpha }^{\prime \prime }+{\alpha }^{\prime 2}=M_p^{-2}{e}^{2\alpha }({\rho }_v-3p_v),\end{array}$$

(3.5)

where $\alpha \left(\eta \right)=loga\left(\eta \right)$ and the reduced Planck mass $M_p=\sqrt{{\displaystyle \frac{3}{4\pi G}}}$ is implied. According to five-vectors theory of gravity [11], the first Friedmann equation (3.4) is satisfied up to some constant, and the main parts of the vacuum energy density and pressure

$$\begin{array}{c}\hfill {\rho }_v\approx N_{all}{\displaystyle \frac{k_{max}^4}{16{\pi }^2{a}^4}},\end{array}$$

(3.6)

$$\begin{array}{c}\hfill p_v={\displaystyle \frac{1}{3}}{\rho }_v\end{array}$$

(3.7)

do not contribute to the universe expansion. In the formula (3.6), the number $N_{all}$ of all degrees of freedom of the quantum fields in nature appears because the zero-point stress-energy tensor is an additive quantity [21]. The momentum ultraviolet cut-off [13,19,20]

$$k_{max}\approx {\displaystyle \frac{12M_p}{\sqrt{2+N_{sc}}}}$$

(3.8)

is proportional to the Planck mass and includes the number of minimally coupled scalar fields $N_{sc}$ plus two, because the gravitational waves possess two additional degrees of freedom [19].

Without including a real matter, and if the constant in Eq. (3.4) compensates a vacuum energy (3.6) exactly, one comes to the static Minkowski space-time. Further, we will consider the perturbations under this background and set $\alpha \left(\eta \right)=0$ in (3.2).

Generally, a vacuum can be considered as some fluid, i.e., “ether” [20], but with some stochastic properties among the elastic ones. Let us return to the stress-energy tensor (3.3) and introduce other variables

$$\begin{array}{c}\hfill \wp (\eta ,\mathit{x})=a^4(\eta ,\mathit{x})\rho (\eta ,\mathit{x}),\end{array}$$

(3.9)

$$\begin{array}{c}\hfill \mathsf{\Pi }(\eta ,\mathit{x})=a^4(\eta ,\mathit{x})p(\eta ,\mathit{x})\end{array}$$

(3.10)

for the reasons which will be explained below. The perturbations around the uniform values can be written now as $\wp (\eta ,\mathit{x})={\rho }_v+\delta \wp (\eta ,\mathit{x})$, $\mathsf{\Pi }(\eta ,\mathit{x})=p_v+\delta \mathsf{\Pi }(\eta ,\mathit{x})$. The vacuum-ether 4-velocity u is represented in the form of

$$u^{\mu }=\{(1-\mathsf{\Phi }(\eta ,\mathit{x})),\nabla {\displaystyle \frac{v(\eta ,\mathit{x})}{\wp (\eta ,\mathit{x})+\mathsf{\Pi }(\eta ,\mathit{x})}}\}\approx \{(1-\mathsf{\Phi }(\eta ,\mathit{x})),\nabla {\displaystyle \frac{v(\eta ,\mathit{x})}{{\rho }_v+p_v}}\},$$

(3.11)

where $v(\eta ,\mathit{x})$ is a scalar function. Expanding all perturbations into the Fourier series $\delta \wp (\eta ,\mathit{x})={\sum }_{\mathit{k}}\delta {\wp }_{\mathit{k}}\left(\eta \right)e^{i\mathit{k}\mathit{x}},\mathsf{\Phi }(\eta ,\mathit{x})={\sum }_{\mathit{k}}{\mathsf{\Phi }}_{\mathit{k}}\left(\eta \right)e^{i\mathit{k}\mathit{x}}...$ etc. results in the equations for the perturbations:

$$\begin{array}{c}\hfill -6{\widehat{\mathsf{\Phi }}}_{\mathit{k}}^{\prime }+k^2{\widehat{F}}_{\mathit{k}}^{\prime }+{\displaystyle \frac{18}{M_p^2}}{\widehat{v}}_{\mathit{k}}=0,\end{array}$$

(3.12)

$$\begin{array}{c}\hfill -6k^2{\widehat{\mathsf{\Phi }}}_{\mathit{k}}+k^4{\widehat{F}}_{\mathit{k}}+{\displaystyle \frac{18}{M_p^2}}\delta {\widehat{\wp }}_{\mathit{k}}=0,\end{array}$$

(3.13)

$$\begin{array}{c}\hfill -12{\widehat{\mathsf{\Phi }}}_{\mathit{k}}-3{\widehat{F}}_{\mathit{k}}^{\prime \prime }+k^2{\widehat{F}}_{\mathit{k}}=0,\end{array}$$

(3.14)

$$\begin{array}{c}\hfill -9{\widehat{\mathsf{\Phi }}}_{\mathit{k}}^{\prime \prime }-9k^2{\widehat{\mathsf{\Phi }}}_{\mathit{k}}+k^4{\widehat{F}}_{\mathit{k}}-{\displaystyle \frac{9}{M_p^2}}\left(3\delta {\widehat{\mathsf{\Pi }}}_{\mathit{k}}-\delta {\widehat{\wp }}_{\mathit{k}}\right)=0,\end{array}$$

(3.15)

$$\begin{array}{c}\hfill -\delta {\widehat{\wp }}_{\mathit{k}}^{\prime }+k^2{\widehat{v}}_{\mathit{k}}=0,\end{array}$$

(3.16)

$$\begin{array}{c}\hfill \delta {\widehat{\mathsf{\Pi }}}_{\mathit{k}}+{\widehat{v}}_{\mathit{k}}^{\prime }=0.\end{array}$$

(3.17)

It is remarkable that the choice of the variables (3.9), (3.10), (3.11) means that the values ${\rho }_v$ and $p_v$ do not appear in the system (3.12)-(3.17). The second point is that the continuity and Newton’s second law equations (3.16), (3.17) do not contain metric perturbation.

From now on, we will begin to consider the perturbation in Eqs. (3.12)-(3.17) as operators by writing a “hat" under every quantity. Here, we do not suppose the strong nonlinearity [22] and assume a smallness of the quantum fluctuations of space-time in this particular conformally-unimodular metric. Let us emphasize that the system (3.12)-(3.17) for a perturbation evolution is exact in the first order on perturbations. However, it is not closed. To obtain a closed system, one needs, for instance, to specify the equation of state for a perturbation of pressure. Still, as an approximation, we could calculate pressure and energy density strictly by using the quantum field theory under the unperturbed Minkowski space-time. Expressing $F_{\mathit{k}}$ from Eq. (3.13) and substituting it into Eq. (3.15) leads to

$${\widehat{\mathsf{\Phi }}}_{\mathit{k}}^{\prime \prime }+{\displaystyle \frac{1}{3}}{k}^2{\widehat{\mathsf{\Phi }}}_{\mathit{k}}+{\displaystyle \frac{1}{M_p^2}}\left(3\delta {\widehat{\mathsf{\Pi }}}_{\mathit{k}}+\delta {\widehat{\wp }}_{\mathit{k}}\right)=0.$$

(3.18)

Below, we will approximately consider an operator $3\delta {\widehat{\mathsf{\Pi }}}_{\mathit{k}}+\delta {\widehat{\wp }}_{\mathit{k}}$ by using the creation and annihilation operators under the Minkowski space-time background. Such an approximation allows closing the system (3.12)-(3.17). It is important that, on the whole, the general relativity, permitting an arbitrary curved space-time background, does not allow a well-defined and invariant vacuum state [23].

3.1. Quantum Fields as a Source for Energy Density and Pressure Perturbation

Let us consider a single massless scalar field as an example of a quantum field. Energy density and pressure of the scalar field in the pure Minkowski space-time (without metric perturbation) has the form [21]

$$\begin{array}{c}\hfill \widehat{p}(\eta ,\mathit{x})={\displaystyle \frac{{\widehat{\phi }}^{\prime 2}}{2}}-{\displaystyle \frac{{(\nabla \widehat{\phi })}^2}{6}},\end{array}$$

(2.39)

$$\begin{array}{c}\hfill \widehat{\rho }(\eta ,\mathit{x})={\displaystyle \frac{{\widehat{\phi }}^{\prime 2}}{2}}+{\displaystyle \frac{{(\nabla \widehat{\phi })}^2}{2}}\end{array}$$

(2.40)

All the quantities may be expanded into the Fourier series $\widehat{\phi }(\eta ,\mathit{x})={\sum }_{\mathit{k}}{\widehat{\varphi }}_{\mathit{k}}\left(\eta \right)e^{i\mathit{k}\mathit{x}}$, $\widehat{p}(\eta ,\mathit{x})={\sum }_{\mathit{k}}{\widehat{p}}_{\mathit{k}}\left(\eta \right)e^{i\mathit{k}\mathit{x}}$ etc. For $\mathit{k}\ne 0$, the approximate identifying $\delta {\widehat{\mathsf{\Pi }}}_{\mathit{k}}={\widehat{p}}_{\mathit{k}}$ and $\delta {\widehat{\wp }}_{\mathit{k}}={\widehat{\rho }}_{\mathit{k}}$ results in

$$\begin{array}{c}\hfill \delta {\widehat{\mathsf{\Pi }}}_{\mathit{k}}=\sum _{\mathit{q}}{\displaystyle \frac{1}{2}}{\widehat{\varphi }}_{\mathit{q}}^{+\prime }{\widehat{\varphi }}_{\mathit{q}+\mathit{k}}^{\prime }-{\displaystyle \frac{1}{6}}(\mathit{q}+\mathit{k})\mathit{q}{\widehat{\varphi }}_{\mathit{q}}^{+}{\widehat{\varphi }}_{\mathit{q}+\mathit{k}},\end{array}$$

(2.41)

$$\begin{array}{c}\hfill \delta {\widehat{\wp }}_{\mathit{k}}=\sum _{\mathit{q}}{\displaystyle \frac{1}{2}}{\widehat{\varphi }}_{\mathit{q}}^{+\prime }{\widehat{\varphi }}_{\mathit{q}+\mathit{k}}^{\prime }+{\displaystyle \frac{1}{2}}(\mathit{q}+\mathit{k})\mathit{q}{\widehat{\varphi }}_{\mathit{q}}^{+}{\widehat{\varphi }}_{\mathit{q}+\mathit{k}},\end{array}$$

(2.42)

so that the quantity $3\delta {\widehat{\mathsf{\Pi }}}_{\mathit{k}}+\delta {\widehat{\wp }}_{\mathit{k}}$ from Eq. (3.18) is reduced to

$$3\delta {\widehat{\mathsf{\Pi }}}_{\mathit{k}}+\delta {\widehat{\wp }}_{\mathit{k}}=2\sum _{\mathit{q}}{\widehat{\varphi }}_{\mathit{q}}^{+\prime }{\widehat{\varphi }}_{\mathit{q}+\mathit{k}}^{\prime }.$$

(3.23)

Writing the quantized field explicitly with creation and annihilation operators [23]

$${\widehat{\varphi }}_{\mathit{k}}\left(\eta \right)={\displaystyle \frac{1}{\sqrt{2{\omega }_k}}}\left({\widehat{a}}_{-\mathit{k}}^{+}{e}^{i{\omega }_k\eta }+{\widehat{a}}_{\mathit{k}}{e}^{-i{\omega }_k\eta }\right),$$

(3.24)

allows obtaining from Eqs. (3.23) and (3.24)

$$\begin{array}{c}\hfill 3\delta {\widehat{\mathsf{\Pi }}}_{\mathit{k}}+\delta {\widehat{\wp }}_{\mathit{k}}=\sum _{\mathit{q}}\sqrt{{\omega }_q{\omega }_{|\mathit{q}+\mathit{k}|}}({\widehat{a}}_{-\mathit{q}}{\widehat{a}}_{-\mathit{q}-\mathit{k}}^{+}{e}^{i({\omega }_{|\mathit{q}+\mathit{k}|}-{\omega }_q)\eta }+{\widehat{a}}_{\mathit{q}}^{+}{\widehat{a}}_{\mathit{q}+\mathit{k}}{e}^{i({\omega }_q-{\omega }_{|\mathit{q}+\mathit{k}|})\eta }\\ \hfill -{\widehat{a}}_{-\mathit{q}}{\widehat{a}}_{\mathit{q}+\mathit{k}}{e}^{-i({\omega }_{|\mathit{q}+\mathit{k}|}+{\omega }_q)\eta }-{\widehat{a}}_{\mathit{q}}^{+}{\widehat{a}}_{-\mathit{q}-\mathit{k}}^{+}{e}^{i({\omega }_{|\mathit{q}+\mathit{k}|}+{\omega }_q)\eta }),\end{array}$$

(3.25)

where ${\omega }_{\mathit{k}}=\left|\mathit{k}\right|$ for a massless scalar field. As is seen from Eq. (3.25), the perturbations have the general form:

$$3\delta {\widehat{\mathsf{\Pi }}}_{\mathit{k}}+\delta {\widehat{\wp }}_{\mathit{k}}=\sum _m{\widehat{\mathcal{P}}}_{m\mathit{k}}{e}^{i{\mathsf{\Omega }}_{m\mathit{k}}\eta },$$

(3.26)

where the frequencies ${\mathsf{\Omega }}_{m\mathit{k}}$ take the values of ${\omega }_q-{\omega }_{|\mathit{q}+\mathit{k}|}$, $-{\omega }_q+{\omega }_{|\mathit{q}+\mathit{k}|}$, ${\omega }_q+{\omega }_{|\mathit{q}+\mathit{k}|}$ and $-{\omega }_q-{\omega }_{|\mathit{q}+\mathit{k}|}$. That allows finding the solution of Eq. (3.18) as

$${\widehat{\mathsf{\Phi }}}_{\mathit{k}}\left(\eta \right)=-{\displaystyle \frac{1}{M_p^2}}\sum _m{\displaystyle \frac{{\widehat{\mathcal{P}}}_{m\mathit{k}}{e}^{i{\mathsf{\Omega }}_{m\mathit{k}}\eta }}{{\mathsf{\Omega }}_{m\mathit{k}}^2-k^2/3}}.$$

(3.27)

Using Eqs. (3.14) and (3.27), one comes to

$${\widehat{F}}_{\mathit{k}}\left(\eta \right)=-{\displaystyle \frac{4}{M_p^2}}\sum _m{\displaystyle \frac{{\widehat{\mathcal{P}}}_{m\mathit{k}}{e}^{i{\mathsf{\Omega }}_{m\mathit{k}}\eta }}{{\mathsf{\Omega }}_{m\mathit{k}}^4-k^4/9}}.$$

(3.28)

Under Eqs. (3.25) and (3.27), the final expression for the metric perturbation ${\widehat{\mathsf{\Phi }}}_{\mathit{k}}\left(\eta \right)$ acquires the form

$$\begin{array}{c}\hfill {\widehat{\mathsf{\Phi }}}_{\mathit{k}}\left(\eta \right)={\displaystyle \frac{1}{M_p^2}}\sum _{\mathit{q}}\sqrt{{\omega }_q{\omega }_{|\mathit{q}+\mathit{k}|}}\left({\displaystyle \frac{1}{{({\omega }_{|\mathit{q}+\mathit{k}|}+{\omega }_q)}^2-k^2/3}}\right({\widehat{a}}_{-\mathit{q}}{\widehat{a}}_{\mathit{q}+\mathit{k}}{e}^{-i({\omega }_{|\mathit{q}+\mathit{k}|}+{\omega }_q)\eta }+\\ \hfill {\widehat{a}}_{\mathit{q}}^{+}{\widehat{a}}_{-\mathit{q}-\mathit{k}}^{+}{e}^{i({\omega }_{|\mathit{q}+\mathit{k}|}+{\omega }_q)\eta })-{\displaystyle \frac{1}{{({\omega }_{|\mathit{q}+\mathit{k}|}-{\omega }_q)}^2-k^2/3}}({\widehat{a}}_{-\mathit{q}}{\widehat{a}}_{-\mathit{q}-\mathit{k}}^{+}{e}^{i({\omega }_{|\mathit{q}+\mathit{k}|}-{\omega }_q)\eta }+\end{array}$$

(3.29)

$$\begin{array}{c}\hfill {\widehat{a}}_{\mathit{q}}^{+}{\widehat{a}}_{\mathit{q}+\mathit{k}}{e}^{i({\omega }_q-{\omega }_{|\mathit{q}+\mathit{k}|})\eta }\left)\right),\\ \hfill {\widehat{F}}_{\mathit{k}}\left(\eta \right)={\displaystyle \frac{4}{M_p^2}}\sum _{\mathit{q}}\sqrt{{\omega }_q{\omega }_{|\mathit{q}+\mathit{k}|}}\left({\displaystyle \frac{1}{{({\omega }_{|\mathit{q}+\mathit{k}|}+{\omega }_q)}^4-k^4/9}}\right({\widehat{a}}_{-\mathit{q}}{\widehat{a}}_{\mathit{q}+\mathit{k}}{e}^{-i({\omega }_{|\mathit{q}+\mathit{k}|}+{\omega }_q)\eta }+\\ \hfill {\widehat{a}}_{\mathit{q}}^{+}{\widehat{a}}_{-\mathit{q}-\mathit{k}}^{+}{e}^{i({\omega }_{|\mathit{q}+\mathit{k}|}+{\omega }_q)\eta })-{\displaystyle \frac{1}{{({\omega }_{|\mathit{q}+\mathit{k}|}-{\omega }_q)}^4-k^4/9}}({\widehat{a}}_{-\mathit{q}}{\widehat{a}}_{-\mathit{q}-\mathit{k}}^{+}{e}^{i({\omega }_{|\mathit{q}+\mathit{k}|}-{\omega }_q)\eta }+\end{array}$$

$$\begin{array}{c}\hfill {\widehat{a}}_{\mathit{q}}^{+}{\widehat{a}}_{\mathit{q}+\mathit{k}}{e}^{i({\omega }_q-{\omega }_{|\mathit{q}+\mathit{k}|})\eta }\left)\right),\end{array}$$

(3.30)

Expressions (3.29), (3.30) allows calculating the correlators

$$\begin{array}{c}\hfill S_{\mathsf{\Phi }}(\tau -\eta ,\mathit{k})=<0\left|{\widehat{\mathsf{\Phi }}}_{\mathit{k}}^{+}\left(\eta \right){\widehat{\mathsf{\Phi }}}_{\mathit{k}}\left(\tau \right)\right|0>=\end{array}$$

$$\begin{array}{c}\hfill {\displaystyle \frac{18}{M_p^4}}\sum _{\mathit{q}}{\displaystyle \frac{e^{i(\tau -\eta )({\omega }_{\mathit{q}}+{\omega }_{\mathit{q}+\mathit{k}})}{\omega }_{\mathit{q}}{\omega }_{\mathit{k}+\mathit{q}}}{{\left(k^2-3{({\omega }_{\mathit{q}}+{\omega }_{\mathit{k}+\mathit{q}})}^2\right)}^2}}={\displaystyle \frac{18}{{\left(2\pi \right)}^3{M}_p^4}}\int {\displaystyle \frac{e^{i(\tau -\eta )({\omega }_{\mathit{q}}+{\omega }_{\mathit{q}+\mathit{k}})}{\omega }_{\mathit{q}}{\omega }_{\mathit{k}+\mathit{q}}{d}^3\mathit{q}}{{\left(k^2-3{({\omega }_{\mathit{q}}+{\omega }_{\mathit{k}+\mathit{q}})}^2\right)}^2}},\\ S_F(\tau -\eta ,\mathit{k})=<0\left|{\widehat{F}}_{\mathit{k}}^{+}\left(\eta \right){\widehat{F}}_{\mathit{k}}\left(\tau \right)\right|0>=\end{array}$$

(3.31)

$$\begin{array}{c}\hfill {\displaystyle \frac{32*81}{M_p^4}}\sum _{\mathit{q}}{\displaystyle \frac{e^{i(\tau -\eta )({\omega }_{\mathit{q}}+{\omega }_{\mathit{q}+\mathit{k}})}{\omega }_{\mathit{q}}{\omega }_{\mathit{k}+\mathit{q}}}{{\left(k^4-9{({\omega }_{\mathit{q}}+{\omega }_{\mathit{k}+\mathit{q}})}^4\right)}^2}}={\displaystyle \frac{32*81}{{\left(2\pi \right)}^3{M}_p^4}}\int {\displaystyle \frac{e^{i(\tau -\eta )({\omega }_{\mathit{q}}+{\omega }_{\mathit{q}+\mathit{k}})}{\omega }_{\mathit{q}}{\omega }_{\mathit{k}+\mathit{q}}{d}^3\mathit{q}}{{\left(k^4-9{({\omega }_{\mathit{q}}+{\omega }_{\mathit{k}+\mathit{q}})}^4\right)}^2}},\end{array}$$

(3.32)

which are related to the space-time correlators

$$\begin{array}{c}\hfill <0|\widehat{F}(\eta ,\mathit{x})\widehat{F}(\tau ,{\mathit{x}}^{\prime })|0>=\sum _{\mathit{k}}{S}_F(\tau -\eta ,\mathit{k})e^{i\mathit{k}(\mathit{x}-{\mathit{x}}^{\prime })}={\displaystyle \frac{1}{{\left(2\pi \right)}^3}}\int S_F(\tau -\eta ,\mathit{k})e^{i\mathit{k}(\mathit{x}-{\mathit{x}}^{\prime })}{d}^3\mathit{k},\\ \hfill <0|\widehat{\mathsf{\Phi }}(\eta ,\mathit{x})\widehat{\mathsf{\Phi }}(\tau ,{\mathit{x}}^{\prime })|0>=\sum _{\mathit{k}}{S}_{\mathsf{\Phi }}(\tau -\eta ,\mathit{k})e^{i\mathit{k}(\mathit{x}-{\mathit{x}}^{\prime })}={\displaystyle \frac{1}{{\left(2\pi \right)}^3}}\int S_{\mathsf{\Phi }}(\tau -\eta ,\mathit{k})e^{i\mathit{k}(\mathit{x}-{\mathit{x}}^{\prime })}{d}^3\mathit{k}.\end{array}$$

We will also need the correlators at different continuous values of $\mathit{k}$, ${\mathit{k}}^{\prime }$ :

$$\begin{array}{c}\hfill <0|{\widehat{F}}_{\mathit{k}}\left(\eta \right){\widehat{F}}_{{\mathit{k}}^{\prime }}\left(\tau \right)|0>=\int <0|\widehat{F}(\eta ,\mathit{x})\widehat{F}(\tau ,{\mathit{x}}^{\prime })|0>e^{-i\mathit{k}\mathit{x}-i{\mathit{k}}^{\prime }{\mathit{x}}^{\prime }}{d}^3\mathit{x}{d}^3{\mathit{x}}^{\prime }=\\ \hfill {\left(2\pi \right)}^3{S}_F(\tau -\eta ,\mathit{k})\delta (\mathit{k}+{\mathit{k}}^{\prime }),\end{array}$$

(3.33)

where $\delta \left(\mathit{k}\right)$ is the Dirac delta function. Explicit calculation gives

$$\begin{array}{c}\hfill {\tilde{S}}_{\mathsf{\Phi }}(\omega ,k)={\displaystyle \frac{1}{2\pi }}\int S_{\mathsf{\Phi }}(\eta ,k)e^{-i\omega \eta }d\eta ={\displaystyle \frac{18}{{\left(2\pi \right)}^3{M}_p^4}}\int {\displaystyle \frac{\delta ({\omega }_{\mathit{q}}+{\omega }_{\mathit{q}+\mathit{k}}-\omega ){\omega }_{\mathit{q}}{\omega }_{\mathit{k}+\mathit{q}}{d}^3\mathit{q}}{{\left(k^2-3{({\omega }_{\mathit{q}}+{\omega }_{\mathit{k}+\mathit{q}})}^2\right)}^2}}=\end{array}$$

$$\begin{array}{c}\hfill \left\{\begin{array}{c}{\displaystyle \frac{1}{160{\pi }^2{M_p}^4}}\left(5+{\displaystyle \frac{4k^4}{{\left(k^2-3{\omega }^2\right)}^2}}\right),k<\omega \\ 0,otherwise\end{array}\right.,\end{array}$$

(3.34)

$$\begin{array}{c}\hfill {\tilde{S}}_F(\omega ,k)={\displaystyle \frac{32*81}{{\left(2\pi \right)}^3{M}_p^4}}\int {\displaystyle \frac{\delta ({\omega }_{\mathit{q}}+{\omega }_{\mathit{q}+\mathit{k}}-\omega ){\omega }_{\mathit{q}}{\omega }_{\mathit{k}+\mathit{q}}{d}^3\mathit{q}}{{\left(k^4-9{({\omega }_{\mathit{q}}+{\omega }_{\mathit{k}+\mathit{q}})}^4\right)}^2}}=\left\{\begin{array}{c}{\displaystyle \frac{27\left(15{\omega }^4-10k^2{\omega }^2+3k^4\right)}{10{\pi }^2{M}_p^4{\left(k^4-9{\omega }^4\right)}^2}},k<\omega \\ 0,\mathrm{otherwise}\end{array}\right..\end{array}$$

(3.35)

4. Migdal Equation for Photon Density Matrix Evolution

The kernel of a photon density matrix is defined by $\rho (\mathit{r},{\mathit{r}}^{\prime },\eta )=\psi (\mathit{r},\eta ){\psi }^{+}({\mathit{r}}^{\prime },\eta )$, but we begin with consideration of a density matrix as an operator, and write the equation for its evolution in the standard form using the Hamiltonian $H_0$ and interaction (2.20):

$$i{\partial }_{\eta }\widehat{\rho }=[{\widehat{H}}_0+\widehat{V},\widehat{\rho }].$$

(4.1)

A formal solution of Eq. (4.1) could be written as

$$\widehat{\rho }\left(\eta \right)=-i{\int }_{-\infty }^{\eta }{e}^{i{\widehat{H}}_0(\tau -\eta )}[\widehat{V}\left(\tau \right),\rho \left(\tau \right)]e^{-i{\widehat{H}}_0(\tau -\eta )}d\tau .$$

(4.2)

This expression can be substituted back into Eq. (4.1), and one comes to

$$i{\partial }_{\eta }\widehat{\rho }=[{\widehat{H}}_0,\widehat{\rho }]-i{\int }_{-\infty }^{\eta }[\widehat{V}\left(\eta \right),e^{i{\widehat{H}}_0(\tau -\eta )}[\widehat{V}\left(\tau \right),\widehat{\rho }\left(\tau \right)]e^{-i{\widehat{H}}_0(\tau -\eta )}]d\tau .$$

(4.3)

Change of the time variable $\tau \to \tau +\eta$ in the integral leads to

$$i{\partial }_{\eta }\widehat{\rho }=[{\widehat{H}}_0,\widehat{\rho }]-i{\int }_{-\infty }^0[\widehat{V}\left(\eta \right),e^{i{\widehat{H}}_0\tau }[\widehat{V}(\eta +\tau ),\widehat{\rho }(\eta +\tau )]e^{-i{\widehat{H}}_0\tau }]d\tau .$$

(4.4)

Further approximation is to write $\widehat{\rho }(\tau +\eta )\approx e^{-i{\widehat{H}}_0\tau }\widehat{\rho }\left(\eta \right)e^{i{\widehat{H}}_0\tau }$ on the right-hand side of (4.4), and obtain in the second order on the interaction [24] :

$$\begin{array}{c}\hfill i{\partial }_{\eta }\widehat{\rho }=[{\widehat{H}}_0,\widehat{\rho }]-i{\int }_{-\infty }^0[\widehat{V}\left(\eta \right),e^{i{\widehat{H}}_0\tau }[\widehat{V}(\tau +\eta ),e^{-i{\widehat{H}}_0\tau }\widehat{\rho }\left(\eta \right)e^{i{\widehat{H}}_0\tau }]e^{-i{\widehat{H}}_0\tau }]d\tau =[{\widehat{H}}_0,\widehat{\rho }]-\\ \hfill i{\int }_{-\infty }^0(\widehat{V}\left(\eta \right)e^{i{\widehat{H}}_0\tau }\widehat{V}(\eta +\tau )e^{-i{\widehat{H}}_0\tau }\widehat{\rho }\left(\eta \right)-\widehat{V}\left(\eta \right)\widehat{\rho }\left(\eta \right)e^{i{\widehat{H}}_0\tau }\widehat{V}(\eta +\tau )e^{-i{\widehat{H}}_0\tau }-\\ \hfill e^{i{\widehat{H}}_0\tau }\widehat{V}(\tau +\eta )e^{-i{\widehat{H}}_0\tau }\widehat{\rho }\left(\eta \right)\widehat{V}\left(\eta \right)+\widehat{\rho }\left(\eta \right)e^{i{\widehat{H}}_0\tau }\widehat{V}(\tau +\eta )e^{-i{\widehat{H}}_0\tau }\widehat{V}\left(\eta \right))d\tau .\end{array}$$

(4.5)

As a result, the equation for ${\rho }_{\mathit{p}{\mathit{p}}^{\prime }}=\int e^{-i\mathit{p}\mathit{r}}\widehat{\rho }{e}^{i{\mathit{p}}^{\prime }\mathit{r}}{d}^3\mathit{r}/{\left(2\pi \right)}^3$ acquires the form

$$\begin{array}{c}\hfill i{\partial }_{\eta }{\rho }_{\mathit{p}{\mathit{p}}^{\prime }}=\mathit{\alpha }\mathit{p}{\rho }_{\mathit{p}{\mathit{p}}^{\prime }}-{\rho }_{\mathit{p}{\mathit{p}}^{\prime }}\mathit{\alpha }{\mathit{p}}^{\prime }-{\displaystyle \frac{i}{{\left(2\pi \right)}^6}}\int {\int }_{-\infty }^0(V_{\mathit{p}\mathit{q}}\left(\eta \right)e^{i\mathit{\alpha }\mathit{q}\tau }{V}_{\mathit{q}{\mathit{q}}^{\prime }}(\eta +\tau )e^{-i\mathit{\alpha }{\mathit{q}}^{\prime }\tau }{\rho }_{{\mathit{q}}^{\prime }{\mathit{p}}^{\prime }}\left(\eta \right)-\\ \hfill V_{\mathit{p}\mathit{q}}\left(\eta \right){\rho }_{\mathit{q}{\mathit{q}}^{\prime }}\left(\eta \right)e^{i\mathit{\alpha }{\mathit{q}}^{\prime }\tau }{V}_{{\mathit{q}}^{\prime }{\mathit{p}}^{\prime }}(\eta +\tau )e^{-i\mathit{\alpha }{\mathit{p}}^{\prime }\tau }-e^{i\mathit{\alpha }\mathit{p}\tau }{V}_{\mathit{p}\mathit{q}}(\eta +\tau )e^{-i\mathit{\alpha }\mathit{q}\tau }{\rho }_{\mathit{q}{\mathit{q}}^{\prime }}\left(\eta \right)V_{{\mathit{q}}^{\prime }{\mathit{p}}^{\prime }}\left(\eta \right)+\\ \hfill {\rho }_{\mathit{p}\mathit{q}}\left(\eta \right)e^{i\alpha \mathit{q}\tau }{V}_{\mathit{q}{\mathit{q}}^{\prime }}(\eta +\tau )e^{-i\mathit{\alpha }{\mathit{q}}^{\prime }\tau }{V}_{{\mathit{q}}^{\prime }{\mathit{p}}^{\prime }}\left(\eta \right))d\tau d^3\mathit{q}{d}^3{\mathit{q}}^{\prime }.\end{array}$$

(4.6)

Let us remind that $V_{\mathit{q}{\mathit{q}}^{\prime }}=\int e^{-i\mathit{q}\mathit{x}}\widehat{V}{e}^{i{\mathit{q}}^{\prime }\mathit{x}}d\mathit{x}$ is a function of $\mathit{q},{\mathit{q}}^{\prime }$ and simultaneously $6\times 6$ matrix. According to (2.13), (2.20), (3.1)

$$V_{\mathit{q}{\mathit{q}}^{\prime }}\left(\tau \right)={\displaystyle \frac{1}{2}}\left({\mathcal{D}}_{\mathit{q}-{\mathit{q}}^{\prime }}\mathit{\alpha }{\mathit{q}}^{\prime }+\mathit{\alpha }\mathit{q}{\mathcal{D}}_{\mathit{q}-{\mathit{q}}^{\prime }}\right)F_{\mathit{q}-{\mathit{q}}^{\prime }}\left(\tau \right)\equiv {\mathcal{V}}_{\mathit{q}{\mathit{q}}^{\prime }}{F}_{\mathit{q}-{\mathit{q}}^{\prime }}\left(\tau \right),$$

(4.7)

where ${\mathcal{V}}_{\mathit{q}{\mathit{q}}^{\prime }}$ denotes matrix part of interaction and scalar function $F_{\mathit{q}}\left(\tau \right)=\int F(\mathit{x},\tau )e^{-i\mathit{q}\mathit{x}}{d}^3\mathit{x}$, corresponds to $F(\mathit{x},\tau )={\sum }_{\mathit{k}}{F}_{\mathit{k}}\left(\tau \right)e^{i\mathit{k}\mathit{x}}={\displaystyle \frac{1}{{\left(2\pi \right)}^3}}\int F_{\mathit{k}}\left(\tau \right)e^{i\mathit{k}\mathit{x}}{d}^3\mathit{k}$. Under Eq. (3.1), the six-dimensional matrix ${\mathcal{D}}_{\mathit{q}}$ is represented in the form of four three-dimensional blocks

$${\mathcal{D}}_{\mathit{q}}=\left(\begin{array}{cc}\mathit{q}\otimes \mathit{q}-{\displaystyle \frac{q^2}{3}}I& 0\\ 0& \mathit{q}\otimes \mathit{q}-{\displaystyle \frac{q^2}{3}}I\end{array}\right).$$

(4.48)

The gravitational field is considered as a random field, originating from the zero-point fluctuations. That means that $F_{\mathit{q}-{\mathit{q}}^{\prime }}$ is a random quantity as well as the density matrix. After averaging, one comes to the approximate equation for a mean value of the density matrix $<{\rho }_{\mathit{p}{\mathit{p}}^{\prime }}>$, but to avoid introducing a new designation, we will denote the averaged density matrix by the same symbol and write:

$$\begin{array}{c}\hfill i{\partial }_{\eta }{\rho }_{\mathit{p}{\mathit{p}}^{\prime }}=\mathit{\alpha }\mathit{p}{\rho }_{\mathit{p}{\mathit{p}}^{\prime }}-{\rho }_{\mathit{p}{\mathit{p}}^{\prime }}\mathit{\alpha }{\mathit{p}}^{\prime }-\\ \hfill {\displaystyle \frac{i}{{\left(2\pi \right)}^6}}\int {\int }_{-\infty }^0({\mathcal{V}}_{\mathit{p}\mathit{q}}{e}^{i\mathit{\alpha }\mathit{q}\tau }{\mathcal{V}}_{\mathit{q}{\mathit{q}}^{\prime }}{e}^{-i\mathit{\alpha }{\mathit{q}}^{\prime }\tau }{\rho }_{{\mathit{q}}^{\prime }{\mathit{p}}^{\prime }}\left(\eta \right)<F_{\mathit{p}-\mathit{q}}\left(\eta \right)F_{\mathit{q}-{\mathit{q}}^{\prime }}(\eta +\tau )>-\\ \hfill {\mathcal{V}}_{\mathit{p}\mathit{q}}{\rho }_{\mathit{q}{\mathit{q}}^{\prime }}{e}^{i\mathit{\alpha }{\mathit{q}}^{\prime }\tau }{\mathcal{V}}_{{\mathit{q}}^{\prime }{\mathit{p}}^{\prime }}{e}^{-i\mathit{\alpha }{\mathit{p}}^{\prime }\tau }<F_{\mathit{p}-\mathit{q}}\left(\eta \right)F_{{\mathit{q}}^{\prime }-{\mathit{p}}^{\prime }}(\eta +\tau )>-\\ \hfill e^{i\mathit{\alpha }\mathit{p}\tau }{\mathcal{V}}_{\mathit{p}\mathit{q}}{e}^{-i\mathit{\alpha }\mathit{q}\tau }{\rho }_{\mathit{q}{\mathit{q}}^{\prime }}\left(\eta \right){\mathcal{V}}_{{\mathit{q}}^{\prime }{\mathit{p}}^{\prime }}<F_{\mathit{p}-\mathit{q}}(\eta +\tau )F_{{\mathit{q}}^{\prime }-{\mathit{p}}^{\prime }}\left(\eta \right)>+\\ \hfill {\rho }_{\mathit{p}\mathit{q}}\left(\eta \right)e^{i\alpha \mathit{q}\tau }{\mathcal{V}}_{\mathit{q}{\mathit{q}}^{\prime }}{e}^{-i\mathit{\alpha }{\mathit{q}}^{\prime }\tau }{\mathcal{V}}_{{\mathit{q}}^{\prime }{\mathit{p}}^{\prime }}<F_{\mathit{q}-{\mathit{q}}^{\prime }}(\eta +\tau )F_{{\mathit{q}}^{\prime }-{\mathit{p}}^{\prime }}\left(\eta \right)>)d\tau d^3\mathit{q}{d}^3{\mathit{q}}^{\prime }.\end{array}$$

(4.9)

Correlators in Eq. (4.9) are taken from (3.33) and contain delta functions. That allows us to perform one integration and reduce Eq. (4.9) to

$$\begin{array}{c}\hfill i{\partial }_{\eta }{\rho }_{\mathit{p}{\mathit{p}}^{\prime }}=\mathit{\alpha }\mathit{p}{\rho }_{\mathit{p}{\mathit{p}}^{\prime }}-{\rho }_{\mathit{p}{\mathit{p}}^{\prime }}\mathit{\alpha }{\mathit{p}}^{\prime }-{\displaystyle \frac{i}{{\left(2\pi \right)}^3}}\int {\int }_{-\infty }^0({\mathcal{V}}_{\mathit{p}\mathit{q}}{e}^{i\mathit{\alpha }\mathit{q}\tau }{\mathcal{V}}_{\mathit{q}\mathit{p}}{e}^{-i\mathit{\alpha }\mathit{p}\tau }{\rho }_{\mathit{p}{\mathit{p}}^{\prime }}\left(\eta \right)S_F(\tau ,\mathit{q}-\mathit{p})-\\ \hfill {\mathcal{V}}_{\mathit{p}\mathit{q}}{\rho }_{\mathit{q},\mathit{q}-\mathit{p}+{\mathit{p}}^{\prime }}\left(\eta \right)e^{i\mathit{\alpha }(\mathit{q}-\mathit{p}+{\mathit{p}}^{\prime })\tau }{\mathcal{V}}_{\mathit{q}-\mathit{p}+{\mathit{p}}^{\prime },{\mathit{p}}^{\prime }}{e}^{-i\mathit{\alpha }{\mathit{p}}^{\prime }\tau }{S}_F(\tau ,\mathit{q}-\mathit{p})-\\ \hfill e^{i\mathit{\alpha }\mathit{p}\tau }{\mathcal{V}}_{\mathit{p}\mathit{q}}{e}^{-i\mathit{\alpha }\mathit{q}\tau }{\rho }_{\mathit{q},\mathit{q}-\mathit{p}+{\mathit{p}}^{\prime }}\left(\eta \right){\mathcal{V}}_{\mathit{q}-\mathit{p}+{\mathit{p}}^{\prime },{\mathit{p}}^{\prime }}{S}_{(}-\tau ,\mathit{q}-\mathit{p})+\\ \hfill {\rho }_{\mathit{p}{\mathit{p}}^{\prime }}\left(\eta \right)e^{i\alpha {\mathit{p}}^{\prime }\tau }{\mathcal{V}}_{{\mathit{p}}^{\prime }\mathit{q}}{e}^{-i\mathit{\alpha }\mathit{q}\tau }{\mathcal{V}}_{\mathit{q}{\mathit{p}}^{\prime }}{S}_F(-\tau ,\mathit{q}-{\mathit{p}}^{\prime }))d\tau d^3\mathit{q}.\end{array}$$

(4.10)

As a measure of the purity of a state, a quantity could be introduced

$$\mathcal{C}=\int Tr{\rho }_{\mathit{p}{\mathit{p}}^{\prime }}{\rho }_{{\mathit{p}}^{\prime }\mathit{p}}{d}^3{\mathit{p}}^{\prime }{d}^3\mathit{p}=\int Tr{\rho }_{\mathit{p}\mathit{p}}{d}^3\mathit{p}=1,$$

(4.11)

which equals unity for a completely pure state. The quantity $\mathcal{C}$ measures decoherence, and the symbol $Tr$ denote track of a $6\times 6$ matrix. A time evolution of $\mathcal{C}$ is given by

$$\begin{array}{c}\hfill {\displaystyle \frac{\partial \mathcal{C}}{\partial \eta }}=\int Tr\left({\displaystyle \frac{\partial {\rho }_{\mathit{p}{\mathit{p}}^{\prime }}}{\partial \eta }}{\rho }_{{\mathit{p}}^{\prime }\mathit{p}}+{\rho }_{\mathit{p}{\mathit{p}}^{\prime }}{\displaystyle \frac{\partial {\rho }_{{\mathit{p}}^{\prime }\mathit{p}}}{\partial \eta }}\right)d^3{\mathit{p}}^{\prime }{d}^3\mathit{p}=\\ \hfill -{\displaystyle \frac{2}{{\left(2\pi \right)}^3}}\int {\int }_{-\infty }^{0}Tr\left(\right({\mathcal{V}}_{\mathit{p}\mathit{q}}{e}^{i\mathit{\alpha }\mathit{q}\tau }{\mathcal{V}}_{\mathit{q}\mathit{p}}{e}^{-i\mathit{\alpha }\mathit{p}\tau }{\rho }_{\mathit{p}{\mathit{p}}^{\prime }}\left(\eta \right){\rho }_{{\mathit{p}}^{\prime }\mathit{p}}\left(\eta \right)-\\ \hfill {\mathcal{V}}_{\mathit{p}\mathit{q}}{\rho }_{\mathit{q},\mathit{q}-\mathit{p}+{\mathit{p}}^{\prime }}\left(\eta \right)e^{i\mathit{\alpha }(\mathit{q}-\mathit{p}+{\mathit{p}}^{\prime })\tau }{\mathcal{V}}_{\mathit{q}-\mathit{p}+{\mathit{p}}^{\prime },{\mathit{p}}^{\prime }}{e}^{-i\mathit{\alpha }{\mathit{p}}^{\prime }\tau }{\rho }_{{\mathit{p}}^{\prime }\mathit{p}}\left(\eta \right))S(\tau ,\mathit{q}-\mathit{p})-\\ \hfill (e^{i\mathit{\alpha }\mathit{p}\tau }{\mathcal{V}}_{\mathit{p}\mathit{q}}{e}^{-i\mathit{\alpha }\mathit{q}\tau }{\rho }_{\mathit{q},\mathit{q}-\mathit{p}+{\mathit{p}}^{\prime }}\left(\eta \right){\mathcal{V}}_{\mathit{q}-\mathit{p}+{\mathit{p}}^{\prime },{\mathit{p}}^{\prime }}{\rho }_{{\mathit{p}}^{\prime }\mathit{p}}\left(\eta \right)-\\ \hfill {\rho }_{\mathit{p}{\mathit{p}}^{\prime }}\left(\eta \right){\rho }_{{\mathit{p}}^{\prime }\mathit{p}}\left(\eta \right)e^{i\alpha \mathit{p}\tau }{\mathcal{V}}_{\mathit{p}\mathit{q}}{e}^{-i\mathit{\alpha }\mathit{q}\tau }{\mathcal{V}}_{\mathit{q}\mathit{p}}\left)S(\tau ,\mathit{q}-\mathit{p})\right)d\tau d^3{\mathit{p}}^{\prime }{d}^3\mathit{p}{d}^3\mathit{q}.\end{array}$$

(4.12)

One could be usable to expand the matrix $\mathit{\alpha }\mathit{p}$ over eagenmodes

$$\mathit{\alpha }\mathit{p}|m,\mathit{p}>={\epsilon }_m\left(\mathit{p}\right)|m,\mathit{p}>.$$

(4.13)

There are two longitudinal modes with zero energy, two modes with positive energy, and two modes with negative energy, as shown in the Appendix. Expansion over modes reduces the equation (4.14) to

$$\begin{array}{c}\hfill {\displaystyle \frac{\partial \mathcal{C}}{\partial \eta }}=-{\displaystyle \frac{2}{{\left(2\pi \right)}^3}}\int \sum _{m,n=1}^6(<n,\mathit{p}|{\rho }_{\mathit{p}{\mathit{p}}^{\prime }}\left(\eta \right){\rho }_{{\mathit{p}}^{\prime }\mathit{p}}\left(\eta \right){\mathcal{V}}_{\mathit{p},\mathit{p}+\mathit{q}}|m,\mathit{p}+\mathit{q}><m,\mathit{p}+\mathit{q}|{\mathcal{V}}_{\mathit{q}+\mathit{p},\mathit{p}}|n,\mathit{p}>\\ \hfill \Delta ({\epsilon }_m(\mathit{p}+\mathit{q})-{\epsilon }_n\left(\mathit{p}\right)+\omega )-\\ \hfill <n,\mathit{p}|{\rho }_{\mathit{p}{\mathit{p}}^{\prime }}\left(\eta \right){\mathcal{V}}_{{\mathit{p}}^{\prime },\mathit{q}+{\mathit{p}}^{\prime }}{\rho }_{\mathit{q}+{\mathit{p}}^{\prime },\mathit{q}+\mathit{p}}\left(\eta \right)|m,\mathit{p}+\mathit{q}><m,\mathit{p}+\mathit{q}|{\mathcal{V}}_{\mathit{q}+\mathit{p},\mathit{p}}|n,\mathit{p}>\\ \hfill \Delta ({\epsilon }_m(\mathit{p}+\mathit{q})-{\epsilon }_n\left(\mathit{p}\right)+\omega )-\\ \hfill <n,\mathit{p}|{\mathcal{V}}_{\mathit{p},\mathit{q}+\mathit{p}}|m,\mathit{p}+\mathit{q}><m,\mathit{p}+\mathit{q}|{\rho }_{\mathit{q}+\mathit{p},\mathit{q}+{\mathit{p}}^{\prime }}\left(\eta \right){\mathcal{V}}_{\mathit{q}+{\mathit{p}}^{\prime },{\mathit{p}}^{\prime }}{\rho }_{{\mathit{p}}^{\prime }\mathit{p}}\left(\eta \right)|n,\mathit{p}>\\ \hfill \Delta (-{\epsilon }_m(\mathit{p}+\mathit{q})+{\epsilon }_n\left(\mathit{p}\right)-\omega )+\\ \hfill <m,\mathit{p}+\mathit{q}|{\mathcal{V}}_{\mathit{q}+\mathit{p},\mathit{p}}{\rho }_{\mathit{p}{\mathit{p}}^{\prime }}\left(\eta \right){\rho }_{{\mathit{p}}^{\prime }\mathit{p}}\left(\eta \right)|n,\mathit{p}><n,\mathit{p}|{\mathcal{V}}_{\mathit{p},\mathit{q}+\mathit{p}}|m,\mathit{p}+\mathit{q}>\\ \hfill \Delta (-{\epsilon }_m(\mathit{p}+\mathit{q})+{\epsilon }_n\left(\mathit{p}\right)-\omega )){\tilde{S}}_F(\omega ,q)d\omega d^3{\mathit{p}}^{\prime }{d}^3\mathit{p}{d}^3\mathit{q},\end{array}$$

(4.14)

where $\tilde{S}(\omega ,k)={\displaystyle \frac{1}{2\pi }}\int S(\eta ,k)e^{-i\omega \eta }d\eta$ and

$$\Delta \left(\omega \right)={\int }_{-\infty }^0{e}^{i\omega \tau }d\tau =\pi \delta \left(\omega \right)-i\mathcal{P}{\displaystyle \frac{1}{\omega }}.$$

(4.15)

We could suggest the following program of calculations: one substitutes the density matrix of a wave packet of free electromagnetic waves into the right-hand side of (4.14) and obtains an estimate for ${\displaystyle \frac{\partial \mathcal{C}}{\partial \eta }}$, which allows extracting the typical decoherence time. Further calculations are performed for the Gaussian wave packet

$$\begin{array}{c}\hfill {\mathcal{D}}_{\mathit{p}}\left(\eta \right)={\displaystyle \frac{2{(\Delta p)}^{9/2}}{{\pi }^{3/2}}}{\displaystyle \frac{\mathit{e}\times \mathit{p}}{|\mathit{e}\times \mathit{p}|}}{e}^{-i{\epsilon }_{\mathit{p}}\eta -{(\mathit{p}-{\mathit{p}}_0)}^2/{(\Delta p)}^2},\\ \hfill {\mathcal{B}}_{\mathit{p}}\left(\eta \right)={\displaystyle \frac{1}{{\epsilon }_{\mathit{p}}}}\mathcal{D}(\mathit{p},\eta ),\end{array}$$

(4.16)

where ${\epsilon }_{\mathit{p}}=\left|\mathit{p}\right|=p$, $\Delta p$ is the width of the packet in momentum space, $\mathit{e}$ is some vector characterizing the wave polarization. The $6\times 6$ density matrix of the free electromagnetic field

$${\rho }_{\mathit{p}{\mathit{p}}^{\prime }}\left(\eta \right)={\psi }_{\mathit{p}}\left(\eta \right){\psi }_{{\mathit{p}}^{\prime }}^{+}\left(\eta \right)=\left(\begin{array}{c}{\mathcal{D}}_{\mathit{p}}\\ i{\mathcal{B}}_{\mathit{p}}\end{array}\right)\otimes \left(\begin{array}{cc}{\mathcal{D}}_{{\mathit{p}}^{\prime }}^{*}& -i{\mathcal{B}}_{{\mathit{p}}^{\prime }}^{*}\end{array}\right)$$

(4.17)

can be constructed.

The first and last terms in the brackets on the right-hand side of Eq. (4.14) do not contribute to the quantity ${\displaystyle \frac{\partial \mathcal{C}}{\partial \eta }}$. Second and third terms contain ${\rho }_{\mathit{q}+\mathit{p},\mathit{q}+{\mathit{p}}^{\prime }}$, ${\rho }_{\mathit{p},{\mathit{p}}^{\prime }}$, which makes the integral over $d^3\mathit{p}{d}^3{\mathit{p}}^{\prime }{d}^3\mathit{q}$ convergent due to the finiteness of a momentum wave packet. The integral over $\omega$ is also convergent, because ${\tilde{S}}_F(q,\omega )\sim 1/{\omega }^4$ according to (3.35). It turns out to be that the terms with the Dirac delta function in (4.15) do not contribute to ${\displaystyle \frac{\partial \mathcal{C}}{\partial \eta }}$ due to energy conservation, because $\epsilon (\mathit{p}+\mathit{q})+q-\epsilon \left(\mathit{p}\right)>0$. Thus, the decoherence effect is a purely off-shell effect, which is not related to the real on-shell scattering. After the calculation of tracks the integral (4.14) takes the form

$$\begin{array}{c}\hfill {\displaystyle \frac{\partial \mathcal{C}}{\partial \eta }}={\int }_{\omega >q}g(\mathit{p},{\mathit{p}}^{\prime },\mathit{q},\eta )exp\left(-{\displaystyle \frac{{(\mathit{p}-{\mathit{p}}_0)}^2+{({\mathit{p}}^{\prime }-{\mathit{p}}_0)}^2+{({\mathit{p}}^{\prime }+\mathit{q}-{\mathit{p}}_0)}^2+{(\mathit{p}+\mathit{q}-{\mathit{p}}_0)}^2}{\Delta p^2}}\right)\\ \hfill {\displaystyle \frac{{\tilde{S}}_F(q,\omega )}{\epsilon (\mathit{p}+\mathit{q})+\omega -\epsilon \left(\mathit{p}\right)}}{d}^3\mathit{p}{d}^3{\mathit{p}}^{\prime }{d}^3\mathit{q}d\omega ,\end{array}$$

(4.18)

where $\epsilon \left(\mathit{p}\right)=p=\sqrt{p_x^2+p_y^2+p_z^2}$, ${\mathit{p}}_0=\{0,0,p_0\}$, and $g(\mathit{p},{\mathit{p}}^{\prime },\mathit{q},\eta )$ is a some function originating from the calculation of traces. Integral (4.18) is complicated. Simplification consists in an expansion of the function $g(\mathit{p},{\mathit{p}}^{\prime },\mathit{q},\eta )$ into the Taylor series in the vicinity $\mathit{p}={\mathit{p}}_0,{\mathit{p}}^{\prime }={\mathit{p}}_0$ and leads to estimation, which finally looks as follows

$${\displaystyle \frac{\partial \mathcal{C}}{\partial \eta }}\approx -{\displaystyle \frac{2.5\times {10}^{-6}{(\Delta p)}^8}{p_0^2{M}_p^4}}\eta ,$$

(4.19)

and, after integration over time $\eta$, we come to

$$\mathcal{C}\approx 1-{\displaystyle \frac{1.25\times {10}^{-6}{(\Delta p)}^8}{p_0^2{M}_p^4}}{(\eta -{\eta }_i)}^2,$$

(4.20)

where ${\eta }_i$ is the initial conformal time at which a pure wave packet was emitted. For an expanding universe, the conformal time is related to the redshift $z=1/a-1$ as $dz=-{\displaystyle \frac{1}{a^2}}{\displaystyle \frac{da}{d\eta }}d\eta$, that gives $\eta \left(z\right)={\eta }_i-{\int }_{z_i}^z{\displaystyle \frac{dz}{H\left(z\right)}}\approx {\eta }_i+{\displaystyle \frac{z_i-z}{H_0}}$, where $H_0$ is the Hubble constant. In the terms of redshift, the expression (4.20) is rewritten as

$$\mathcal{C}\approx 1-{\displaystyle \frac{1.25\times {10}^{-6}{(\Delta p)}^8}{p_0^2{M}_p^4{H}_0^2}}{(z_i-z)}^2,$$

(4.21)

A condition that $\mathcal{C}$ turns to zero at present time ${\eta }_0$, when $z=0$ gives an estimation of $z_i$:

$$z_i\sim 900{\displaystyle \frac{H_0{p}_0{M}_p^2}{{(\Delta p)}^4}}.$$

(4.22)

For instance, at $p_0\sim 30GeV$ and $\Delta p/p_0\sim 0.05$ we have $z_i\sim 0.25$. That means that radiation emitted at this $z_i$ must be fully decoherent today. Thus, the effect seems rather observable, although the decoherence depends strongly on $\Delta p/p_0$.

5. Decoherence and Depolarization

A decrease of $\mathcal{C}$ suggests that spatial and polarization coherence is lost. Let’s say spatial coherence is already lost, and only polarization decoherence remains. How is the polarization decoherence related to the polarization of radiation? The simplest way to clarify this is to consider the usual two-dimensional photon polarization matrix, corresponding to two possible photon polarizations. Let’s say this density matrix has the form

$$\rho =\left(\begin{array}{cc}\alpha & \beta \\ {\beta }^{*}& 1-\alpha \end{array}\right).$$

(5.1)

$$Tr\left({\rho }^2\right)-1=2\left({\alpha }^2-\alpha +{\left|\beta \right|}^2\right)$$

(5.2)

The minimum of this quantity, corresponding to maximal decoherence, is reached at $\alpha =1/2$ and $\beta =0$ if to consider $\alpha$ and $\beta$ as independent variables. In this case, radiation will be unpolarized. It seems that decoherence also destroys the polarization of radiation. Indeed, this question needs more careful investigation because it is related to the conservation of orbital and total angular momenta [10].

6. Conclusion

We have considered the decoherence of a Gaussian wave packet in vacuum. The attractive feature of the resulting formula is that it does not depend on the ultraviolet cut-off because the corresponding correlator $S_F(q,\omega )$ of the metric fluctuations decreases in a sufficient degree with a frequency increase.

The next stage should be a discussion of possible experiments to observe the incoherent radiation from distant astrophysical sources.