Intrinsic Entropy of Squeezed Quantum Fields and Nonequilibrium Quantum Dynamics of Cosmological Perturbations

Abstract

Density contrasts in the universe are governed by scalar cosmological perturbations which, when expressed in terms of gauge-invariant variables, contain a classical component from scalar metric perturbations and a quantum component from inflaton field fluctuations. It has long been known that the effect of cosmological expansion on a quantum field amounts to squeezing. Thus, the entropy of cosmological perturbations can be studied by treating them in the framework of squeezed quantum systems. Entropy of a free quantum field is a seemingly simple yet subtle issue. In this paper, different from previous treatments, we tackle this issue with a fully developed nonequilibrium quantum field theory formalism for such systems. We compute the covariance matrix elements of the parametric quantum field and solve for the evolution of the density matrix elements and the Wigner functions, and, from them, derive the von Neumann entropy. We then show explicitly why the entropy for the squeezed yet closed system is zero, but is proportional to the particle number produced upon coarse-graining out the correlation between the particle pairs. We also construct the bridge between our quantum field-theoretic results and those using the probability distribution of classical stochastic fields by earlier authors, preserving some important quantum properties, such as entanglement and coherence, of the quantum field.

1. Introduction

Entropy of quantum cosmological perturbations is an important topic which has by now clocked almost three decades of investigations [1,2,3]. In terms of its theoretical foundation, it is built upon the bigger issue of (I) the entropy of quantum fields, where investigation started in the mid-1980s [4,5,6,7] continuing on into the 1990s [8,9,10,11,12], 2000s [13,14,15,16,17] and has been rekindled in recent years [18,19,20,21,22]. There are two aspects in this theme: (1) the quantum field theory component depicting particle creation from the vacuum and quantum cosmological perturbations; (2) the nonequilibrium statistical mechanics aspect describing the evolutionary dynamics of a quantum many-body system—the fluctuations in, and the dissipation of, a quantum field. These two components when combined make up the quantum field theory of nonequilibrium systems [23,24,25], the two major paradigms being the Boltzmann correlation hierarchy and the Langevin open systems [26]. These two components or aspects enter into many theoretical issues of fundamental interest such as quantum (de)coherence, quantum correlations and quantum entanglement. Concepts in open quantum systems [27,28,29] have been introduced, and advanced techniques in quantum field theory [7,23,24,25] have been utilized for these purposes [26]. Against this backdrop we wish to investigate the entropy of quantum fields and cosmological perturbations in a dynamical setting, such as in the early universe. Saving the discussions of the two theoretical factors for the next section, we focus here on the cosmological aspect; (II) The cosmological factor also has two components. They are: (1) gravitational perturbation theory which describes how weak classical perturbations of scalar (density contrast) vector (vorticity) and tensor (gravitational waves) types evolve in a dynamical spacetime; and (2) quantum matter field processes such as particle creation from vacuum fluctuations amplified by the expansion of the universe and their consequences.

1.1. The Cosmological Aspect: Gravitational Perturbations and Quantum Fluctuations

1.1.1. Classical Gravitational Perturbations: Scalar, Vector and Tensor Components

Gravitational perturbation theory has been a well established subject in cosmological structure formation since the 1946 seminal paper of Lifshitz based on the amplification of density contrasts related to the (scalar sector of the) metric perturbations [30,31,32,33,34]. What seeded the structures in the classical gravitational perturbation theory was assumed to be from white noises. After the advent of the inflationary universe [35,36,37] where the vacuum energy density of a quantum scalar field, the inflaton, is believed to have driven the cosmos to (near-) exponential expansion for a certain duration in the early universe, one needs to take into account how the quantum scalar field and its fluctuations are coupled to the (scalar sector of the) metric perturbations controlling the density contrasts. If we regard Lifshitz’s 1946 paper as opening the first (classical) stage of cosmological structure formation investigation based on gravitational perturbation theory; this second (quantum) stage involving the inflaton’s fluctuations began in 1982 [38,39,40,41,42,43].

With the attention now focused on everything we see today in the universe as originating from quantum field mediated-gravitational perturbations, it is worth mentioning another important process involving fluctuations of quantum fields in the early universe, specifically, cosmological particle creation from the vacuum. The former class of activities placed in the inflationary universe context could take place as early as the GUT time (${10}^{-35}$ s) while the latter happened even earlier, predominantly at the Planck time (${10}^{-43}$ s).

1.1.2. Quantum Field Processes Involving Vacuum Fluctuations

Quantum field processes in curved spacetime such as Hawking radiation from black holes or cosmological particle creation in the very early universe are described in great detail in many monographs, for example [44,45,46,47,48]. The amplitude of a classical wave mode can be parametrically amplified by a time-dependent drive [49]. It is the same way for vacuum fluctuations in a quantum field, resulting in the creation of particle pairs: the expansion of the universe acting like a drive, parametrically amplifying the quantum noise, giving rise to spontaneous particle production [49,50]. If there were particles already present in an initial state they will be amplified with the same amplification factor in stimulated production. This is like the spontaneous and stimulated emission of atoms in quantum optics. In fact, quantum field processes in a time-dependent background, be it in a cosmological spacetime or in an external laboratory field [51], can be captured by the ‘squeezing’ of quantum states [52,53,54,55]. The vacuum is ‘squeezed’ in the evolutionary history while particles are produced [56]. A summary description of cosmological particle creation in terms of squeezing can be found in, for example, [11,57,58].

We mention quantum cosmological perturbations and cosmological particle creation together because they are subjected to the same mechanism, which amplifies these perturbations or fluctuations, namely, as in parametric oscillators, where the frequencies of the normal modes are time-dependent. While parametrically amplified quantum fluctuations engender particle creation, parametrically amplified gravitational perturbations engender galaxies and structures, either classically with seeds of white noise or by the inflaton field’s quantum fluctuations.

1.1.3. Distinguish Classical Perturbations from Quantum Fluctuations

It is of theoretical significance to make the distinction between classical linear gravitational perturbations, which are believed to be the progenitors of galaxies and structures we see today, and vacuum fluctuations of a quantum field, which engender the spontaneous creation of particle pairs, a subject fundamental to quantum field theory in curved spacetime. When we talk about the cosmological density contrasts, the isocurvature perturbations, the vorticity and the primordial gravitational waves, we are referring to quantities derived from the scalar, vector and tensor perturbations of the background spacetime. Density contrasts are derived from the scalar sector of the metric perturbations, a subject well explored in classical general relativity, culminating in Bardeen’s gauge invariant quantities [34,42]. In inflationary cosmology they are coupled to a quantum scalar field, the inflaton. To get one compact equation of motion for the density contrast, mixed metric perturbations + scalar field variables are used, such as the gauge invariant Mukhanov-Sasaki variable. This is all standard and fine. However, when one moves to the quantum theory of cosmological perturbations in inflationary universe, one has to be careful what this means, especially when dealing with decoherence and quantum to classical transition issues. The scalar (inflaton) field $\widehat{\varphi }(\mathit{x},t)$ is intrinsically quantum in nature and comes with its quantum fluctuations. Often a background field expansion $\widehat{\varphi }(\mathit{x},t)=\overline{\varphi }(\mathit{x},t)+\widehat{\delta }\varphi (\mathit{x},t)$ is performed, assuming the background field $\overline{\varphi }(\mathit{x},t)$ is classical and the quantum character shows through its fluctuations $\widehat{\delta }\varphi (\mathit{x},t)$. However, the gravitational perturbation component in this mixed variable is of classical origin. The scalar sector of the metric perturbations related to the Newtonian potential is a constraint, not a dynamical degree of freedom (the tensor modes, the gravitational waves, are). Its nature is determined by (or ‘slaved’ to [59,60]) the matter source. In general relativity when the matter source is classical, this scalar sector of the metric perturbation is classical. In inflationary cosmology, what determines the density contrast comes from both the classical scalar metric perturbations and the quantum inflaton field. When one says, ‘quantize the mixed variable’, one should bear in mind that the intrinsic nature of metric perturbations remains classical. In an extreme case, one may even conjure up situations where the quantum fluctuations of the scalar field are made to vanish, such as “choosing a (co-moving) gauge which for scalar perturbations makes the velocity perturbation vanish. For single field inflation, this means that the time coordinate is defined so that at any given time the scalar field equals its unperturbed value” (one is riding up and down with the scalar field’s fluctuations), “with all perturbations relegated to components of the metric” ([61], Section 5.3D). This does not mean that gravity has become quantum, only that the scalar perturbations now acquire a quantum nature by virtue of the presence of the inflaton field. Put it in another way, when there is no inflaton, one returns to purely classical general relativity. The Newtonian force is slaved to the source which is classical. There is no way for the gravitational perturbations to become quantum.

Now about the tensor sector. Gravity’s dynamical (or propagating) degrees of freedom reside only in the tensor sector, that is, the gravitational waves. Primordial gravitational waves are described by the tensor component of gravitational metric perturbations. They have also been studied at the classical level since 1946. One can consider quantizing the linearized tensor perturbations, whence they become the primordial gravitons (Note gravitational waves are weak metric perturbations, like sound waves, barring the differences between transverse and longitudinal waves. Gravitons are quantized linear perturbations of spacetime, like phonons, in the nature of collective excitations. This is a far cry from quantum gravity, defined as theories for the microscopic structures of spacetime at the Planck scale [62]). These gravitons still obey deterministic equations of motion. They are not stochastic intrinsically. Only when one considers a large number of primordial gravitons and use statistical means to describe their distributions would they become ‘stochastic’ (Probing further into their connection in cosmological perturbation theory Roura and Verdaguer [63] show that in the theoretical framework of stochastic gravity [45], at the Gaussian level, the stochastic variables give an equivalent description as the quantized linear perturbations [42]).

1.2. Related Issues: Decoherence and Entanglement

Before we delve into a full discussion of the entropy of cosmological perturbations and fluctuations we want to mention two issues in quantum information related to entropy, namely, decoherence and entanglement. This would be brief because each topic would merit a separate paper or two to describe.

(1) Decoherence of cosmological perturbations. Again, many questions are asked and different approaches have been suggested. To give two early examples, in Guth and Pi [64], by using an inverted harmonic oscillator model as an example and invoking the uncertainty relation they showed that the inflaton will, at late times, behave classically. Note that this is for a closed system thus, strictly speaking, it is at most dephasing, rather than decoherence. In Starobinsky’s stochastic inflation model [65], there are two key points: (a) the long wavelength perturbations are assumed to behave classically and (b) the short wavelength perturbations are treated as white noise. Under what conditions could (a) be realized and (b) be implemented? How could it be that a free field partitioned into two sectors in an expanding universe that the short wavelength sector would behave like white noise? One should also clarify the conditions this noise can effectively decohere the long wavelength sector by its backreaction. These questions were asked and challenged by, for example, [12,66]. This area of research begun in earnest in the early 1990s [12,67,68,69,70,71,72] is still being pursued vibrantly today. For a sample of recent work, see, for example, [21,22], where other earlier references can be found. We shall discuss this topic with a detailed background summary in a companion paper under preparation [73];

(2) Entanglement and entropy. The entanglement between particle pairs, one with momentum $\mathit{k}$, the other with $-\mathit{k}$, such as studied in [19,74] is relevant to our present consideration of entropy associated with particle creation. Another related topic is entanglement entropy. The seminal papers [75,76] explore the entropy of free quantum fields in Minkowski space with a partition, providing a more general and basic statistical mechanical way to understand black hole entropy. It was the precursor of the by-now familiar topic of entanglement entropy [77,78] hotly pursued in the last two decades. In cosmology when the spacetime is dynamical the natural ‘partition’ which separates modes into two sectors is the Hubble horizon defined by the inverse of the expansion rate. Entanglement between two Unruh–DeWitt detectors has been studied concerning the effects of spacetime curvature [79,80,81,82] and topology [83] and under other conditions in the emergent field of relativistic quantum information [84], but that is not the main concern in our present study. We now turn to the two aspects mentioned in the beginning, that of quantum field theory and nonequilibrium dynamics.

2. Entropy of Quantum Fields

Defining and understanding entropy for quantum fields is a task of fundamental significance in both quantum field theory and nonequilibrium statistical mechanics. For free fields in a dynamical setting (driven by some external source or in a dynamical spacetime as in cosmology) with no spatial boundary or event horizon present (thus not involving entanglement entropy) the authors of [4] provided a first answer to the following question: Is there entropy production in particle creation from the vacuum? Adhering to quantum field theory, because the vacuum is a pure state and particle creation is a unitary process, the answer should be no. However, even in textbooks one sees that an entropy, say, of photons, can be assigned proportionally to the number of particles present. These two seemingly contradictory answers each seem to stand on its firm ground. How does one reason out their differences? This paradox was what lured one of us into exploring this issue and came out with some interesting discoveries. The first observation is: the first answer seems to focus on the initial state while the second answer on the final outcome. For the second answer to make sense some information must have been lost. Two possibilities come to one’s mind. Either (a) some essential information is ignored, implicitly when one defines the field entropy in terms of particle numbers, or (b) some kind of coarse-graining measure is introduced explicitly which curtails some information of the system.

For (a), when one argues that the entropy of a free quantum field is proportional to the particles created, one often implicitly adopts a Fock space representation to capture the particle number, while ignoring completely the information about the correlation and the coherence in the particle pairs. This was pointed out in [4]. Note that the entropy of a closed system is defined by taking the full trace, which is independent of the basis or representation, and there is no priori definition or requirement that the entropy be a function of the particle number. This point is illustrated explicitly in [85] where an equation governing the quantum coherence is presented alongside that governing the particle numbers, from which the entropy is calculated. For (b), in a realistic setting, after particles are created from the vacuum, they interact, and information about their interaction could be lost, e.g., if one focuses on the scattering cross section but ignores some of their correlations. This second point is illustrated in an interacting harmonic oscillator model of field theory in [5] and later with greater depth in [14], which proves an H-theorem in such coarse-grained systems.

2.1. Intrinsic Entropy of Free Quantum Fields

Entropy of quantum fields is a fundamental, seemingly simple, yet somewhat tricky issue. Early conceptual inquiries [4,5] made explicit the specific underlying conditions which allow for the entropy of free quantum fields to be related to particle creation numbers (for boson fields). Namely, by adopting a Fock space representation and making statements only in terms of number operators, one implicitly ignores all quantum phase information which determines the coherence (See [85] for the interplay of both quantities). From this, it is easy to understand why the entropy associated with particle production is proportional to the degree of squeezing, in cosmology by the expansion of the universe [3,8,9,10,86], or more down to earth, the moving mirror in the dynamical Casimir effect [87].

2.2. Entropy of Interacting Quantum Fields

We have alluded to this class of theories above, the simplest exemplified by the anharmonic oscillator models in quantum mechanics. A more systematic analysis is based on the Boltzmann paradigm applied to the BBGKY hierarchy, namely, the ‘slaving’ of the higher order correlation functions to the lower ones, with Boltzmann equation being the lowest order, describing the dynamics of the one point distribution function with source from the two particle collision integral. The assumption of a causal factorizable initial condition on the two-particle correlation function gives rise to a dissipative dynamics in the one particle distribution function. For interacting quantum field, it is the concept of correlation entropy defined in the Schwinger–Dyson hierarchy for a certain nth order correlation functions (see [13,14,88,89]). Notice the fundamental difference between the ‘truncation’ and the ‘slaving’ procedures, the former giving rise to a unitary equation such as for the mean field in the Vlasov equation, the latter a (nonunitary) dissipative Boltzmann equation.

2.3. Entropy Measures

As we stated in the beginning and as is well known, entropy arises from a loss of information, either by choice or by necessity. Thus, any mention of entropy measure must be accompanied by the specification of what information is being dropped, lost, or operationally inaccessible. That, in turn, depends on one’s coarse-graining choices or one’s level of ignorance. We give as an example two commonly used coarse-graining measures.

2.3.1. Coarse-Graining by Discrepancies of Scales

If the relevant physical variables in a system have clear discrepancies in scales, be it slow vs. fast, low vs. high frequency, light vs. heavy mass, one can decide based on physical conditions which scale is of the most physical interest and proceed to coarse-grain the other variables, one level at a time in an orderly nested way. Famous examples are the slow variables of van Hove, the adiabatic invariants for a parametric process, the Born-Oppenheimer approximation wherein the heavy mass (e.g., nucleus) motion is considered as effectively static with respect to the light mass (e.g., electron), etc. Formally, these situations are best described by the projection operator formalism of Zwanzig-Mori-Nakajima [90]. In quantum systems, one measure often adopted is the coarse-graining of the phase information since the phase often varies rapidly. This is behind the random phase approximation. The rotating wave approximation (RWA) is one of the two preambles of quantum statistical mechanics, where the density matrix of a quantum system is assumed to be diagonal. By eliminating all phase information one can consider only probabilities, which is why quantum statistical mechanics taught from textbooks misses all the quantum coherence and entanglement issues.

2.3.2. Coarse-Graining by Partitioning

If there is a clear cut partition in the system one is interested in, one can use it to separate the variables of special interest from other variables perhaps of lesser interest. In the stochastic inflation example, the horizon is the divider between the long wavelength sector and the short wavelenght sector. The physical issues of concern are accordingly different: for the former, decoherence—how could they be viewed as classical? For the latter, noise—do they really constitute white noise? For interacting quantum fields, if there is only one field and one has chosen the appropriate partition one can invoke the coarse-grained effective action (CGEA) [91,92,93,94] to describe the backreaction effects of the coarse-grained sector viewed as the environment to the system. The ‘in–in’, Schwinger–Keldysh, or ‘closed-time-path’ version of CGEA is equivalent to the Feynman–Vernon influence functional formalism. This approach was used to describe cosmological decoherence in [12,68,69,72]. If there is more than one field, it is even easier, because the system field and the environment field are distinct and their coupling does not change in time. One can decide which field should be considered as one’s system and coarse-grain the others [20,95,96,97]. One does not need to keep track of the time-varying partition in an expanding universe.

In cosmology, when the spacetime is dynamical, the natural ‘partition’ is between the modes of wavelengths greater or less than the horizon, defined by the inverse of the Hubble expansion rate—between the super and sub horizon sectors. This is referred to as entanglement entropy in cosmology. Recent work includes [22,98,99,100,101], where earlier references can be found.

2.4. Approaches Taken Here in Relation to Earlier Work

From the pioneering work in the 1980s and 1990s mentioned above, it is encouraging to see that the conceptual framework of open quantum systems and the technical tools of nonequilibrium quantum field theory are more widely adopted in current theoretical cosmological research. For recent work see, for example, [20,21,22,97,100,101,102] and references therein. We now focus on the specific main points on the entropy of cosmological perturbations, analyzed here from a squeezed open quantum system perspective using techniques from nonequilibrium quantum field theory and relate them to earlier work.

(1) The theoretical framework useful for our purpose is known as ‘squeezed open quantum systems’ [57,58,86] with cosmology being an important exemplary category [11]. For the entropy measure we shall follow:

(A) the vein of [4,85] where, using a Fock basis measuring the number of particles created, or the cosmological perturbations, we show explicitly when the quantum correlations between the particle pairs are kept, there is no entropy generation, whereas when the correlation is ignored, we obtain the well known expression for the entropy proportional to the number of particles. Related to [85] is the work of [17] where the entropy is calculated from the von Neumann entropy constructed from the reduced density matrix defined with respect to a second adiabatic order number state (Of interest is their conclusion that there is no decoherence as claimed in [103], despite their ‘cute’ way of explaining why there is decoherence);

(B) Our demonstration that there is entropy generated by coarse-graining the $-\mathit{k}$ component of the particle pair corroborates the results obtained by [74] via a different pathway, by showing the destruction of the entanglement between the particle pair.

(2) The formal structure. Many authors, including those of the major work [1,2] (BMP) and [19] (CP) on the entropy of cosmological perturbations, correctly start with the formal structure of an interacting quantum field theory. Treating interacting fields in the framework of the classical BBGKY or the quantum Schwinger–Dyson hierarchy—the Boltzmann paradigm—is a serious challenge. It was taken up in earnest in the work of [12,13,14,88], developing the concepts of correlation histories, correlation dynamics and correlation noise. Their master effective action and the way how the higher order correlators should be ‘slaved’ to a lower order correlator provide the basis for formulating effective field theories for open quantum systems [104], an important theme which saw fruitful applications to cosmological issues in [20,97]. In the same vein, entropy from non-Gaussian states is pursued by [18,105]. However, BMP and CP quickly assume a Gaussian truncation of the BBGKY or the Schwinger–Dyson hierarchy, thus restricting their attention only to the lowest two-point functions. This is enough to treat linear perturbations and effectively places their theories on the same footing as free fields, where considerations in [4,74] would suffice. CP refer to the reduced density matrices defined by truncating the hierarchy at the first level as Gaussian and homogeneous density matrices (GHDM). Coarse-graining the phase information is also used in two major works: In [1,2], their coarse-grained entropy is obtained by averaging the two point correlation function over the phases. In [19] the growth of entropy associated with the coarse-graining is described by the set of Gaussian two-mode reduced density matrices ${\rho }_{(+\mathit{k},-\mathit{k})}^{\mathrm{red}}$, each characterizes the loss of entanglement between the two modes in the presence of interactions.

This is also the starting point of the way BMP defines entropy for the gravitational field. However, instead of treating a quantized field, they consider a classical stochastic field. Thus their considerations are predicated upon the quantum perturbations having effectively decohered: “Provided there is decoherence and $sinh{r}_{\mathit{k}}\gg 1$ (where $r_{\mathit{k}}$ is the squeeze parameter of the $\mathit{k}$th mode), then we can take the classical limit where the classical correlation functions can be identified with quantum expectation values, i.e., $\langle \varphi \left(x\right)\varphi \left(y\right)\rangle =\langle 0|\widehat{\varphi }\left(x\right)\widehat{\varphi }\left(y\right)|0\rangle$”. Note the angular brackets on the lefthand side denotes a statistical average. The decoherence procedure is clearly stated in [3] as dropping the off-diagonal components of the density matrix. At the time these papers were written decoherence as an important theoretical issue was not that well understood. See more recent papers, such as [18,22,73,105], where the developments in the last 30 years can be traced back. We have devoted a subsection at the end of this paper to clarify the theoretical differences between a quantum field theory treatment and a classical stochastic field theory treatment albeit the two can serve some common purposes.

2.5. Our Findings and Organization

Our narrative takes on four steps, showing: (1) that the cosmological expansion results in the squeezing of the quantum field; (2) the nonequilibrium dynamics of a squeezed quantum field as a closed system, seen through the evolution of the density matrix operator in terms of the covariance matrix elements; (3) that the von Neumann entropy is zero in this closed system, as expected, and proportional to the particle numbers when the quantum correlation between the particle pair is coarse-grained away. We derive (4) the Wigner function of the quantum Gaussian field driven by a parametric process, which gives a clear picture of how the coarse-grained entropy may emerge, and whose phase-space description leads to the quantum thermodynamics of this system. With this, we also discuss (5) the differences and commonalities between our quantum field-theoretical methodology and approaches in the literature using the probability distribution of classical stochastic fields, especially on the ‘uniquely quantum’ features such as decoherence and entanglement issues. While facts in (1) and (3) are largely known before, we try to capture a global perspective and point out the subtleties in the meaning of entropy and how it is defined in closed quantum systems. Points (2), (4), and (5) form the technical core and the results we obtained are new.

This paper is organized as follows: in Section 3, to be self-contained, we provide a quick overview of the essential ingredients of first-order classical cosmological perturbation theory, and link the cosmological perturbations to the classical parametric field via the Mukhanov–Sasaki variables. In Section 4 we formally compute the covariance matrix elements of the parametric field as the basic building blocks of a Gaussian system. Then in Section 5, we decompose the field into a collection of parametric oscillators to take advantage of the oft-used nonequilibrium formalism to construct the density matrix operator and the density matrix elements. Doing so equips us to tackle the full dynamical evolution of the physical observables and derive the nonequilibrium thermodynamics of the system under study. We present the time evolution of the density matrix and its elements in Section 6, where we combine the language of the Bogoliubov transformation and the aforementioned nonequilibrium formalism for a succinct description of the full dynamics. In Section 7, we show entropy production associated with the particle pairs in the parametric development of the cosmological perturbations. In Section 8, we construct the Wigner function of the quantum field and show how it evolves in time. As a bonus we also show the relation between our quantum field theoretic results and those using probability distribution of classical stochastic fields by earlier authors. We point out that stochastic cannot be a complete replacement of quantum, especially over the core quantum issues like entanglement and quantum coherence of the field.

3. Cosmological Perturbation Theory in a Nutshell

In this short section we give a brief summary of cosmological perturbation theory [42,106] to make our presentation self-contained and to define the notations, which have evolved to become quite standardized over the years. Readers familiar with this topic can skip to the next section.

3.1. Gauge Transformations

The linear perturbations $\delta g_{\mu \nu }$ off a background spacetime with metric ${\overline{g}}_{\mu \nu }$. In the case of a Friedmann-Lemaître-Robertson-Walker (FLRW) universe background with scale factor a we can write the metric as:

$$g_{\mu \nu }={\overline{g}}_{\mu \nu }+\delta g_{\mu \nu }=a^2\left({\gamma }_{\mu \nu }+h_{\mu \nu }\right),$$

(1)

where $h_{\mu \nu }$ and its derivatives are assumed to be small compared to the background metric ${\gamma }_{\mu \nu }$, which is used for the raising and lowering of the indices for the metric perturbations. Namely,

$$\begin{array}{cccc}\hfill h^{\mu }{}_{\nu }& ={\gamma }^{\mu \alpha }{h}_{\alpha \nu },\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hfill & \hfill h^{\mu \nu }& ={\gamma }^{\mu \alpha }{\gamma }^{\nu \beta }{h}_{\alpha \beta },\hfill \end{array}$$

(2)

such that the inverse metric of the perturbed universe is:

$$g^{\mu \nu }=a^{-2}\left({\gamma }^{\mu \nu }-h^{\mu \nu }\right).$$

(3)

The metric perturbations $h_{\mu \nu }$ can be decomposed into the scalar, vector and tensor components, but the scalar component of the perturbations is essential for the consideration of cosmological density contrasts. The line element can therefore be parametrized as:

$$d{s}^2=a^2\left(\eta \right)\left\{-\left(1+2\mathcal{A}\right)d{\eta }^2-2{\mathcal{B}}_{,i}d\eta d{x}^i+\left[\left(1-2\psi \right){\delta }_{ij}+2{\partial }_i{\partial }_j\mathcal{E}\right]d{x}^{i}d{x}^j\right\},$$

(4)

by four scalars $\mathcal{A}$, $\mathcal{B}$, $\psi$ and $\mathcal{E}$. Here, $\eta$ is the conformal time coordinate and $x^i$ are the spatial coordinates and we assume the background universe is spatially flat.

Now consider the choice of gauge and the construction of gauge-invariant variables. The infinitesimal coordinate transformation between two coordinate systems, say $\left\{x^{\alpha }\right\}$ and $\left\{{\tilde{x}}^{\alpha }\right\}$, in the perturbed spacetime is given by:

$${\tilde{x}}^{\alpha }=x^{\alpha }+{\xi }^{\alpha },$$

(5)

where ${\xi }^{\alpha }$ and their derivatives are first-order small in h. It then can be shown that the four scalars in the line element will transform as:

$$\begin{array}{cccccc}\hfill \tilde{\mathcal{A}}& =\mathcal{A}-{\xi }^0{}^{\prime }-\frac{a^{\prime }}{a}{\xi }^0,\hfill & \hfill \tilde{\psi }& =\psi +\frac{a^{\prime }}{a}{\xi }^0,\hfill & \hfill \tilde{\mathcal{B}}-{\tilde{\mathcal{E}}}^{\prime }& =\mathcal{B}-{\mathcal{E}}^{\prime }+{\xi }^0,\hfill \end{array}$$

(6)

where the prime denotes the derivative with respect to the conformal time, and $\xi$ is the scalar component of the spatial vector ${\xi }^i$. Among them, we can form two gauge-invariant Bardeen potentials:

$$\begin{array}{cccc}\hfill \Phi & =\mathcal{A}+\frac{a^{\prime }}{a}\left(\mathcal{B}-{\mathcal{E}}^{\prime }\right)+{\left(\mathcal{B}-{\mathcal{E}}^{\prime }\right)}^{\prime },\hspace{1em}\hfill & \hfill \Psi & =\psi -\frac{a^{\prime }}{a}\left(\mathcal{B}-{\mathcal{E}}^{\prime }\right),\hfill \end{array}$$

(7)

which are invariant under the coordinate transformation (5), and $\tilde{\Phi }=\Phi$ and $\tilde{\Psi }=\Psi$. If we choose $\xi =-\mathcal{E}$ and ${\xi }^0=-(\mathcal{B}-{\mathcal{E}}^{\prime })$, then we have $\tilde{\mathcal{B}}=0$ and $\tilde{\mathcal{E}}=0$. This is the Newtonian gauge, in which $\tilde{\mathcal{A}}=\tilde{\Phi }$, and $\tilde{\psi }=\tilde{\Psi }$ become gauge-invariant. In terms of the gauge-invariant quantities the line element reduces to:

$$d{s}^2=a^2\left\{-\left(1+2\Phi \right)d{\eta }^2+\left(1-2\Psi \right){\delta }_{ij}d{x}^{i}d{x}^j\right\}.$$

(8)

The Bardeen potentials will be equal if the stress tensor of matter has no anisotropic perturbations.

3.2. Scalar Field Perturbations

Next we turn to the perturbation of a real scalar field (We may consider non-minimal coupling: adding a term proportional to the scalar curvature for a spatially-flat spacetime would amount to an additional term in the time-dependent frequency of the modes. This is discussed in e.g., Section 2.1.4 of [45]. If the ‘non-minimal’ coupling is conformal coupling, then under a conformal transformation of the field, in the conformal time, one gets an equation of the same form as a quantum field in flat space with time-dependent ‘effective mass’, as is studied in (21)).

$$S_{\varphi }=-\int d^{4}x\sqrt{-g}\left\{\frac{1}{2}{g}^{\mu \nu }{\partial }_{\mu }\varphi {\partial }_{\nu }\varphi +V\left(\varphi \right)\right\},$$

(9)

where V is its potential. We shall use a background field expansion $\varphi =\overline{\varphi }+\delta \varphi$ around a background field $\overline{\varphi }$ and focus on the perturbations $\delta \varphi$. We will assume that the unperturbed spacetime is homogeneous and isotropic, and that the background field $\overline{\varphi }\left(\eta \right)$ in this spacetime is a function of time only. The scalar field perturbations $\delta \varphi (\mathit{x},\eta )$ obey the wave equation:

$$\delta {\varphi }^{″}-{\partial }^2\delta \varphi +2\mathcal{H}\delta {\varphi }^{\prime }+a^2\frac{{\partial }^{2}V\left(\overline{\varphi }\right)}{\partial {\overline{\varphi }}^2}\delta \varphi =\left({\mathcal{A}}^{\prime }+3{\mathcal{D}}^{\prime }+{\partial }^2\mathcal{B}\right){\overline{\varphi }}^{\prime }-2a^2\mathcal{A}\frac{\partial V\left(\overline{\varphi }\right)}{\partial \overline{\varphi }},$$

(10)

where $\mathcal{H}=a^{\prime }/a$ is the Hubble expansion rate in conformal time, and ${\partial }^2$ denotes ${\sum }_i{\partial }_i^2$. We see that the scalar field perturbations couple only to scalar metric perturbations.

After some algebra, one can derive the Einstein equations for the scalar perturbations given by:

$$\begin{array}{cc}\hfill {\partial }^2\Phi -3\mathcal{H}{\Phi }^{\prime }-3{\mathcal{H}}^2\Phi & =4\pi a^2\left[\frac{1}{a^2}\left({\overline{\varphi }}^{\prime }\delta {\varphi }^{\prime }-{\overline{\varphi }}^{\prime 2}\Phi \right)+\frac{\partial V\left(\overline{\varphi }\right)}{\partial \overline{\varphi }}\delta \varphi \right],\hfill \end{array}$$

(11)

$$\begin{array}{cc}\hfill {\Phi }^{\prime }+\mathcal{H}\Phi & =4\pi {\overline{\varphi }}^{\prime }\delta \varphi ,\hfill \end{array}$$

(12)

$$\begin{array}{cc}\hfill {\Phi }^{″}+3\mathcal{H}{\Phi }^{\prime }+\left(2{\mathcal{H}}^{\prime }+{\mathcal{H}}^2\right)\Phi & =4\pi a^2\left[\frac{1}{a^2}\left({\overline{\varphi }}^{\prime }\delta {\varphi }^{\prime }-{\overline{\varphi }}^{\prime 2}\Phi \right)-\frac{\partial V\left(\overline{\varphi }\right)}{\partial \overline{\varphi }}\delta \varphi \right],\hfill \end{array}$$

(13)

with the background field $\overline{\varphi }$ satisfying:

$${\overline{\varphi }}^{″}+2\mathcal{H}{\overline{\varphi }}^{\prime }+a^2\frac{\partial V\left(\overline{\varphi }\right)}{\partial \overline{\varphi }}=0.$$

(14)

Recall that $\Phi$ is the Bardeen potential. Taking a suitable superposition of the metric perturbation and the scalar field perturbation, Mukhanov and Sasaki define the gauge-invariant variable u in the conformal Newtonian gauge as:

$$u=\frac{a^2}{a^{\prime }}{\overline{\varphi }}^{\prime }\left(\Phi +\frac{a^{\prime }}{a{\overline{\varphi }}^{\prime }}\delta \varphi \right)$$

(15)

and rewrite the Einstein equation in the conformal Newtonian gauge into the Mukhanov–Sasaki equation:

$$\begin{array}{cccccc}\hfill u^{″}-{\partial }^{2}u-\frac{z^{″}}{z}u& =0,\hfill & \hfill & \mathrm{with}\hspace{1em}\hspace{1em}\hspace{1em}\hfill & \hfill z& =\frac{a^2}{a^{\prime }}{\overline{\varphi }}^{\prime }.\hfill \end{array}$$

(16)

With the scalar field playing the role of the inflaton, in the slow-roll limit, the spacetime is nearly de Sitter and (16) will assume the form

$$u^{″}-{\partial }^{2}u-{\mathcal{H}}^2\left(2+5ϵ-3\eta \right)u=0,$$

(17)

where the slow-roll parameters $ϵ$ and $\eta$ are given by:

$$\begin{array}{cccc}\hfill ϵ& =\frac{m_{\mathrm{p}}^2}{2}\frac{V^{\prime 2}}{V^2},\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hfill & \hfill \eta & =m_{\mathrm{p}}^2\frac{V^{″}}{V},\hfill \end{array}$$

(18)

with the Planck mass $m_{\mathrm{p}}^2=1/8\pi$. Alternatively, we may rewrite (16) to obtain the wave equation of the curvature perturbation:

$${\mathcal{R}}^{″}-{\partial }^2\mathcal{R}+2\mathcal{H}\left(1+2ϵ-\eta \right){\mathcal{R}}^{\prime }=0,$$

(19)

in the slow-roll limit, where the curvature perturbation $\mathcal{R}$ is related to the Mukhanov-Sasaki variable u by $u=-z\mathcal{R}$. We observe both Equation (16) and (19) are in the form of an equation of motion for a parametric field $u(\mathit{x},\eta )$ or $\mathcal{R}(\mathit{x},\eta )$ in flat space. The gravitational influence is wrapped in the function of $z\left(\eta \right)$ which, according to (16), contains the effects of the scale factor or the inflaton potential $V\left(\overline{\varphi }\right)$. Thus, for our purpose here, it suffices to focus on this equation and discuss the dynamics of a generic parametrically-driven quantum scalar field in flat space as we shall do in the next sections.

4. Dynamics of a Parametrically-Driven Quantum Field

We consider a real scalar parametric field $\varphi (\mathit{x},t)$ in a flat spacetime, whose Lagrangian density is given by:

$$\mathcal{L}=-\frac{1}{2}\left[{\eta }^{\mu \nu }{\partial }_{\mu }\varphi {\partial }_{\nu }\varphi +m^2\left(t\right){\varphi }^2\right],$$

(20)

with the Minkowski metric ${\eta }^{\mu \nu }$. We use the convention with signature $(-1,+1,+1,+1)$. The time-dependent function $m^2\left(t\right)$ will effectively account for the parametric processes from the external agent. The corresponding equation of motion is:

$$-{\partial }_{\mu }\left({\eta }^{\mu \nu }{\partial }_{\nu }\varphi \right)+m^2\left(t\right)\varphi =0.$$

(21)

Expanding the field variable $\varphi$ in terms of the mode function ${\varphi }_{\mathit{k}}\left(t\right)$ of the form

$$\begin{array}{cccccc}\hfill \varphi (\mathit{x},t)& =\sum _{\mathit{k}}{e}^{+i\mathit{k}·\mathit{x}}{\varphi }_{\mathit{k}}\left(t\right),\hfill & \hfill & \mathrm{with}\hspace{1em}\hspace{1em}\hspace{1em}\hfill & \hfill {\varphi }_{+\mathit{k}}^{\ast }\left(t\right)& ={\varphi }_{-\mathit{k}}\left(t\right),\hfill \end{array}$$

(22)

we find the action given by:

$$\begin{array}{c}\hfill S=-\frac{1}{2}\int dt\sum _{\mathit{k}}\left[-{\dot{\varphi }}_{\mathit{k}}^{\phantom{\ast }}{\dot{\varphi }}_{\mathit{k}}^{\ast }+\left({\mathit{k}}^2+m^2\right){\varphi }_{\mathit{k}}^{\phantom{\ast }}{\varphi }_{\mathit{k}}^{\ast }\right].\end{array}$$

(23)

Since ${\dot{\varphi }}_{+\mathit{k}}^{\phantom{\ast }}{\dot{\varphi }}_{+\mathit{k}}^{\ast }={\dot{\varphi }}_{-\mathit{k}}^{\phantom{\ast }}{\dot{\varphi }}_{-\mathit{k}}^{\ast }$ and the summand in (23) takes the same form for $\pm \mathit{k}$, we can write the action only in terms of the $\mathit{k}>0$ modes:

$$S=\int dtL({\varphi }_{\mathit{k}},{\dot{\varphi }}_{\mathit{k}})=-\int dt\sum _{\mathit{k}>0}\left[-{\dot{\varphi }}_{\mathit{k}}^{\phantom{\ast }}{\dot{\varphi }}_{\mathit{k}}^{\ast }+\left({\mathit{k}}^2+m^2\right){\varphi }_{\mathit{k}}^{\phantom{\ast }}{\varphi }_{\mathit{k}}^{\ast }\right].$$

(24)

With the momentum conjugate to ${\varphi }_{\mathit{k}}$ defined by:

$${\pi }_{\mathit{k}}=\frac{\partial L}{\partial {\dot{\varphi }}_{\mathit{k}}}={\dot{\varphi }}_{\mathit{k}}^{\ast },$$

(25)

the Euler–Lagrange equation gives the equation of motion for ${\varphi }_{\mathit{k}}\left(t\right)$

$${\ddot{\varphi }}_{\mathit{k}}+\left({\mathit{k}}^2+m^2\right){\varphi }_{\mathit{k}}=0.$$

(26)

For the treatment of quantum fields we promote the field to be and operator. In the Heisenberg picture the expectation value will be computed with respect to the initial state. Then for each mode $\mathit{k}$ we have:

$$\begin{array}{cccc}\hfill {\widehat{\varphi }}_{+\mathit{k}}^{†}\left(t\right)& ={\widehat{\varphi }}_{-\mathit{k}}^{\phantom{†}}\left(t\right),\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hfill & \hfill {\widehat{\pi }}_{\mathit{k}}& ={\dot{\widehat{\varphi }}}_{\mathit{k}}^{†}={\dot{\widehat{\varphi }}}_{-\mathit{k}}^{\phantom{†}}.\hfill \end{array}$$

(27)

On the other hand, since the momentum $\widehat{\pi }$ conjugate to the field operator $\widehat{\varphi }$ is given by $\widehat{\pi }=\dot{\widehat{\varphi }}$, we have:

$$\begin{array}{c}\hfill \widehat{\pi }=\dot{\widehat{\varphi }}=\sum _{\mathit{k}}{\dot{\widehat{\varphi }}}_{\mathit{k}}{e}^{+i\mathit{k}·\mathit{x}}=\sum _{\mathit{k}}{\dot{\widehat{\varphi }}}_{-\mathit{k}}^{\phantom{†}}{e}^{-i\mathit{k}·\mathit{x}}=\sum _{\mathit{k}}{\dot{\widehat{\varphi }}}_{\mathit{k}}^{†}{e}^{-i\mathit{k}·\mathit{x}}=\sum _{\mathit{k}}{\widehat{\pi }}_{\mathit{k}}{e}^{-i\mathit{k}·\mathit{x}}.\end{array}$$

(28)

Thus, the Fourier expansion of the momentum operator $\widehat{\pi }$ takes a different form from that of $\widehat{\varphi }$ in (22). This is necessary to ensure the standard form of the equal-time canonical commutation relation:

$$\left[\widehat{\varphi }(\mathit{x},t),\widehat{\pi }({\mathit{x}}^{\prime },t)\right]=\sum _{\mathit{k},{\mathit{k}}^{\prime }}\left[{\widehat{\varphi }}_{\mathit{k}}\left(t\right),{\widehat{\pi }}_{{\mathit{k}}^{\prime }}\left(t\right)\right]e^{i\mathit{k}·\mathit{x}-i{\mathit{k}}^{\prime }·{\mathit{x}}^{\prime }}=\sum _{\mathit{k}}{e}^{i\mathit{k}·(\mathit{x}-{\mathit{x}}^{\prime })}=i{\delta }^{\left(3\right)}(\mathit{x}-{\mathit{x}}^{\prime }),$$

(29)

if we require:

$$\left[{\widehat{\varphi }}_{\mathit{k}}\left(t\right),{\widehat{\pi }}_{{\mathit{k}}^{\prime }}\left(t\right)\right]=i{\delta }_{\mathit{k},{\mathit{k}}^{\prime }}.$$

(30)

In the general case, (27) allows an expansion:

$$\begin{array}{cccccc}\hfill {\widehat{\varphi }}_{\mathit{k}}\left(t\right)& ={\widehat{a}}_{+\mathit{k}}^{\phantom{†}}{u}_{\omega }^{\phantom{\ast }}\left(t\right)+{\widehat{a}}_{-\mathit{k}}^{†}{u}_{\omega }^{\ast }\left(t\right),\hfill & \hfill & \mathrm{and}\hspace{1em}\hspace{1em}\hfill & \hfill {\widehat{\pi }}_{\mathit{k}}\left(t\right)& ={\widehat{a}}_{-\mathit{k}}^{\phantom{†}}{\dot{u}}_{\omega }^{\phantom{\ast }}\left(t\right)+{\widehat{a}}_{+\mathit{k}}^{†}{\dot{u}}_{\omega }^{\ast }\left(t\right).\hfill \end{array}$$

(31)

Here, $u_{\omega }\left(t\right)$ is a solution of (26), satisfying the positive frequency condition at a specified initial time with ${\omega }^2\left(t\right)={\mathit{k}}^2+m^2\left(t\right)$, and ${\widehat{a}}_{\mathit{k}}^{\phantom{†}}$, ${\widehat{a}}_{\mathit{k}}^{†}$ are the annihilation and creation operators associated with $u_{\omega }$ at the initial time, satisfying the usual commutation relations:

$$\begin{array}{cc}\hfill \left[{\widehat{a}}_{\mathit{k}}^{\phantom{†}},{\widehat{a}}_{{\mathit{k}}^{\prime }}^{†}\right]& ={\delta }_{\mathit{k}{\mathit{k}}^{\prime }},\hfill \end{array}$$

(32)

and zero otherwise. Then the commutation relation (30) is consistently obeyed due to the Wronskian condition of $u_{\omega }$:

$$u_{\omega }^{\phantom{\ast }}\left(t\right){\dot{u}}_{\omega }^{\ast }\left(t\right)-u_{\omega }^{\ast }\left(t\right){\dot{u}}_{\omega }^{\phantom{\ast }}\left(t\right)=i.$$

(33)

We may invert (31) to express ${\widehat{a}}_{\mathit{k}}$ in terms of ${\widehat{\varphi }}_{\mathit{k}}$ and ${\widehat{\pi }}_{\mathit{k}}$,

$${\widehat{a}}_{\mathit{k}}=-i\left({\dot{u}}_{\omega }^{\ast }{\widehat{\varphi }}_{+\mathit{k}}-u_{\omega }^{\ast }{\widehat{\pi }}_{-\mathit{k}}^{\phantom{\ast }}\right),$$

(34)

and confirm the commutation relations (32) for the creation and the annihilation operators.

From the expansions (31), we may compute the corresponding correlation or covariance matrix elements among various modes [107]:

$$\begin{array}{cc}\hfill \frac{1}{2}\langle \left\{{\widehat{\varphi }}_{\mathit{k}},{\widehat{\varphi }}_{{\mathit{k}}^{\prime }}\right\}\rangle & =\frac{1}{2}\langle \left\{{\widehat{a}}_{+\mathit{k}}^{\phantom{†}},{\widehat{a}}_{+{\mathit{k}}^{\prime }}^{\phantom{†}}\right\}\rangle u_{\omega }^{\phantom{\ast }}{u}_{{\omega }^{\prime }}^{\phantom{\ast }}+\frac{1}{2}\langle \left\{{\widehat{a}}_{-\mathit{k}}^{†},{\widehat{a}}_{+{\mathit{k}}^{\prime }}^{\phantom{†}}\right\}\rangle u_{\omega }^{\ast }{u}_{{\omega }^{\prime }}^{\phantom{\ast }}+\frac{1}{2}\langle \left\{{\widehat{a}}_{+\mathit{k}}^{\phantom{†}},{\widehat{a}}_{-{\mathit{k}}^{\prime }}^{†}\right\}\rangle u_{\omega }^{\phantom{\ast }}{u}_{{\omega }^{\prime }}^{\ast }\hfill \\ \hfill & +\frac{1}{2}\langle \left\{{\widehat{a}}_{-\mathit{k}}^{†},{\widehat{a}}_{-{\mathit{k}}^{\prime }}^{†}\right\}\rangle u_{\omega }^{\ast }{u}_{{\omega }^{\prime }}^{\ast },\hfill \end{array}$$

(35)

$$\begin{array}{cc}\hfill \frac{1}{2}\langle \left\{{\widehat{\varphi }}_{\mathit{k}}^{\phantom{†}},{\widehat{\varphi }}_{{\mathit{k}}^{\prime }}^{†}\right\}\rangle & =\frac{1}{2}\langle \left\{{\widehat{a}}_{+\mathit{k}}^{\phantom{†}},{\widehat{a}}_{+{\mathit{k}}^{\prime }}^{†}\right\}\rangle u_{\omega }^{\phantom{\ast }}{u}_{{\omega }^{\prime }}^{\ast }+\frac{1}{2}\langle \left\{{\widehat{a}}_{-\mathit{k}}^{†},{\widehat{a}}_{+{\mathit{k}}^{\prime }}^{†}\right\}\rangle u_{\omega }^{\ast }{u}_{{\omega }^{\prime }}^{\ast }+\frac{1}{2}\langle \left\{{\widehat{a}}_{+\mathit{k}}^{\phantom{†}},{\widehat{a}}_{-{\mathit{k}}^{\prime }}^{\phantom{†}}\right\}\rangle u_{\omega }^{\phantom{\ast }}{u}_{{\omega }^{\prime }}^{\phantom{\ast }}\hfill \\ \hfill & +\frac{1}{2}\langle \left\{{\widehat{a}}_{-\mathit{k}}^{†},{\widehat{a}}_{-{\mathit{k}}^{\prime }}^{\phantom{†}}\right\}\rangle u_{\omega }^{\ast }{u}_{{\omega }^{\prime }}^{\phantom{\ast }},\hfill \end{array}$$

(36)

$$\begin{array}{cc}\hfill \frac{1}{2}\langle \left\{{\widehat{\varphi }}_{\mathit{k}}^{\phantom{†}},{\widehat{\pi }}_{{\mathit{k}}^{\prime }}^{\phantom{†}}\right\}\rangle & =\frac{1}{2}\langle \left\{{\widehat{a}}_{+\mathit{k}}^{\phantom{†}},{\widehat{a}}_{+{\mathit{k}}^{\prime }}^{†}\right\}\rangle u_{\omega }^{\phantom{\ast }}{\dot{u}}_{{\omega }^{\prime }}^{\ast }+\frac{1}{2}\langle \left\{{\widehat{a}}_{-\mathit{k}}^{†},{\widehat{a}}_{+{\mathit{k}}^{\prime }}^{†}\right\}\rangle u_{\omega }^{\ast }{\dot{u}}_{{\omega }^{\prime }}^{\ast }+\frac{1}{2}\langle \left\{{\widehat{a}}_{+\mathit{k}}^{\phantom{†}},{\widehat{a}}_{-{\mathit{k}}^{\prime }}^{\phantom{†}}\right\}\rangle u_{\omega }^{\phantom{\ast }}{\dot{u}}_{{\omega }^{\prime }}^{\phantom{\ast }}\hfill \\ \hfill & +\frac{1}{2}\langle \left\{{\widehat{a}}_{-\mathit{k}}^{†},{\widehat{a}}_{-{\mathit{k}}^{\prime }}^{\phantom{†}}\right\}\rangle u_{\omega }^{\ast }{\dot{u}}_{{\omega }^{\prime }}^{\phantom{\ast }}.\hfill \end{array}$$

(37)

The results for $\langle \left\{{\widehat{\pi }}_{\mathit{k}},{\widehat{\pi }}_{{\mathit{k}}^{\prime }}\right\}\rangle$ will be obtained by replacing $u_{\omega }$ in (35) by ${\dot{u}}_{\omega }$. These elements are extremely useful in constructing physical observables of a Gaussian system in that, for each mode, we can simplify calculations involving an infinite-dimensional density matrix operator into those based a $2\times 2$ covariance matrix.

If the state used to evaluate the expectation value is a stationary state, that is, $\langle {\widehat{a}}_{\mathit{k}}^2\rangle =0$, then (35)–(37) reduce to:

$$\begin{array}{cc}\hfill \frac{1}{2}\langle \left\{{\widehat{\varphi }}_{\mathit{k}},{\widehat{\varphi }}_{{\mathit{k}}^{\prime }}\right\}\rangle & ={\delta }_{+\mathit{k},-{\mathit{k}}^{\prime }}\left[\left(\langle {\widehat{N}}_{+\mathit{k}}\rangle +\frac{1}{2}\right)+\left(\langle {\widehat{N}}_{-\mathit{k}}\rangle +\frac{1}{2}\right)\right]u_{\omega }^{\phantom{\ast }}{u}_{\omega }^{\ast },\hfill \end{array}$$

(38)

$$\begin{array}{cc}\hfill \frac{1}{2}\langle \left\{{\widehat{\varphi }}_{\mathit{k}}^{\phantom{†}},{\widehat{\varphi }}_{{\mathit{k}}^{\prime }}^{†}\right\}\rangle & ={\delta }_{+\mathit{k},+{\mathit{k}}^{\prime }}\left[\left(\langle {\widehat{N}}_{+\mathit{k}}\rangle +\frac{1}{2}\right)+\left(\langle {\widehat{N}}_{-\mathit{k}}\rangle +\frac{1}{2}\right)\right]u_{\omega }^{\phantom{\ast }}{u}_{\omega }^{\ast },\hfill \end{array}$$

(39)

$$\begin{array}{cc}\hfill \frac{1}{2}\langle \left\{{\widehat{\varphi }}_{\mathit{k}}^{\phantom{†}},{\widehat{\pi }}_{{\mathit{k}}^{\prime }}^{\phantom{†}}\right\}\rangle & ={\delta }_{+\mathit{k},+{\mathit{k}}^{\prime }}\left[\left(\langle {\widehat{N}}_{+\mathit{k}}\rangle +\frac{1}{2}\right)u_{\omega }^{\phantom{\ast }}{\dot{u}}_{\omega }^{\ast }+\left(\langle {\widehat{N}}_{-\mathit{k}}\rangle +\frac{1}{2}\right){\dot{u}}_{\omega }^{\phantom{\ast }}{u}_{\omega }^{\ast }\right],\hfill \end{array}$$

(40)

where ${\widehat{N}}_{\mathit{k}}={\widehat{a}}_{\mathit{k}}^{†}{\widehat{a}}_{\mathit{k}}^{\phantom{†}}$ is the number operator of mode $\mathit{k}$. Here we see that both $\pm \mathit{k}$ modes contribute. Equations (38)–(40) imply that

$$\begin{array}{cccc}\hfill \langle {\widehat{\varphi }}_{\pm \mathit{k}}^2\left(t\right)\rangle & =0,\hfill & \hfill \langle {\widehat{\varphi }}_{\pm \mathit{k}}^{\phantom{†}}\left(t\right){\widehat{\pi }}_{\pm \mathit{k}}^{†}\left(t\right)\rangle & =0.\hfill \end{array}$$

(41)

These ingredients will allow us to construct the density matrix operator of the field, which will be used to compute the expectation values or observables of the field.

5. Nonequilibrium Evolution of the Density Matrix

We now proceed to construct the density matrix operator of the quantum parametric scalar field using the formalisms developed in [108]. We begin by writing the field mode operators in the form of harmonic oscillators.

As the operators ${\widehat{\varphi }}_{\mathit{k}}$, ${\widehat{\pi }}_{\mathit{k}}$ are not Hermitian, we decompose them into the real and the imaginary parts by:

$$\begin{array}{cc}\hfill {\widehat{\varphi }}_{\mathit{k}}& =\frac{1}{\sqrt{2}}\left({\widehat{\varphi }}_{\mathit{k}}^{\left(1\right)}+i{\widehat{\varphi }}_{\mathit{k}}^{\left(2\right)}\right).\hfill \end{array}$$

(42)

Hence the action (24) becomes:

$$\begin{array}{c}\hfill S=\frac{1}{2}\int dt\sum _{\mathit{k}>0}\left[\left({\dot{\varphi }}_{\mathit{k}}^{\left(1\right)2}+{\dot{\varphi }}_{\mathit{k}}^{\left(2\right)2}\right)-{\omega }^2\left(t\right)\left({\varphi }_{\mathit{k}}^{\left(1\right)2}+{\varphi }_{\mathit{k}}^{\left(2\right)2}\right)\right].\end{array}$$

(43)

with the real and the imaginary parts acting as uncoupled parametric oscillators and ${\omega }^2\left(t\right)={\mathit{k}}^2+m^2\left(t\right)$.

The corresponding conjugate momenta are given by:

$$\begin{array}{cccc}\hfill p_{\mathit{k}}^{\left(1\right)}& =\frac{\partial L}{\partial {\dot{\varphi }}_{\mathit{k}}^{\left(1\right)}}={\dot{\varphi }}_{\mathit{k}}^{\left(1\right)},\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hfill & \hfill p_{\mathit{k}}^{\left(2\right)}& =\frac{\partial L}{\partial {\dot{\varphi }}_{\mathit{k}}^{\left(2\right)}}={\dot{\varphi }}_{\mathit{k}}^{\left(2\right)},\hfill \end{array}$$

(44)

which implies

$${\widehat{\pi }}_{\mathit{k}}={\dot{\widehat{\varphi }}}_{\mathit{k}}^{†}=\frac{1}{\sqrt{2}}\left(p_{\mathit{k}}^{\left(1\right)}-i{p}_{\mathit{k}}^{\left(2\right)}\right).$$

(45)

Note the sign before the imaginary part of ${\widehat{\pi }}_{\mathit{k}}$. This then ensures that Equation (30) is obeyed:

$$\left[{\widehat{\varphi }}_{\mathit{k}},{\widehat{\pi }}_{{\mathit{k}}^{\prime }}\right]=\frac{1}{2}\left[{\widehat{\varphi }}_{\mathit{k}}^{\left(1\right)}+i{\widehat{\varphi }}_{\mathit{k}}^{\left(2\right)},p_{\mathit{k}}^{\left(1\right)}-i{p}_{\mathit{k}}^{\left(2\right)}\right]=i{\delta }_{\mathit{k},{\mathit{k}}^{\prime }},$$

(46)

if $\left[{\widehat{\varphi }}_{\mathit{k}}^{\left(i\right)},p_{{\mathit{k}}^{\prime }}^{\left(j\right)}\right]=i{\delta }_{ij}{\delta }_{\mathit{k},{\mathit{k}}^{\prime }}$. Furthermore, since ${\widehat{\varphi }}_{+\mathit{k}}^{†}={\widehat{\varphi }}_{-\mathit{k}}^{\phantom{†}}$, we find:

$$\begin{array}{cccccc}\hfill {\widehat{\varphi }}_{+\mathit{k}}^{\left(1\right)}& =+{\widehat{\varphi }}_{-\mathit{k}}^{\left(1\right)},\hfill & \hfill & \mathrm{and}\hspace{1em}\hspace{1em}\hspace{1em}\hfill & \hfill {\widehat{\varphi }}_{+\mathit{k}}^{\left(2\right)}& =-{\widehat{\varphi }}_{-\mathit{k}}^{\left(2\right)}.\hfill \end{array}$$

(47)

That is, ${\widehat{\varphi }}_{\mathit{k}}^{\left(1\right)}$ is an even function of $\mathit{k}$, but ${\widehat{\varphi }}_{\mathit{k}}^{\left(2\right)}$ is an odd function of $\mathit{k}$. Here one may already expect that neither the operator ${\widehat{\varphi }}_{\mathit{k}}$ nor ${\widehat{\varphi }}_{\mathit{k}}^{\left(i\right)}$ can offer an unambiguous specification of the $\pm \mathit{k}$ modes.

To construct the density matrix operator, we need the covariance matrix elements of ${\widehat{\varphi }}_{\mathit{k}}^{\left(i\right)}$. This can be done by first inverting the decompositions (42) and (45)

$$\begin{array}{cccc}\hfill {\widehat{\varphi }}_{\mathit{k}}^{\left(1\right)}& =\frac{1}{\sqrt{2}}\left({\widehat{\varphi }}_{\mathit{k}}^{\phantom{†}}+{\widehat{\varphi }}_{\mathit{k}}^{†}\right),\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hfill & \hfill {\widehat{\varphi }}_{\mathit{k}}^{\left(2\right)}& =\frac{1}{i\sqrt{2}}\left({\widehat{\varphi }}_{\mathit{k}}^{\phantom{†}}-{\widehat{\varphi }}_{\mathit{k}}^{†}\right),\hfill \end{array}$$

(48)

$$\begin{array}{cccc}\hfill {\widehat{p}}_{\mathit{k}}^{\left(1\right)}& =\frac{1}{\sqrt{2}}\left({\widehat{\pi }}_{\mathit{k}}^{†}+{\widehat{\pi }}_{\mathit{k}}^{\phantom{†}}\right),\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hfill & \hfill {\widehat{p}}_{\mathit{k}}^{\left(2\right)}& =\frac{1}{i\sqrt{2}}\left({\widehat{\pi }}_{\mathit{k}}^{†}-{\widehat{\varphi }}_{\mathit{k}}^{\phantom{†}}\right).\hfill \end{array}$$

(49)

From (38)–(40), we obtain the covariance matrix elements of ${\widehat{\varphi }}_{\mathit{k}}^{\left(i\right)}$ with respect to the initial stationary state given by:

$$\begin{array}{cc}\hfill \langle {\widehat{\varphi }}_{\mathit{k}}^{\left(1\right)2}\rangle & =\left[\left(\langle {\widehat{N}}_{+\mathit{k}}\rangle +\frac{1}{2}\right)+\left(\langle {\widehat{N}}_{-\mathit{k}}\rangle +\frac{1}{2}\right)\right]u_{\omega }^{\phantom{\ast }}{u}_{\omega }^{\ast }=\langle {\widehat{\varphi }}_{-\mathit{k}}^{\left(1\right)2}\rangle ,\hfill \end{array}$$

(50)

$$\begin{array}{cc}\hfill \frac{1}{2}\langle \left\{{\widehat{\varphi }}_{\mathit{k}}^{\left(1\right)},{\widehat{p}}_{\mathit{k}}^{\left(1\right)}\right\}\rangle & =\frac{1}{2}\left[\left(\langle {\widehat{N}}_{+\mathit{k}}\rangle +\frac{1}{2}\right)+\left(\langle {\widehat{N}}_{-\mathit{k}}\rangle +\frac{1}{2}\right)\right]\left(u_{\omega }^{\phantom{\ast }}{\dot{u}}_{\omega }^{\ast }+{\dot{u}}_{\omega }^{\phantom{\ast }}{u}_{\omega }^{\ast }\right)=\frac{1}{2}\langle \left\{{\widehat{\varphi }}_{-\mathit{k}}^{\left(1\right)},{\widehat{p}}_{-\mathit{k}}^{\left(1\right)}\right\}\rangle ,\hfill \end{array}$$

(51)

and then

$$\begin{array}{cccc}\hfill \langle {\widehat{\varphi }}_{\mathit{k}}^{\left(2\right)2}\rangle & =\langle {\widehat{\varphi }}_{\mathit{k}}^{\left(1\right)2}\rangle ,\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hfill & \hfill \frac{1}{2}\langle \left\{{\widehat{\varphi }}_{\mathit{k}}^{\left(2\right)},{\widehat{p}}_{\mathit{k}}^{\left(2\right)}\right\}\rangle & =\frac{1}{2}\langle \left\{{\widehat{\varphi }}_{\mathit{k}}^{\left(1\right)},{\widehat{p}}_{\mathit{k}}^{\left(1\right)}\right\}\rangle .\hfill \end{array}$$

(52)

Finally, we find:

$$\begin{array}{cccc}\hfill \frac{1}{2}\langle \left\{{\widehat{\varphi }}_{\mathit{k}}^{\left(1\right)},{\widehat{\varphi }}_{\mathit{k}}^{\left(2\right)}\right\}\rangle & =0,\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hfill & \hfill \frac{1}{2}\langle \left\{{\widehat{\varphi }}_{\mathit{k}}^{\left(1\right)},{\widehat{p}}_{\mathit{k}}^{\left(2\right)}\right\}\rangle & =0,\hfill \end{array}$$

(53)

for the initial stationary state.

Having seen how the field modes are expressed in terms of harmonic oscillators, we can now construct the density matrix operators of the field modes.

5.1. Density Matrix Operator

Let us quickly review the general form of the density matrix operator of a Gaussian state [108]. Consider the one-mode case, and let $(\widehat{q},\widehat{p})$ be the conjugated canonical operators of a Gaussian system. The density operator of the Gaussian state in general takes the form:

$$\begin{array}{c}\hfill \widehat{\rho }(\widehat{q},\widehat{p})=\frac{2}{\sqrt{e^{2\chi }-1}}exp\left\{-e^{-\chi }{cosh}^{-1}(coth\chi )\left[\mathsf{a}{\widehat{q}}^2+\mathsf{b}{\widehat{p}}^2-\mathsf{c}\left\{\widehat{q},\widehat{p}\right\}\right]\right\},\end{array}$$

(54)

where the coefficients can be expressed in terms of the corresponding time-dependent covariance matrix elements:

$$\begin{array}{cccccccccc}\hfill \mathsf{a}& =\langle {\widehat{p}}^2\rangle ,\hfill & \hfill \mathsf{b}& =\langle {\widehat{q}}^2\rangle ,\hfill & \hfill \mathsf{c}& =\frac{1}{2}\langle \left\{\widehat{q},\widehat{p}\right\}\rangle ,\hfill & \hfill & \mathrm{and}\hspace{1em}\hfill & \hfill e^{2\chi }& =4\left(\mathsf{a}\mathsf{b}-{\mathsf{c}}^2\right).\hfill \end{array}$$

(55)

By means of the creation and annihilation operators, the density matrix operator (54) becomes:

$$\begin{array}{c}\hfill \widehat{\rho }(\widehat{a},{\widehat{a}}^{†})=\frac{2}{\sqrt{e^{2\chi }-1}}exp\left\{-2\mathsf{\Xi }{e}^{-\chi }{cosh}^{-1}(coth\chi )\left[\kappa {\widehat{a}}^2+{\kappa }^{\ast }{\widehat{a}}^{†2}+\frac{\tau }{2}\left\{\widehat{a},{\widehat{a}}^{†}\right\}\right]\right\},\end{array}$$

(56)

where

$$\begin{array}{cccc}\hfill \kappa & =\frac{1}{2}sinh2\eta e^{-i\psi },\hspace{1em}\hspace{1em}\hspace{1em}\hfill & \hfill \tau & =cosh2\eta ,\hfill \end{array}$$

(57)

$$\begin{array}{cccc}\hfill \mathsf{\Xi }& =coth\frac{\beta }{2},\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hfill & \hfill e^{2\chi }& ={\mathsf{\Xi }}^2\left({\tau }^2-4{\left|\kappa \right|}^2\right)={\mathsf{\Xi }}^2={coth}^2\frac{\beta }{2},\hfill \end{array}$$

(58)

and $\beta$, $\zeta =\eta e^{i\psi }$ are inverse-temperature-like and squeeze-like parameters. Hence we have:

$$\begin{array}{cccc}\hfill 2\mathsf{\Xi }{e}^{-\chi }{cosh}^{-1}(coth\chi )& =2\beta ,\hfill & \hfill \frac{2}{\sqrt{e^{2\chi }-1}}& =2sinh\frac{\beta }{2}.\hfill \end{array}$$

(59)

The second expression in (59) is essentially the inverse of the nonequilibrium partition function [108] associated with the density matrix operator (54) or (56). Here, we have assumed that the first moments vanish for the state considered; otherwise, we may let $\widehat{a}\to \widehat{a}-\langle \widehat{a}\rangle$ or $\widehat{q}\to \widehat{q}-\langle \widehat{q}\rangle$ for the more general cases.

Observe that the factor:

$$\kappa {\widehat{a}}^2+{\kappa }^{\ast }{\widehat{a}}^{†2}+\frac{\tau }{2}\left\{\widehat{a},{\widehat{a}}^{†}\right\}$$

(60)

can be ‘rotated’ into the form:

$$\frac{1}{2}\left\{\widehat{b},{\widehat{b}}^{†}\right\}.$$

(61)

Suppose such a $\widehat{b}$ is related to $\widehat{a}$ and ${\widehat{a}}^{†}$ by $\widehat{b}=\mu \widehat{a}+\nu {\widehat{a}}^{†}$, with ${\left|\mu \right|}^2-{\left|\nu \right|}^2=1$, we have:

$$\begin{array}{cc}\hfill \frac{1}{2}\left\{\widehat{b},{\widehat{b}}^{†}\right\}& =\frac{\mu {\nu }^{\ast }}{2}\left\{\widehat{a},\widehat{a}\right\}+\frac{{\mu }^{\ast }\nu }{2}\left\{{\widehat{a}}^{†},{\widehat{a}}^{†}\right\}+\frac{{\left|\mu \right|}^2+{\left|\nu \right|}^2}{2}\left\{\widehat{a},{\widehat{a}}^{†}\right\}.\hfill \end{array}$$

(62)

Comparing with (60), we obtain:

$$\begin{array}{cccc}\hfill \mu {\nu }^{\ast }& =\kappa =\frac{1}{2}sinh2\eta e^{-i\psi },\hspace{1em}\hspace{1em}\hfill & \hfill {\left|\mu \right|}^2+{\left|\nu \right|}^2& =2\left({\left|\nu \right|}^2+\frac{1}{2}\right)=\tau =cosh2\eta .\hfill \end{array}$$

(63)

If we assume $\mu \in \mathbb{R}$, then a useful choice of $(\mu ,\nu )$ is

$$\begin{array}{cccc}\hfill \mu & =cosh\eta ,\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hfill & \hfill \nu & =sinh\eta e^{-i\psi },\hfill \end{array}$$

(64)

that is, a squeeze transformation. A complex $\mu$ will correspond to an additional rotation, so that (64) will contain the phase factor that accounts for the rotation.

Since from (43), the real part and the imaginary parts of ${\widehat{\varphi }}_{\mathit{k}}$ behave like independent harmonic oscillators, we expect that the corresponding density operator for each $\mathit{k}>0$ mode will be:

$$\begin{array}{cc}\hfill \widehat{\rho }({\widehat{\varphi }}_{\mathit{k}}^{\left(1\right)},{\widehat{\varphi }}_{\mathit{k}}^{\left(2\right)})& ={\left(\frac{2}{\sqrt{e^{2{\chi }_{\mathit{k}}}-1}}\right)}^{2}exp\left\{-e^{-{\chi }_{\mathit{k}}}{cosh}^{-1}(coth{\chi }_{\mathit{k}})\left[{\mathsf{a}}_{\mathit{k}}^{\left(1\right)}{\widehat{\varphi }}_{\mathit{k}}^{\left(1\right)2}+{\mathsf{b}}_{\mathit{k}}^{\left(1\right)}{\widehat{p}}_{\mathit{k}}^{\left(1\right)2}\right.\right.\hfill \\ \hfill & -\left.\left.{\mathsf{c}}_{\mathit{k}}^{\left(1\right)}\left\{{\widehat{\varphi }}_{\mathit{k}}^{\left(1\right)},{\widehat{p}}_{\mathit{k}}^{\left(1\right)}\right\}+{\mathsf{a}}_{\mathit{k}}^{\left(2\right)}{\widehat{\varphi }}_{\mathit{k}}^{\left(2\right)2}+{\mathsf{b}}_{\mathit{k}}^{\left(2\right)}{\widehat{p}}_{\mathit{k}}^{\left(2\right)2}-{\mathsf{c}}_{\mathit{k}}^{\left(2\right)}\left\{{\widehat{\varphi }}_{\mathit{k}}^{\left(2\right)},{\widehat{p}}_{\mathit{k}}^{\left(2\right)}\right\}\right]\right\},\hfill \end{array}$$

(65)

where we have assumed that ${\chi }_{\mathit{k}}$ is the same for the real and the imaginary components, to be justified later. Equations (50)–(53) imply that:

$$\begin{array}{cccccc}\hfill {\mathsf{a}}_{+\mathit{k}}^{\left(2\right)}& ={\mathsf{a}}_{+\mathit{k}}^{\left(1\right)}={\mathsf{a}}_{-\mathit{k}}^{\left(1\right)}={\mathsf{a}}_{-\mathit{k}}^{\left(2\right)},\hfill & \hfill {\mathsf{b}}_{+\mathit{k}}^{\left(2\right)}& ={\mathsf{b}}_{+\mathit{k}}^{\left(1\right)}={\mathsf{b}}_{-\mathit{k}}^{\left(1\right)}={\mathsf{b}}_{-\mathit{k}}^{\left(2\right)},\hfill & \hfill {\mathsf{c}}_{+\mathit{k}}^{\left(2\right)}& ={\mathsf{c}}_{+\mathit{k}}^{\left(1\right)}={\mathsf{c}}_{-\mathit{k}}^{\left(1\right)}={\mathsf{c}}_{-\mathit{k}}^{\left(2\right)},\hfill \end{array}$$

(66)

so let us assign

$$\begin{array}{cccccc}\hfill {\mathsf{a}}_{\mathit{k}}& \equiv {\mathsf{a}}_{\mathit{k}}^{\left(1\right)}={\mathsf{a}}_{\mathit{k}}^{\left(2\right)},\hspace{1em}\hspace{1em}\hspace{1em}\hfill & \hfill {\mathsf{b}}_{\mathit{k}}& \equiv {\mathsf{b}}_{\mathit{k}}^{\left(1\right)}={\mathsf{b}}_{\mathit{k}}^{\left(2\right)},\hspace{1em}\hspace{1em}\hspace{1em}\hfill & \hfill {\mathsf{c}}_{\mathit{k}}& \equiv {\mathsf{c}}_{\mathit{k}}^{\left(1\right)}={\mathsf{c}}_{\mathit{k}}^{\left(2\right)},\hfill \end{array}$$

(67)

and write Equation (65) as:

$$\begin{array}{cc}\hfill \widehat{\rho }({\widehat{\varphi }}_{\mathit{k}}^{\left(1\right)},{\widehat{\varphi }}_{\mathit{k}}^{\left(2\right)})& ={\left(\frac{2}{\sqrt{e^{2{\chi }_{\mathit{k}}}-1}}\right)}^2\hfill \\ \hfill & \times exp\left\{-e^{-{\chi }_{\mathit{k}}}{cosh}^{-1}(coth{\chi }_{\mathit{k}})\left[{\mathsf{a}}_{\mathit{k}}\left({\widehat{\varphi }}_{\mathit{k}}^{\left(1\right)2}+{\widehat{\varphi }}_{\mathit{k}}^{\left(2\right)2}\right)+{\mathsf{b}}_{\mathit{k}}\left({\widehat{p}}_{\mathit{k}}^{\left(1\right)2}+{\widehat{p}}_{\mathit{k}}^{\left(2\right)2}\right)\right.\right.\hfill \\ \hfill & -\left.\left.{\mathsf{c}}_{\mathit{k}}\left(\left\{{\widehat{\varphi }}_{\mathit{k}}^{\left(1\right)},{\widehat{p}}_{\mathit{k}}^{\left(1\right)}\right\}+\left\{{\widehat{\varphi }}_{\mathit{k}}^{\left(2\right)},{\widehat{p}}_{\mathit{k}}^{\left(2\right)}\right\}\right)\right]\right\},\hfill \end{array}$$

(68)

with $\mathit{k}>0$.

When we take into account all modes, the density matrix operator of the field modes $\widehat{\rho }(\left\{{\widehat{\varphi }}_{\mathit{k}}\right\},\left\{{\widehat{\pi }}_{\mathit{k}}\right\})={\prod }_{\mathit{k}>0}\widehat{\rho }({\widehat{\varphi }}_{\mathit{k}}^{\left(1\right)},{\widehat{\varphi }}_{\mathit{k}}^{\left(2\right)})$ corresponding to (43) can be found to be:

$$\begin{array}{cc}\hfill \widehat{\rho }(\left\{{\widehat{\varphi }}_{\mathit{k}}\right\},\left\{{\widehat{\pi }}_{\mathit{k}}\right\})& =\left(\prod _{\mathit{k}}\frac{2}{\sqrt{e^{2{\chi }_{\mathit{k}}}-1}}\right)exp\left\{-\sum _{\mathit{k}}{e}^{-{\chi }_{\mathit{k}}}{cosh}^{-1}(coth{\chi }_{\mathit{k}})\left[{\mathsf{a}}_{\mathit{k}}{\widehat{\varphi }}_{+\mathit{k}}^{\phantom{†}}{\widehat{\varphi }}_{+\mathit{k}}^{†}\right.\right.\hfill \\ \hfill & +\left.\left.{\mathsf{b}}_{\mathit{k}}{\widehat{\pi }}_{+\mathit{k}}^{\phantom{†}}{\widehat{\pi }}_{+\mathit{k}}^{†}-{\mathsf{c}}_{\mathit{k}}\left\{{\widehat{\varphi }}_{\mathit{k}},{\widehat{\pi }}_{\mathit{k}}\right\}\right]\right\},\hfill \end{array}$$

(69)

where we have assumed that the parameter ${\chi }_{+\mathit{k}}={\chi }_{-\mathit{k}}$, to be justified later. What remains is to express the coefficients ${\mathsf{b}}_{\mathit{k}}$, ${\mathsf{a}}_{\mathit{k}}$ and ${\mathsf{c}}_{\mathit{k}}$ by the original field mode operators ${\widehat{\varphi }}_{\mathit{k}}$ and ${\widehat{\pi }}_{\mathit{k}}$. Since ${\mathsf{b}}_{\mathit{k}}=\langle {\widehat{\varphi }}_{\mathit{k}}^{\left(1\right)2}\rangle =\langle {\widehat{\varphi }}_{\mathit{k}}^{\left(2\right)2}\rangle$, we find:

$$\begin{array}{c}\hfill {\mathsf{b}}_{\mathit{k}}=\frac{1}{2}\left(\langle {\widehat{\varphi }}_{\mathit{k}}^{\left(1\right)2}\rangle +\langle {\widehat{\varphi }}_{\mathit{k}}^{\left(2\right)2}\rangle \right)=\frac{1}{2}\left(\langle {\widehat{\varphi }}_{\mathit{k}}^{\phantom{†}}{\widehat{\varphi }}_{\mathit{k}}^{†}\rangle +\langle {\widehat{\varphi }}_{-\mathit{k}}^{\phantom{†}}{\widehat{\varphi }}_{-\mathit{k}}^{†}\rangle \right)=\frac{1}{2}\langle \left\{{\widehat{\varphi }}_{\mathit{k}}^{\phantom{†}},{\widehat{\varphi }}_{\mathit{k}}^{†}\right\}\rangle .\end{array}$$

(70)

Likewise, from:

$$\begin{array}{c}\hfill {\mathsf{c}}_{\mathit{k}}=\frac{1}{2}\langle \left\{{\widehat{\varphi }}_{\mathit{k}}^{\left(1\right)},{\widehat{p}}_{\mathit{k}}^{\left(1\right)}\right\}\rangle =\frac{1}{2}\langle \left\{{\widehat{\varphi }}_{\mathit{k}}^{\left(2\right)},{\widehat{p}}_{\mathit{k}}^{\left(2\right)}\right\}\rangle ,\end{array}$$

(71)

we obtain

$$\begin{array}{c}\hfill {\mathsf{c}}_{\mathit{k}}=\frac{1}{2}\left[\frac{1}{2}\langle \left\{{\widehat{\varphi }}_{\mathit{k}},{\widehat{\pi }}_{\mathit{k}}\right\}\rangle +\frac{1}{2}\langle \left\{{\widehat{\varphi }}_{\mathit{k}}^{†},{\widehat{\pi }}_{\mathit{k}}^{†}\right\}\rangle \right]=\frac{1}{2}\left[\frac{1}{2}\langle \left\{{\widehat{\varphi }}_{\mathit{k}},{\widehat{\pi }}_{\mathit{k}}\right\}\rangle +\frac{1}{2}\langle \left\{{\widehat{\varphi }}_{-\mathit{k}},{\widehat{\pi }}_{-\mathit{k}}\right\}\rangle \right].\end{array}$$

(72)

At the first sight, the result for ${\mathsf{c}}_{\mathit{k}}$ is not quite satisfactory because we cannot clearly separate contributions from the $\pm \mathit{k}$ modes. However, later we will see this is the form needed to construct the Wigner function of the whole field, instead of for each mode. This shows explicitly that there is an inseparable correlation between the $\pm \mathit{k}$ modes.

In terms of the mode expansion (31), the elements ${\mathsf{b}}_{\mathit{k}}$, ${\mathsf{a}}_{\mathit{k}}$ and ${\mathsf{c}}_{\mathit{k}}$ take the forms:

$$\begin{array}{cc}\hfill {\mathsf{b}}_{\mathit{k}}& =\left[\left(\langle {\widehat{N}}_{+\mathit{k}}\rangle +\frac{1}{2}\right)+\left(\langle {\widehat{N}}_{-\mathit{k}}\rangle +\frac{1}{2}\right)\right]u_{\omega }^{\phantom{\ast }}{u}_{\omega }^{\ast },\hfill \end{array}$$

(73)

$$\begin{array}{cc}\hfill {\mathsf{a}}_{\mathit{k}}& =\left[\left(\langle {\widehat{N}}_{+\mathit{k}}\rangle +\frac{1}{2}\right)+\left(\langle {\widehat{N}}_{-\mathit{k}}\rangle +\frac{1}{2}\right)\right]{\dot{u}}_{\omega }^{\phantom{\ast }}{\dot{u}}_{\omega }^{\ast },\hfill \end{array}$$

(74)

$$\begin{array}{cc}\hfill {\mathsf{c}}_{\mathit{k}}& =\frac{1}{2}\left[\left(\langle {\widehat{N}}_{+\mathit{k}}\rangle +\frac{1}{2}\right)+\left(\langle {\widehat{N}}_{-\mathit{k}}\rangle +\frac{1}{2}\right)\right]\left(u_{\omega }^{\phantom{\ast }}{\dot{u}}_{\omega }^{\ast }+{\dot{u}}_{\omega }^{\phantom{\ast }}{u}_{\omega }^{\ast }\right),\hfill \end{array}$$

(75)

with ${\widehat{N}}_{\mathit{k}}^{\phantom{†}}={\widehat{a}}_{\mathit{k}}^{†}{\widehat{a}}_{\mathit{k}}^{\phantom{†}}$, if the initial state is stationary. We again observe that they are the same for $\pm \mathit{k}$ modes and thus the corresponding parameters ${\chi }_{\mathit{k}}$, ${\beta }_{\mathit{k}}$, ${\zeta }_{\mathit{k}}$, ${\kappa }_{\mathit{k}}$ and ${\tau }_{\mathit{k}}$ in (57) and (58) will also be the same for $\pm \mathit{k}$ modes, because they are functions of the covariance matrix elements ${\mathsf{b}}_{\mathit{k}}$, ${\mathsf{a}}_{\mathit{k}}$ and ${\mathsf{c}}_{\mathit{k}}$. Note that from [108], we observe that the density matrix operator thus constructed is fairly general for the configuration under study. In addition, it is worth noticing that, since the second expression in (59) can be identified as the inverse of the nonequilibrium partition function for the canonical variables $({\widehat{\varphi }}_{\mathit{k}},{\widehat{\pi }}_{\mathit{k}})$, we find that the prefactor before the exponential in (69) turns out to be the inverse of the nonequilibrium partition function of the quantum field, from which we can extract the nonequilibrium free energy of the field. Since these quantities are nonequilibrium in nature, they are all functions of time. Moreover, since the field is driven by the parametric process, the state of the field at any moment is a general Gaussian state, given an initial Gaussian state. It means that the nonequilibrium partition function or the nonequilibrium free energy for such a Gaussian system is introduced in a context much more general than conventional equilibrium quantum thermodynamics.

5.2. Density Matrix Elements

The matrix elements of the density operator of a Gaussian state in the ${\varphi }_{\mathit{k}}^{\left(i\right)}$ basis can be written as:

$$\begin{array}{cc}\hfill \varrho ({\varphi }_{\mathit{k}}^{\left(i\right)},{\varphi }_{\mathit{k}}^{\prime \left(i\right)})& ={\mathcal{N}}_{\varrho }exp\left(-{\mathcal{A}}_{\mathit{k}}^{\phantom{2}}{\mathsf{\Delta }}_{\mathit{k}}^{\left(i\right)2}-i2{\mathcal{B}}_{\mathit{k}}^{\phantom{2}}{\mathsf{\Delta }}_{\mathit{k}}^{\left(i\right)}{\Sigma }_{\mathit{k}}^{\left(i\right)}-{\mathcal{C}}_{\mathit{k}}^{\phantom{2}}{\Sigma }_{\mathit{k}}^{\left(i\right)2}\right),\hfill \end{array}$$

(76)

where $i=1$, 2, and the relative coordinate ${\mathsf{\Delta }}_{\mathit{k}}^{\left(i\right)}$, center-of-mass coordinate ${\Sigma }_{\mathit{k}}^{\left(i\right)}$, and normalization constant ${\mathcal{N}}_{\varrho }$ are respectively given by:

$$\begin{array}{cccccc}\hfill {\mathsf{\Delta }}_{\mathit{k}}^{\left(i\right)}& ={\varphi }_{\mathit{k}}^{\left(i\right)}-{\varphi }_{\mathit{k}}^{\prime \left(i\right)},\hfill & \hfill {\Sigma }_{\mathit{k}}^{\left(i\right)}& =\frac{{\varphi }_{\mathit{k}}^{\left(i\right)}+{\varphi }_{\mathit{k}}^{\prime \left(i\right)}}{2},\hfill & \hfill {\mathcal{N}}_{\varrho }& =\frac{\sqrt{{\mathcal{C}}_{\mathit{k}}}}{\sqrt{\pi }},\hfill \end{array}$$

(77)

and

$$\begin{array}{cccccc}\hfill {\mathcal{A}}_{\mathit{k}}& =\frac{{\mathsf{a}}_{\mathit{k}}^{\phantom{2}}{\mathsf{b}}_{\mathit{k}}^{\phantom{2}}-{\mathsf{c}}_{\mathit{k}}^2}{2{\mathsf{b}}_{\mathit{k}}},\hfill & \hfill {\mathcal{B}}_{\mathit{k}}& =-\frac{{\mathsf{c}}_{\mathit{k}}}{2{\mathsf{b}}_{\mathit{k}}},\hfill & \hfill {\mathcal{C}}_{\mathit{k}}& =\frac{1}{2{\mathsf{b}}_{\mathit{k}}},\hfill \end{array}$$

(78)

with

$$\begin{array}{cccccc}\hfill {\mathsf{b}}_{\mathit{k}}& =\langle {\varphi }_{\mathit{k}}^{\left(i\right)2}\rangle ,\hfill & \hfill {\mathsf{a}}_{\mathit{k}}& =\langle p_{\mathit{k}}^{\left(i\right)2}\rangle ,\hfill & \hfill {\mathsf{c}}_{\mathit{k}}& =\frac{1}{2}\langle \left\{{\varphi }_{\mathit{k}}^{\left(i\right)},p_{\mathit{k}}^{\left(i\right)}\right\}\rangle .\hfill \end{array}$$

(79)

If the initial state of the density matrix is a vacuum state, then we always have:

$${\mathsf{b}}_{\mathit{k}}^{\phantom{2}}{\mathsf{a}}_{\mathit{k}}^{\phantom{2}}-{\mathsf{c}}_{\mathit{k}}^2=\frac{1}{4},$$

(80)

so that (76) reduces to:

$$\begin{array}{cc}\hfill \varrho ({\varphi }_{\mathit{k}}^{\left(i\right)},{\varphi }_{\mathit{k}}^{\prime \left(i\right)})& ={\mathcal{N}}_{\varrho }exp\left[-\frac{1}{4b_{\mathit{k}}}\left({\varphi }_{\mathit{k}}^{\left(i\right)2}+{\varphi }_{\mathit{k}}^{\prime \left(i\right)2}\right)+i\frac{c_{\mathit{k}}}{2b_{\mathit{k}}}\left({\varphi }_{\mathit{k}}^{\left(i\right)2}-{\varphi }_{\mathit{k}}^{\prime \left(i\right)2}\right)\right].\hfill \end{array}$$

(81)

6. Evolution of Mode Functions: Bogoliubov Transformation

The dynamics of the field is fully encoded in its density matrix operator or its elements. The time evolution of the density matrix can be described by the unitary time-evolution operator $\widehat{U}(t_f,t_i)$ associated with the Lagrangian (24),

$$\widehat{\rho }\left(t_f\right)=\widehat{U}(t_f,t_i)\widehat{\rho }\left(t_i\right){\widehat{U}}^{†}(t_f,t_i).$$

(82)

A more succinct description of the evolution of a linear quantum system is by using the language of Bogoliubov transformations in the Heisenberg picture. We begin with a brief overview of the Bogoliubov transformation.

6.1. Active and Passive Views of the Bogoliubov Transformation

Suppose the quantum scalar field $\widehat{\varphi }$ is expanded with respect to two different sets of mode functions $u_{\mathit{k}}^{\left(\mathrm{IN}\right)}\left(t\right)$ and $u_{\mathit{k}}^{\left(\mathrm{OUT}\right)}\left(t\right)$. Then according to (31) the expansion takes the form:

$$\begin{array}{cc}\hfill \widehat{\varphi }(\mathit{x},t)& =\sum _{\mathit{k}}\left[{\widehat{a}}_{\mathit{k}}^{\phantom{†}}{u}_{\mathit{k}}^{\left(\mathrm{IN}\right)}(\mathit{x},t)+{\widehat{a}}_{\mathit{k}}^{†}{u}_{\mathit{k}}^{\left(\mathrm{IN}\right)\ast }(\mathit{x},t)\right],\hfill \end{array}$$

(83)

or, in terms of the out-mode functions, in the form

$$\begin{array}{cc}\hfill \widehat{\varphi }(\mathit{x},t)& =\sum _{\mathit{k}}\left[{\widehat{b}}_{\mathit{k}}^{\phantom{†}}{u}_{\mathit{k}}^{\left(\mathrm{OUT}\right)}(\mathit{x},t)+{\widehat{b}}_{\mathit{k}}^{†}{u}_{\mathit{k}}^{\left(\mathrm{OUT}\right)\ast }(\mathit{x},t)\right],\hfill \end{array}$$

(84)

where ${\widehat{a}}_{\mathit{k}}$ and ${\widehat{b}}_{\mathit{k}}$ are the annihilation operators associated $u_{\mathit{k}}^{\left(\mathrm{IN}\right)}\left(t\right)$ and $u_{\mathit{k}}^{\left(\mathrm{OUT}\right)}\left(t\right)$ respectively. If both sets of mode functions are complete, then they will be related by a linear transformation:

$$\begin{array}{cccccc}\hfill u_{\mathit{k}}^{\left(\mathrm{IN}\right)}(\mathit{x},t)& ={\alpha }_{\mathit{k}}^{\phantom{\left(\mathrm{OUT}\right)}}{u}_{\mathit{k}}^{\left(\mathrm{OUT}\right)}(\mathit{x},t)+{\beta }_{\mathit{k}}^{\phantom{\left(\mathrm{OUT}\right)}}{u}_{-\mathit{k}}^{\left(\mathrm{OUT}\right)\ast }(\mathit{x},t),\hfill & \hfill & \mathrm{with} \hfill & \hfill u_{\mathit{k}}^{\left(\mathrm{IN}\right)}(\mathit{x},t)& =e^{i\mathit{k}·\mathit{x}}{u}_{\omega }^{\left(\mathrm{IN}\right)}\left(t\right),\hfill \end{array}$$

(85)

where ${\alpha }_{\mathit{k}}$, ${\beta }_{\mathit{k}}\in \mathbb{C}$ are called the Bogoliubov coefficients [11,109,110], satisfying:

$$\begin{array}{c}\hfill {\left|{\alpha }_{\mathit{k}}\right|}^2-{\left|{\beta }_{\mathit{k}}\right|}^2=1,\end{array}$$

(86)

for the configuration of our interest, and $u_{\omega }^{\left(\mathrm{IN}\right)}\left(t\right)$ is the correspondence of $u_{\omega }\left(t\right)$ in (31) for the IN modes. In fact, Equation (85) is equivalent to:

$$u_{\omega }^{\left(\mathrm{IN}\right)}\left(t\right)={\alpha }_{\mathit{k}}^{\phantom{\left(\mathrm{OUT}\right)}}{u}_{\omega }^{\left(\mathrm{OUT}\right)}\left(t\right)+{\beta }_{\mathit{k}}^{\phantom{\left(\mathrm{OUT}\right)}}{u}_{\omega }^{\left(\mathrm{OUT}\right)\ast }\left(t\right),$$

(87)

after we factor out $e^{i\mathit{k}·\mathit{x}}$. Thus, for an isotropic, spatially-flat space, the Bogoliubov coefficients depend only on $\left|\mathit{k}\right|$ or $\omega$.

The annihilation operator ${\widehat{b}}_{\mathit{k}}$ associated with the OUT mode function $u_{\mathit{k}}^{\left(\mathrm{OUT}\right)}$ will be related to the operators $({\widehat{a}}_{\mathit{k}}^{\phantom{†}},{\widehat{a}}_{\mathit{k}}^{†})$ by:

$${\widehat{b}}_{\mathit{k}}^{\phantom{†}}={\alpha }_{\mathit{k}}^{\phantom{†}}{\widehat{a}}_{\mathit{k}}^{\phantom{†}}+{\beta }_{-\mathit{k}}^{\phantom{†}\ast }{\widehat{a}}_{-\mathit{k}}^{†}.$$

(88)

This implies that in terms of ${\widehat{b}}_{\mathit{k}}$ we have nonzero particle number $N_{\mathit{k}}$ and quantum coherence $C_{\mathit{k}}$

$$\begin{array}{cc}\hfill N_{\mathit{k}}^{\left(\mathrm{OUT}\right)}& =\langle {\widehat{b}}_{+\mathit{k}}^{†}{\widehat{b}}_{+\mathit{k}}^{\phantom{†}}\rangle ={\left|{\beta }_{\mathit{k}}\right|}^2,\hfill \end{array}$$

(89)

$$\begin{array}{cc}\hfill C_{\mathit{k}}^{\left(\mathrm{OUT}\right)}& =\langle {\widehat{b}}_{+\mathit{k}}{\widehat{b}}_{-\mathit{k}}\rangle ={\alpha }_{\mathit{k}}^{\phantom{\ast }}{\beta }_{\mathit{k}}^{\ast },\hfill \end{array}$$

(90)

in the vacuum state associated with ${\widehat{a}}_{\mathit{k}}$ if ${\beta }_{\mathit{k}}\ne 0$, while in comparison, in terms of ${\widehat{a}}_{\mathit{k}}$, they are:

$$\begin{array}{cc}\hfill N_{\mathit{k}}^{\left(\mathrm{IN}\right)}& =\langle {\widehat{a}}_{+\mathit{k}}^{†}{\widehat{a}}_{+\mathit{k}}^{\phantom{†}}\rangle =0,\hfill \end{array}$$

(91)

$$\begin{array}{cc}\hfill C_{\mathit{k}}^{\left(\mathrm{IN}\right)}& =\langle {\widehat{a}}_{+\mathit{k}}{\widehat{a}}_{-\mathit{k}}\rangle =0.\hfill \end{array}$$

(92)

The physical meaning of these results will be explained later.

We may choose ${\alpha }_{\mathit{k}}$ in (86) to be real if we factor out the additional rotation in the transformation (85). This implies that generically the Bogoliubov coefficients can be parametrized by

$$\begin{array}{cccc}\hfill {\alpha }_{\mathit{k}}& =cosh{\eta }_{\mathit{k}},\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hfill & \hfill {\beta }_{\mathit{k}}& =-sinh{\eta }_{\mathit{k}}{e}^{-i{\theta }_{\mathit{k}}},\hfill \end{array}$$

(93)

with the time- and mode-dependent parameters ${\eta }_{\mathit{k}}>0$ and $0\le {\theta }_{\mathit{k}}<2\pi$. These expressions remind us the squeeze transformation of the creation and annihilation operators in quantum optics. Actually, we can write the linear transformation (88) as the two-mode squeezing of ${\widehat{a}}_{\mathit{k}}$ by the two-mode squeeze operator:

$${\widehat{S}}_2^{\phantom{†}}\left({\zeta }_{\mathit{k}}^{\phantom{†}}\right)=exp\left[{\zeta }_{\mathit{k}}^{\ast \phantom{†}}{\widehat{a}}_{\mathit{k}}^{\phantom{†}}{\widehat{a}}_{-\mathit{k}}^{\phantom{†}}-{\zeta }_{\mathit{k}}^{\phantom{†}}{\widehat{a}}_{\mathit{k}}^{†}{\widehat{a}}_{-\mathit{k}}^{†}\right],$$

(94)

with the squeeze parameter ${\zeta }_{\mathit{k}}={\eta }_{\mathit{k}}{e}^{i{\theta }_{\mathit{k}}}$ such that:

$$\begin{array}{cc}\hfill {\widehat{b}}_{\mathit{k}}^{\phantom{†}}& ={\widehat{S}}_2^{†}\left({\zeta }_{\mathit{k}}^{\phantom{†}}\right){\widehat{a}}_{\mathit{k}}^{\phantom{†}}{\widehat{S}}_2^{\phantom{†}}\left({\zeta }_{\mathit{k}}^{\phantom{†}}\right)=cosh{\eta }_{\mathit{k}}^{\phantom{†}}{\widehat{a}}_{\mathit{k}}^{\phantom{†}}-e^{+i{\theta }_{\mathit{k}}}sinh{\eta }_{\mathit{k}}^{\phantom{†}}{\widehat{a}}_{-\mathit{k}}^{†}.\hfill \end{array}$$

(95)

When it applies to the vacuum state, it gives a two-mode squeezed vacuum state:

$$|{\zeta }_{\mathit{k}}^{\left(2\right)}\rangle ={\widehat{S}}_2^{\phantom{†}}\left({\zeta }_{\mathit{k}}^{\phantom{†}}\right)|0_{+\mathit{k}},0_{-\mathit{k}}\rangle =\frac{1}{cosh{\eta }_{\mathit{k}}}\sum _{n=0}^{\infty }{\left(-tanh{\eta }_{\mathit{k}}{e}^{+i{\theta }_{\mathit{k}}}\right)}^n|n_{+\mathit{k}},n_{-\mathit{k}}\rangle ,$$

(96)

in terms of the in-particle number states, and the corresponding density matrix operator is given by:

$$\begin{array}{c}\hfill {\widehat{\rho }}_{\mathit{k}}=|{\zeta }_{\mathit{k}}^{\left(2\right)}\rangle \langle {\zeta }_{\mathit{k}}^{\left(2\right)}|,\end{array}$$

(97)

still a pure state, but in general the off-diagonal terms $n\ne m$ of the density matrix do not vanish for each mode $\mathit{k}$. Equation (96) explicitly shows that the resulting state contains particles in pairs, with opposite momenta. Note that Equation (94) implies that ${\zeta }_{\mathit{k}}={\zeta }_{-\mathit{k}}$, so that we use $+\mathit{k}$ instead of $-\mathit{k}$ for the squeeze parameters in (93). The identification (96) is compatible with the fact that the most general Gaussian state described by the density matrix operator (54) is a squeezed thermal state.

Relating the Bogoliubov coefficients to the squeeze operator in (95) provides an active view of the Bogoliubov transformation. It allows a mapping of an operator in the Heisenberg picture at one time to its counterpart at another time, that is, under time evolution. Or, alternatively, in the Schrödinger picture it maps the time evolution of the state. This is in contrast to the implementation of the Bogoliubov coefficients in (85), where they are used to relate different bases (mode expansion) of the same operator. This is the passive view. Hereafter, we will adopt the active view to investigate the time evolution of the parametric field.

Now from this active view, we see that the parametric process in general will drive the field from its initial vacuum state to a squeezed vacuum state, creating particles in pairs. It is also interesting to note that these $N_{\mathit{k}}$ and $C_{\mathit{k}}$ account for three parameters, needed to fully specify a Gaussian state, as ${\mathsf{a}}_{\mathit{k}}$, ${\mathsf{b}}_{\mathit{k}}$, and ${\mathsf{c}}_{\mathit{k}}$ do, and thus they contain complete information about the linear quantum parametric field in a Gaussian state via the numbers of particle it creates and the correlation between the produced particle pairs.

The particles, as seen by ${\widehat{b}}_{\mathit{k}}$, are correlated. According to the active view, Equation (89) only accounts for spontaneous creation since we have assumed that the initial state is a vacuum. If the initial state is not a vacuum state, but with $\langle {\widehat{N}}_{\mathit{k}}^{\left(\mathrm{IN}\right)}\rangle$ particles present we can show that:

$$\begin{array}{c}\hfill \sum _{\mathit{k}>0}\langle {\widehat{N}}_{+\mathit{k}}^{\left(\mathrm{OUT}\right)}\rangle +\langle {\widehat{N}}_{-\mathit{k}}^{\left(\mathrm{OUT}\right)}\rangle =\sum _{\mathit{k}}2\left({\left|{\beta }_{\mathit{k}}^{\phantom{†}}\right|}^2+\frac{1}{2}\right)\left(\langle {\widehat{N}}_{\mathit{k}}^{\left(\mathrm{IN}\right)}\rangle +\frac{1}{2}\right)-\frac{1}{2},\end{array}$$

(98)

where $\langle \cdots \rangle$ are taken with respect to the initial state. Thus, in addition to spontaneous particle creation, there is also stimulated creation due to the presence of particles in the initial state.

So far we have given a discussion of the generic kinematical properties of particle creation in a time-dependent background and made identification with quantum squeezing by the language of Bogoliubov transformations. We now show their connection with dynamics by using a nonequilibrium quantum dynamics formulation.

6.2. Nonequilibrium Evolution of Modes in Time

In the Heisenberg picture, the time evolution of the parametric field is accounted for by the (field) operators, and for the linear field, this information is encoded in the mode function $u_{\omega }\left(t\right)$. In other words, the time evolution of an operator can be translated into the time evolution of the mode function that span the operator. As an initial-value problem, in general, the mode function $u_{\omega }\left(t\right)$ at any time can be expanded by its initial conditions at $t=0$ by:

$$u_{\omega }\left(t\right)=d_{\omega }^{\left(1\right)}\left(t\right)u_{\omega }\left(0\right)+d_{\omega }^{\left(2\right)}\left(t\right){\dot{u}}_{\omega }\left(0\right),$$

(99)

where the overdot represents the time derivative. The two functions $d_{\omega }^{\left(1\right)}\left(t\right)$ and $d_{\omega }^{\left(2\right)}\left(t\right)$ are a special set of homogeneous solutions to the wave Equation (26), satisfying:

$$\begin{array}{cccccccc}\hfill d_{\omega }^{\left(1\right)}\left(0\right)& =1,\hfill & \hfill {\dot{d}}_{\omega }^{\left(1\right)}\left(0\right)& =0,\hfill & \hfill d_{\omega }^{\left(2\right)}\left(0\right)& =0,\hfill & \hfill {\dot{d}}_{\omega }^{\left(2\right)}\left(0\right)& =1.\hfill \end{array}$$

(100)

This approach is convenient in the sense that, for a linear system, almost all physical observables can be packaged into these two functions plus some initial conditions, so it would be computationally efficient once we get hold of the functions $d_{\omega }^{\left(1\right)}\left(t\right)$ and $d_{\omega }^{\left(2\right)}\left(t\right)$.

For convenience, suppose that the in-mode functions have the form:

$$\begin{array}{cccccc}\hfill u_{\omega }^{\left(\mathrm{IN}\right)}\left(t\right)& =\frac{1}{\sqrt{2{\omega }_i}}{e}^{-i{\omega }_{i}t},\hfill & \hfill & \mathrm{and}\hspace{1em}\hspace{1em}\hfill & \hfill {\dot{u}}_{\omega }^{\left(\mathrm{IN}\right)}\left(t\right)& =-i\sqrt{\frac{{\omega }_i}{2}}{e}^{-i{\omega }_{i}t}.\hfill \end{array}$$

(101)

In terms for the fundamental solutions $d_{\omega }^{\left(i\right)}$, we can write the out-mode functions as:

$$u_{\omega }^{\left(\mathrm{OUT}\right)}\left(t\right)=\frac{1}{\sqrt{2{\omega }_i}}\left[d_{\omega }^{\left(1\right)}\left(t\right)-i{\omega }_i{d}_{\omega }^{\left(2\right)}\left(t\right)\right],$$

(102)

according to (99), where ${\omega }_i=\omega \left(0\right)$ is the initial value of the frequency modulation. Thus, the field operator $\widehat{\varphi }(\mathit{x},t)$ at time t will be given by:

$$\begin{array}{c}\hfill \widehat{\varphi }(\mathit{x},t)=\int \frac{d^{3}k}{{\left(2\pi \right)}^{\frac{3}{2}}}\frac{1}{\sqrt{2{\omega }_i}}\left\{{\widehat{a}}_{\mathit{k}}^{\phantom{†}}{e}^{+i\mathit{k}·\mathit{x}}\left[d_{\omega }^{\left(1\right)}\left(t\right)-i{\omega }_i{d}_{\omega }^{\left(2\right)}\left(t\right)\right]+{\widehat{a}}_{\mathit{k}}^{†}{e}^{-i\mathit{k}·\mathit{x}}\left[d_{\omega }^{\left(1\right)}\left(t\right)+i{\omega }_i{d}_{\omega }^{\left(2\right)}\left(t\right)\right]\right\}\end{array}$$

(103)

in contrast to the in-field before the parametric process starts at $t=0$.

$$\begin{array}{c}\hfill {\widehat{\varphi }}_{\mathrm{IN}}(\mathit{x},t)=\int \frac{d^{3}k}{{\left(2\pi \right)}^{\frac{3}{2}}}\frac{1}{\sqrt{2{\omega }_i}}\left\{{\widehat{a}}_{\mathit{k}}^{\phantom{†}}{e}^{+i\mathit{k}·\mathit{x}}{e}^{-i{\omega }_{i}t}+{\widehat{a}}_{\mathit{k}}^{†}{e}^{-i\mathit{k}·\mathit{x}}{e}^{+i{\omega }_{i}t}\right\}.\end{array}$$

(104)

When we apply the two-mode squeezing to this in-field, we can relate them by:

$$\begin{array}{cc}\hfill \widehat{\varphi }(\mathit{x},t)& ={\widehat{S}}_2^{†}\left(\left\{{\zeta }_{\mathit{k}}\right\}\right){\widehat{\varphi }}_{\mathrm{IN}}(\mathit{x},t){\widehat{S}}_2\left(\left\{{\zeta }_{\mathit{k}}\right\}\right)\hfill \\ \hfill & =\int \frac{d^3\mathit{k}}{{\left(2\pi \right)}^{\frac{3}{2}}}\frac{1}{\sqrt{2{\omega }_i}}\left\{\left[{\alpha }_{\mathit{k}}^{\phantom{\ast }}{e}^{-i{\omega }_{i}t}+{\beta }_{\mathit{k}}^{\phantom{\ast }}{e}^{+i{\omega }_{i}t}\right]{\widehat{a}}_{\mathit{k}}^{\phantom{†}}{e}^{+i\mathit{k}·\mathit{x}}+\left[{\beta }_{\mathit{k}}^{\ast }{e}^{-i{\omega }_{i}t}+{\alpha }_{\mathit{k}}^{\ast }{e}^{+i{\omega }_{i}t}\right]{\widehat{a}}_{\mathit{k}}^{†}{e}^{-i\mathit{k}·\mathit{x}}\right\},\hfill \end{array}$$

(105)

where ${\widehat{S}}_2\left(\left\{{\zeta }_{\mathit{k}}\right\}\right)$ is understood as a collection of two-mode squeeze operators for all $\mathit{k}>0$. Compare this expression with (103) and we find:

$$\begin{array}{cc}\hfill u_{\omega }^{\left(\mathrm{OUT}\right)}\left(t\right)& ={\alpha }_{\mathit{k}}^{\phantom{\ast }}\left(t\right)u_{\omega }^{\left(\mathrm{IN}\right)}\left(t\right)+{\beta }_{\mathit{k}}^{\phantom{\ast }}\left(t\right)u_{\omega }^{\left(\mathrm{IN}\right)\ast }\left(t\right),\hfill \\ \hfill \Rightarrow d_{\omega }^{\left(1\right)}\left(t\right)-i{\omega }_i{d}_{\omega }^{\left(2\right)}\left(t\right)& =e^{-i{\omega }_{i}t}{\alpha }_{\mathit{k}}^{\phantom{\ast }}\left(t\right)+e^{+i{\omega }_{i}t}{\beta }_{\mathit{k}}^{\phantom{\ast }}\left(t\right),\hfill \end{array}$$

(106)

for $t>0$. Applying similar arguments to the conjugate momentum $\widehat{\pi }(\mathit{x},t)$, we will obtain:

$$\begin{array}{cc}\hfill {\dot{u}}_{\omega }^{\left(\mathrm{OUT}\right)}\left(t\right)& ={\alpha }_{\mathit{k}}^{\phantom{\ast }}\left(t\right){\dot{u}}_{\omega }^{\left(\mathrm{IN}\right)}\left(t\right)+{\beta }_{\mathit{k}}^{\phantom{\ast }}\left(t\right){\dot{u}}_{\omega }^{\left(\mathrm{IN}\right)\ast }\left(t\right),\hfill \\ \hfill \Rightarrow {\dot{d}}_{\omega }^{\left(1\right)}\left(t\right)-i{\omega }_i{\dot{d}}_{\omega }^{\left(2\right)}\left(t\right)& =-i{\omega }_i{e}^{-i{\omega }_{i}t}{\alpha }_{\mathit{k}}\left(t\right)+i{\omega }_i{e}^{+i{\omega }_{i}t}{\beta }_{\mathit{k}}\left(t\right).\hfill \end{array}$$

(107)

Solving (106) and (107) simultaneously we arrive at:

$$\begin{array}{cc}\hfill {\alpha }_{\mathit{k}}\left(t\right)& =\frac{1}{2{\omega }_i}{e}^{+i{\omega }_{i}t}\left[{\omega }_i{d}_{\omega }^{\left(1\right)}\left(t\right)+i{\dot{d}}_{\omega }^{\left(1\right)}\left(t\right)-i{\omega }_i^2{d}_{\omega }^{\left(2\right)}\left(t\right)+{\omega }_i{\dot{d}}_{\omega }^{\left(2\right)}\left(t\right)\right],\hfill \end{array}$$

(108)

$$\begin{array}{cc}\hfill {\beta }_{\mathit{k}}\left(t\right)& =\frac{1}{2{\omega }_i}{e}^{-i{\omega }_{i}t}\left[{\omega }_i{d}_{\omega }^{\left(1\right)}\left(t\right)-i{\dot{d}}_{\omega }^{\left(1\right)}\left(t\right)-i{\omega }_i^2{d}_{\omega }^{\left(2\right)}\left(t\right)-{\omega }_i{\dot{d}}_{\omega }^{\left(2\right)}\left(t\right)\right].\hfill \end{array}$$

(109)

Thus we have expressed the Bogoliubov coefficients of all modes at any time by the corresponding fundamental solutions $d_{\omega }^{\left(i\right)}\left(t\right)$. Then from (93), we can write the squeeze parameters by the same set of fundamental solutions. For example, we can show:

$${\left|{\alpha }_{\mathit{k}}\right|}^2+{\left|{\beta }_{\mathit{k}}\right|}^2=\frac{1}{2{\omega }_i^2}\left[{\omega }_i^2{d}_{\omega }^{\left(1\right)2}\left(t\right)+{\dot{d}}_{\omega }^{\left(1\right)2}\left(t\right)+{\omega }_i^4{d}_{\omega }^{\left(2\right)2}\left(t\right)+{\omega }_i^2{\dot{d}}_{\omega }^{\left(2\right)2}\left(t\right)\right]=cosh2{\eta }_{\mathit{k}}.$$

(110)

This also implies that the squeeze parameters ${\eta }_{\mathit{k}}$, ${\theta }_{\mathit{k}}$ are actually functions of $\left|\mathit{k}\right|$.

The covariance matrix elements ${\mathsf{b}}_{\mathit{k}}$, ${\mathsf{a}}_{\mathit{k}}$ and ${\mathsf{c}}_{\mathit{k}}$ can also be expressed by these fundamental solutions via the Bogoliubov coefficients. For example, from (48) and (31), we have:

$$\begin{array}{c}\hfill {\widehat{\varphi }}_{\mathit{k}}^{\left(1\right)}\left(t\right)=\frac{1}{\sqrt{2}}\left[{\widehat{a}}_{+\mathit{k}}^{\phantom{†}}{u}_{\omega }^{\phantom{\ast }(\mathrm{OUT}}\left(t\right)+{\widehat{a}}_{-\mathit{k}}^{†}{u}_{\omega }^{\left(\mathrm{OUT}\right)\ast }\left(t\right)+{\widehat{a}}_{-\mathit{k}}^{\phantom{†}}{u}_{\omega }^{\left(\mathrm{OUT}\right)\phantom{\ast }}\left(t\right)+{\widehat{a}}_{+\mathit{k}}^{†}{u}_{\omega }^{\left(\mathrm{OUT}\right)\ast }\left(t\right)\right],\end{array}$$

(111)

and thus at any time:

$$\begin{array}{cc}\hfill {\mathsf{b}}_{\mathit{k}}& =\langle {\widehat{\varphi }}_{\mathit{k}}^{\left(1\right)2}\rangle =2\left\{{\left|{\beta }_{\mathit{k}}\right|}^2+\frac{1}{2}\right\}{u}_{\omega }^{\left(\mathrm{IN}\right)}{u}_{\omega }^{\left(\mathrm{IN}\right)\ast }+{\alpha }_{\mathit{k}}^{\phantom{\ast }}{\beta }_{\mathit{k}}^{\ast }{u}_{\omega }^{\left(\mathrm{IN}\right)2}+{\alpha }_{\mathit{k}}^{\ast }{\beta }_{\mathit{k}}^{\phantom{\ast }}{u}_{\omega }^{\left(\mathrm{IN}\right)\ast 2},\hfill \end{array}$$

(112)

$$\begin{array}{cc}\hfill {\mathsf{a}}_{\mathit{k}}& =\langle p_{\mathit{k}}^{\left(1\right)2}\rangle =2\left\{{\left|{\beta }_{\mathit{k}}\right|}^2+\frac{1}{2}\right\}{\dot{u}}_{\omega }^{\left(\mathrm{IN}\right)}{\dot{u}}_{\omega }^{\left(\mathrm{IN}\right)\ast }+{\alpha }_{\mathit{k}}^{\phantom{\ast }}{\beta }_{\mathit{k}}^{\ast }{\dot{u}}_{\omega }^{\left(\mathrm{IN}\right)2}+{\alpha }_{\mathit{k}}^{\ast }{\beta }_{\mathit{k}}^{\phantom{\ast }}{\dot{u}}_{\omega }^{\left(\mathrm{IN}\right)\ast 2},\hfill \end{array}$$

(113)

$$\begin{array}{cc}\hfill {\mathsf{c}}_{\mathit{k}}& =\frac{1}{2}\langle \left\{{\varphi }_{\mathit{k}}^{\left(1\right)},p_{\mathit{k}}^{\left(1\right)}\right\}\rangle \hfill \\ \hfill & =\left\{{\left|{\beta }_{\mathit{k}}\right|}^2+\frac{1}{2}\right\}\frac{d}{dt}\left[u_{\omega }^{\left(\mathrm{IN}\right)}{u}_{\omega }^{\left(\mathrm{IN}\right)\ast }\right]+\frac{1}{2}{\alpha }_{\mathit{k}}^{\phantom{\ast }}{\beta }_{\mathit{k}}^{\ast }\frac{d}{dt}{u}_{\omega }^{\left(\mathrm{IN}\right)2}+\frac{1}{2}{\alpha }_{\mathit{k}}^{\ast }{\beta }_{\mathit{k}}^{\phantom{\ast }}\frac{d}{dt}{u}_{\omega }^{\left(\mathrm{IN}\right)\ast 2}.\hfill \end{array}$$

(114)

Then we can further use (108) and (109) or simply use only (102) and (107). If the in-mode function is of the form (101), we can easily relate the covariance matrix elements with the particle number (89) and coherence (90) in the out-state by:

$$\begin{array}{cc}\hfill {\mathsf{b}}_{\mathit{k}}\left(t\right)& =d_{\omega }^{\left(1\right)2}\left(t\right)\frac{1}{2{\omega }_i}+d_{\omega }^{\left(2\right)2}\left(t\right)\frac{{\omega }_i}{2}\hfill \end{array}$$

(115)

$$\begin{array}{cc}\hfill & =\left(N_{\mathit{k}}^{\left(\mathrm{OUT}\right)}\left(t\right)+\frac{1}{2}\right)\frac{1}{{\omega }_i}+\frac{1}{2{\omega }_i}{C}_{\mathit{k}}^{\left(\mathrm{OUT}\right)\phantom{\ast }}\left(t\right)e^{-i2{\omega }_{i}t}+\frac{1}{2{\omega }_i}{C}_{\mathit{k}}^{\left(\mathrm{OUT}\right)\ast }\left(t\right)e^{+i2{\omega }_{i}t},\hfill \end{array}$$

(116)

$$\begin{array}{cc}\hfill {\mathsf{a}}_{\mathit{k}}\left(t\right)& ={\dot{d}}_{\omega }^{\left(1\right)2}\left(t\right)\frac{1}{2{\omega }_i}+{\dot{d}}_{\omega }^{\left(2\right)2}\left(t\right)\frac{{\omega }_i}{2}\hfill \end{array}$$

(117)

$$\begin{array}{cc}\hfill & =\left(N_{\mathit{k}}^{\left(\mathrm{OUT}\right)}\left(t\right)+\frac{1}{2}\right){\omega }_i-\frac{{\omega }_i}{2}{C}_{\mathit{k}}^{\left(\mathrm{OUT}\right)\phantom{\ast }}\left(t\right)e^{-i2{\omega }_{i}t}-\frac{{\omega }_i}{2}{C}_{\mathit{k}}^{\left(\mathrm{OUT}\right)\ast }\left(t\right)e^{+i2{\omega }_{i}t},\hfill \end{array}$$

(118)

$$\begin{array}{cc}\hfill {\mathsf{c}}_{\mathit{k}}\left(t\right)& =d_{\mathit{k}}^{\left(1\right)}\left(t\right){\dot{d}}_{\omega }^{\left(1\right)}\left(t\right)\frac{1}{2{\omega }_i}+d_{\omega }^{\left(2\right)}\left(t\right){\dot{d}}_{\mathit{k}}^{\left(2\right)}\left(t\right)\frac{{\omega }_i}{2}\hfill \end{array}$$

(119)

$$\begin{array}{cc}\hfill & =-\frac{i}{2}{C}_{\mathit{k}}^{\left(\mathrm{OUT}\right)\phantom{\ast }}\left(t\right)e^{-i2{\omega }_{i}t}+\frac{i}{2}{C}_{\mathit{k}}^{\left(\mathrm{OUT}\right)\ast }\left(t\right)e^{+i2{\omega }_{i}t},\hfill \end{array}$$

(120)

or vice versa, write the particle number and coherence in terms of the fundamental solutions, attaining in principle all detailed information of their time evolution.

$$\begin{array}{cc}\hfill N_{\mathit{k}}^{\left(\mathrm{OUT}\right)}\left(t\right)& ={sinh}^2{\eta }_{\mathit{k}}\left(t\right)=\frac{1}{4}\left\{{\left[d_{\omega }^{\left(1\right)}\left(t\right)-{\dot{d}}_{\omega }^{\left(2\right)}\left(t\right)\right]}^2+{\left[{\omega }_i{d}_{\omega }^{\left(2\right)}\left(t\right)+\frac{1}{{\omega }_i}{\dot{d}}_{\omega }^{\left(1\right)}\left(t\right)\right]}^2\right\}\hfill \\ \hfill & =\frac{1}{4}\left[d_{\omega }^{\left(1\right)2}\left(t\right)+\frac{1}{{\omega }_i^2}{\dot{d}}_{\omega }^{\left(1\right)2}\left(t\right)+{\omega }_i^2{d}_{\omega }^{\left(2\right)2}\left(t\right)+{\dot{d}}_{\omega }^{\left(2\right)2}\left(t\right)\right]-\frac{1}{2},\hfill \end{array}$$

(121)

and

$$\begin{array}{cc}\hfill C_{\mathit{k}}^{\left(\mathrm{OUT}\right)}\left(t\right)& =-\frac{1}{2}sinh2{\eta }_{\mathit{k}}\left(t\right)e^{+i{\theta }_{\mathit{k}}\left(t\right)}\hfill \\ \hfill & =\frac{e^{i2{\omega }_{i}t}}{4}\left\{{\left[d_{\mathit{k}}^{\left(1\right)}\left(t\right)+\frac{i}{{\omega }_i}{\dot{d}}_{\mathit{k}}^{\left(1\right)}\left(t\right)\right]}^2+{\omega }_i^2{\left[d_{\mathit{k}}^{\left(2\right)}\left(t\right)+\frac{i}{{\omega }_i}{\dot{d}}_{\mathit{k}}^{\left(2\right)}\left(t\right)\right]}^2\right\}.\hfill \end{array}$$

(122)

We observe that, comparing with $N_{\mathit{k}}^{\left(\mathrm{OUT}\right)}\left(t\right)$, the coherence $C_{\mathit{k}}^{\left(\mathrm{OUT}\right)}\left(t\right)$ tends to oscillate more rapidly with time (but not necessarily always so), even though in principle both particle number and coherence are oscillatory, because $d_{\omega }^{\left(i\right)}$ are solutions of the parametric oscillator. To explicitly illustrate this point, we provide a few examples. First consider a parametric oscillator with the time-dependent frequency:

$${\omega }^2\left(t\right)={\mathfrak{R}}^2{cos}^2\varpi t,$$

(123)

where $\varpi$ is the frequency of modulation and $\mathfrak{R}$ is the amplitude. From Figure 1, we find that both $N\left(t\right)$ and $C\left(t\right)$ are oscillatory more or less at the same tempo except that $N\left(t\right)$ is always nonnegative. Here we will suppress the superscript $\left(\mathrm{OUT}\right)$ for the moment because we are not comparing the results at two special instances of time. A contrasting example is the case with frequency.

$${\omega }^2\left(t\right)={\mathfrak{R}}^2\left(1-e^{-2\varpi t}\right),$$

(124)

where the frequency modulation monotonically rises from zero to a constant $\mathfrak{R}$. In this case, both $N\left(t\right)$ and $C\left(t\right)$ barely oscillate and quickly settle down to constant values, as shown in Figure 2. However, in the case of cosmological perturbation, if we consider (17) in the leading order and ignore the slow-roll parameter for the moment, we find that $N_{\mathit{k}}^{\left(\mathrm{OUT}\right)}$ happens to be non-oscillating and monotonically increases like $4^{-1}{\left(k\eta \right)}^{-4}$ as the conformal time $\eta \to 0^{-}$, but the coherence still oscillates with an increasing amplitude over the evolution toward the end of inflation, and the rate of oscillations depends on the mode of perturbations. The oscillatory behavior of the coherence will come to stop at about $k\eta \simeq 1$, that is, the horizon crossing time, after which it evolves monotonically:

$$\underset{k\eta \ll 1}{lim}{C}_{\mathit{k}}\left(\eta \right)\simeq -\frac{1}{4k^4{\eta }^4}+i\frac{1}{k^3{\eta }^3}.$$

(125)

Thus, for the superhorizon modes, both the particle number density and the coherence are non-oscillating. These can be seen in Figure 3. These examples illustrate our point that the contributions of quantum coherence need be taken more seriously, and arguments to justify its cancellation attributed to the rapid oscillations of coherence in comparison to particle number may be too simplistic. This observation is essential in the sense that our subsequent discussion on (possible) entropy production associated with particle creation applies to each mode of the field, and it does not involve the summation over a certain band of the field spectrum nor a measurement over a finite time duration.

For a Gaussian system, the physical quantities we are interested in can be built from the covariance matrix elements. We shall now focus on the von Neumann entropy.

7. Von Neumann Entropy of the Parametric Quantum Field

The von Neumann entropy ${\mathcal{S}}_{\mathit{k}}$ for $({\varphi }_{\mathit{k}}^{\left(1\right)},p_{\mathit{k}}^{\left(1\right)})$

$${\mathcal{S}}_{\mathit{k}}=-Tr\left\{\widehat{\rho }({\widehat{\varphi }}_{\mathit{k}}^{\left(1\right)},{\widehat{\pi }}_{\mathit{k}}^{\left(1\right)})ln\widehat{\rho }({\widehat{\varphi }}_{\mathit{k}}^{\left(1\right)},{\widehat{\pi }}_{\mathit{k}}^{\left(1\right)})\right\}$$

(126)

can be given by the symplectic eigenvalue of the covariance matrix $\mathit{\sigma }$ associated with the canonical variables $({\varphi }_{\mathit{k}}^{\left(1\right)},p_{\mathit{k}}^{\left(1\right)})$

$${\mathcal{S}}_{\mathit{k}}=({\lambda }_{\mathit{k}}+\frac{1}{2})ln({\lambda }_{\mathit{k}}+\frac{1}{2})-({\lambda }_{\mathit{k}}-\frac{1}{2})ln({\lambda }_{\mathit{k}}-\frac{1}{2}),$$

(127)

where the symplectic eigenvalue ${\lambda }_{\mathit{k}}$ is defined by:

$${\lambda }_{\mathit{k}}^2={\mathsf{a}}_{\mathit{k}}^{\phantom{\ast }}{\mathsf{b}}_{\mathit{k}}^{\phantom{\ast }}-{\mathsf{c}}_{\mathit{k}}^2=\langle {\varphi }_{\mathit{k}}^{\left(1\right)2}\rangle \langle p_{\mathit{k}}^{\left(1\right)2}\rangle -\left[\frac{1}{2}\langle \{{\varphi }_{\mathit{k}}^{\left(1\right)},p_{\mathit{k}}^{\left(1\right)}\}\rangle \right]{}^2.$$

(128)

As a reminder, the symplectic eigenvalue $\lambda$ of a $2N$-dimensional covariance matrix $\mathit{\sigma }$ can be found by solving the eigenvalue problem of the form

$$\mathit{\sigma }·\mathit{v}=\lambda \Sigma ·\mathit{v},\mathrm{or}\mathrm{equivalently}\left(\Sigma ·\mathit{\sigma }\right)·\mathit{v}=\lambda \mathit{v},$$

where the matrix $\Sigma$ is the fundamental symplectic matrix

$$\Sigma =\underset{k=1}{\overset{N}{⨁}}\left(\begin{array}{cc}0& +i\\ -i& 0\end{array}\right).$$

From the covariance matrix elements (112)–(114), if the in-mode functions satisfy the Wronskian conditions, then we find the symplectic eigenvalues can be expressed by $N_{\mathit{k}}$ and $C_{\mathit{k}}$ as well,

$$\begin{array}{cc}\hfill {\lambda }_{\mathit{k}}^2={\mathsf{a}}_{\mathit{k}}^{\phantom{\ast }}{\mathsf{b}}_{\mathit{k}}^{\phantom{\ast }}-{\mathsf{c}}_{\mathit{k}}^2& =-\left[{\left(N_{\mathit{k}}^{\left(\mathrm{OUT}\right)}+\frac{1}{2}\right)}^2-{\left|C_{\mathit{k}}^{\left(\mathrm{OUT}\right)}\right|}^2\right]\left[u_k^{\left(\mathrm{IN}\right)}{\dot{u}}_k^{\left(\mathrm{IN}\right)\ast }-{\dot{u}}_k^{\left(\mathrm{IN}\right)}{u}_k^{\left(\mathrm{IN}\right)\ast }\right]{}^2\hfill \\ \hfill & =(N_{\mathit{k}}^{\left(\mathrm{OUT}\right)}+\frac{1}{2})2-{\left|C_{\mathit{k}}^{\left(\mathrm{OUT}\right)}\right|}^2,\hfill \end{array}$$

(129)

independent of the phase of $C_{\mathit{k}}^{\left(\mathrm{OUT}\right)}$ and independent of the explicit expressions of the in-mode functions. The latter merely reflects the invariance of the von Neumann entropy of the system during its unitary time evolution. The former implies that the phase of the coherence plays no role in evaluating the symplectic eigenvalue, and thus the von Neumann entropy. It further suggests we should not assume that a rapidly oscillating coherence can be canceled out and use this to smear off its contribution to the von Neumann entropy for each field mode.

Now we see that the entropy thus generated is formally somewhat different from that of a thermal state of the harmonic oscillator, or of a photon mode, we often come across:

$$\begin{array}{cccc}\hfill {\mathcal{S}}_{\mathit{k}}& =\left(N_{\mathit{k}}+1\right)ln\left(N_{\mathit{k}}+1\right)-N_{\mathit{k}}ln{N}_{\mathit{k}},\hspace{1em}\hspace{1em}\hspace{1em}\hfill & \hfill N_{\mathit{k}}& =\frac{1}{e^{\beta {\omega }_i}-1},\hfill \end{array}$$

(130)

by an additional coherence ${\left|C_{\mathit{k}}\right|}^2$. The expressions of $N_{\mathit{k}}$ is different too. At first sight, it is rather baffling why there exists any entropy (127) for a parametric process from a vacuum state, since the final state of each mode is a pure two-mode squeezed state. It is irrelevant to the basis we use in calculating the trace in (126) because the entropy is defined in terms of the trace and the eigenvalues ${\lambda }_{\mathit{k}}$ are invariant under the linear symplectic transformations. It is of importance to reconcile the paradoxical coexistence of the time-reversal of unitary evolution and the time irreversible entropy production.

The entropy (127), together with (128), seems to have a nonzero value, but the squeeze transformation alone does not change the value of the uncertainty function, that is, ${\lambda }_{\mathit{k}}^2$, so if the initial state before the parametric process is a vacuum state, whose entropy vanishes, then the resulting squeezed vacuum state due to the parametric process will also have zero entropy. This is consistent with the fact that both the vacuum state and the squeezed vacuum state are pure states, and thus their entropies are zero. Hence we expect (127) should vanish. Alternatively, since the density matrix of a closed system is governed by the Liouville equation:

$$\frac{\partial \widehat{\rho }\left(t\right)}{\partial t}=-i\left[\widehat{H}\left(t\right),\widehat{\rho }\left(t\right)\right],$$

(131)

we can explicitly show that the von Neumann entropy is a constant of motion:

$$\begin{array}{c}\hfill \frac{d}{dt}\mathcal{S}\left(t\right)=-\frac{d}{dt}Tr\widehat{\rho }\left(t\right)ln\widehat{\rho }\left(t\right)=iTr\left\{\left[\widehat{H}\left(t\right),\widehat{\rho }\left(t\right)\right]ln\widehat{\rho }\left(t\right)\right\}=0,\end{array}$$

(132)

where we have used the cyclic property of the trace and

$$\begin{array}{cccccc}\hfill Tr\widehat{\rho }\left(t\right)& =1,\hfill & \hfill & \Rightarrow \hfill & \hfill Tr\frac{\partial \widehat{\rho }\left(t\right)}{\partial t}& =0.\hfill \end{array}$$

(133)

That is, the entropy is indeed conserved under the unitary evolution. It implies that even the number of particles is changing, the entropy $\mathcal{S}$ remains the same.

On the other hand from (121) and (122) we can explicitly verify that indeed (129) is:

$${\lambda }_{\mathit{k}}^2={\left(N_{\mathit{k}}^{\left(\mathrm{OUT}\right)}+\frac{1}{2}\right)}^2-{\left|C_{\mathit{k}}^{\left(\mathrm{OUT}\right)}\right|}^2=\frac{1}{4}{\left[d_{\mathit{k}}^{\left(1\right)}\left(t\right){\dot{d}}_{\mathit{k}}^{\left(2\right)}\left(t\right)-{\dot{d}}_{\mathit{k}}^{\left(1\right)}\left(t\right)d_{\mathit{k}}^{\left(2\right)}\left(t\right)\right]}^2=\frac{1}{4},$$

(134)

so that the corresponding entropy (127) of each field mode due to the parametric process does give zero, consistent with the earlier arguments. Thus we definitely will not have any entropy production for the parametric process evolving from the initial vacuum state.

At this point we have assumed that the field evolve unitarily, so the field evolves from a pure state into another one. We also assume that we can have full access to every bit of information the field possesses. However, in reality, a physical system is never perfectly closed, but coupled to some environment. When some information of the system is lost to its environment, or if some coarse-graining measure is introduced in the measurement protocols, entropy is observed to be generated in the reduced system. When a system undergoes nonequilibrium evolution some information is likely to be lost to the environment possessing a much larger number of degrees of freedom. The key point is coarse-graining, resulting in certain information becoming irretrievable. If we coarse-grain the information in the field we shall see nonzero entropy generated in the evolution of a pure state.

8. Evolution of the Wigner Function of the Quantum Field

In the context of creation of particle pairs, the Wigner function of the field offers the clearest picture of how the coarse-grained entropy may emerge, and a phase-space description offers an intimate view of its connection with quantum thermodynamics.

8.1. Wigner Function of the Field Mode

The Wigner function can be related to the density matrix elements (76) by a Fourier transformation,

$${\mathcal{W}}_{\mathit{k}}^{\left(i\right)}({\Sigma }_{\mathit{k}}^{\left(i\right)},P_{\mathit{k}}^{\left(i\right)};t)=\int \frac{du}{2\pi }{e}^{-i{P}_{\mathit{k}}^{\left(i\right)}{\mathsf{\Delta }}_{\mathit{k}}^{\left(i\right)}}\rho ({\Sigma }_{\mathit{k}}^{\left(i\right)}+\frac{{\mathsf{\Delta }}_{\mathit{k}}^{\left(i\right)}}{2},{\Sigma }_{\mathit{k}}^{\left(i\right)}-\frac{{\mathsf{\Delta }}_{\mathit{k}}^{\left(i\right)}}{2};t),$$

(135)

and we obtain:

$${\mathcal{W}}_{\mathit{k}}^{\left(i\right)}({\Sigma }_{\mathit{k}}^{\left(i\right)},P_{\mathit{k}}^{\left(i\right)};t)=\frac{1}{2\pi \sqrt{{\mathsf{a}}_{\mathit{k}}^{\phantom{2}}{\mathsf{b}}_{\mathit{k}}^{\phantom{2}}-{\mathsf{c}}_{\mathit{k}}^2}}exp\left[-\frac{{\mathsf{b}}_{\mathit{k}}^{\phantom{2}}{P}_{\mathit{k}}^{\left(i\right)2}-2{\mathsf{c}}_{\mathit{k}}^{\phantom{2}}{P}_{\mathit{k}}^{\left(i\right)}{\Sigma }_{\mathit{k}}^{\left(i\right)}+{\mathsf{a}}_{\mathit{k}}^{\phantom{2}}{\Sigma }_{\mathit{k}}^{\left(i\right)2}}{2({\mathsf{a}}_{\mathit{k}}^{\phantom{2}}{\mathsf{b}}_{\mathit{k}}^{\phantom{2}}-{\mathsf{c}}_{\mathit{k}}^2)}\right].$$

(136)

The parameters $N_{\mathit{k}}$ and $C_{\mathit{k}}$ include the information we need about the Gaussian state of the linear parametric field. If we make the substitutions:

$$\begin{array}{cc}\hfill {\mathsf{b}}_{\mathit{k}}& =2\left(N_{\mathit{k}}^{\left(\mathrm{OUT}\right)}+\frac{1}{2}\right)u_{\omega }^{\left(\mathrm{IN}\right)}{u}_{\omega }^{\left(\mathrm{IN}\right)\ast }+C_{\mathit{k}}^{\left(\mathrm{OUT}\right)\phantom{\ast }}{u}_{\omega }^{\left(\mathrm{IN}\right)2}+C_{\mathit{k}}^{\left(\mathrm{OUT}\right)\ast }{u}_{\omega }^{\left(\mathrm{IN}\right)\ast 2},\hfill \end{array}$$

(137)

$$\begin{array}{cc}\hfill {\mathsf{a}}_{\mathit{k}}& =2\left(N_{\mathit{k}}^{\left(\mathrm{OUT}\right)}+\frac{1}{2}\right){\dot{u}}_{\omega }^{\left(\mathrm{IN}\right)}{\dot{u}}_{\omega }^{\left(\mathrm{IN}\right)\ast }+C_{\mathit{k}}^{\left(\mathrm{OUT}\right)\phantom{\ast }}{\dot{u}}_{\omega }^{\left(\mathrm{IN}\right)2}+C_{\mathit{k}}^{\left(\mathrm{OUT}\right)\ast }{\dot{u}}_{\omega }^{\left(\mathrm{IN}\right)\ast 2},\hfill \end{array}$$

(138)

$$\begin{array}{cc}\hfill {\mathsf{c}}_{\mathit{k}}& =\left(N_{\mathit{k}}^{\left(\mathrm{OUT}\right)}+\frac{1}{2}\right)\left({\dot{u}}_{\omega }^{\left(\mathrm{IN}\right)}{u}_{\omega }^{\left(\mathrm{IN}\right)\ast }+u_{\omega }^{\left(\mathrm{IN}\right)}{\dot{u}}_{\omega }^{\left(\mathrm{IN}\right)\ast }\right)+C_{\mathit{k}}^{\left(\mathrm{OUT}\right)\phantom{\ast }}{u}_{\omega }^{\left(\mathrm{IN}\right)}{\dot{u}}_{\omega }^{\left(\mathrm{IN}\right)}+C_{\mathit{k}}^{\left(\mathrm{OUT}\right)\ast }{u}_{\omega }^{\left(\mathrm{IN}\right)\ast }{\dot{u}}_{\omega }^{\left(\mathrm{IN}\right)\ast },\hfill \end{array}$$

(139)

and use the facts:

$$\begin{array}{cccc}\hfill u_{\omega }^{\left(\mathrm{IN}\right)\phantom{\ast }}\left(t\right){\dot{u}}_{\omega }^{\left(\mathrm{IN}\right)\ast }\left(t\right)-u_{\omega }^{\left(\mathrm{IN}\right)\ast }\left(t\right){\dot{u}}_{\omega }^{\left(\mathrm{IN}\right)\phantom{\ast }}\left(t\right)& =i,\hfill & \hfill {\left(N_{\mathit{k}}^{\left(\mathrm{OUT}\right)}+\frac{1}{2}\right)}^2-{\left|C_{\mathit{k}}^{\left(\mathrm{OUT}\right)}\right|}^2& =\frac{1}{4}\hfill \end{array}$$

(140)

for the current setting, we find

$${\mathcal{W}}_{\mathit{k}}^{\left(i\right)}({\Sigma }_{\mathit{k}}^{\left(i\right)},P_{\mathit{k}}^{\left(i\right)};t)=\frac{1}{\pi }exp\left\{-2\left[2\left(N_{\mathit{k}}^{\left(\mathrm{OUT}\right)}+\frac{1}{2}\right)D_{\mathit{k}}^{\left(i\right)}{D}_{\mathit{k}}^{\left(i\right)\ast }+C_{\mathit{k}}^{\left(\mathrm{OUT}\right)\phantom{\ast }}{D}_{\mathit{k}}^{\left(i\right)2}+C_{\mathit{k}}^{\left(\mathrm{OUT}\right)\ast }{D}_{\mathit{k}}^{\left(i\right)\ast 2}\right]\right\},$$

(141)

where $D_{\mathit{k}}^{\left(i\right)}=P_{\mathit{k}}^{\left(i\right)}{u}_{\omega }^{\left(\mathrm{IN}\right)}-{\Sigma }_{\mathit{k}}^{\left(i\right)}{\dot{u}}_{\omega }^{\left(\mathrm{IN}\right)}$. We have expressed the coefficients in the exponent of the Gaussian-state Wigner function in terms of $N_{\mathit{k}}$ and $C_{\mathit{k}}$.

Equation (141) is not quite what we want because it is written in terms of the real and the imaginary parts of the field modes. Rather we would like to write it as a function of the field modes ${\varphi }_{\mathit{k}}$ and the conjugate momenta ${\pi }_{\mathit{k}}$. Inspired by (48) and (49)

$$\begin{array}{cccc}\hfill {\widehat{\varphi }}_{\mathit{k}}^{\left(1\right)}& =\frac{1}{\sqrt{2}}\left({\widehat{\varphi }}_{\mathit{k}}^{\phantom{†}}+{\widehat{\varphi }}_{\mathit{k}}^{†}\right),\hfill & \hfill {\widehat{\varphi }}_{\mathit{k}}^{\left(2\right)}& =\frac{1}{i\sqrt{2}}\left({\widehat{\varphi }}_{\mathit{k}}^{\phantom{†}}-{\widehat{\varphi }}_{\mathit{k}}^{†}\right),\hfill \end{array}$$

(142)

$$\begin{array}{cccc}\hfill {\widehat{p}}_{\mathit{k}}^{\left(1\right)}& =\frac{1}{\sqrt{2}}\left({\widehat{\pi }}_{\mathit{k}}^{†}+{\widehat{\pi }}_{\mathit{k}}^{\phantom{†}}\right),\hfill & \hfill {\widehat{p}}_{\mathit{k}}^{\left(2\right)}& =\frac{1}{i\sqrt{2}}\left({\widehat{\pi }}_{\mathit{k}}^{†}-{\widehat{\varphi }}_{\mathit{k}}^{\phantom{†}}\right),\hfill \end{array}$$

(143)

we will combine the Wigner functions ${\mathcal{W}}_{\mathit{k}}^{\left(1\right)}({\Sigma }_{\mathit{k}}^{\left(1\right)},P_{\mathit{k}}^{\left(1\right)};t)$ and ${\mathcal{W}}_{\mathit{k}}^{\left(2\right)}({\Sigma }_{\mathit{k}}^{\left(2\right)},P_{\mathit{k}}^{\left(2\right)};t)$ to form the Wigner function for the modes $\pm \mathit{k}$:

$$\begin{array}{cc}\hfill & {\mathcal{W}}_{\pm \mathit{k}}({\Sigma }_{+\mathit{k}},{\mathsf{\Pi }}_{+\mathit{k}};{\Sigma }_{-\mathit{k}},{\mathsf{\Pi }}_{-\mathit{k}};t)={\mathcal{W}}_{\mathit{k}}^{\left(1\right)}({\Sigma }_{\mathit{k}}^{\left(1\right)},P_{\mathit{k}}^{\left(1\right)};t)\times {\mathcal{W}}_{\mathit{k}}^{\left(2\right)}({\Sigma }_{\mathit{k}}^{\left(2\right)},P_{\mathit{k}}^{\left(2\right)};t)\hfill \\ \hfill & ={\left(\frac{1}{\pi }\right)}^{2}exp\left\{-2\left[2\left(N_{\mathit{k}}^{\left(\mathrm{OUT}\right)}+\frac{1}{2}\right)\left(D_{\mathit{k}}^{\left(1\right)}{D}_{\mathit{k}}^{\left(1\right)\ast }+D_{\mathit{k}}^{\left(2\right)}{D}_{\mathit{k}}^{\left(2\right)\ast }\right)\right.\right.\hfill \\ \hfill & +\left.\left.C_{\mathit{k}}^{\left(\mathrm{OUT}\right)\phantom{\ast }}\left(D_{\mathit{k}}^{\left(1\right)2}+D_{\mathit{k}}^{\left(2\right)2}\right)+C_{\mathit{k}}^{\left(\mathrm{OUT}\right)\ast }\left(D_{\mathit{k}}^{\left(1\right)\ast 2}+D_{\mathit{k}}^{\left(2\right)\ast 2}\right)\right]\right\}.\hfill \end{array}$$

(144)

Observe that we can write $D_{\mathit{k}}^{\left(1\right)}{D}_{\mathit{k}}^{\left(1\right)\ast }+D_{\mathit{k}}^{\left(2\right)}{D}_{\mathit{k}}^{\left(2\right)\ast }$ into

$$\begin{array}{c}\hfill D_{\mathit{k}}^{\left(1\right)}{D}_{\mathit{k}}^{\left(1\right)\ast }+D_{\mathit{k}}^{\left(2\right)}{D}_{\mathit{k}}^{\left(2\right)\ast }=D_{+\mathit{k}}^{\phantom{\ast }}{D}_{+\mathit{k}}^{\ast }+D_{-\mathit{k}}^{\phantom{\ast }}{D}_{-\mathit{k}}^{\ast },\end{array}$$

in which we introduce a new set of variables $D_{\pm \mathit{k}}$,

$$\begin{array}{cc}\hfill \frac{1}{\sqrt{2}}\left(D_{\mathit{k}}^{\left(1\right)}+i{D}_{\mathit{k}}^{\left(2\right)}\right)& ={\mathsf{\Pi }}_{-\mathit{k}}{u}_{\omega }^{\left(\mathrm{IN}\right)\phantom{\ast }}-{\Sigma }_{+\mathit{k}}{\dot{u}}_{\omega }^{\left(\mathrm{IN}\right)\phantom{\ast }}\equiv D_{+\mathit{k}},\hfill \end{array}$$

(145)

$$\begin{array}{cc}\hfill \frac{1}{\sqrt{2}}\left(D_{\mathit{k}}^{\left(1\right)}-i{D}_{\mathit{k}}^{\left(2\right)}\right)& ={\mathsf{\Pi }}_{+\mathit{k}}{u}_{\omega }^{\left(\mathrm{IN}\right)\phantom{\ast }}-{\Sigma }_{-\mathit{k}}{\dot{u}}_{\omega }^{\left(\mathrm{IN}\right)\phantom{\ast }}\equiv D_{-\mathit{k}},\hfill \end{array}$$

(146)

with

$$\begin{array}{cccc}\hfill {\mathsf{\Pi }}_{\mathit{k}}& =\frac{1}{\sqrt{2}}\left(P_{\mathit{k}}^{\left(1\right)}-i{P}_{\mathit{k}}^{\left(2\right)}\right),\hfill & \hfill {\Sigma }_{\mathit{k}}& =\frac{1}{\sqrt{2}}\left({\Sigma }_{\mathit{k}}^{\left(1\right)}+i{\Sigma }_{\mathit{k}}^{\left(2\right)}\right).\hfill \end{array}$$

(147)

Note that in general $D_{+\mathit{k}}^{\ast }\ne D_{-\mathit{k}}^{\phantom{\ast }}$. Similarly we find:

$$\begin{array}{cccc}\hfill D_{\mathit{k}}^{\left(1\right)2}+D_{\mathit{k}}^{\left(2\right)2}& =2D_{+\mathit{k}}^{\phantom{\ast }}{D}_{-\mathit{k}}^{\phantom{\ast }},\hfill & \hfill D_{\mathit{k}}^{\left(1\right)\ast 2}+D_{\mathit{k}}^{\left(2\right)\ast 2}& =2D_{+\mathit{k}}^{\ast }{D}_{-\mathit{k}}^{\ast }.\hfill \end{array}$$

(148)

Therefore Equation (144) becomes

$$\begin{array}{cc}\hfill {\mathcal{W}}_{\pm \mathit{k}}({\Sigma }_{+\mathit{k}},{\mathsf{\Pi }}_{+\mathit{k}};{\Sigma }_{-\mathit{k}},{\mathsf{\Pi }}_{-\mathit{k}};t)& ={\left(\frac{1}{\pi }\right)}^{2}exp\left\{-4\left[\left(N_{\mathit{k}}^{\left(\mathrm{OUT}\right)}+\frac{1}{2}\right)\left(D_{+\mathit{k}}^{\phantom{\ast }}{D}_{+\mathit{k}}^{\ast }+D_{-\mathit{k}}^{\phantom{\ast }}{D}_{-\mathit{k}}^{\ast }\right)\right.\right.\hfill \\ \hfill & +\left.\left.C_{\mathit{k}}^{\left(\mathrm{OUT}\right)\phantom{\ast }}{D}_{+\mathit{k}}^{\phantom{\ast }}{D}_{-\mathit{k}}^{\phantom{\ast }}+C_{\mathit{k}}^{\left(\mathrm{OUT}\right)\ast }{D}_{+\mathit{k}}^{\ast }{D}_{-\mathit{k}}^{\ast }\right]\right\}.\hfill \end{array}$$

(149)

We arrive at a Winger function in terms of a new set of variables $D_{\pm \mathit{k}}$, a superposition of the field mode ${\Sigma }_{\mathit{k}}$ and its conjugate momentum ${\mathsf{\Pi }}_{\mathit{k}}$, but they have an unambiguous sense of parity due to $D_{+\mathit{k}}^{\ast }\ne D_{-\mathit{k}}^{\phantom{\ast }}$, unlike ${\Sigma }_{\mathit{k}}$ and ${\mathsf{\Pi }}_{\mathit{k}}$. In contrast, the latter set of variables in mode $+\mathit{k}$ can be related to their counterparts in the $-\mathit{k}$ mode by complex conjugation.

This nice feature allows us to unambiguously identify the effect of coarse-graining the field. Suppose we lose track of the information regarding the $-\mathit{k}$ modes of the created particles, so we coarse-grain out the contribution of the $-\mathit{k}$ mode in the Wigner Equation (149). Then we obtain:

$$\begin{array}{cc}\hfill {\mathcal{W}}_{\mathit{k}}({\Sigma }_{\mathit{k}},{\mathsf{\Pi }}_{\mathit{k}};t)& =\int d{D}_{-\mathit{k}}^{\phantom{\ast }}d{D}_{-\mathit{k}}^{\ast }{\mathcal{W}}_{\pm \mathit{k}}({\Sigma }_{+\mathit{k}},{\mathsf{\Pi }}_{+\mathit{k}};{\Sigma }_{-\mathit{k}},{\mathsf{\Pi }}_{-\mathit{k}};t)\hfill \\ \hfill & =\frac{1}{2\pi (N_{\mathit{k}}^{\left(\mathrm{OUT}\right)}+1/2)}exp\left\{-\frac{1}{(N_{\mathit{k}}^{\left(\mathrm{OUT}\right)}+1/2)}{D}_{\mathit{k}}^{\phantom{\ast }}{D}_{\mathit{k}}^{\ast }\right\}.\hfill \end{array}$$

(150)

The resulting reduced Wigner function is independent of $C_{\mathit{k}}^{\left(\mathrm{OUT}\right)}$. It implies that the correlation between the particles of $\pm \mathit{k}$ is completely lost. Thus we expect the corresponding entropy will be a function of $N_{\mathit{k}}^{\left(\mathrm{OUT}\right)}$ only. In principle, a similar reduced-system argument can be carried out by taking the appropriate partial trace of the density matrix elements of the field (76). The eigenvalues of the resulting reduced density matrix elements can then be used to compute the associated von Neumann entropy. Both approaches will provide an equivalent description of the Gaussian system since the Wigner function and the density matrix elements are related by a Fourier transformation.

To better understand what (150) tells, and more specifically, the meaning of the variables $D_{\mathit{k}}$, we first note that (150) is essentially (141) with $C_{\mathit{k}}^{\left(\mathrm{OUT}\right)}=0$ and a rescaling of the variable $D_{\mathit{k}}^{\left(i\right)}$ by:

$$D_{\mathit{k}}^{\left(i\right)}\mapsto \frac{1}{2N_{\mathit{k}}^{\left(\mathrm{OUT}\right)}+1}{D}_{\mathit{k}}^{\left(i\right)},$$

(151)

in (141), so we expect the corresponding entropy will be given by (127) with $C_{\mathit{k}}^{\left(\mathrm{OUT}\right)}=0$ in ${\lambda }_{\mathit{k}}$, that is,

$$\begin{array}{cc}\hfill {\mathcal{S}}_{\mathit{k}}& =\left(N_{\mathit{k}}^{\left(\mathrm{OUT}\right)}+1\right)ln\left(N_{\mathit{k}}^{\left(\mathrm{OUT}\right)}+1\right)-N_{\mathit{k}}^{\left(\mathrm{OUT}\right)}ln{N}_{\mathit{k}}^{\left(\mathrm{OUT}\right)},\hfill \end{array}$$

(152)

with a corresponding particle number $N_{\mathit{k}}^{\left(\mathrm{OUT}\right)}$ for each $\mathit{k}>0$ mode. Thus the previous arguments lead to the conclusion that the field will have entropy production associated with the creation of particle pairs only when we lose track of some information embodied in the field—in this case, all the information about one partner of the pair. The nonzero entropy arises from the lack of complete information of the field, not from the creation of field quanta. Moreover, since the state of each field mode before coarse-graining is pure, the nonzero value of the von Neumann entropy (152) also indicates that the particle pair is entangled [74]. This may be alternatively seen if we decompose a two-mode squeezed state in terms of the particle number states of two modes as in (96). Though not obvious, the righthand side of (96) cannot be separated into a product of the states of $\pm \mathit{k}$ modes like $|{\Psi }_{+\mathit{k}}\rangle \otimes |{\Phi }_{-\mathit{k}}\rangle$. The von Neumann entropy provides an easy measure about this nonseparability. Thus now it is clear that what is deeper inside correlation or coherence between the particle pair is indeed quantum entanglement between them, as explicated in [74].

To find out the meaning of $D_{\mathit{k}}$, we first write the field by the out-mode functions,

$$\begin{array}{cccc}\hfill {\Sigma }_{\mathit{k}}& ={\mathfrak{a}}_{+\mathit{k}}^{\phantom{\ast }}{u}_{\omega }^{\left(\mathrm{OUT}\right)}+{\mathfrak{a}}_{-\mathit{k}}^{\ast }{u}_{\omega }^{\left(\mathrm{OUT}\right)\ast },\hfill & \hfill {\mathsf{\Pi }}_{\mathit{k}}& ={\mathfrak{a}}_{-\mathit{k}}^{\phantom{\ast }}{\dot{u}}_{\omega }^{\left(\mathrm{OUT}\right)}+{\mathfrak{a}}_{+\mathit{k}}^{\ast }{\dot{u}}_{\omega }^{\left(\mathrm{OUT}\right)\ast },\hfill \end{array}$$

(153)

where ${\mathfrak{a}}_{\mathit{k}}$ are the expansion coefficient, the c-number counterpart of the in-mode operator ${\widehat{a}}_{\mathit{k}}$. Then, since the in-mode functions and the out-mode functions are related by:

$$u_{\omega }^{\left(\mathrm{OUT}\right)}\left(t\right)={\alpha }_{\mathit{k}}{u}_{\omega }^{\left(\mathrm{IN}\right)}+{\beta }_{\mathit{k}}{u}_{\omega }^{\left(\mathrm{IN}\right)\ast },$$

(154)

we find

$$\begin{array}{cc}\hfill D_{\mathit{k}}& ={\mathsf{\Pi }}_{\mathit{k}}^{\ast }{u}_{\omega }^{\left(\mathrm{IN}\right)\phantom{\ast }}-{\Sigma }_{\mathit{k}}^{\phantom{\ast }}{\dot{u}}_{\omega }^{\left(\mathrm{IN}\right)\phantom{\ast }}=\left({\alpha }_{\mathit{k}}^{\ast }{\mathfrak{a}}_{-\mathit{k}}^{\ast }+{\beta }_{\mathit{k}}^{\phantom{\ast }}{\mathfrak{a}}_{+\mathit{k}}^{\phantom{\ast }}\right)\left(u_{\omega }^{\left(\mathrm{IN}\right)}{\dot{u}}_{\omega }^{\left(\mathrm{IN}\right)\ast }-u_{\omega }^{\left(\mathrm{IN}\right)\ast }{\dot{u}}_{\omega }^{\left(\mathrm{IN}\right)}\right).\hfill \end{array}$$

(155)

The second pair of parentheses give an i due to the Wronskian condition of the mode functions, while the expressions inside the first pair of parentheses is worth further investigation. From (88), we have:

$${\mathfrak{b}}_{-\mathit{k}}^{\ast }={\alpha }_{-\mathit{k}}^{\ast }{\mathfrak{a}}_{-\mathit{k}}^{\ast }+{\beta }_{+\mathit{k}}^{\phantom{\ast }}{\mathfrak{a}}_{\mathit{k}}^{\phantom{\ast }}={\alpha }_{\omega }^{\ast }{\mathfrak{a}}_{-\mathit{k}}^{\ast }+{\beta }_{\omega }^{\phantom{\ast }}{\mathfrak{a}}_{\mathit{k}}^{\phantom{\ast }},$$

(156)

because for our configuration the Bogoliubov coefficients are only functions of $\left|\mathit{k}\right|$. Thus we find that the variable $D_{\mathit{k}}^{\phantom{\ast }}$ actually is

$$D_{\mathit{k}}^{\phantom{\ast }}=i{\mathfrak{b}}_{-\mathit{k}}^{\ast },$$

(157)

that is, expansion coefficient of ${\Sigma }_{\mathit{k}}$ with respect to the IN mode, so the Wigner function (150) now becomes

$${\mathcal{W}}_{\mathit{k}}\left(t\right)=\frac{1}{2\pi (N_{\mathit{k}}^{\left(\mathrm{OUT}\right)}+1/2)}exp\left\{-\frac{1}{(N_{\mathit{k}}^{\left(\mathrm{OUT}\right)}+1/2)}{\mathfrak{b}}_{-\mathit{k}}^{\phantom{\ast }}{\mathfrak{b}}_{-\mathit{k}}^{\ast }\right\}.$$

(158)

The labeling of the Wigner function can be a little awkward because now it is expanded by the variables of the $-\mathit{k}$ mode.

It is interesting to note that, once we sever the correlation between the $+\mathit{k}$ and the $-\mathit{k}$ modes, the resulting Wigner function seems to be “instantaneously” stationary in the out-variables. However, it actually is not, because the associated coefficients are all time-dependent.

In summary, for the parametric processes of interest here, the field evolves unitarily, so even though particles are produced, there is no entropy production. Quantum entanglement exists between the created particle pairs. However, if some information of the field is lost resulting in some degree of nonunitarity, then the entropy of the field can increase even if the field is initially in a pure state. In particular, an entropy of the form (152) is obtained if information is lost associated with one of the particle pairs produced for each mode.

8.2. Wigner Function of the Quantum Field

It has been been shown that the statistical properties of the classical stochastic field can be described by a classical distribution functional of the field variable and its conjugate momentum over an infinite-dimensional phase space spanned by the canonical pair [1,2]. However, it is not explained in that context how a classical stochastic field can act as a full delegate of a quantum field. It will be more desirable to have a full quantum description for the statistical process of the field. Here, we provide such a formulation based on (136).

We will construct the Wigner function in terms of the full field $\varphi (\mathit{x},t)$, instead of its modes ${\varphi }_{\mathit{k}}\left(t\right)$. From (136), we have the Wigner function of all modes given by:

$$\begin{array}{cc}\hfill \prod _{\mathit{k},i}{\mathcal{W}}_{\mathit{k}}^{\left(i\right)}({\varphi }_{\mathit{k}}^{\left(i\right)},p_{\mathit{k}}^{\left(i\right)};t)& =\left(\prod _{\mathit{k}>0}{\mathcal{N}}_{\mathit{k}}^{\left(w\right)2}\right)exp\left\{-\frac{1}{2}\sum _{\mathit{k}>0}\left[{\mathsf{A}}_{\mathit{k}}^{\phantom{2}}\left({\varphi }_{\mathit{k}}^{\left(1\right)2}+{\varphi }_{\mathit{k}}^{\left(2\right)2}\right)+{\mathsf{B}}_{\mathit{k}}^{\phantom{2}}\left(p_{\mathit{k}}^{\left(1\right)2}+p_{\mathit{k}}^{\left(2\right)2}\right)\right.\right.\hfill \\ \hfill & +\left.\left.2{\mathsf{C}}_{\mathit{k}}^{\phantom{2}}\left({\varphi }_{\mathit{k}}^{\left(1\right)}{p}_{\mathit{k}}^{\left(1\right)}+{\varphi }_{\mathit{k}}^{\left(2\right)}{p}_{\mathit{k}}^{\left(2\right)}\right)\right]\right\},\hfill \end{array}$$

(159)

where we note that ${\varphi }_{\mathit{k}}^{\left(i\right)}$ and $p_{\mathit{k}}^{\left(i\right)}$ are c-number, and

$$\begin{array}{cccc}\hfill {\mathcal{N}}_{\mathit{k}}^{\left(w\right)}& =\frac{1}{2\pi \sqrt{{\mathsf{a}}_{\mathit{k}}^{\phantom{2}}{\mathsf{b}}_{\mathit{k}}^{\phantom{2}}-{\mathsf{c}}_{\mathit{k}}^2}},\hfill & \hfill {\mathsf{A}}_{\mathit{k}}^{\phantom{2}}& =\frac{{\mathsf{a}}_{\mathit{k}}}{{\mathsf{a}}_{\mathit{k}}^{\phantom{2}}{\mathsf{b}}_{\mathit{k}}^{\phantom{2}}-{\mathsf{c}}_{\mathit{k}}^2},\hfill \end{array}$$

(160)

$$\begin{array}{cccc}\hfill {\mathsf{B}}_{\mathit{k}}^{\phantom{2}}& =\frac{{\mathsf{b}}_{\mathit{k}}}{{\mathsf{a}}_{\mathit{k}}^{\phantom{2}}{\mathsf{b}}_{\mathit{k}}^{\phantom{2}}-{\mathsf{c}}_{\mathit{k}}^2},\hfill & \hfill {\mathsf{C}}_{\mathit{k}}^{\phantom{2}}& =-\frac{{\mathsf{c}}_{\mathit{k}}}{{\mathsf{a}}_{\mathit{k}}^{\phantom{2}}{\mathsf{b}}_{\mathit{k}}^{\phantom{2}}-{\mathsf{c}}_{\mathit{k}}^2}.\hfill \end{array}$$

(161)

Since following (48) and (49), we immediately have the c-number counterparts of (142) and (143) given by:

$$\begin{array}{cccc}\hfill {\varphi }_{\mathit{k}}^{\left(1\right)}& =\frac{1}{\sqrt{2}}\left({\varphi }_{\mathit{k}}^{\phantom{\ast }}+{\varphi }_{\mathit{k}}^{\ast }\right),\hfill & \hfill {\varphi }_{\mathit{k}}^{\left(2\right)}& =\frac{1}{i\sqrt{2}}\left({\varphi }_{\mathit{k}}^{\phantom{\ast }}-{\varphi }_{\mathit{k}}^{\ast }\right),\hfill \end{array}$$

(162)

$$\begin{array}{cccc}\hfill p_{\mathit{k}}^{\left(1\right)}& =\frac{1}{\sqrt{2}}\left({\pi }_{\mathit{k}}^{\ast }+{\pi }_{\mathit{k}}^{\phantom{\ast }}\right),\hfill & \hfill p_{\mathit{k}}^{\left(2\right)}& =\frac{1}{i\sqrt{2}}\left({\pi }_{\mathit{k}}^{\ast }-{\varphi }_{\mathit{k}}^{\phantom{\ast }}\right),\hfill \end{array}$$

(163)

so that we obtain

$$\begin{array}{cc}\hfill {\varphi }_{\mathit{k}}^{\left(1\right)2}+{\varphi }_{\mathit{k}}^{\left(2\right)2}& =2{\varphi }_{\mathit{k}}^{\phantom{\ast }}{\varphi }_{\mathit{k}}^{\ast }=2{\varphi }_{-\mathit{k}}^{\phantom{\ast }}{\varphi }_{-\mathit{k}}^{\ast }=2{\varphi }_{\mathit{k}}^{\phantom{\ast }}{\varphi }_{-\mathit{k}}^{\phantom{\ast }}={\varphi }_{\mathit{k}}^{\phantom{\ast }}{\varphi }_{\mathit{k}}^{\ast }+{\varphi }_{-\mathit{k}}^{\phantom{\ast }}{\varphi }_{-\mathit{k}}^{\ast },\hfill \end{array}$$

(164)

$$\begin{array}{cc}\hfill p_{\mathit{k}}^{\left(1\right)2}+p_{\mathit{k}}^{\left(2\right)2}& =2{\pi }_{\mathit{k}}^{\phantom{\ast }}{\pi }_{\mathit{k}}^{\ast },\hfill \end{array}$$

(165)

$$\begin{array}{cc}\hfill {\varphi }_{\mathit{k}}^{\left(1\right)}{p}_{\mathit{k}}^{\left(1\right)}+{\varphi }_{\mathit{k}}^{\left(2\right)}{p}_{\mathit{k}}^{\left(2\right)}& ={\varphi }_{\mathit{k}}^{\phantom{\ast }}{\pi }_{\mathit{k}}^{\phantom{\ast }}+{\varphi }_{\mathit{k}}^{\ast }{\pi }_{\mathit{k}}^{\ast }={\varphi }_{\mathit{k}}^{\phantom{\ast }}{\pi }_{\mathit{k}}^{\phantom{\ast }}+{\varphi }_{-\mathit{k}}^{\phantom{\ast }}{\pi }_{-\mathit{k}}^{\phantom{\ast }}.\hfill \end{array}$$

(166)

If the initial state is stationary, then from (66) and (67) we find:

$$\begin{array}{cccccc}\hfill {\mathsf{A}}_{+\mathit{k}}^{\phantom{2}}& ={\mathsf{A}}_{-\mathit{k}}^{\phantom{2}},\hfill & \hfill {\mathsf{B}}_{+\mathit{k}}^{\phantom{2}}& ={\mathsf{B}}_{-\mathit{k}}^{\phantom{2}},\hfill & \hfill {\mathsf{C}}_{+\mathit{k}}^{\phantom{2}}& ={\mathsf{C}}_{-\mathit{k}}^{\phantom{2}}.\hfill \end{array}$$

(167)

and then we obtain

$$\begin{array}{cc}\hfill \prod _{\mathit{k},i}{\mathcal{W}}_{\mathit{k}}^{\left(i\right)}({\varphi }_{\mathit{k}}^{\left(i\right)},p_{\mathit{k}}^{\left(i\right)};t)& =\left(\prod _{\mathit{k}}{\mathcal{N}}_{\mathit{k}}^{\left(w\right)}\right)exp\left\{-\frac{1}{2}\sum _{\mathit{k}}\left[{\mathsf{A}}_{\mathit{k}}^{\phantom{2}}{\varphi }_{\mathit{k}}^{\phantom{\ast }}{\varphi }_{\mathit{k}}^{\ast }+{\mathsf{B}}_{\mathit{k}}^{\phantom{2}}{\pi }_{\mathit{k}}^{\phantom{\ast }}{\pi }_{\mathit{k}}^{\ast }+2{\mathsf{C}}_{\mathit{k}}^{\phantom{2}}{\varphi }_{\mathit{k}}^{\phantom{\ast }}{\pi }_{\mathit{k}}^{\phantom{\ast }}\right]\right\}.\hfill \end{array}$$

(168)

Following some algebraic manipulations similar to that in [1,2], we observe that the functional “product” of the stationary functions,

$$f(\mathit{x}-\mathit{z})={\int }_{-\infty }^{\infty }{d}^{3}yg(\mathit{x}-\mathit{y})h(\mathit{y}-\mathit{z}),$$

(169)

can be written as

$$\begin{array}{cc}\hfill f(\mathit{x}-\mathit{z})={\int }_{-\infty }^{\infty }{d}^{3}yg(\mathit{x}-\mathit{y})h(\mathit{y}-\mathit{z})& =\int \frac{d^{3}k}{{\left(2\pi \right)}^3}{g}_{\mathit{k}}{h}_{\mathit{k}}{e}^{+i\mathit{k}·(\mathit{x}-\mathit{z})},\hfill \end{array}$$

(170)

and thus $f_{\mathit{k}}=g_{\mathit{k}}{h}_{\mathit{k}}$, where the spatial Fourier transformation is defined by

$$f\left(\mathit{x}\right)=\int \frac{d^{3}k}{{\left(2\pi \right)}^3}{f}_{\mathit{k}}{e}^{+i\mathit{k}·\mathit{x}}.$$

(171)

It implies that the inverse in this sense is given by

$$\begin{array}{cccccccccc}\hfill {\delta }^{\left(3\right)}(\mathit{x}-\mathit{z})& ={\int }_{-\infty }^{\infty }{d}^{3}yf(\mathit{x}-\mathit{y})g(\mathit{y}-\mathit{z}),\hfill & \hfill & \Rightarrow \hspace{1em}\hspace{1em}\hfill & \hfill 1& =f_{\mathit{k}}{g}_{\mathit{k}},\hspace{1em}\hfill & \hfill & \Rightarrow \hspace{1em}\hspace{1em}\hfill & \hfill f_{\mathit{k}}^{-1}& =\frac{1}{g_{\mathit{k}}}.\hfill \end{array}$$

(172)

That is, in this context, the notion of the inverse $g(\mathit{x}-\mathit{y})=f^{-1}(\mathit{x}-\mathit{y})$ of a function $f(\mathit{x}-\mathit{y})$ does not mean the usual inverse function, nor $1/f(\mathit{x}-\mathit{y})$. Instead, their Fourier transforms obey the standard rule of the multiplication.

Therefore, the convolution relation for $\varphi (\mathit{x},t)$ reduces to

$$\begin{array}{cc}\hfill {\int }_{-\infty }^{\infty }{d}^{3}x{\int }_{-\infty }^{\infty }{d}^{3}y\varphi (\mathit{x},t)\mathsf{A}(\mathit{x}-\mathit{y};t)\varphi (\mathit{y},t)& =\int \frac{d^{3}k}{{\left(2\pi \right)}^3}{\varphi }_{\mathit{k}}^{\phantom{\ast }}\left(t\right){\mathsf{A}}_{\mathit{k}}\left(t\right){\varphi }_{\mathit{k}}^{\ast }\left(t\right),\hfill \end{array}$$

(173)

On the other hand, the convolution relations are different for $\pi (\mathit{x},t)$ due to (28). For example we find

$$\begin{array}{cc}\hfill {\int }_{-\infty }^{\infty }{d}^{3}x{\int }_{-\infty }^{\infty }{d}^{3}y\varphi (\mathit{x},t)\mathsf{C}(\mathit{x}-\mathit{y};t)\pi (\mathit{y},t)& =\int \frac{d^{3}k}{{\left(2\pi \right)}^3}{\varphi }_{\mathit{k}}^{\phantom{\ast }}\left(t\right){\mathsf{C}}_{-\mathit{k}}^{\phantom{\ast }}\left(t\right){\pi }_{\mathit{k}}^{\phantom{\ast }}\left(t\right)\hfill \\ \hfill & =\int \frac{d^{3}k}{{\left(2\pi \right)}^3}{\varphi }_{\mathit{k}}^{\phantom{\ast }}\left(t\right){\mathsf{C}}_{\mathit{k}}^{\phantom{\ast }}\left(t\right){\pi }_{\mathit{k}}^{\phantom{\ast }}\left(t\right),\hfill \end{array}$$

(174)

and

$$\begin{array}{cc}\hfill {\int }_{-\infty }^{\infty }{d}^{3}x{\int }_{-\infty }^{\infty }{d}^{3}y\pi (\mathit{x},t)\mathsf{B}(\mathit{x}-\mathit{y};t)\pi (\mathit{y},t)& =\int \frac{d^{3}k}{{\left(2\pi \right)}^3}{\pi }_{\mathit{k}}^{\phantom{\ast }}\left(t\right){\mathsf{B}}_{\mathit{k}}^{\phantom{\ast }}\left(t\right){\pi }_{\mathit{k}}^{\ast }\left(t\right).\hfill \end{array}$$

(175)

This implies that the Wigner function of the full field will be given by

$$\begin{array}{cc}\hfill \mathcal{W}(\varphi ,\pi ;t)& =\prod _{\mathit{k},i}{\mathcal{W}}_{\mathit{k}}^{\left(i\right)}({\varphi }_{\mathit{k}}^{\left(i\right)},p_{\mathit{k}}^{\left(i\right)};t)\hfill \\ \hfill & ={\mathcal{N}}^{\left(w\right)}exp\left\{-\frac{1}{2}{\int }_{-\infty }^{\infty }{d}^{3}x{\int }_{-\infty }^{\infty }{d}^{3}y\left[\varphi (\mathit{x},t)\mathsf{A}(\mathit{x}-\mathit{y};t)\varphi (\mathit{y},t)\right.\right.\hfill \\ \hfill & +\left.\left.\pi (\mathit{x},t)\mathsf{B}(\mathit{x}-\mathit{y};t)\pi (\mathit{y},t)+2\varphi (\mathit{x},t)\mathsf{C}(\mathit{x}-\mathit{y};t)\pi (\mathit{y},t)\right]\right\},\hfill \end{array}$$

(176)

where we have used the fact that the covariance matrix elements

$$\begin{array}{cccccc}\hfill {\mathsf{b}}_{\mathit{k}}& =\frac{1}{2}\langle \left\{{\widehat{\varphi }}_{\mathit{k}}^{\phantom{†}},{\widehat{\varphi }}_{\mathit{k}}^{†}\right\}\rangle ,\hfill & \hfill {\mathsf{a}}_{\mathit{k}}& =\frac{1}{2}\langle \left\{{\widehat{\pi }}_{\mathit{k}}^{\phantom{†}},{\widehat{\pi }}_{\mathit{k}}^{†}\right\}\rangle ,\hfill & \hfill {\mathsf{c}}_{\mathit{k}}& =\frac{1}{2}\left[\frac{1}{2}\langle \left\{{\widehat{\varphi }}_{\mathit{k}},{\widehat{\pi }}_{\mathit{k}}\right\}\rangle +\frac{1}{2}\langle \left\{{\widehat{\varphi }}_{\mathit{k}}^{†},{\widehat{\pi }}_{\mathit{k}}^{†}\right\}\rangle \right],\hfill \end{array}$$

are all real. The normalization constant ${\mathcal{N}}^{\left(w\right)}$ will be determined by the functional integral,

$$\int \mathcal{D}\varphi \mathcal{D}\pi \mathcal{W}(\varphi ,\pi ;t)=1.$$

(177)

The “inverse” of the coefficient functions $\mathsf{A}(\mathit{x}-\mathit{y};t)$, $\mathsf{B}(\mathit{x}-\mathit{y};t)$ and $\mathsf{C}(\mathit{x}-\mathit{y};t)$ take simpler forms than the coefficient functions themselves because

$$\begin{array}{cccccc}\hfill {\mathsf{A}}_{\mathit{k}}^{-1}& ={\mathsf{b}}_{\mathit{k}}-{\mathsf{c}}_{\mathit{k}}^{\phantom{-1}}{\mathsf{a}}_{\mathit{k}}^{-1}{\mathsf{c}}_{\mathit{k}}^{\phantom{-1}},\hfill & \hfill {\mathsf{B}}_{\mathit{k}}^{-1}& ={\mathsf{a}}_{\mathit{k}}-{\mathsf{c}}_{\mathit{k}}^{\phantom{-1}}{\mathsf{b}}_{\mathit{k}}^{-1}{\mathsf{c}}_{\mathit{k}}^{\phantom{-1}},\hfill & \hfill {\mathsf{C}}_{\mathit{k}}^{-1}& ={\mathsf{b}}_{\mathit{k}}-{\mathsf{a}}_{\mathit{k}}^{\phantom{-1}}{\mathsf{c}}_{\mathit{k}}^{-1}{\mathsf{b}}_{\mathit{k}}^{\phantom{-1}}.\hfill \end{array}$$

(178)

Then, we find

$$\begin{array}{cc}\hfill {\mathsf{A}}^{-1}(\mathit{x}-\mathit{y};t)& =\int \frac{d^{3}k}{{\left(2\pi \right)}^3}{\mathsf{A}}_{\mathit{k}}^{-1}{e}^{+i\mathit{k}·(\mathit{x}-\mathit{y})}=\int \frac{d^{3}k}{{\left(2\pi \right)}^3}\left[{\mathsf{b}}_{\mathit{k}}-{\mathsf{c}}_{\mathit{k}}^{\phantom{1}}{\mathsf{a}}_{\mathit{k}}^{-1}{\mathsf{c}}_{\mathit{k}}^{\phantom{1}}\right]e^{+i\mathit{k}·(\mathit{x}-\mathit{y})}\hfill \end{array}$$

$$\begin{array}{cc}\hfill & =\mathsf{b}(\mathit{x}-\mathit{y};t)-\int d^{3}u{d}^{3}v\mathsf{c}(\mathit{x}-\mathit{u};t){\mathsf{a}}^{-1}(\mathit{u}-\mathit{v};t)\mathsf{c}(\mathit{v}-\mathit{y};t),\hfill \end{array}$$

(179)

$$\begin{array}{cc}\hfill {\mathsf{B}}^{-1}(\mathit{x}-\mathit{y};t)& =\mathsf{a}(\mathit{x}-\mathit{y};t)-\int d^{3}u{d}^{3}v\mathsf{c}(\mathit{x}-\mathit{u};t){\mathsf{b}}^{-1}(\mathit{u}-\mathit{v};t)\mathsf{c}(\mathit{v}-\mathit{y};t),\hfill \end{array}$$

(180)

$$\begin{array}{cc}\hfill {\mathsf{C}}^{-1}(\mathit{x}-\mathit{y};t)& =\mathsf{c}(\mathit{x}-\mathit{y};t)-\int d^{3}u{d}^{3}v\mathsf{a}(\mathit{x}-\mathit{u};t){\mathsf{c}}^{-1}(\mathit{u}-\mathit{v};t)\mathsf{b}(\mathit{v}-\mathit{y};t),\hfill \end{array}$$

(181)

in the sense of functional multiplication. In deriving the above results we have used

$$\begin{array}{cc}\hfill \mathsf{b}(\mathit{x}-\mathit{y};t)=\int \frac{d^{3}k}{{\left(2\pi \right)}^3}{\mathsf{b}}_{\mathit{k}}\left(t\right)e^{+i\mathit{k}·(\mathit{x}-\mathit{y})}& =\int \frac{d^{3}k}{{\left(2\pi \right)}^3}\frac{1}{2}\langle \left\{{\widehat{\varphi }}_{\mathit{k}}^{\phantom{†}}\left(t\right),{\widehat{\varphi }}_{\mathit{k}}^{†}\left(t\right)\right\}\rangle e^{+i\mathit{k}·(\mathit{x}-\mathit{y})}\hfill \\ \hfill & =\frac{1}{2}\langle \left\{\widehat{\varphi }(\mathit{x},t),\widehat{\varphi }(\mathit{y},t)\right\}\rangle ,\hfill \end{array}$$

(182)

due to (38), and similarly

$$\mathsf{a}(\mathit{x}-\mathit{y};t)=\int \frac{d^{3}k}{{\left(2\pi \right)}^3}{\mathsf{a}}_{\mathit{k}}\left(t\right)e^{+i\mathit{k}·(\mathit{x}-\mathit{y})}=\frac{1}{2}\langle \left\{\widehat{\pi }(\mathit{x},t),\widehat{\pi }(\mathit{y},t)\right\}\rangle ,$$

and

$$\begin{array}{cc}\hfill \mathsf{c}(\mathit{x}-\mathit{y};t)=\int \frac{d^{3}k}{{\left(2\pi \right)}^3}{\mathsf{c}}_{\mathit{k}}\left(t\right)e^{+i\mathit{k}·(\mathit{x}-\mathit{y})}& =\int \frac{d^{3}k}{{\left(2\pi \right)}^3}\frac{1}{2}\left[\frac{1}{2}\langle \left\{{\widehat{\varphi }}_{\mathit{k}},{\widehat{\pi }}_{\mathit{k}}\right\}\rangle +\frac{1}{2}\langle \left\{{\widehat{\varphi }}_{\mathit{k}}^{†},{\widehat{\pi }}_{\mathit{k}}^{†}\right\}\rangle \right]e^{+i\mathit{k}·(\mathit{x}-\mathit{y})}\hfill \\ \hfill & =\frac{1}{2}\langle \left\{\widehat{\varphi }(\mathit{x},t),\widehat{\pi }(\mathit{y},t)\right\}\rangle ,\hfill \end{array}$$

(183)

owing to (40). Again we remind that ${\mathsf{A}}^{-1}(\mathit{x}-\mathit{y};t)$ does not represent the reciprocal of $\mathsf{A}(\mathit{x}-\mathit{y};t)$ per se.

We therefore see that the distribution function of the classical stochastic field is interpreted as the Wigner function (176) of the quantum scalar field. This may not come as a surprise. For the Gaussian system we consider here, the Wigner function is positive definite because the covariance matrix elements ${\mathsf{a}}_{\mathit{k}}$ and ${\mathsf{b}}_{\mathit{k}}$ are positive and they satisfy ${\mathsf{a}}_{\mathit{k}}{\mathsf{b}}_{\mathit{k}}-{\mathsf{c}}_{\mathit{k}}^2>0$. Thus, it can serve as a classical probability distribution. Although this seems to justify the results obtained from classical stochastic field theory as in [1,2], the latter theory is not quite self-contained. It still needs inputs from the two-point functions of the quantum field to yield the stochastic counterparts. This bridging protocol is the source of persistent confusion in the literature so we to provide a pure quantum field theory approach. Thus, even though both look identical formally, the frameworks they are base upon are drastically distinct. Other than the lurid difference of quantum non-commutativity (In the classical stochastic field approach, different distribution functions may be needed to account for different operator orderings in quantum theory. This may be seen from the fact that in quantum theory the Wigner function corresponds to a fully symmetrized ordering [111,112,113]). Since in the current case, the quantities of interest can be expressed in terms of the covariance matrix elements, which are the expectation values of the canonical operators in symmetrized ordering, that may explain why the stochastic field approach may work, there are a few additional subtleties. The Wigner functions for the more general non-Gaussian configurations are not positive-definite, they carry more information than what is in the corresponding classical probability distributions. Classical field theory cannot account for quantum entanglement, in which the nature of quantum state plays an important role. Moreover, the Shannon entropy associated with classical distributions has the property of monotonicity. That is, given the combined systems A and B, we have ${\mathcal{S}}_s[A+B]\ge max\{{\mathcal{S}}_s\left[A\right],{\mathcal{S}}_s\left[B\right]\}$ for the Shannon entropy, where ${\mathcal{S}}_s\left[A\right]$ denotes the Shannon entropy of the subsystem A. On the other hand, for the von Neumann entropy $\mathcal{S}$, the closest property to the classical monotonicity is the theorem of Araki–Lieb triangle inequality:

$$\left|\mathcal{S}\right[A]-\mathcal{S}[B\left]\right|\le \mathcal{S}[A+B]\le \mathcal{S}\left[A\right]+\mathcal{S}\left[B\right].$$

(184)

The second inequalities is known as the subadditivity inequality for von Neumann entropy, and holds with equality if and only if systems A and B are uncorrelated, that is, ${\varrho }^{AB}={\varrho }^A\otimes {\varrho }^B$. In the case we discussed earlier, let the $+\mathit{k}$ mode be the subsystem A, and the $-\mathit{k}$ mode the subsystem B. Then we have shown $\mathcal{S}[A,B]=0$ but $\mathcal{S}\left[A\right]=\mathcal{S}\left[B\right]>0$. The entropy of a subsystem is larger than the entropy of the combined system. These common properties are now often invoked in the discussions of entanglement entropy.

9. Discussion and Summary

In this concluding section we begin with a short description of how to apply what we have learned from the analysis of the nonequilibrium dynamics of squeezed quantum systems to the cosmological perturbation problem. We can only provide some useful pointers here, leaving a fuller implementation to subsequent work. A summary of this paper follows, addressing the objective, the methodology and our findings in some key issues.

9.1. Return to Cosmological Perturbations

A good portion of our effort in this work is devoted to the analysis of the evolution of a quantum field under a time-dependent parametric drive (squeezing) and how these results can be applied to the study of particle production and entropy generation as well as addressing the issues of quantum correlations, coherence and entanglement in flat space. At first sight it seems too remotely related to those of the cosmological perturbations in an inflationary universe. The connection between the parametric scalar field in flat space and the cosmological perturbations in (near) de Sitter space lies in (16), where, in terms of the Mukhanov–Sasaki variable, the cosmological perturbations follow an equation of motion of the exact form as for a parametric field in flat space. (Not even a conformal transformation is needed, and in flat space there is no distinction between conformal time and ordinary time). The gravitational effects of the cosmological spacetime is packaged in the function $z\left(\eta \right)$, as can be seen in (16). Thus our results can be readily used for the Mukhanov–Sasaki variable. Since neither the scalar field perturbation $\delta \varphi$ and scalar component of the metric perturbation $\Phi$ is gauge or coordinate independent, the Mukhanov–Sasaki variable, a linear combination of them according to (15), provides a gauge-invariant description for the cosmological perturbations. This allows us to take a more direct route in applying the present results to the cosmological perturbation problem. We may choose a convenient gauge, such as the spatially-flat gauge [42,61], in which the scalar sector of the metric perturbations is completely turned off. In this gauge the Mukhanov–Sasaki variable for the cosmological perturbations becomes proportional to the field perturbation only, $u(\mathit{x},\eta )=a\left(\eta \right)\delta \varphi (\mathit{x},\eta )$. Thus all the results we have presented on entropy generation and quantum entanglement can be readily applied to both $u(\mathit{x},\eta )$ and $\delta \varphi (\mathit{x},\eta )$, with the same physical contents in the context of cosmological perturbations. The minor difference only lies in the functional form of the frequency modulation.

9.2. Objective: Entropy of Quantum Fields and Cosmological Perturbations

We now provide a summary of the objective, our method and the issues explored.

A central subject of interest in cosmology is structure formation. The classical theory began with Lifshitz’s gravitational perturbation theory in 1946. Another important subject is particle creation from the quantum vacuum, believed to be abundant at the Planck time, first explored by Parker in 1966. Both aspects enter in theories of cosmological perturbations from inflationary cosmology, in 1982, where a quantum field is supposed to drive the universe into inflation and its fluctuations engender structures.

The entropy budget of the universe generated from various sources is also an essential concern in cosmology. Main focus was placed on particle processes at various stages of the cosmic evolution. Entropy associated with cosmological perturbations met with a new level of challenge after inflation took center stage. A necessary ingredient is to understand the entropy of quantum field processes, for free [4] and interacting fields [5]. That began in 1986. Incorporating both the quantum theory of gravitational perturbations [42] and the ideas about the entropy of quantum fields into a theory for the entropy of cosmological perturbations took shape in 1993 in the work of [1,2].

Our present work attempts to provide a synthesis of all the essential elements on this subject matter since that time, constructed from a more rigorous and comprehensive theoretical framework based on the nonequilibrium quantum field theory of squeezed quantum systems. We highlight its advantages as follows.

9.3. Methodology: Nonequilibrium Dynamics of Squeezed Quantum Systems

The technical core of this theory as applied to the present problem was presented in Section 5, Section 6 and Section 8 while the results meeting our objective, in Section 7. The systematic theoretical development of the formalism can be found in [26,108,109,114]. Applying this theory for squeezed open quantum systems, as long as they are Gaussian, enables us to obtain exact solutions for the fully nonequilibrium time-evolution of the cosmological perturbations, up to the Gaussian level.

We mention a few specific advantages: Essentially all quantities of interest can be expressed in terms of two fundamental solutions of the equation of motion of the system plus initial conditions. This make it computationally efficient, and easier to identify the structures of the models, compared with some earlier approaches, say [11]. With an exact formulation one can leave the necessary mathematical approximations or physical procedures, such as the regularization of the stress tensor associated produced particles’ density to the very end of the calculations where the actual expression is finally need.

A very useful tool special to Gaussian systems is the covariance matrix, whose power lies in the fact that the calculations involving the infinite-dimensional density matrix elements can be reduced to those of a handful of covariance matrix elements. The authors of [19] have used them to compute the von Neumann entropy of the field, but they can be also be used to form the symplectic invariants [115] and the observables of the Gaussian systems. These physical quantities are very useful for various applications in the context of open quantum systems like dynamical equilibration [116,117], quantum entanglement [118], nonequilibrium quantum thermodynamics [114,119].

In this work, we follow the same spirit. We first express the quantities of interest, not just limited to the von Neumann entropy, in terms of the covariance matrix elements. We then will express these elements in terms of the aforementioned fundamental solutions. One can break down a complex calculation to modular forms, which makes the comparisons of results from different theories easier to identify, and be able to go further as needed, sailing along guided by the formalism.

9.4. Issues: Quantum Correlations, Coherence, Phases

Entropy generation from quantum cosmological perturbations, as explained above, is weighted upon the issue of entropy of quantum fields associated with vacuum particle creation. Cosmological expansion having the same effect on a quantum state as squeezing, the issue becomes that of entropy generation from a squeezed quantum system. This is a good entry point to dissect the evolving explanations of entropy generation from quantum cosmological perturbations. We also see this as a good intellectual exercise which can reveal the subtleties of this issue, as it bears on the concepts of quantum correlations, coherence and what (not) to do with the phase.

Intuitively one may easily attribute entropy generation to the particle number from pair creation, as a result of the parametric amplification of the vacuum. However, an actual calculation shows that, starting in a vacuum state, the field will end up in the two-mode squeezed vacuum stage after the parametric process, which is a pure state and no entropy change. The authors of [4] pointed out that if one chooses to describe the particle creation process in a Fock representation and totally ignores their phase relations one would come up with this conclusion. Entropy is ‘generated’ because one choose to ignore certain information of the field which could be just as important. How one deals with the phase information is a different story. A common argument is to appeal to the random phase approximation (RHA) which assumes that rapid phase oscillation tends to smear out its own footprints. However the relevant phases are rather illusive. Let us examine this issue more closely.

Phase information can be placed in two categories: (1) the phase acquired during the evolution residing in the changes of the rotation phase plus the squeeze phase read off from the Bogoliubov coefficients; (2) initial phase of the squeeze parameter or the state. The former is fixed and dynamically determined, but the latter is largely unknown and in principle cannot be canceled as such. Thus the information of the latter is mostly inaccessible. The authors of [85] concretized this by writing down the quantum Vlasov equations of the particle number density and the coherence of the created particles, in an attempt to establish the connection between quantum coherence and the particle creation process which is non-Markovian. They argued that a nonzero entropy production associated with particle creation can emerge, if the variables associated with rapid phases are either plainly discarded or averaged out. This strategy is also adopted in [1,2,3]. However, this way to resolve the issue is not completely satisfactory in the sense that: (1) both the particle number density and coherence may be oscillatory in time for the parametric process, so even though the coherence may oscillate more rapidly than the particle number does, the difference is not always dramatic enough to justify discarding the coherence; (2) The authors of [74] show an example that coherence does not always oscillate. In particular, in the inflation setting, the coherence actually becomes non-oscillatory after the modes of the cosmological perturbations crossed the Hubble horizon; and (3) the low momentum modes have slow oscillation to make cancellation less effective. In addition, resorting to phase cancellation over all modes does not apply here because the emergent entropy occurs in each mode.

There is no denying that coherence tends to cancel because it is not positive definite as is the particle number density. Nonetheless it seems more reasonable to go after the essence of the coherence rather than its phases literally. The authors of [74] further argued that behind the coherence between pair-produced particles, more precisely, there is quantum entanglement between them at play. They showed the clear presence of entanglement by the von Neumann entropy of the reduced system after the degree of freedom of one party in the created particle pair is completely traced out. One can treat this as the emergent entropy associated with particle creation once one completely ignores or lose track of the information of one party in the particle pair. This we view as a better explanation of the root cause of emergent entropy associated with particle creation, and by extension, with cosmological perturbations involving quantum fields like in the inflationary universe.

We regained the result in [74] using a totally different formulation, based on the nonequilibrium dynamics of squeezed quantum fields. A quantum description is necessary for addressing the foundational issues of quantum decoherence and quantum entanglement, and a nonequilibrium dynamics formulation is necessary for time-dependent situations such as in a cosmological evolution. In relation to this we showed that a classical stochastic description of the quantum field often used, such as in [1,2], where quantum observables are believed to be obtainable by taking suitable ensemble averages of a classical stochastic field, though attractively simple and works in some aspects (A classical stochastic theory can be used to calculate certain quantum features of a Gaussian system, but it can not describe them in full. For example, once the classical stochastic theory acquires the two-point functions from the quantum field inputs, it can be used to calculate the covariance matrix elements, from which the entanglement measure can be computed. However there is no way for a classical stochastic theory to formulate, say, an entangled state, because it cannot account for the amplitude. In the end one still needs inputs from the quantum theory), proves inadequate to capture important quantum attributes fundamental to quantum information issues. The difficulties of attempting to substitute a quantum field by classical stochastic fields become almost insurmountable when one confronts nonGaussian systems. There, the Wigner function, known to possess the full information contained in a density matrix, is non-positive definite and can no longer be interpreted as a probability distribution. That is where the more intricate quantum features of the system reside.