# Nonequilibrium Schwinger-Keldysh formalism for density matrix states:
analytic properties and implications in cosmology

Andrei O. Barvinsky<EMAIL_ADDRESS>Theory Department, Lebedev Physics
Institute, Leninsky Prospect 53, Moscow 119991, Russia  Institute for
Theoretical and Mathematical Physics, Moscow State University, 119991,
Leninskie Gory, GSP-1, Moscow, Russia  Nikita Kolganov
<EMAIL_ADDRESS>Moscow Institute of Physics and Technology,
141700, Institutskiy pereulok, 9, Dolgoprudny, Russia  Institute for
Theoretical and Mathematical Physics, Moscow State University, 119991,
Leninskie Gory, GSP-1, Moscow, Russia  Theory Department, Lebedev Physics
Institute, Leninsky Prospect 53, Moscow 119991, Russia

###### Abstract

Motivated by cosmological Hartle-Hawking and microcanonical density matrix
prescriptions for the quantum state of the Universe we develop Schwinger-
Keldysh in-in formalism for generic nonequilibrium dynamical systems with the
initial density matrix. We build the generating functional of in-in Green’s
functions and expectation values for a generic density matrix of the Gaussian
type and show that the requirement of particle interpretation selects a
distinguished set of positive/negative frequency basis functions of the wave
operator of the theory, which is determined by the density matrix parameters.
Then we consider a special case of the density matrix determined by the
Euclidean path integral of the theory, which in the cosmological context can
be considered as a generalization of the no-boundary pure state to the case of
the microcanonical ensemble, and show that in view of a special reflection
symmetry its Wightman Green’s functions satisfy Kubo-Martin-Schwinger
periodicity conditions which hold despite the nonequilibrium nature of the
physical setup. Rich analyticity structure in the complex plane of the time
variable reveals the combined Euclidean-Lorentzian evolution of the theory,
which depending on the properties of the initial density matrix can be
interpreted as a decay of a classically forbidden quantum state.

_To the memory of Jim Hartle_

###### Contents

  1. 1 Introduction
  2. 2 Summary of main results
    1. 2.1 Schwinger-Keldysh technique for models with density matrix state
    2. 2.2 Euclidean density matrix
    3. 2.3 Reflection symmetry and analyticity properties
  3. 3 Preliminaries
    1. 3.1 Condensed notations
    2. 3.2 Canonical formalism
    3. 3.3 The solution of Dirichlet and Neumann boundary value problems
    4. 3.4 The relation between Dirichlet and Neumann Green’s functions
    5. 3.5 Canonical quantization
    6. 3.6 Bogoliubov transformations
    7. 3.7 Fock space and the coherent states
  4. 4 Generating functional in the path integral formalism
    1. 4.1 Gaussian states
    2. 4.2 In-in boundary value problem
    3. 4.3 Neumann type basis functions and Green’s function representation
    4. 4.4 Keldysh rotation
    5. 4.5 Special choice of basis functions and particle interpretation
    6. 4.6 Euclidean density matrix state
    7. 4.7 Analytic continuation and KMS condition
  5. 5 Simple applications
    1. 5.1 Harmonic oscillator
    2. 5.2 General one-dimensional system
    3. 5.3 The case of a pure state: vacuum no-boundary wavefunction
  6. 6 Discussion and conclusions
  7. A Inversion of matrices
  8. B Derivation of Eq. (4.32)
  9. C Derivation of Eq. ()
  10. D Properties of Gaussian density matrices

##  1 Introduction

The purpose of this paper is to construct the Schwinger-Keldysh in-in
formalism [1, 2] for expectation values and correlation functions in a rather
generic non-equilibrium system with the initial state in the form of a special
density matrix. This density matrix is itself assumed to be determined by the
dynamical content of the system. The motivation for this construction comes
from the scope of ideas of quantum cosmology suggesting that the initial state
of the Universe should be prescribed not from some ad hoc and freely variable
initial conditions like in a generic Cauchy problem, but rather intrinsically
fixed by the field theory model of the Universe. The pioneering implementation
of these ideas was the prescription of the Harle-Hawking no-boundary
cosmological wavefunction [3, 4], _no-boundary_ connotation indicating the
absence of the notion of the initial Cauchy (boundary) surface of spacetime.
Such a pescription replaces the existence of this surface by the requirement
of regularity of all fields at all spacetime points treated in the past as
regular internal points of spacetime manifold.

Applied to a wide class of spatially closed cosmological models this
prescription qualitatively leads to the picture of expanding Friedmann
Universe with the Lorentzian signature spacetime nucleating from the domain of
a Euclidean space with the topology of a 4-dimensional hemisphere, the
Euclidean and Lorentzian metrics being smoothly matched by analytical
continuation in the complex plane of time coordinate. This picture allows one
to avoid initial singularity in the cosmological evolution and, in particular,
serves as initial conditions for inflationary scenarios. This is because it
implies a pure vacuum state of quantum matter perturbations on top of a quasi-
exponentially expanding metric background, both the background and this vacuum
state being generated by tunneling from the classically forbidden
(underbarrier) state of the Universe, described by the Euclidean spacetime
with the imaginary time. Correlation functions of quantum cosmological
perturbations in this vacuum state have a good fit to nearly flat red-tilted
primordial spectrum of the cosmic microwave background radiation (CMBR) [5, 6]
and other features of the observable large scale structure of the Universe [7,
*Planck2].

Limitation of this no-boundary concept consists in the fact that it covers
only the realm of pure quantum states. Moreover, it prescribes a particular
quantum state which in the lowest order of the perturbation theory yields a
special vacuum state. In fact, the idea of Hartle-Hawking no-boundary initial
conditions came from the understanding that the vacuum state wavefunction
$\varPsi[\varphi(\mathbf{x})]$ of a generic free fields model in flat
spacetime can be built by the path integral over the field histories
$\phi(\tau,\mathbf{x})$ on a half-space interpolating between a given
3-dimensional configuration $\varphi(\mathbf{x})$ on the boundary plane of
$\tau=0$ and the vanishing value of these fields at the Euclidean time
$\tau\to-\infty$. Beyond perturbation theory, in the models with a bounded
from below spectrum of their Hamiltonian this procedure yields the lowest
energy eigenstate. Thus, the Hartle-Hawking no-boundary wavefunction is the
generalization of this distinguished state to a special case of curved
spatially closed spacetime, which can be formulated even though the notion of
nontrivially conserved energy does not exist for such a situation.

A natural question arises how to generalize this picture to the physical setup
with the density matrix replacing this distinguished pure state. The attempt
to do this encounters the problem of constructing the set of physical states
$|\psi\rangle$ along with the set of their weights $w_{\psi}$ participating in
the construction of the density matrix
$\hat{\rho}=\sum_{\psi}w_{\psi}|\psi\rangle\langle\psi|$. This problem looks
unmanageable without additional assumptions, but the simplest possible
assumption — universal microcanonical equipartition of all physical states —
allows one to write down the density matrix in a closed form provided one has
a complete set of equations which determine a full set of $|\psi\rangle$.
These are the Wheeler-DeWitt equations $\hat{H}_{\mu}|\psi\rangle=0$ which are
quantum Dirac constraints in gravity theory selecting the physical states [9],
$\mu$ being the label enumerating the full set of Hamiltonian and
diffeomorphism constraints, which includes also a continuous range of spatial
coordinates. The density matrix becomes a formal operator projector on the
subspace of these states, which can be written down as an operator delta
functions

$\hat{\rho}=\frac{1}{Z}\prod_{\mu}\delta(\hat{H}_{\mu}),$ (1.1)

the factor $Z$ being a partition function which provides the normalization
$\operatorname{tr}\hat{\rho}=1$ [10]. Important feature of this formal
projector is that a detailed construction of the delta function of
_noncommuting_ operators $\hat{H}_{\mu}$ (which form an open algebra of first
class constraints) leads to the representation of this projector in terms of
the Batalin-Fradkin-Vilkovisky or Faddeev-Popov path integral of quantum
gravity [10, 11] and makes it tractable within perturbation theory.

In contrast to the Hartle-Hawking prescription formulated exclusively in
Euclidean spacetime this density matrix expression is built within unitary
Lorentzian quantum gravity formalism [12]. Euclidean quantum gravity, however,
arises in this picture at the semiclssical level as a mathematical tool of
perturbative loop expansion. The partition function $Z$ of the density matrix
(its normalization coefficient) should be determined by the above path
integral over closed periodic histories, and the dominant semiclassical
contribution comes from the saddle points — periodic solutions of classical
equations of motion. The practice of applications to concrete cosmological
models shows, however, that such solutions do not exist in spacetime with the
Lorentzian signature, but can be constructed in Euclidean spacetime. The
deformation of the integration contour in the complex plane of both dynamical
variables and their time argument suggests that these Euclidean configurations
can be taken as a ground for a dominant contribution of semiclassical
expansion. This gives rise to the following definition of the Euclidean path
integral density matrix.

Let the classical background have at least two turning points and describe the
periodic (classically forbidden or underbarrier) motion between them in
imaginary Lorentzian time (or real Euclidean time $\tau$). Then the two-point
kernel
$\rho_{E}(\varphi_{+},\varphi_{-})=\langle\varphi_{+}|\mathinner{\hat{\rho}_{E}}|\varphi_{-}\rangle$
of the density matrix in question is defined by

$\rho_{E}(\varphi_{+},\varphi_{-})=\frac{1}{Z}\int
D\phi\,e^{-S_{E}[\phi]}\Bigr{|}_{\phi(\tau_{\pm})=\varphi_{\pm}},$ (1.2)

where $S_{E}[\phi]$ is the Euclidean action of the field perturbations
$\phi(\tau)$ on top of the given background, defined on the period of the
Euclidean time, $\tau_{-}\leq\tau\leq\tau_{+}$, the functional integration
runs over field histories interpolating between their values $\varphi_{\pm}$ —
the arguments of the density matrix kernel. $Z$ is the partition function
given by the path integral over the periodic histories with the period
$\beta=\tau_{+}-\tau_{-}$,

$Z=\int D\phi\,e^{-S_{E}[\phi]}\Bigr{|}_{\phi(\tau_{+})=\phi(\tau_{-})},$
(1.3)

providing the normalization $\operatorname{tr}\hat{\rho}_{E}=1$. Hermiticity
of this density matrix, which in view of its reality reduces to its symmetry
$\rho_{E}(\varphi_{+},\varphi_{-})=\rho_{E}(\varphi_{-},\varphi_{+})$, implies
that the background solution is a bounce that has a reflection symmetry with
respect to the middle turning point at $\tfrac{\tau_{+}+\tau_{-}}{2}$, and the
turning points $\tau_{\pm}$ are in fact identified.

Up to a normalization the expression (1.2) is the evolution operator of the
Schroedinger equation in imaginary time, $t=-i\tau$, with the quantum
Hamiltonian $\hat{H}_{S}(\tau)$ calculated on top of the _non-stationary_
background. The Hamiltonian operator here is written down in the Schroedinger
picture (which is indicated by the subscript $S$) and explicitly depends on
the Euclidean time because of this non-stationarity, so that the evolution
operator is the Dyson chronological $\tau$-ordered exponent

$\rho_{E}(\varphi_{+},\varphi_{-})={\rm
const}\times\langle\varphi_{+}|\mathrm{T}e^{-\int_{\tau_{-}}^{\tau_{+}}d\tau\,\hat{H}_{S}(\tau)}|\varphi_{-}\rangle.$
(1.4)

Because of the properties of the turning points (zero derivatives of the
background field) the Euclidean background can be smoothly matched at
$\tau_{\pm}$ with the classically allowed and _real_ background solution of
equations of motion parameterized by real Lorentzian time $t$. The evolution
of quantum perturbations on this Lorentzian branch of the background is then
driven by the unitary version of the $t$-ordered exponent (1.4)

$\hat{U}(t_{+},t_{-})=\mathrm{T}e^{-i\int^{t_{+}}_{t_{-}}dt\,\hat{H}_{S}(t)}$
(1.5)

with the Hermitian time-dependent Hamiltonian which is evaluated on this
Lorentzian background. In the cosmological context, when the spatial sections
of spacetime of $S^{3}$-topology are represented by circles of a variable
scale factor, the graphical image of the combined Euclidean-Lorentzian
evolution operator $\hat{U}(T,0)\hat{\rho}_{E}\hat{U}^{\dagger}(T,0)$ is
depicted on Fig. 1. It shows the Euclidean spacetime instanton with the
topology $R^{1}\times S^{3}$, $R^{1}=[\tau_{-},\tau_{+}]$, bounded at the
turning points $\tau_{\pm}$ by two minimal surfaces $\varSigma_{\pm}$ with a
vanishing extrinsic curvature. This instanton represents the density matrix
$\hat{\rho}_{E}$ and connects the Lorenzian spacetime branches. These branches
correspond to the unitary and anti-unitary evolution from $\varSigma_{\pm}$ in
some finite interval of the Lorentzian time $0\leq t\leq T$.111Of course, the
second Lorentzian branch could have been attached to the middle turning point
$\tfrac{\tau_{+}+\tau_{-}}{2}$ of the total period, but this reflection
asymmetric setup would correspond to the calculation of the in-out amplitude
of underbarrier tunneling through the Euclidean domain, which is not the goal
of this paper.

Figure 1: Picture of instanton representing the density matrix. Gray lines
depict the Lorentzian Universe nucleating from the instanton at the minimal
surfaces $\varSigma_{-}$ and $\varSigma_{+}$.

The pictorial representation of the cosmological partition function $Z$ in
view of cancellation of unitary evolution factors,
$\operatorname{tr}\bigl{(}\hat{U}(T,0)\hat{\rho}_{E}\hat{U}^{\dagger}(T,0)\bigr{)}=\operatorname{tr}\hat{\rho}_{E}=1$,
contains only the Euclidean part of Fig. 1. It is represented by the closed
cosmological instanton with the identified surfaces
$\varSigma_{+}=\varSigma_{-}$ and their 3-dimensional field configurations
$\varphi_{+}=\varphi_{-}$ (following from the identification of the arguments
in $\operatorname{tr}\hat{\rho}_{E}=\int
d\varphi\,\rho_{E}(\varphi,\varphi)$). The origin of this instanton having a
donut topology $S^{1}\times S^{3}$ is shown on Fig. 2.

The Euclidean space bridge incorporates the density matrix correlations
between the fields on opposite Lorentzian branches, which only vanish for the
density matrix of the pure state factorizable in the product of the
wavefunction $\varPsi(\varphi_{+})$ and its complex conjugated counterpart
$\varPsi^{*}(\varphi_{-})$,
$\rho_{E}(\varphi_{+},\varphi_{-})=\varPsi(\varphi_{+})\,\varPsi^{*}(\varphi_{-})$.
In the cosmological context this situation is depicted on Fig. 3 with two
disconnected Euclidean-Lorentzian manifolds corresponding to these factors.
Each of them corresponds to the Hartle-Hawking state, and the partition
function is based on the instanton with $S^{4}$-topology of Fig. 4. The latter
originates by glueing together two 4-dimensional hemispheres (discs
$D^{4}_{\pm}$) along their common equatorial boundary.

Figure 2: Origin of the partition function instanton from the density matrix
instanton by the procedure of gluing the boundaries $\varSigma_{+}$ and
$\varSigma_{-}$ — tracing the density matrix. Figure 3: Density matrix of the
pure Hartle-Hawking state represented by the union of two no-boundary
instantons.

So the goal of this paper is to construct the generating functional of
expectation values and correlation functions of Heisenberg operators defined
with respect to such a density matrix. Motivated by applications of quantum
cosmology, this is essentially non-equilibrium physical setup, because the
cosmological inflationary background is very non-stationary. Because of this
it raises many questions which for the impure density matrix case go
essentially beyond what is known about the Hartle-Hawking state. In
particular, despite non-equilibrium nature this pure state selects a
distinguished set of positive/negative frequency basis functions of the so-
called Euclidean vacuum which for the de Sitter metric background turns out to
be a special case of the de Sitter invariant vacuum [13, 14, 15]. But for a
density matrix case this distinguished choice is unknown and, moreover, its
reasonable particle interpretation is not granted at all to be possible.

Figure 4: Origin of the partition function instanton from the density matrix
instanton by the procedure of gluing the boundaries $\varSigma_{+}$ and
$\varSigma_{-}$ — tracing the density matrix.

The notion of the Euclidean quantum gravity density matrix was pioneered in
[16]. Then, within the concept of the above type, it was built in a concrete
inflationary model driven by the trace anomaly of Weyl invariant fields [17].
Interpreted as a microcanonical density matrix of spatially closed cosmology
[10]222This interpretation follows from the analogy with the microcanonical
ensemble whose density matrix is a projector on the subspace of fixed
conserved energy. As mentioned above, in the absence of the notion of
conserved energy the role of this projection in closed cosmology is played by
the delta function of Hamiltonian and momentum constraints — the projector on
their conserved zero value. it was later shown to be a very promising
candidate for the initial quantum state of the Universe. In particular, it
includes the existence of the quasi-thermal stage preceding the inflation
[17], provides the origin of the Higgs-type or $R^{2}$-type inflationary
scenario [18] with subplanckian Hubble scale [19] and suppresses the
contribution of Hartle-Hawking instantons to zero. Thus, this model allows one
to circumvent the main difficulty of the Hartle-Hawking prescription —
insufficient amount of inflation in the Hartle-Hawking ensemble of universes
dominated by vanishingly small values of the effective cosmological constant.
Elimination of this infrared catastrophe is, on the one hand, the quantum
effect of the trace anomaly which flips the sign of the Euclidean effective
action and sends it to $+\infty$ [17, 20]. On the other hand, this is the
hill-top nature of inflation starting from the maximum of inflaton potential
rather than from its minimum [21]. Finally, this model suggests that quantum
origin of the Universe is the subplanckian phenomenon subject to semiclassical
$1/N$-perturbation theory in the number of numerous higher-spin conformal
fields [22]. Thus, it sounds reliable even in the absence of currently
unavailable non-perturbative methods of quantum gravity.

All these conclusions have been recently reviewed in [23] including certain
preliminary results on the primordial CMBR spectra, which might even bear
potentially observable thermal imprint of the pre-inflation stage of this
model [24]. However, detailed calculation of this spectrum and of higher order
correlation functions requires the construction of the in-in Schwinger-Keldysh
formalism extended to the setup with the initial density matrix of the above
type.

Schwinger-Keldysh formalism [1, 2] was intensively applied in quantum gravity
and cosmology, and the number of publications on this subject is
overwhelmingly high, so that we briefly mention only their minor part.
Together with early applications [25, 26, 27] and the pioneering calculation
of non-gaussianities in cosmological perturbation spectra [28] these works
include the calculation of cosmological correlation functions [29, 30], the
results on cosmological singularity avoidance due to nonlocal effects [31],
equivalence of the Euclidean and in-in formalisms in de Sitter QFT [32, 33]
and even the analysis of initial conditions within Schwinger-Keldysh formalism
[34]. Among recent results one should mention the development of a special
effective field theory method based on analyticity and unitarity features of
in-in formalism [35], its applications to four-point correlators in
inflationary cosmology [36] and numerous conformal field theory and holography
ramifications of Schwinger-Keldysh technique (see, for example [37, 35] and
references therein). However, the success of these works essentially relies on
working with the model of spatially flat Universe — extension to spatially
closed cosmology with $S^{3}$-sections is likely to invalidate many of these
exact analytical results. At the same time, despite a general belief that
inflation completely washes out details of initial quantum state, learning its
imprint on the Universe requires to go beyond $K=0$ FRW model. Moreover,
recent analysis of the large scale Planck 2018 data, associated with the
Hubble tension problem in modern precision cosmology [38], testifies at more
that 99% confidence level in favor of the closed Universe preferring a
positive curvature with $K=+1$ [39, 40]. Remarkably, the model of
microcanonical initial conditions in early quantum cosmology of [17, 10]
exists only for $K=1$. Therefore, robust observational evidence in favour of a
positive spatial curvature serves as an additional motivation for this model
and justifies our goals.

Having said enough about the motivation coming from cosmology for the density
matrix modification of the in-in formalism coming from cosmology, let us
emphasize that the usefulness of this modification extends to a much wider
area. Note that the expression (1.4) for the case of a static background is
nothing but a well-known density matrix of the equilibrium canonical ensemble
at the inverse temperature $\beta=\tau_{+}-\tau_{-}$,

$\hat{\rho}=\frac{1}{Z}e^{-\beta\hat{H}}.$ (1.6)

Its evolution in time gives rise to Matsubara technique of thermal Green’s
functions [41] and thermofield dynamics [42] which satisfies nontrivial
analyticity properties in the complex plane of time including periodicity in
the direction of imaginary axis — Kubo-Martin-Schwinger (KMS) condition [43,
44]. Many of these properties depend on the condition of equilibrium and
associated with the conservation of energy. What we suggest here is the
generalization of this technique to non-equilibrium situation with the
Hamiltonian explicitly depending on time, which would be important in many
areas of quantum field theory, high energy and condensed matter physics. To
cover as wide scope of models and problems as possible we will try being
maximally generic and use condensed notations applicable in generic dynamical
systems.

In this paper we will basically consider the elements of the diagrammatic
technique for the density matrix in-in formalism. Therefore we restrict
ourselves with the systems having a quadratic action on top of the non-
stationary background subject to reflection symmetry discussed above. The one-
loop preexponential factors of this formalism will be considered elswhere.

The paper is organized as follows. Section 2 contains the summary of notations
and main results. It includes the formulation of in-in generating functional
in the generic non-equilibrium system with a Gaussian type initial density
matrix, the selection of distinguished set of positive/negative frequency
basis functions of the wave operator, determined by the density matrix
parameters, and application of this formalism to a special density matrix
based on the Euclidean path integral, this case demonstrating special
reflection symmetry, analyticity and KMS periodicity properties. Section 3
presents preliminary material of canonical quantization and the technique of
boundary value problems and relevant Green’s functions in a generic dynamical
system. Section 4 contains detailed derivation of all the results. Section 5
is devoted to the demonstration of the formalism on concrete examples, while
Section 6 contains a concluding discussion along with the prospects of future
research. Several appendices give technical details of derivations and present
certain nontrivial properties of Green’s functions and Gaussian type density
matrices.

##  2 Summary of main results

###  2.1 Schwinger-Keldysh technique for models with density matrix state

We consider a generic system with the action $S[\phi]$ quadratic in dynamical
variables $\phi=\phi^{I}(t)$, the index $I$ including both the discrete tensor
labels and in field-theoretical context also the spatial coordinates,

$\displaystyle S[\phi]=\frac{1}{2}\int
dt\Bigl{(}\dot{\phi}^{T}A\dot{\phi}+\dot{\phi}^{T}B\phi+\phi^{T}B^{T}\dot{\phi}+\phi^{T}C\phi\Bigr{)}.$
(2.1)

Here $A=A^{T}\equiv A_{IJ}$, $B\equiv B_{IJ}$ and $C=C^{T}\equiv C_{IJ}$ are
the matrices acting in the vector space of $\phi^{J}$, the superscript T
denoting the transposition, $\phi$ being a column and $\phi^{T}$ — a row (the
use of these canonical condensed notations including also spatial integration
over contracted indices $I$ will be discussed in much detail in Section 3).
What is most important throughout the paper, all these matrices are generic
functions of time $A=A(t)$, $B=B(t)$, $C=C(t)$, reflecting non-equilibrium and
non-stationary physical setup. This action will be considered as a quadratic
part of the full nonlinear action in field perturbations $\phi$ on a certain
background whose possible symmetries will be inherited by these coefficients
as certain restrictions on their time dependence. These restrictions will be
very important for the results of the paper and will be discussed below, but
otherwise this time dependence is supposed to be rather generic.

The prime object of our interest will be the Schwinger-Keldysh generating
functional of the in-in expectation values and correlation functions of
Heisenberg operators in the physical state described by the initial density
matrix $\hat{\rho}$. This is the functional of two sources

$Z[J_{1},J_{2}]=\operatorname{tr}\left[\hat{U}_{J_{1}}(T,0)\mathinner{\hat{\rho}}\hat{U}^{\dagger}_{-J_{2}}(T,0)\right].$
(2.2)

Here the trace is taken over the Hilbert space of the canonically quantized
field $\hat{\phi}$ and $\hat{U}_{J}(T,0)$ is the operator of unitary evolution
from $t=0$ to $t=T$ with the time dependent Hamiltonian corresponding to the
action (2.1) and modified by the source term $-J^{T}(t)\phi(t)\equiv-
J_{I}(t)\phi^{I}(t)$ with the source $J^{T}(t)=J_{I}(t)$. In the Schroedinger
picture (labelled by $S$) it reads as the chronologically ordered operator
$\mathbb{\mathrm{T}}$-exponent

$\hat{U}_{J}(T,0)=\mathbb{\mathrm{T}}e^{-i\int\limits^{T}_{0}dt\,\bigl{(}\hat{H}_{S}(t)-J(t)\hat{\phi}_{S}\bigr{)}}.$
(2.3)

We will consider the class of density matrices whose kernel in the coordinate
representation
$\langle\varphi_{+}|\mathinner{\hat{\rho}}|\varphi_{-}\rangle=\rho(\varphi_{+},\varphi_{-})$
has the following Gaussian form — exponentiated quadratic and linear forms in
$\varphi_{\pm}$,

$\displaystyle\rho({\boldsymbol{\varphi}})=\text{const}\times\exp\left\\{-\frac{1}{2}{\boldsymbol{\varphi}}^{T}{\boldsymbol{\varOmega}}\,{\boldsymbol{\varphi}}+{{\boldsymbol{j}}}^{T}\\!{\boldsymbol{\varphi}}\right\\},$
(2.4) $\displaystyle{\boldsymbol{\varphi}}=\begin{bmatrix}\,\varphi_{+}\\\
\,\varphi_{-}\end{bmatrix},\quad{{\boldsymbol{j}}}=\begin{bmatrix}[l]\,{j}_{+}\\\
\,{j}_{-}\end{bmatrix},$ (2.5)

where we assembled $\varphi_{\pm}$ into the two-component column multiplets
(denoted by boldfaced letters) ${\boldsymbol{\varphi}}$, did the same with the
coefficients ${\boldsymbol{j}}$ of the linear form and introduced the $2\times
2$ block-matrix ${\boldsymbol{\varOmega}}$ acting in the space of such two-
component multiplets

${\boldsymbol{\varOmega}}=\begin{bmatrix}[l]\;R&S\\\
\;S^{*}&R^{*}\end{bmatrix},\quad R=R^{T},\quad S=S^{\dagger}.$ (2.6)

The blocks of this matrix $R=R_{IJ}$, $S=S_{IJ}$ and their complex and
Hermitian conjugated versions, $S^{\dagger}\equiv S^{T*}$, should satisfy
these transposition and conjugation properties in order to provide Hermiticity
of the density matrix. The same concerns the “sources” $j_{\pm}$ in the
definition of ${\boldsymbol{j}}$, $j_{+}=j_{-}^{*}\equiv j$. Transposition
operation above applies also to two-component objects, so that
${\boldsymbol{\varphi}}^{T}=[\varphi_{+}^{T}\,\;\varphi_{-}^{T}]$. Below we
will denote $2\times 2$ block matrices and relevant 2-block component columns
and rows by boldfaced letters.

Such a choice of the density matrix is motivated by the fact that for a block-
diagonal ${\boldsymbol{\varOmega}}$ it reduces to a pure quasi-vacuum state,
its “source” ${\boldsymbol{j}}$ allows one to induce nonzero mean value of the
quantum field and by the differentiation with respect to ${\boldsymbol{j}}$
one can generate a much wider class of density matrices with “interaction”
terms in the exponential. Normalizability of the density matrix of course
implies that the real part of ${\boldsymbol{\varOmega}}$ should be a positive
definite matrix.

The path integral representation for the coordinate kernels of the unitary
evolution operator (2.2) allows one to rewrite the generating functional
$Z[J_{1},J_{2}]$ as the double path integral. For this purpose it is useful to
introduce the two-component notations for the histories $\phi_{1}(t)$ and
$\phi_{2}(t)$ as well as for their sources,

$\displaystyle\phi_{1},\phi_{2}\mapsto{\boldsymbol{\phi}}=\begin{bmatrix}[l]\,\phi_{1}\\\
\,\phi_{2}\end{bmatrix},\quad
J_{1},J_{2}\mapsto{\boldsymbol{J}}=\begin{bmatrix}[l]\,J_{1}\\\
\,J_{2}\end{bmatrix},$ (2.7)

In terms of these notations the generating functional reads

$Z[{\boldsymbol{J}}]=\int
D[{\boldsymbol{\phi}},{\boldsymbol{\varphi}}]\,\exp\Biggl{\\{}i{\boldsymbol{S}}[{\boldsymbol{\phi}}]+i\int_{0}^{T}dt\,{\boldsymbol{J}}^{T}{\boldsymbol{\phi}}\\\
-\frac{1}{2}{\boldsymbol{\varphi}}^{T}{\boldsymbol{\varOmega}}\,{\boldsymbol{\varphi}}+{\boldsymbol{j}}^{T}\\!{\boldsymbol{\varphi}}\Biggr{\\}},$
(2.8)

where the total action is obviously

${\boldsymbol{S}}[{\boldsymbol{\phi}}]=S[\phi_{1}]-S[\phi_{2}]$ (2.9)

with the actions $S[\phi_{1,2}]$ given by (2.1) in the integration range from
$t=0$ to $t=T$ and the total integration measure over ${\boldsymbol{\phi}}$
and ${\boldsymbol{\varphi}}$

$D[{\boldsymbol{\phi}},{\boldsymbol{\varphi}}]=\int
d\varphi_{+}\,d\varphi_{-}\\!\\!\\!\\!\\!\\!\\!\\!\int\limits_{\begin{subarray}{c}{\boldsymbol{\phi}}(0)={\boldsymbol{\varphi}}\\\
\phi_{1}(T)=\phi_{2}(T)\end{subarray}}\\!\\!\\!\\!\\!\\!\\!\\!D\phi_{1}\,D\phi_{2}.$
(2.10)

Here $d\varphi$ and $D\phi$ denote respectively the integration measures over
variables at a given moment of time and the integration measures
$D\phi=\prod_{t}d\phi(t)$ over time histories subject to indicated boundary
conditions.

Calculation of this Gaussian path integral leads to the expression

$Z[{\boldsymbol{J}}]=\text{const}\times\exp\Biggl{\\{}-\frac{i}{2}\int_{0}^{T}dt\,dt^{\prime}\,{\boldsymbol{J}}^{T}\\!(t){\boldsymbol{G}}(t,t^{\prime}){\boldsymbol{J}}(t)\\\
-\int_{0}^{T}dt\,{\boldsymbol{J}}^{T}\\!(t){\boldsymbol{G}}(t,0)\,{\boldsymbol{j}}+\frac{i}{2}{\boldsymbol{j}}^{T}{\boldsymbol{G}}(0,0)\,{\boldsymbol{j}}\Biggr{\\}},$
(2.11)

where we disregard the source-independent prefactor. The bilinear in the full
set of sources exponential is the total action in the integrand of (2.8) at
its saddle point — the point of stationarity of the action with respect to
variations of both the histories ${\boldsymbol{\phi}}(t)$ and their boundary
data ${\boldsymbol{\varphi}}$ at $t=0$. The condition of stationarity
generates the boundary value problem for this saddle point including the
linear second order equation of motion for ${\boldsymbol{\phi}}(t)$ and the
full set of boundary conditions at $t=0$ and $t=T$. This problem is posed and
solved in Section 4 in terms of its Green’s function
${\boldsymbol{G}}(t,t^{\prime})$ subject to homogeneous version of these
boundary conditions. The Green’s function has a block-matrix form typical of
Schwinger-Keldysh in-in formalism composed of the Feynman
$G_{\mathrm{T}}(t,t^{\prime})$, anti-Feynman
$G_{\bar{\mathrm{T}}}(t,t^{\prime})$ and off-diagonal Wightman Green’s
functions blocks (4.3),

${\boldsymbol{G}}(t,t^{\prime})=\begin{bmatrix}[r]G_{\mathrm{T}}(t,t^{\prime})&G_{<}(t,t^{\prime})\\\
G_{>}(t,t^{\prime})&G_{\bar{\mathrm{T}}}(t,t^{\prime})\end{bmatrix},$ (2.12)

which are related to one another by the equalities
$G_{\bar{\mathrm{T}}}(t,t^{\prime})=[G_{\mathrm{T}}(t,t^{\prime})]^{*}$ and
$G_{>}(t,t^{\prime})=G_{<}^{T}(t^{\prime},t)$ and satisfy respectively
inhomogeneous and homogeneous wave equations

$FG_{\mathrm{T},\bar{\mathrm{T}}}(t,t^{\prime})=\delta(t-t^{\prime}),\quad
FG_{\gtrless}(t,t^{\prime})=0$ (2.13)

with the operator $F$ — the Hessian of the action (2.1),
$F\delta(t-t^{\prime})=\delta^{2}S[\phi]/\delta\phi(t)\,\delta\phi(t^{\prime})$,

$F=-\frac{d}{dt}A\frac{d}{dt}-\frac{d}{dt}B+B^{T}\frac{d}{dt}+C.$ (2.14)

The block-matrix Green’s function ${\boldsymbol{G}}(t,t^{\prime})$, as is
usually done in boundary value problems, can be built in terms of the full set
of basis functions ${\boldsymbol{v}}_{\pm}$ of this operator, satisfying the
boundary conditions of the variational problem for the action in (2.8). This
will explicitly be done in Section 4, but in view of the complexity of these
boundary conditions intertwining the $\phi_{1,2}$-branches of the field space
these basis functions do not have a clear particle interpretation, that is
separation into positive and negative frequency parts. This difficulty is
caused by the convolution of problems associated, on the one hand, with the
non-equilibrium nature of a generic background (rather generic dependence of
the operator coefficients $A(t)$, $B(t)$ and $C(t)$ on time) and, on the other
hand, with the in-in physical setup involving a nontrivial density matrix.

Despite these difficulties, there exists a distinguished set of basis
functions for the wave operator which have a clear particle interpretation,
and this is one of the main results of the paper. This set is related by
Bogoliubov transformations to ${\boldsymbol{v}}_{\pm}(t)$ and is uniquely
prescribed by the full set of complex conjugated positive and negative
frequency basis functions of the operator (2.14) $v(t)$ and $v^{*}(t)$ which
satisfy the intial value problem at $t=0$,

$\displaystyle\begin{aligned} &Fv(t)=0,\\\
&(iW-\omega)v(t)\bigr{|}_{t=0}=0,\;\,(iW+\omega^{*})v^{*}(t)\bigr{|}_{t=0}=0,\end{aligned}$
(2.15)

where $W$ is what we call the _Wronskian_ operator

$W=A\frac{d}{dt}+B,$ (2.16)

which participates in the Wronskian relation for the operator $F$, which is
valid for arbitrary two complex fields $\phi_{1,2}(t)$,

$\phi_{2}^{T}F\phi_{1}-(F\phi_{2})^{T}\phi_{1}=-\frac{d}{dt}\left[\phi_{2}^{T}W\phi_{1}-(W\phi_{2})^{T}\phi_{1}\right]$
(2.17)

and, moreover, serves as the definition of the conserved (not positive-
definite) inner product in the space of solutions of the homogeneous wave
equation, $F\phi_{1,2}=0$,

$(\phi_{1},\phi_{2})=i\phi_{1}^{\dagger}(W\phi_{2})-i(W\phi_{1})^{\dagger}\phi_{2}.$
(2.18)

We will call the boundary conditions (LABEL:eq:green_lor_neum_bdy0) and
associated with them Green’s functions the Neumann ones333Strictly speaking,
these are the analogue of Robin boundary conditions, because they contain
together with the derivative transversal to the boundary also the tangential
terms composed of the coefficient $B$ and $\omega$..

Important point of the definition (LABEL:eq:green_lor_neum_bdy0) is that the
frequency matrix $\omega$ (remember that in the generic setup this is a matrix
$\omega_{IJ}$ acting in the vector space of $\phi^{J}$) is not directly
contained in the blocks of the matrix (2.6), but follows from the requirement
of the particle interpretation of the basis functions $v(t)$. This requirement
can be formulated as follows. One defines the creation-annihilation operators
$\hat{a}^{\dagger}$ and $\hat{a}$ in terms of which the Heisenberg operator
$\hat{\phi}(t)$ is decomposed as the sum of positive-negative basis functions
$v(t)$ and $v^{*}(t)$,
$\hat{\phi}(t)=v(t)\,\hat{a}+v^{*}(t)\,\hat{a}^{\dagger}$. Then there exist
non-anomalous and anomalous particle averages with respect to the density
matrix,

$\displaystyle\nu=\operatorname{tr}\bigl{[}\hat{\rho}\,\hat{a}^{\dagger}\hat{a}\bigr{]},\quad\kappa=\operatorname{tr}\bigl{[}\hat{\rho}\,\hat{a}\,\hat{a}\bigr{]},$
(2.19)

and the requirement of vanishing anomalous average $\kappa=0$ allows one to
assign the average $\nu$ the interpretation of the set of occupation numbers
associated with $\hat{\rho}$. This requirement serves as the equation for the
frequency matrix $\omega$ which, as it is shown in Section 4, can be
explicitly solved for a special case of the real matrix
${\boldsymbol{\varOmega}}$. This solution reads

$\omega=R^{1/2}\sqrt{I-\sigma^{2}}R^{1/2},\quad\sigma\equiv R^{-1/2}SR^{-1/2}$
(2.20)

and gives the expression for the occupation number matrix in terms of the
single symmetric matrix $\sigma$ after the orthogonal rotation by the
orthogonal matrix $\varkappa$,

$\displaystyle\nu=\frac{1}{2}\varkappa\left(\sqrt{\frac{I-\sigma}{I+\sigma}}-1\right)\varkappa^{T},$
(2.21)
$\displaystyle\varkappa\equiv\big{[}\omega^{1/2}R^{-1}\omega^{1/2}\big{]}^{1/2}\omega^{-1/2}R^{1/2}=\bigl{(}\varkappa^{T}\bigr{)}^{-1}.$
(2.22)

As shown in Appendix D, the existence of this particle interpretation with a
positive definite matrix $\nu$ fully matches with conditions of
normalizability, boundedness and positivity of the density matrix
incorporating positive definiteness of matrices $I\pm\sigma$ and negative
definiteness of $\sigma$.

With the normalization of these distinguished basis functions $v$ to unity

$(v_{A},v_{B})=-(v^{*}_{A},v^{*}_{B})=\delta_{AB},$ (2.23)

where $A$ is the index enumerating the full set of basis functions, the blocks
of the in-in Green’s function (2.12) take the form

$\displaystyle
iG_{\mathrm{T}}(t,t^{\prime})=v(t)\,v^{\dagger}(t^{\prime})\,\theta(t-t^{\prime})+v^{*}(t)\,v^{T}(t^{\prime})\,\theta(t^{\prime}-t)$
$\displaystyle\qquad\qquad+\,v(t)\,\nu\,v^{\dagger}(t^{\prime})+v^{*}(t)\,\nu\,v^{T}(t^{\prime}),$
(2.24) $\displaystyle
iG_{>}(t,t^{\prime})=v(t)\,\bigl{(}\nu+I\bigr{)}\,v^{\dagger}(t^{\prime})+v^{*}(t)\,\nu\,v^{T}(t^{\prime}).$
(2.25)

Here the terms of the type $v(t)\,v^{\dagger}(t^{\prime})$ should be
understood as the matrix products
$\sum_{A}v^{I}_{A}(t)\,v^{*J}_{A}(t^{\prime})$ (one should bear in mind that
the basis function $v(t)=v^{I}_{A}(t)$ represents the square (but asymmetric)
matrix whose upper indices label the field $\phi^{I}$ components, whereas the
subscript indices $A$ enumerate the basis functions in their full linear
independent set). Correspondingly,
$v(t)\,\nu\,v^{\dagger}(t^{\prime})=\sum_{A,B}v^{I}_{A}(t)\,\nu^{AB}v^{*J}_{B}(t^{\prime})$,
etc.

This form of the Green’s functions is very familiar from thermofield dynamics
for simple equilibrium condensed matter systems, when all the matrices of the
above type become diagonal in the momentum space of field modes labeled by
$A=\mathbf{p}$, $\sum_{A}=\int d^{3}\mathbf{p}/(2\pi)^{3/2}$ and
$\nu_{AB}=\nu_{\mathbf{p},\mathbf{p}^{\prime}}=(\exp(\beta\omega_{\mathbf{p}})-1)^{-1}\delta(\mathbf{p}-\mathbf{p}^{\prime})$
represents expected occupation number for Bose-Einstein statistics at inverse
temperature $\beta$ (detailed consideration of this example is presented in
Section 5). Remarkably, the occupation number picture generalizes to
nonequilibrium systems of a very general type — the function of the single
symmetric matrix in the parentheses of Eq. (2.21) can be diagonalized by extra
orthogonal rotation (additional to that of $\varkappa$), and its eigenvalues
would serve as occupation numbers in the generic nonequilibrium state with the
initial density matrix.

###  2.2 Euclidean density matrix

As discussed in Introduction, in quantum cosmology context the density matrix
itself can be given in terms of the Euclidean path integral and thus
dynamically determined by individual properties of the system including its
action functional. So we consider the path integral expression for the
Euclidean density matrix
$\langle\varphi_{+}|\mathinner{\hat{\rho}_{E}[J_{E}]}|\varphi_{-}\rangle\equiv\rho_{E}(\varphi_{+},\varphi_{-};J_{E}]$,

$\rho_{E}(\varphi_{+},\varphi_{-};J_{E}]=\frac{1}{Z}\\!\\!\\!\\!\\!\\!\\!\int\limits_{\phi(\tau_{\pm})=\varphi_{\pm}}\\!\\!\\!\\!\\!\\!\\!D\phi\,\exp\biggl{\\{}-S_{E}[\phi]\\\
-\int_{\tau_{-}}^{\tau_{+}}d\tau\,J_{E}(\tau)\phi(\tau)\biggr{\\}},$ (2.26)

where integration runs over histories $\phi(\tau)$ in Euclidean time $\tau$ on
the segment $[\tau_{-},\tau_{+}]$, interpolating between the arguments of the
density matrix $\varphi_{\pm}$. In what follows we will assume that
$\tau_{-}=0$ and $\tau_{+}=\beta$. The Euclidean action, supplied by the
Euclidean source $J_{E}(\tau)$ probing the interior of the Euclidean
spacetime, has a structure similar to the Lorentzian action and can be
obtained from (2.1) by the replacement

$t,A(t),B(t),C(t)\,\mapsto\,\tau,A_{E}(\tau),B_{E}(\tau),C_{E}(\tau).$ (2.27)

These functions of the Euclidean time are rather generic, except that they
should not contradict the basic property of the density matrix (with the
source $J_{E}$ switched off) — its Hermiticity. Sufficient conditions
providing this property,
$\rho_{E}(\varphi_{+},\varphi_{-};0]=\rho^{*}_{E}(\varphi_{-},\varphi_{+};0]$,
read

$\displaystyle\begin{aligned} &A_{E}(\beta-\tau)=A_{E}^{*}(\tau),\quad
B_{E}(\beta-\tau)=-B_{E}^{*}(\tau),\\\
&C_{E}(\beta-\tau)=C_{E}^{*}(\tau).\end{aligned}$ (2.28)

For real values of these coefficients these relations reduce to the reflection
symmetry of the action and the whole formalism relative to inversions with
respect to the center of the Euclidean segment $[0,\beta]$. Here we consider
this property as given, but it can be derived from the assumption that the
quadratic Euclidean action is built on top of the Euclidean spacetime
background — the bounce which solves full nonlinear equations of motion of the
theory and represents the periodic (underbarrier) motion of the system between
two turning points. One of these points is associated with the center of the
above Euclidean segment $\tfrac{\tau_{+}+\tau_{-}}{2}=\tfrac{\beta}{2}$, and
the other one corresponds to the (identified) points of nucleation
$\tau_{\pm}$ from the classically forbidden Euclidean regime to the Lorentzian
regime, the latter being described by the two branches of the Schwinger-
Keldysh formalism (labelled above by 1 and 2). For an equilibrium situation
with constant coefficients (2.27) at the inverse temperature
$\beta=\tau_{+}-\tau_{-}$ this setup is even simpler and corresponds to the
density matrix of the thermal canonical ensemble.

The Gaussian integration in (2.26) allows one to express the result in terms
of the Green’s function of the Hessian of the Euclidean action $F_{E}$,
$F_{E}\delta(\tau-\tau^{\prime})=\delta^{2}S_{E}[\phi]/\delta\phi(\tau)\,\delta\phi(\tau^{\prime})$,
which can be obtained from (2.14) by the replacement (2.27), subject to
Dirichlet boundary conditions,

$F_{E}\,G_{D}(\tau,\tau^{\prime})=\delta(\tau-\tau^{\prime}),\quad
G_{D}(\tau_{\pm},\tau^{\prime})=0.$ (2.29)

The resulting density matrix looks like the expression (2.4) amended by the
quadratic form in the Euclidean source,

$\rho_{E}(\varphi_{+},\varphi_{-};J_{E}]={\rm
const}\times\exp\biggl{\\{}-\frac{1}{2}{\boldsymbol{\varphi}}^{T}{\boldsymbol{\varOmega}}_{E}\,{\boldsymbol{\varphi}}+{\boldsymbol{j}}^{T}_{E}\,{\boldsymbol{\varphi}}\\\
+\frac{1}{2}\int_{0}^{\beta}d\tau\,d\tau^{\prime}J_{E}(\tau)\,G_{D}(\tau,\tau^{\prime})J_{E}(\tau^{\prime})\biggr{\\}}$
(2.30)

with the special expressions for the matrix ${\boldsymbol{\varOmega}}_{E}$ and
the source ${\boldsymbol{j}}_{E}$. The matrix ${\boldsymbol{\varOmega}}_{E}$
reads

$\displaystyle{\boldsymbol{\varOmega}}_{E}=\begin{bmatrix}[r]-{W}_{\\!E}\,G_{D}(\beta,\beta)\overset{\leftarrow}{W}_{\\!E}&\overset{\rightarrow}{W}_{\\!E}\,G_{D}(\beta,0)\overset{\leftarrow}{W}_{\\!E}\\\
\overset{\rightarrow}{W}_{\\!E}\,G_{D}(0,\beta)\overset{\leftarrow}{W}_{\\!E}&-\overset{\rightarrow}{W}_{\\!E}\,G_{D}(0,0)\overset{\leftarrow}{W}_{\\!E}\end{bmatrix}$
(2.31)

where we use the arrow to indicate the direction in which the Wronskian
operator is acting on the corresponding first or second time argument of the
Green’s function, so that for the left action the following rule holds
$\phi^{T}\\!(\tau)\overset{\leftarrow}{W}_{\\!E}\equiv(W_{E}\phi(\tau))^{T}$,
and
$\overset{\rightarrow}{W}_{\\!E}\,G_{D}(\beta,0)\overset{\leftarrow}{W}_{\\!E}$
of course implies
$\overset{\rightarrow}{W}_{\\!E}\,G_{D}(\tau,\tau^{\prime})\overset{\leftarrow}{W}_{\\!E}|_{\tau=\beta,\tau^{\prime}=0}$,
etc. The column ${\boldsymbol{j}}_{E}$ is given by the following integral

$\displaystyle{\boldsymbol{j}}_{E}$
$\displaystyle=\int_{0}^{\beta}d\tau^{\prime}\begin{bmatrix}[r]W_{E}\,G_{D}(\beta,\tau^{\prime})\,\\\
-W_{E}\,G_{D}(0,\tau^{\prime})\,\end{bmatrix}J_{E}(\tau^{\prime}).$ (2.32)

Using (2.11) in the generating functional with the Euclidean density matrix
(2.30)

$Z[{\boldsymbol{J}},J_{E}]\equiv\operatorname{tr}\left[\hat{U}_{J_{1}}(T,0)\mathinner{\hat{\rho}_{E}[J_{E}]}\hat{U}^{\dagger}_{-J_{2}}(T,0)\right]$
(2.33)

one directly finds the total Schwinger-Keldysh generating functional with the
full set of sources probing the two Lorentzian branches and the Euclidean
branch of the in-in formalism

$\displaystyle Z[{\boldsymbol{J}},$ $\displaystyle J_{E}]={\rm
const}\times\exp\Biggl{\\{}-\frac{i}{2}\int_{0}^{T}dt\,dt^{\prime}\,{\boldsymbol{J}}^{T}\\!(t){\boldsymbol{G}}(t,t^{\prime}){\boldsymbol{J}}(t^{\prime})$
$\displaystyle-\int_{0}^{T}dt\,{\boldsymbol{J}}^{T}\\!(t)\,{\boldsymbol{G}}(t,0)\,{\boldsymbol{j}}_{E}+\frac{i}{2}\,{\boldsymbol{j}}^{T}_{E}\,{\boldsymbol{G}}(0,0)\,{\boldsymbol{j}}_{E}$
$\displaystyle+\frac{1}{2}\int_{0}^{\beta}d\tau\,d\tau^{\prime}\,J_{E}(\tau)\,G_{D}(\tau,\tau^{\prime})\,J_{E}(\tau^{\prime})\Biggr{\\}}.$
(2.34)

###  2.3 Reflection symmetry and analyticity properties

The obtained expression for $Z[{\boldsymbol{J}},J_{E}]$ features a nontrivial
mixup of the Neumann and Dirichlet Green’s functions of different Lorentzian
$F$ and Euclidean $F_{E}$ wave operators, but it becomes essentially unified
if we assume that the Lorentzian and Euclidean actions are related by the
analytic continuation of the form

$iS[\phi(t)]\bigr{|}_{t=-i\tau}=-S_{E}[\phi_{E}(\tau)],$ (2.35)

where the Lorentzian and Euclidean histories are also related by the same
continuation rule $\phi(t)|_{t=-i\tau}=\phi_{E}(\tau)$. This, in particular,
implies that the coefficients of the operators $F_{E}$ and $F$ are related by

$\displaystyle\begin{aligned} &A_{E}(\tau)=A(-i\tau),\quad
B_{E}(\tau)=-iB(-i\tau),\\\ &C_{E}(\tau)=-C(-i\tau),\end{aligned}$ (2.36)

so that $F|_{t=-i\tau}=-F_{E}$. The origin of these relations, especially in
connection with reality condition for the coefficients of both Lorentzian and
Euclidean operators at their respective _real_ $t$ and $\tau$ arguments can be
traced back to the properties of the full nonlinear action which gives rise to
its quadratic part on top of a special background solution of full equations
of motion. It is assumed that the Euclidean background solution has a turning
point at $\tau=0$ where all real field variables have zero $\tau$-derivatives
and can be smoothly continued to the imaginary axis of $\tau=it$ where they
become again real functions of real $t$. This leads to the above continuation
rule with real $A(t),B(t),C(t)$ and $A_{E}(\tau),B_{E}(\tau),C_{E}(\tau)$ at
real $t$ and $\tau$.

With this analytic continuation rule, the expression (2.2) for
$Z[{\boldsymbol{J}},J_{E}]$ can indeed be uniformly rewritten in terms of the
Lorentzian, Euclidean and mixed Lorentzian-Euclidean Green’s functions, all of
them subject to one and the same set of Neumann type boundary conditions which
select in the Lorentzian branch of the in-in formalism a distinguished set of
positive and negative frequency basis functions. This expression reads

$\displaystyle Z[\mathbb{J}]={\rm
const}\times\exp\Biggl{\\{}\frac{1}{2}\int_{\mathbb{C}}dz\,dz^{\prime}\,\mathbb{J}^{T}(z)\,\mathbb{G}(z,z^{\prime})\;\mathbb{J}(z^{\prime})\Biggr{\\}},$
(2.37)

where the $z$-integration runs respectively over $t$ or $\tau$ in the domain
$\mathbb{C}=[0\leq t\leq T]\cup[0\leq\tau\leq\beta]$ depending on which of
these Lorentzian or Euclidean time variables is in the argument of the
following block matrix Green’s function $\mathbb{G}(z,z^{\prime})$ and the
corresponding source $\mathbb{J}(z)$,

$\displaystyle\mathbb{G}(z,z^{\prime})\\!=\\!\\!\begin{bmatrix}-i{\boldsymbol{G}}(t,t^{\prime})&\;{\boldsymbol{G}}_{LE}^{<}(t,\tau^{\prime})\,\\\
\;{\boldsymbol{G}}^{>}_{LE}(\tau,t^{\prime})&\;G_{E}(\tau,\tau^{\prime})\,\end{bmatrix}\\!,\;\mathbb{J}(z)\\!=\\!\begin{bmatrix}[c]\,{\boldsymbol{J}}(t)\,\\\
\,J_{E}(\tau)\,\end{bmatrix}\\!.$ (2.38)

Here the Euclidean and Lorentzian-Euclidean blocks of the total Green’s
function

$\displaystyle
G_{E}(\tau,\tau^{\prime})=G_{E}^{>}(\tau,\tau^{\prime})\,\theta(\tau-\tau^{\prime})$
$\displaystyle\qquad\qquad\qquad\qquad+G_{E}^{<}(\tau,\tau^{\prime})\,\theta(\tau^{\prime}-\tau),$
(2.39)
$\displaystyle{\boldsymbol{G}}_{LE}^{<}(t,\tau)=\begin{bmatrix}\,G^{1}_{LE}(t,\tau)\;\\\
G^{2}_{LE}(t,\tau)\end{bmatrix}=\begin{bmatrix}&I&\\\
&I&\end{bmatrix}G_{LE}^{<}(t,\tau),$ (2.40)
$\displaystyle{\boldsymbol{G}}^{>}_{LE}(\tau,t)=\big{[}{\boldsymbol{G}}^{<}_{LE}(t,\tau)\big{]}^{T},$
(2.41)

express in terms of the relevant Euclidean and Lorentzian-Euclidean Wightman
functions

$\displaystyle G_{E}^{>}(\tau,\tau^{\prime})$
$\displaystyle=u_{+}(\tau)(\nu+I)u_{-}^{T}(\tau^{\prime})\\!+\\!u_{-}(\tau)\,\nu\,u_{+}^{T}(\tau^{\prime}),$
(2.42) $\displaystyle G_{E}^{<}(\tau,\tau^{\prime})$
$\displaystyle=\bigl{[}G_{E}^{>}(\tau^{\prime},\tau)\bigl{]}^{T},$ (2.43)
$\displaystyle G_{LE}^{<}(t,\tau)$
$\displaystyle=v(t)\,(\nu+I)\,u_{-}^{T}(\tau)+v^{*}(t)\,\nu\,u_{+}^{T}(\tau).$
(2.44)

In their turn these Green’s functions, as one can see, are built according to
one and the same universal pattern out of the full set of Lorentzian $v(t)$
and $v^{*}(t)$ and Euclidean $u_{\pm}(\tau)$ basis functions. All these
functions are subject to Neumann boundary conditions
(LABEL:eq:green_lor_neum_bdy0) and

$(W_{E}+\omega)u_{+}|_{\tau=\beta}=0,\qquad(W_{E}-\omega)u_{-}|_{\tau=0}=0.$
(2.45)

For $\omega$ fixed by the above condition of particle interpretation, leading
to the expressions (2.20)–(2.22), the Euclidean basis functions $u_{\pm}$ have
a remarkable property. They satisfy at opposite ends of the Euclidean segment
$\tau_{\mp}$ the same boundary conditions444In fact, the requirement of
$\kappa=0$ in (2.19) turns out to be the necessary and sufficient condition
for this property of Euclidean basis functions, that is the coincidence of
boundary condition for $u_{\pm}(\tau)$ at both ends of the time segment leads
to $\kappa=0$.

$(W_{E}+\omega)u_{+}|_{\tau=0}=0,\quad(W_{E}-\omega)u_{-}|_{\tau=\beta}=0.$
(2.46)

If one smoothly continues the operator $F_{E}$ beyond the segment
$\tau_{-}\leq\tau\leq\tau_{+}$, then it becomes periodic with the period
$\beta$ (which is possible because $\tau_{\pm}$ are assumed to be the turning
points of the background solution on top of which the Hessian of the nonlinear
action of the theory is built). This means that the basis functions $u_{\pm}$
of this operator become quasi-periodic — $u_{\pm}(\tau+\beta)$ expresses as a
linear combination of the same basis functions $u_{\pm}(\tau)$ (no mixing
between $u_{-}$ and $u_{+}$ occurs in their monodromy matrix). As shown in
Section 4, with the normalization $u_{\pm}(0)=1/\sqrt{2\omega}$ this quasi-
periodicity property reads in terms of the occupation number matrix (2.21)

$\displaystyle u_{-}(\tau+\beta)$
$\displaystyle=u_{-}(\tau)\frac{\nu+I}{\nu},$ (2.47) $\displaystyle
u_{+}(\tau+\beta)$ $\displaystyle=u_{+}(\tau)\frac{\nu}{\nu+I}.$

Together with the reflection symmetry relative to the middle point of the
Euclidean time segment (2.28) the periodicity of the operator $F_{E}$ implies
its reflection symmetry with respect to the point $\tau=0$

$\displaystyle\begin{aligned} &A_{E}(\tau)=A_{E}(-\tau),\quad
B_{E}(\tau)=-B_{E}(-\tau),\\\ &C_{E}(\tau)=C_{E}(-\tau).\end{aligned}$ (2.48)

Therefore, similarly to quasi-periodicity the basis functions
$u_{\pm}(\pm\tau)$ are also related by the analogue of the anti-diagonal
monodromy matrix $L$, $u_{+}(\tau)=u_{-}(-\tau)\,L$, which is trivial in view
of the normalization $u_{\pm}(0)=1/\sqrt{2\omega}$,

$u_{+}(\tau)=u_{-}(-\tau).$ (2.49)

The above relations introduce the analytic structure which allows one to
express all basis and Green’s functions on the Euclidean-Lorentzian domain
$\mathbb{C}$ in terms of one analytic function $V(z)$ of the complexified time
variable $z=t-i\tau$. This follows from the fact, mentioned above, that the
Lorentzian wave operator $F$ can be regarded as the analytic continuation of
the Euclidean operator $F_{E}$ into the complex plane of time at the point
$z=0$, $F\equiv F(t,d/dt)=-F_{E}|_{\tau=it}$. As a consequence its basis
function $v(t)$ in view of its boundary conditions and boundary conditions
(2.46) for the Euclidean function $u_{+}(\tau)$ also turns out to be the
analytic continuation of the latter,

$v(t)=u_{+}(it).$ (2.50)

Therefore, the operators $F$ and $-F_{E}$ as well as the full set of their
basis functions $v(t)$ and $u_{\pm}(\tau)$ can be represented respectively as
the boundary values at the real and imaginary axes of the complex $z$-plane of
the complex operator $F_{\mathbb{C}}$ and the solution $V(z)$ of its
homogeneous wave equation,

$\displaystyle
F_{\mathbb{C}}V(z)\equiv\biggl{[}-\frac{d}{dz}A(z)\frac{d}{dz}-\frac{d}{dz}B(z)+B^{T}(z)\frac{d}{dz}$
$\displaystyle\qquad\qquad\quad+C(z)\biggr{]}V(z)=0,\quad z=t-i\tau,$ (2.51)
$\displaystyle\bigl{(}iW_{\mathbb{C}}-\omega\bigr{)}V(z)\bigr{|}_{z=0}=0,\;\,W_{\mathbb{C}}\equiv
A(z)\frac{d}{dz}+B(z).$ (2.52)

The function $V(z)$ gives rise to basis functions as

$\displaystyle v(t)=V(z)\bigr{|}_{z=t},\quad u_{\pm}(\tau)=V(z)\bigr{|}_{z=\mp
i\tau},$ (2.53)

and thus can be used in (2.38) for the construction of all Green’s functions
of the Schwinger-Keldysh in-in formalism. Conversely $V(z)$ can be obtained by
analytic continuation of the single Euclidean function $u_{+}(\tau)$ from the
imaginary axes $z=-i\tau$,

$\displaystyle V(z)=u_{+}(iz)=u_{+}(\tau+it),$ (2.54)

and in view of reality of $u_{+}(\tau)$ for real $\tau$ it has the property
$[V(z)]^{*}=V(-z^{*})$.

Important corollary of these analiticity properties is that in view of the
monodromy relations for Euclidean basis functions (2.47) the Lorenzian basis
functions become quasi-periodic in the imaginary time

$\displaystyle
v(t-i\beta)\\!=\\!v(t)\frac{\nu+I}{\nu},\;\,v^{*}(t-i\beta)\\!=\\!v^{*}(t)\frac{\nu}{\nu+I}.$
(2.55)

Due to inverse matrix factors of positive and negative basis functions here
the Lorentzian Wightman functions $G_{>}(t,t^{\prime})$ (given by the
expression (2.25)) and $G_{<}(t,t^{\prime})=G^{T}_{>}(t^{\prime},t)$ satisfy
the relation

$G_{>}(t-i\beta,t^{\prime})=G_{<}(t,t^{\prime}),$ (2.56)

which is nothing but Kubo-Martin-Schwinger condition [43, 44]. It is important
that this condition is satisfied in the generic non-equilibrium system with
the special Euclidean density matrix (2.26) even despite the fact that no
notion of conserved energy can be formulated in such a physical setup.

Figure 5: Euclidean-Lorentzian contour $\mathbb{C}$ on the Riemann surface of
complex time $z=t-i\tau$. Wightman functions are periodic in imaginary
(Euclidean) time direction with a period $\beta$, whereas the basis function
$v(z)$ suffers a jump at the cut denoted by the horizontal dashed line, the
two Lorentzian time branches running along the shores of this cut.

The tubular Riemann surface of complex time $z=t-i\tau$ whose main sheet is
compactified in $\tau$ to the circle of curcumference $\beta$ is shown on Fig.
5. The boundaries of the main sheet of this surface form two shores of the cut
depicted by dashed line, along which two branches of Lorentzian evolution are
running. This rich analytic structure of Euclidean-Lorentzian evolution
suggests that the equivalence of the Euclidean and Lorentzian formalisms
proven beyond tree level for interacting QFT on top of the de Sitter spacetime
[32, 33] might be extended to a generic reflection-symmetry background
underlying our definition of the Euclidean density matrix.

##  3 Preliminaries

To derive the aforementioned results we dwell here in more detail on the
introduced above notations and develop a canonical formalism and quantization
of the underlying theory. In particular, we pose rather generic initial value
and boundary value problems for equations of motion and discuss the properties
of the related Green’s functions.

###  3.1 Condensed notations

The elements of the field space will be denoted as $\phi^{I}(t)$, where the
index $I$ is, in fact, a multi-index, and contains both the dependence on
spatial coordinates denoted as $\mathbf{x}$ and discrete spin-tensor labels
$i$, $I=(\mathbf{x},i)$. Thus, we can equivalently write the fields in the
form, emphasizing its dependence on the spatial coordinates
$\phi^{I}(t)=\phi^{i}(t,\mathbf{x})$.

Assuming that equations of motion are of the second order in time derivatives
one has the most general quadratic action of the theory of the form (2.1)
where we explicitly specify the initial and final moments of time range
$t_{\pm}$,

$\displaystyle
S[\phi]=\frac{1}{2}\int_{t_{-}}^{t_{+}}dt\,\Bigl{(}\dot{\phi}^{T}A\dot{\phi}+\dot{\phi}^{T}B\phi+\phi^{T}B^{T}\dot{\phi}+\phi^{T}C\phi\Bigr{)}.$
(3.1)

Here dots denote the derivatives with respect to time $t$, and $A$, $B$ and
$C$ are the time-dependent real bilinear forms in the space of fields.
Moreover, $A$ and $C$ are assumed to be symmetric. The explicit action of
these bilinear forms on the fields, e.g. for $A$ reads

$(A\phi)_{I}(t)=A_{IJ}(t)\phi^{J}(t)\equiv\sum_{j}\int
d\mathbf{x}^{\prime}\,A_{ij}(t,\mathbf{x},\mathbf{x}^{\prime})\,\phi^{j}(t,\mathbf{x}^{\prime}),$
(3.2)

where $A_{ij}(t,\mathbf{x},\mathbf{x}^{\prime})$ is the kernel of the
operator. Thus, the first term in (3.1) has the following explicit structure

$\displaystyle\dot{\phi}^{T}A\dot{\phi}$
$\displaystyle=\dot{\phi}^{I}(A\dot{\phi})_{I}$ $\displaystyle=\sum_{ij}\int
d\mathbf{x}\,d\mathbf{x}^{\prime}\,\dot{\phi}^{i}(t,\mathbf{x})A_{ij}(t,\mathbf{x},\mathbf{x}^{\prime})\dot{\phi}^{j}(t,\mathbf{x}).$
(3.3)

The superscript $T$ applied to the bilinear form denotes the functional matrix
transposition operation which implies the transposition of discrete and
spatial labels of the corresponding kernel, but does not touch the time
variable

$(B^{T})_{ij}(t,\mathbf{x},\mathbf{x}^{\prime})=B_{ji}(t,\mathbf{x}^{\prime},\mathbf{x}).$
(3.4)

Consequently, the second and the third terms in (3.1) are the same. However,
we will keep them separate for symmetry reasons.

In local non-gauge theories the kernels of the above coefficients are
represented by delta functions of spatial coordinates and their finite order
derivatives. For local gauge theories treated within reduction to the physical
sector in certain gauges these coefficients can become nonlocal in space, but
locality in time derivatives within canonical quantization should be strictly
observed.

The equations of motion, obtained by varying the action (3.1) with respect to
$\phi$ have the form

$F\phi(t)=0,\quad
F\equiv-\frac{d}{dt}A\frac{d}{dt}-\frac{d}{dt}B+B^{T}\frac{d}{dt}+C,$ (3.5)

where the wave operator $F$, or the Hessian of the action (3.1), has already
been defined above by Eq. (2.14). Another form of this operator, obtained by
integration by parts and involving both left and right time derivatives, the
direction of their action being indicated by arrows

$\overset{\leftrightarrow}{F}\equiv\overset{\leftarrow}{\frac{d}{dt}}A\overset{\rightarrow}{\frac{d}{dt}}+\overset{\leftarrow}{\frac{d}{dt}}B+B^{T}\overset{\rightarrow}{\frac{d}{dt}}+C,$
(3.6)

allows one to rewrite the quadratic action (3.1) in even more condensed form

$\displaystyle S[\phi]$
$\displaystyle=\frac{1}{2}\int_{t_{-}}^{t_{+}}dt\,\phi^{T}\overset{\leftrightarrow}{F}\,\phi$
$\displaystyle=\frac{1}{2}\int_{t_{-}}^{t_{+}}dt\,\phi^{T}(F\phi)+\frac{1}{2}\phi^{T}(W\phi)\Bigr{|}_{t_{-}}^{t_{+}}.$
(3.7)

Here the Wronskian operator $W$ is defined by (2.16) and the origin of the
boundary term at $t_{\pm}$ is the result of integration by parts, which is
also associated with the Wronskian relation (2.17).

###  3.2 Canonical formalism

The Hamiltonian formalism of the theory with the action (3.1), which is the
first step to the canonical quantization begins with the determination of the
momentum $\pi$ canonically conjugated to the field $\phi$

$\pi=\frac{\partial L}{\partial\dot{\phi}}=A\dot{\phi}+B\phi=W\phi,\quad
W=A\frac{d}{dt}+B,$ (3.8)

where $L$ is the Lagrangian of the action (3.1). The corresponding Hamiltonian
has the form

$\displaystyle
H=\pi^{T}\dot{\phi}-L=\frac{1}{2}(\pi-B\phi)^{T}A^{-1}(\pi-B\phi)-\frac{1}{2}\phi^{T}C\phi.$
(3.9)

Together with the Poisson bracket $\\{\phi^{I},\pi_{J}\\}=\delta^{I}_{J}$ it
defines the dynamics of the system. The Hamiltonian equations of motion read

$\displaystyle\dot{\phi}$ $\displaystyle=\\{\phi,H\\}=A^{-1}(\pi-B\phi)$
(3.10a) $\displaystyle\dot{\pi}$
$\displaystyle=\\{\pi,H\\}=B^{T}A^{-1}(\pi-B\phi)+C\phi.$ (3.10b)

Transition to the Lagrangian formalism by expressing $\pi$ in terms of $\phi$
and $\dot{\phi}$ obviously leads to equations of motion (3.5) following from
the variation of the action (3.1).

Let us denote the basis of the independent solutions to (3.5) as
$v_{\pm}{}^{I}_{A}(t)$, where the multi-index $A$ enumerates the number of the
particular solution and has the same range as the index $I$. The general
solution in terms of basis functions reads

$\phi^{I}(t)=v_{+}{}_{A}^{I}(t)\,\alpha^{+A}+v_{-}{}_{A}^{I}(t)\,\alpha^{-A}$
(3.11)

and can be rewritten in shortened notations as

$\phi(t)=v_{+}(t)\,\alpha^{+}+v_{-}(t)\,\alpha^{-}.$ (3.12)

Here $\alpha^{\pm A}$ constitute a set of constants, specifying particular
initial conditions. Using (3.8), we find the corresponding solution for the
momentum

$\pi(t)=Wv_{+}(t)\,\alpha^{+}+Wv_{-}(t)\,\alpha^{-},$ (3.13)

so, the evolution of phase space variables can be rewritten in the joint form
as

$\begin{bmatrix}\,\phi(t)\,\\\
\,\pi(t)\end{bmatrix}=\mathcal{M}(t)\begin{bmatrix}[l]\,\alpha^{+}\\\
\,\alpha^{-}\end{bmatrix},\;\mathcal{M}(t)=\begin{bmatrix}\,v_{+}(t)&v_{-}(t)\,\\\
\,Wv_{+}(t)&Wv_{-}(t)\,\end{bmatrix}.$ (3.14)

Now, we can equip the space of initial conditions, consisting of
$\alpha^{\pm}$, with the Poisson bracket structure inherited from the Poisson
brackets of $\phi$ and $\pi$. Substituting (3.14) into the left hand side of

$\displaystyle\left\\{\begin{pmatrix}\phi^{I}\\\
\pi_{J}\end{pmatrix},\begin{pmatrix}\phi^{I^{\prime}}&\pi_{J^{\prime}}\end{pmatrix}\right\\}=\begin{bmatrix}\,0&\delta^{I}_{J^{\prime}}\;\\\
\;-\delta_{J}^{I^{\prime}}&0\;\end{bmatrix},$ (3.15)

we have in condensed notations

$\displaystyle\mathcal{M}(t)\begin{bmatrix}\;\\{\alpha^{+},\alpha^{+}\\}&\\{\alpha^{+},\alpha^{-}\\}\;\\\
\;\\{\alpha^{-},\alpha^{+}\\}&\\{\alpha^{-},\alpha^{-}\\}\;\end{bmatrix}\mathcal{M}^{T}(t)=\begin{bmatrix}[r]\,0&&I&\\\
\,-I&&0&\end{bmatrix},$ (3.16)

where $I$ denotes the identity matrix. The identity above fixes the pairwise
Poisson brackets of $\alpha^{\pm}$. Let us denote the right hand side of this
equality, playing the role of the Poisson bivector in the Darboux coordinates,
as

$\mathcal{P}\equiv\begin{bmatrix}[r]0&&I&\\\ -I&&0&\end{bmatrix}.$ (3.17)

Introducing also the matrix $\mathcal{D}$ as inverse to the matrix of the
pairwise Poisson brackets

$\mathcal{D}=\begin{bmatrix}\;\Delta_{++}&&\Delta_{+-}\\\
\;\Delta_{-+}&&\Delta_{--}\end{bmatrix}\equiv-\begin{bmatrix}\,\\{\alpha^{+},\alpha^{+}\\}&\\{\alpha^{+},\alpha^{-}\\}\,\\\
\,\\{\alpha^{-},\alpha^{+}\\}&\\{\alpha^{-},\alpha^{-}\\}\,\end{bmatrix}^{-1}$
(3.18)

where the matrices $\Delta$ denote the corresponding block-elements of
$\mathcal{D}$, we can invert the equality (3.16) as

$\mathcal{M}^{T}(t)\,\mathcal{P}\,\mathcal{M}(t)=\mathcal{D}.$ (3.19)

Thus, one can express the inverse of $\mathcal{M}(t)$ in terms of its
transpose, namely

$\mathcal{M}^{-1}(t)=\mathcal{D}^{-1}\,\mathcal{M}^{T}(t)\,\mathcal{P}.$
(3.20)

Before proceeding further, let us show explicitly that the right hand side
(3.19) is indeed independent of time $t$. To demonstrate this, we contract
l.h.s of the equation (3.5) where field $\phi=\phi_{1}$, with another field
$\phi_{2}$, and subtract the same quantity, but with $F$, acting on $\phi_{2}$
($\phi_{1,2}$ are not necessarily solve e.o.m.). The result can be written as

$\phi_{2}^{T}F\phi_{1}-(F\phi_{2})^{T}\phi_{1}=-\frac{d}{dt}\left[\phi_{2}^{T}W\phi_{1}-(W\phi_{2})^{T}\phi_{1}\right].$
(3.21)

Thus, for $\phi_{1,2}$ — solutions of (3.5) l.h.s. vanishes, so we have

$\phi_{2}^{T}W\phi_{1}-(W\phi_{2})^{T}\phi_{1}=\text{const}.$ (3.22)

It is easy to see, that each element of (3.18) has the form (3.22) as above,
where the role of solutions $\phi_{1}$, $\phi_{2}$ is played by the basis
functions $v^{+}$, $v^{-}$. Applying the matrix transposition operator to both
sides of (3.16), we obtain that the matrix $\mathcal{D}$ is skew-symmetric,
since $\mathcal{P}^{T}=-\mathcal{P}$. In terms of the block elements of
$\mathcal{D}$ this means that

$\Delta_{+-}^{T}=-\Delta_{-+},\;\Delta_{++}^{T}=-\Delta_{++},\;\Delta_{--}^{T}=-\Delta_{--}.$
(3.23)

Moreover, using the fact that the coefficient matrices $A$, $B$, and $C$ in
(3.1) are real, we conclude that basis functions $v_{+}$, $v_{-}$ can also be
chosen to be real. Thus, the matrix $\mathcal{D}$ is real skew-symmetric, so
there is a time-independent linear transformation $\mathcal{S}$, bringing it
to the canonical form, i.e.
$\mathcal{S}^{T}\mathcal{D}\mathcal{S}=\mathcal{P}$. Without the loss of
generality one set $\mathcal{D}=\mathcal{P}$ by default555This choice allows
to give an additional interpretation for the equation (3.19), which becomes
$\mathcal{M}^{T}(t)\mathcal{P}\mathcal{M}(t)=\mathcal{P}$. Namely, the matrix
$\mathcal{M}(t)$ performs a time-dependent symplectomorphism of the Poisson
bivector $\mathcal{P}$.. However, for the reasons which will become clear soon
(see equation (4.21) below), we will assume that $\mathcal{D}$ has the
following more general form

$\mathcal{D}=\begin{bmatrix}[l]\;\;\,0&&\Delta_{+-}\\\
\;\;\Delta_{-+}&&\,0\end{bmatrix},$ (3.24)

where

$\Delta_{+-}=-\Delta_{-+}^{T}=v_{+}^{T}Wv_{-}-(Wv_{+})^{T}v_{-}.$ (3.25)

In terms of the basis functions, the vanishing of the diagonal blocks of
$\mathcal{D}$ implies that $v_{+}$, $v_{-}$ are chosen such that

$\displaystyle\begin{aligned}
\Delta_{++}&=v_{+}^{T}Wv_{+}-(Wv_{+})^{T}v_{+}=0,\\\
\Delta_{--}&=v_{-}^{T}Wv_{-}-(Wv_{-})^{T}v_{-}=0.\end{aligned}$ (3.26)

This can always be done by an appropriate transformation of the basis
functions, possibly mixing $v_{+}$ and $v_{-}$. Consequently, the pairwise
Poisson brackets of $\alpha^{+}$ and $\alpha^{-}$ take the form

$\displaystyle\\{\alpha^{+},\alpha^{-}\\}=-\\{\alpha^{+},\alpha^{-}\\}=-\Delta_{-+}^{-1},$
(3.27a)
$\displaystyle\\{\alpha^{+},\alpha^{+}\\}=\\{\alpha^{-},\alpha^{-}\\}=0.$
(3.27b)

As noted above, one can go further, and set $\Delta_{+-}=-\Delta_{-+}^{T}=I$.

Now, let us modify the Hamiltonian by introducing time-dependent sources
$J_{\phi}$, $J_{\pi}$ for the field and its conjugate momentum

$H\;\mapsto\;H+J_{\phi}^{T}\phi+J_{\pi}^{T}\pi.$ (3.28)

The modified equations of motion can be written as

$\displaystyle\frac{d}{dt}\begin{bmatrix}\phi_{J}(t)\\\
\pi_{J}(t)\end{bmatrix}=\mathcal{A}(t)\begin{bmatrix}\phi_{J}(t)\\\
\pi_{J}(t)\end{bmatrix}+\mathcal{P}\begin{bmatrix}[r]J_{\phi}(t)\\\
J_{\pi}(t)\end{bmatrix},$ (3.29)
$\displaystyle\mathcal{A}(t)\equiv\begin{bmatrix}-A^{-1}B&&&A^{-1}\\\
-B^{T}A^{-1}B+C&&&B^{T}A^{-1}\end{bmatrix},$

where the subscript $J$ of $\phi$, $\pi$ emphasizes the presence of the
sources in equations of motion. We will find a solution to modified equations
of motion using the constant variation method. Namely, we start with the
solution (3.14) to equations of motion with vanishing sources, but make the
integration constants $\alpha^{+}$, $\alpha^{-}$ in its definition time-
dependent

$\begin{bmatrix}\phi_{J}(t)\\\
\pi_{J}(t)\end{bmatrix}=\mathcal{M}(t)\begin{bmatrix}\alpha^{+}(t)\\\
\alpha^{-}(t)\end{bmatrix},$ (3.30)

Then, we substitute the result to the modified e.o.m. and obtain

$\mathcal{M}(t)\frac{d}{dt}\\!\begin{bmatrix}\alpha^{+}(t)\\\
\alpha^{-}(t)\end{bmatrix}=\mathcal{P}\begin{bmatrix}[r]J_{\phi}(t)\\\
J_{\pi}(t)\end{bmatrix},$ (3.31)

where we exploit the fact that $\mathcal{M}(t)$ satisfies the system (3.10).
Using the equality (3.20) for the inverse of the matrix $\mathcal{M}(t)$ and
integrating the equation on $\alpha^{+}(t)$ and $\alpha^{-}(t)$ we obtain

$\begin{bmatrix}\alpha^{+}(t)\\\
\alpha^{-}(t)\end{bmatrix}=\begin{bmatrix}\alpha^{+}_{0}\\\
\alpha^{-}_{0}\end{bmatrix}-\int_{t_{-}}^{t}dt^{\prime}\;\mathcal{D}^{-1}\mathcal{M}^{T}(t^{\prime})\begin{bmatrix}[r]J_{\phi}(t^{\prime})\\\
J_{\pi}(t^{\prime})\end{bmatrix},$ (3.32)

where $\alpha^{+}_{0}$ and $\alpha^{-}_{0}$ are integration constants.
Substitution back to (3.30) gives the solution to the equations (3.29)

$\begin{bmatrix}\phi_{J}(t)\\\
\pi_{J}(t)\end{bmatrix}=\begin{bmatrix}\phi_{0}(t)\\\
\pi_{0}(t)\end{bmatrix}-\int_{t_{-}}^{t}dt^{\prime}\;\mathcal{M}(t)\,\mathcal{D}^{-1}\mathcal{M}^{T}(t^{\prime})\begin{bmatrix}[r]J_{\phi}(t^{\prime})\\\
J_{\pi}(t^{\prime})\end{bmatrix},$ (3.33)

where the initial conditions $\phi_{0}(t)$, $\pi_{0}(t)$ are related to
constants of integration by

$\begin{bmatrix}\phi_{0}(t)\\\
\pi_{0}(t)\end{bmatrix}=\mathcal{M}(t)\begin{bmatrix}\alpha^{+}_{0}\\\
\alpha^{-}_{0}\end{bmatrix},$ (3.34)

and represent the solution to homogeneous equation, i.e. for vanishing sources
$J_{\phi}$ and $J_{\pi}$.

Now, let us focus on the case of vanishing momentum source and also redefine
the field source for the convenience

$J_{\pi}(t)=0,\qquad J(t)\equiv-J_{\phi}(t).$ (3.35)

The corresponding e.o.m. in the Lagrange form reads

$F\phi_{J}(t)+J(t)=0.$ (3.36)

From (3.33), one obtains the explicit form of the solution for $\phi(t)$,
which is

$\phi_{J}(t)=\phi_{0}(t)-\int_{t_{-}}^{t_{+}}dt^{\prime}\,G_{R}(t,t^{\prime})J(t^{\prime}),$
(3.37)

where $G_{R}(t,t^{\prime})$ is called the retarded Green’s function and
expressed through the top-left block of the matrix
$\mathcal{M}(t)\,\mathcal{D}^{-1}\mathcal{M}^{T}(t^{\prime})$, specifically

$G_{R}(t,t^{\prime})\\!=\\!-\Bigl{(}v_{+}(t)\Delta_{-+}^{-1}v_{-}^{T}(t^{\prime})+v_{-}(t)\Delta_{+-}^{-1}v_{+}^{T}(t^{\prime})\Bigr{)}\,\theta(t-t^{\prime}).$
(3.38)

The fact that $\Delta_{++}=\Delta_{--}=0$ is crucial in obtaining this simple
expression for $G_{R}$. From (3.37) we find that $G_{R}$ satisfies the
equation

$FG_{R}(t,t^{\prime})=I\,\delta(t-t^{\prime})$ (3.39)

and is uniquely determined by the condition

$G_{R}(t,t^{\prime})=0,\qquad t<t^{\prime}.$ (3.40)

The latter fact follows, in particular, from the fact that any two Green’s
functions of the same differential operator differ by the solution of the
homogeneous equation. Once some Green’s function, satisfying the condition
(3.40) is found, a shift by a solution to homogeneous equation will violate
this condition. Alternatively, $G_{R}$ can be defined via initial value
problem

$G_{R}(t,t^{\prime})\bigr{|}_{t^{\prime}=t+0}=0,\;WG_{R}(t,t^{\prime})\bigr{|}_{t^{\prime}=t+0}=-I.$
(3.41)

The fact that solution (3.37) is expressed through the retarded Green’s
function means that $\phi(t)$ is subject to the following initial (rather than
boundary) value problem

$\phi_{J}(t_{-})=\phi_{0}(t_{-}),\;W\phi_{J}(t_{-})=W\phi_{0}(t_{-})\equiv\pi_{0}(t_{-}).$
(3.42)

###  3.3 The solution of Dirichlet and Neumann boundary value problems

The Green’s functions, solving the boundary problems, can be obtained from the
retarded Green’s function by shifting it by the solution of the homogeneous
equation (3.5). In particular, one constructs the so-called symmetric Green’s
function as

$\displaystyle G_{S}(t,t^{\prime})$
$\displaystyle=G_{R}(t,t^{\prime})+v_{+}(t)\Delta_{-+}^{-1}v_{-}^{T}(t^{\prime})$
$\displaystyle=-v_{+}(t)\Delta_{-+}^{-1}v_{-}^{T}(t^{\prime})\,\theta(t-t^{\prime})$
$\displaystyle\quad+v_{-}(t)\Delta_{+-}^{-1}v_{+}^{T}(t^{\prime})\,\theta(t^{\prime}-t).$
(3.43)

It is symmetric under the simultaneous transposition and exchange of the time
arguments, i.e. $G_{S}^{T}(t,t^{\prime})=G_{S}(t^{\prime},t)$. Unlike the
retarded Green’s function it is defined non-uniquely and the concrete boundary
conditions should be specified. These are in one-to-one correspondence to the
boundary conditions, satisfied by the basis functions $v_{+}$ and $v_{-}$ at
the higher and lower time limits $t=t_{+}$ and $t=t_{-}$, respectively.

In particular, to solve the inhomogeneous equation (3.36) supplemented with
the vanishing Dirichlet boundary conditions

$\phi_{J}(t_{\pm})=0,$ (3.44)

one can use the Dirichlet Green’s function subject to the same boudary
conditions

$G_{D}(t_{\pm},t^{\prime})=0\qquad\leftrightarrow\qquad v_{\pm}(t_{\pm})=0,$
(3.45)

so that the solution reads

$\phi_{J}(t)=-\int_{t_{-}}^{t_{+}}dt^{\prime}\,G_{D}(t,t^{\prime})\,J(t^{\prime}).$
(3.46)

Similarly, in solving Neumann boundary problem

$(iW\mp\omega_{\pm})\phi_{J}(t_{\pm})=0,$ (3.47)

one defines the corresponding Neumann Green’s function demanding

$(iW\mp\omega_{\pm})G_{N}(t_{\pm},t^{\prime})=0\leftrightarrow(iW\mp\omega_{\pm})v_{\pm}(t_{\pm})=0,$
(3.48)

and obtains the solution as

$\phi_{J}(t)=-\int_{t_{-}}^{t_{+}}dt^{\prime}\,G_{N}(t,t^{\prime})\,J(t^{\prime}).$
(3.49)

Notably, the Dirichlet and Neumann Green’s functions, which are subject to
homogeneous boundary conditions, allows one to solve the modified boundary
problems, namely with inhomogeneous boundary conditions. Namely, the solutions
can be obtained as follow. First, we exploit the equality (3.21) and perform
in it the substitutions $\phi_{2}\mapsto\phi(t^{\prime})$, $\phi_{1}\mapsto
G(t^{\prime},t)$, where $\phi(t^{\prime})$ solves (3.36) and $G(t^{\prime},t)$
is some Green’s function, solving $FG(t^{\prime},t)=\delta(t-t^{\prime})$.
Next, integrating both sides of the equality over $t^{\prime}$ from $t_{-}$ to
$t_{+}$, we obtain

$\displaystyle\phi_{J}(t)={}$
$\displaystyle-\int_{t_{-}}^{t_{+}}dt^{\prime}\,G(t,t^{\prime})\,J(t^{\prime})$
$\displaystyle+(WG(t_{+},t))^{T}\phi(t_{+})-(WG(t_{+},t))^{T}\phi_{J}(t_{+})$
$\displaystyle-G^{T}(t_{+},t)\,W\phi_{J}(t_{+})+G^{T}(t_{-},t)\,W\phi_{J}(t_{-})$
(3.50)

Now, suppose we are to solve (3.36) supplemented by inhomogeneous boundary
conditions (in contrast to homogeneous ones (3.44))

$\phi_{J}(t_{\pm})=\varphi_{\pm},$ (3.51)

for some constants $\varphi_{+}$, $\varphi_{-}$. Substituting these conditions
to (3.3) together with Dirichlet Green’s function $G\mapsto G_{D}$, satisfying
(3.45), we observe that the third line vanishes, so we get

$\phi_{J}(t)=-{\boldsymbol{w}}^{T}(t)\begin{bmatrix}\;\varphi_{+}\\\
\;\varphi_{-}\end{bmatrix}-\int_{t_{-}}^{t_{+}}dt^{\prime}\,G_{D}(t,t^{\prime})\,J(t^{\prime}),$
(3.52)

where we introduce the notation for the two-component row as the transposition
of the newly introduced column

$\displaystyle{\boldsymbol{w}}^{T}(t)$
$\displaystyle\equiv\begin{bmatrix}\;G_{D}(t,t_{+})\overset{\leftarrow}{W}&&-G_{D}(t,t_{-})\overset{\leftarrow}{W}\;\end{bmatrix}$
$\displaystyle=[{\boldsymbol{w}}(t)]^{T},$ (3.53)
$\displaystyle{\boldsymbol{w}}(t)$
$\displaystyle\equiv\begin{bmatrix}[r]\overset{\rightarrow}{W}G_{D}(t_{+},t)\\\
-\overset{\rightarrow}{W}G_{D}(t_{-},t)\end{bmatrix},$ (3.54)

and $\overset{\leftarrow}{W}$ denotes the Wronskian operator (2.16) acting
from the right on the second argument of $G_{D}(t,t^{\prime})$ at the total
boundary of the time domain at $t_{\pm}$ (the sign taking into account the
outward pointing time derivative in $W$) — the notation used above in (2.31).
The transposition law here, of course, takes into account the symmetry of
Dirichlet Green’s function,

$\displaystyle[G_{D}(t,t_{+})\overset{\leftarrow}{W}]^{T}$
$\displaystyle=\bigl{(}A(t_{+})\frac{d}{dt_{+}}+B(t_{+})\bigr{)}G_{D}^{T}(t,t_{+})$
$\displaystyle=\overset{\rightarrow}{W}G_{D}(t_{+},t).$ (3.55)

The quantity ${\boldsymbol{w}}(t)$ introduced above has the following
important property. Namely, evaluating both sides of (3.52) at the boundary
points $t=t_{\pm}$, and using (3.45) we observe that

$\displaystyle{\boldsymbol{w}}^{T}(t_{+})=\begin{bmatrix}\,-I&&0\,\end{bmatrix},\qquad{\boldsymbol{w}}^{T}(t_{-})$
$\displaystyle=\begin{bmatrix}\,0&&-I\,\end{bmatrix}.$ (3.56)

Similarly, one can consider inhomogeneous Neumann boundary conditions

$\bigl{(}\pm iW-\omega_{\pm}\bigr{)}\phi_{J}(t_{\pm})=j_{\pm},$ (3.57)

with some boundary sources $j_{+}$ and $j_{-}$. Substitution of this condition
and Neumann Green’s function $G\mapsto G_{N}$, satisfying (3.48), to (3.3)
gives the solution to (3.36) with the boundary conditions above

$\phi_{J}(t)=-i\,{\boldsymbol{g}}_{N}^{T}(t)\begin{bmatrix}\;j_{+}\\\
\;j_{-}\end{bmatrix}-\int_{t_{-}}^{t_{+}}dt^{\prime}\,G_{N}(t,t^{\prime})J(t^{\prime}).$
(3.58)

Here ${\boldsymbol{g}}_{N}^{T}(t)$ is the notation analogous to (3.53) — the
row built in terms of the Neumann Green’s function kernels with the second
argument located at the total 2-point boundary of the time domain (points
$t_{-}$ and $t_{+}$),

${\boldsymbol{g}}_{N}^{T}(t)\equiv\begin{bmatrix}\;G_{N}(t,t_{+})&&G_{N}(t,t_{-})\;\end{bmatrix}.$
(3.59)

###  3.4 The relation between Dirichlet and Neumann Green’s functions

There is important explicit connection between Dirichlet and Neumann Green’s
functions, which can be derived in the following way. The idea is to consider
the problem with homogeneous Neumann boundary conditions (3.47) as the
Dirichlet problem with some nontrivial boundary values $\varphi_{\pm}$.
Substituting the solution of this problem (3.52) into (3.47) one obtains a
linear equation on $\varphi_{\pm}$, which can be solved as

$\begin{bmatrix}\;\varphi_{+}\\\
\;\varphi_{-}\end{bmatrix}=(i{\boldsymbol{\omega}}+{\boldsymbol{\varOmega}})^{-1}\int_{t_{-}}^{t_{+}}dt\,{\boldsymbol{w}}(t)\,J(t),$
(3.60)

where the matrices ${\boldsymbol{\omega}}$ and ${\boldsymbol{\varOmega}}$ read

$\displaystyle\qquad\qquad\qquad{\boldsymbol{\omega}}\equiv\begin{bmatrix}\;\omega_{+}&0\\\
0&\;\omega_{-}\end{bmatrix},$ (3.61)
$\displaystyle{\boldsymbol{\varOmega}}\equiv\begin{bmatrix}[r]-\overset{\rightarrow}{W}G_{D}(t_{+},t_{+})\overset{\leftarrow}{W}&\overset{\rightarrow}{W}G_{D}(t_{+},t_{-})\overset{\leftarrow}{W}\;\\\
\overset{\rightarrow}{W}G_{D}(t_{-},t_{+})\overset{\leftarrow}{W}&-\overset{\rightarrow}{W}G_{D}(t_{-},t_{-})\overset{\leftarrow}{W}\;\end{bmatrix}.$
(3.62)

Substituting these $\varphi_{\pm}$ back into (3.52) gives

$\displaystyle\phi_{J}(t)$
$\displaystyle=-\int_{t_{-}}^{t_{+}}dt^{\prime}\,\Bigl{[}G_{D}(t,t^{\prime})$
$\displaystyle+{\boldsymbol{w}}^{T}(t)\,(i{\boldsymbol{\omega}}+{\boldsymbol{\varOmega}})^{-1}{\boldsymbol{w}}(t^{\prime})\Bigr{]}J(t^{\prime}),$
(3.63)

which implies, after comparing with (3.49), the following expression for the
Neumann Green’s functrion

$G_{N}(t,t^{\prime})=G_{D}(t,t^{\prime})+{\boldsymbol{w}}^{T}(t)\,(i{\boldsymbol{\omega}}+{\boldsymbol{\varOmega}})^{-1}{\boldsymbol{w}}(t^{\prime}).$
(3.64)

Here we use the notations (3.53)–(3.54) introduced above. Substituting
$t=t_{\pm}$ to the both sides of the equality and using (3.56), we get the
equality

${\boldsymbol{g}}_{N}^{T}(t)=-{\boldsymbol{w}}^{T}(t)\,(i{\boldsymbol{\omega}}+{\boldsymbol{\varOmega}})^{-1},$
(3.65)

that allows us to express the Dirichlet Green’s function from (3.64) via the
Neumann one as

$G_{D}(t,t^{\prime})=G_{N}(t,t^{\prime})-{\boldsymbol{g}}_{N}^{T}(t)\,(i{\boldsymbol{\omega}}+{\boldsymbol{\varOmega}})\,{\boldsymbol{g}}_{N}(t^{\prime}),$
(3.66)

where we use the notation (3.59) for the row
${\boldsymbol{g}}_{N}(t)=[G_{N}(t,t_{+})\;\,G_{N}(t,t_{-})]$ and its
transpose. Using (3.56) once again, we can write down the expression for the
block matrix of boundary values of the Neumann function ${\boldsymbol{g}}_{N}$
at both ends of the time segment (double bar denoting the restriction of both
arguments to $t_{\pm}$)

${\boldsymbol{G}}_{N}\|=\\!\begin{bmatrix}[r]\;G_{N}(t_{+},t_{+})&G_{N}(t_{+},t_{-})\;\\\
\;G_{N}(t_{-},t_{+})&G_{N}(t_{-},t_{-})\;\end{bmatrix}\\!=\\!(i{\boldsymbol{\omega}}+{\boldsymbol{\varOmega}})^{-1}.$
(3.67)

###  3.5 Canonical quantization

Before proceeding to the canonical quantization of the theory (3.1), whose
Hamiltonian formalism was constructed in the previous subsection, let us make
a more specific choice of basis functions, which is more convenient for the
quantization purposes. We first choose the basis functions $v_{\pm}(t)$ real,
and such that the matrix $\mathcal{D}$ defined by (3.24) has a canonical form,
$\mathcal{D}=\mathcal{P}$. Together with the reality of $\phi(t)$ this implies
also the reality of the corresponding integration constants $\alpha^{\pm}$.
Next, we combine these basis functions and integration constants into the
following complex conjugated pairs

$\displaystyle\begin{bmatrix}\,v_{+}&&v_{-}\end{bmatrix}\mapsto\begin{bmatrix}\,v&&v^{*}\end{bmatrix}=\frac{1}{\sqrt{2}}\begin{bmatrix}\,v_{+}&&v_{-}\end{bmatrix}\begin{bmatrix}[c]I&I&\\\
-iI&iI&\end{bmatrix},$ (3.68) $\displaystyle\begin{bmatrix}\,\alpha^{+}\\\
\,\alpha^{-}\end{bmatrix}\;\mapsto\;\begin{bmatrix}[l]&\alpha\\\
&\alpha^{*}\end{bmatrix}=\frac{1}{\sqrt{2}}\begin{bmatrix}[r]\,I&iI\,\\\
\,I&-iI\,\end{bmatrix}\begin{bmatrix}\,\alpha^{+}\\\
\,\alpha^{-}\end{bmatrix}.$ (3.69)

After this change of basis, the matrix $\mathcal{D}$ becomes

$\mathcal{D}\mapsto i\mathcal{P}=\begin{bmatrix}[c]\phantom{-}0&&iI&\\\
-iI&&0&\end{bmatrix}.$ (3.70)

According to (3.27a), this implies the following pairwise Poisson brackets of
$\alpha$, $\alpha^{*}$

$\\{\alpha,\alpha^{*}\\}=-\\{\alpha^{*},\alpha\\}=-iI,\;\\{\alpha,\alpha\\}=\\{\alpha^{*},\alpha^{*}\\}=0.$
(3.71)

In terms of the new basis functions, the evolution law (3.14) of the field and
the canonical momentum becomes

$\begin{bmatrix}\phi(t)\\\ \pi(t)\end{bmatrix}=\begin{bmatrix}v(t)&v^{*}(t)\\\
Wv(t)&Wv^{*}(t)\end{bmatrix}\begin{bmatrix}[l]&\alpha\\\
&\alpha^{*}\end{bmatrix}.$ (3.72)

The equation (3.20) takes the form

$\mathcal{M}^{-1}(t)=i\,\mathcal{P}\,\mathcal{M}^{T}\\!(t)\,\mathcal{P},\quad\mathcal{M}(t)=\begin{bmatrix}v(t)&v^{*}(t)\\\
Wv(t)&Wv^{*}(t)\end{bmatrix}$ (3.73)

and allows to invert the equality (3.72) as

$\begin{bmatrix}[l]\;\alpha\\\
\;\alpha^{*}\end{bmatrix}=i\,\mathcal{P}\,\mathcal{M}^{T}\\!(t)\,\mathcal{P}\begin{bmatrix}\,\phi(t)\,\\\
\,\pi(t)\,\end{bmatrix}.$ (3.74)

Evaluating at $t=t_{-}$ and substituting back to (3.72), one obtains evolving
phase space variables in terms of the basis functions $v(t)$, $v^{*}(t)$ and
initial data,

$\begin{bmatrix}\,\phi(t)\,\\\
\,\pi(t)\,\end{bmatrix}=i\,\mathcal{M}(t)\,\mathcal{P}\,\mathcal{M}^{T}\\!(t_{-})\,\mathcal{P}\begin{bmatrix}\phi(t_{-})\\\
\pi(t_{-})\end{bmatrix}.$ (3.75)

Now, we are ready perform the canonical quantization of the system under
consideration, whose Hamiltonian form was obtained in the previous subsection.
We will quantize it in the Heisenberg picture. Thus, we map the solutions of
the Hamiltonian equations to the corresponding Heisenberg operators, i.e.
$\phi(t),\pi(t)\mapsto\hat{\phi}(t),\hat{\pi}(t)$, whereas the Poisson bracket
is replaced by the commutator times the factor $i$, so that we obtain the
equal-time quantum commutation relations
$[\hat{\phi}(t),\hat{\pi}(t)]=i\hat{I}$, where $\hat{I}$ is the identity
operator in the Hilbert space. Thus, the Hamiltonian equations (3.10) are
mapped to the corresponding Heisenberg equations, defining the evolution of
the operators

$\displaystyle\frac{d}{dt}\hat{\phi}(t)$
$\displaystyle=-i[\hat{\phi}(t),\hat{H}(t)]=A^{-1}(\hat{\pi}(t)-B\hat{\phi}(t))$
(3.76a) $\displaystyle\frac{d}{dt}\hat{\pi}(t)$
$\displaystyle=-i[\hat{\pi}(t),\hat{H}(t)]$
$\displaystyle=B^{T}A^{-1}\bigl{(}\hat{\pi}(t)-B\hat{\phi}(t)\bigr{)}+C\hat{\phi}(t).$
(3.76b)

Here $\hat{H}(t)$ is the classical Hamiltonian (3.9) where the field and the
momentum are replaced by the corresponding Heisenberg operators.

Linearity of the system obviously implies that the classical Hamiltonian and
the Heisenberg equations formally coincide and their solutions are in one-to-
one correspondence. In particular, the relation (3.8) between the field $\phi$
and its conjugate momentum $\pi$ is literally the same at classical and
quantum levels $\hat{\pi}(t)=W\hat{\phi}(t)$. Formal coincidence and linearity
of the Hamiltonian and Heisenberg equations allow one to obtain the solution
of the latter ones from classical equations (3.75)

$\begin{bmatrix}\,\hat{\phi}(t)\,\\\
\,\hat{\pi}(t)\,\end{bmatrix}=i\,\mathcal{M}(t)\,\mathcal{P}\,\mathcal{M}^{T}\\!(t_{-})\,\mathcal{P}\begin{bmatrix}\,\hat{\phi}(t_{-})\,\\\
\hat{\pi}(t_{-})\end{bmatrix}.$ (3.77)

Similarly, our quantization procedure implies that the integration constants
$\alpha$, $\alpha^{*}$ are in one-to-one correspondence to the
creation/annihilation operators $\hat{a}$, $\hat{a}^{\dagger}$. According to
(3.72) the operators $\hat{\phi}(t)$, $\hat{\pi}(t)$ are decomposed in the
creation/annihilation operators as

$\begin{bmatrix}\,\hat{\phi}(t)\,\\\
\hat{\pi}(t)\end{bmatrix}=\begin{bmatrix}v(t)&v^{*}(t)\\\
\,Wv(t)&Wv^{*}(t)\,\end{bmatrix}\begin{bmatrix}[l]\;\hat{a}\\\
\;\hat{a}^{\dagger}\end{bmatrix},$ (3.78)

that can be inverted similar to (3.74) as

$\begin{bmatrix}[l]\;\hat{a}\\\
\;\hat{a}^{\dagger}\end{bmatrix}=i\,\mathcal{P}\,\mathcal{M}^{T}\\!(t)\,\mathcal{P}\begin{bmatrix}\,\hat{\phi}(t)\,\\\
\hat{\pi}(t)\end{bmatrix}.$ (3.79)

The fact that $\hat{a}$ and $\hat{a}^{\dagger}$ are indeed Hermitian conjugate
to each other immediately follows from the Hermicity of $\hat{\phi}(t)$.
Indeed, comparing $\hat{\phi}(t)$ to it’s conjugate

$\displaystyle\begin{aligned}
\hat{\phi}(t)&=v(t)\,\hat{a}+v^{*}(t)\,\hat{a}^{\dagger},\\\
\hat{\phi}^{\dagger}(t)&=\bigl{(}v(t)\,\hat{a}+v^{*}(t)\,\hat{a}^{\dagger}\bigr{)}^{\dagger}=v^{*}(t)\,\hat{a}^{\dagger}+v(t)\,\hat{a},\end{aligned}$
(3.80)

we find the coincidence, for which the choice (3.68) of two complex conjugated
basis functions is crucial. The commutation relations of the
creation/annihilation operators are inherited from the Poisson brackets
(3.71), namely

$\displaystyle\begin{aligned} &[\hat{a}^{A},\hat{a}^{\dagger
B}]=-[\hat{a}^{\dagger B},\hat{a}^{B}]=\delta^{AB}\hat{I},\\\
&[\hat{a}^{A},\hat{a}^{B}]=[\hat{a}^{\dagger A},\hat{a}^{\dagger
B}]=0.\end{aligned}$ (3.81)

Though the explicit solution to the Heisenberg equations (3.77) is obtained,
we still have no the expression for the evolution operator in a closed form.
The latter solves the Schroedinger equation in the form of the chronological
ordering,

$i\frac{d}{dt}\hat{U}(t,t^{\prime})=\hat{H}_{S}(t)\hat{U}(t,t^{\prime}),\;\hat{U}(t_{+},t_{-})=\mathrm{T}e^{-i\int_{t_{-}}^{t_{+}}dt\,\hat{H}_{S}(t)},$
(3.82)

where $\hat{H}_{S}(t)$ is the Hamiltonian in the Schroedinger picture, so that
its time dependence is only due to time-dependent coefficients $A$, $B$, and
$C$. The operators $\hat{\phi}$, $\hat{\pi}$ in the Schroedinger picture are
identified with the Heisenberg ones, evaluated at the initial time

$\hat{\phi}=\hat{\phi}(t_{-}),\qquad\hat{\pi}=\hat{\pi}(t_{-}).$ (3.83)

In the presence of the source, $H\mapsto H-J^{T}\phi$, the solution (3.78) to
the Heisenberg equation generalizes to

$\hat{\phi}_{J}(t)=\hat{\phi}(t)-\int_{t_{-}}^{t_{+}}dt^{\prime}\,G_{R}(t,t^{\prime})J(t^{\prime}),$
(3.84)

that can be easily derived from (3.37). Here $\hat{\phi}(t)$ is the solution
(3.78) to the sourceless Heisenberg equation. The Schroedinger equation for
the evolution operator in the presence of the source

$i\frac{d}{dt}\hat{U}_{J}(t,t^{\prime})=\bigl{(}\hat{H}_{S}(t)-J^{T}(t)\hat{\phi}\bigr{)}\,\hat{U}_{J}(t,t^{\prime}),\;\hat{U}_{J}(t,t)=\hat{I}$
(3.85)

can be solved by the chronological Dyson $t$-exponent (cf. Eq. (2.3)),

$\hat{U}_{J}(t_{+},t_{-})=\mathbb{\mathrm{T}}\,e^{-i\int^{t_{+}}_{t_{-}}dt\,\bigl{(}\hat{H}_{S}(t)-J(t)\hat{\phi}_{S}\bigr{)}}.$
(3.86)

Another representation for the evolution operator follows from functional
integral formalism. If one introduces the coordinate representation,
associated with the Scroedinger operators (3.83),

$\hat{\phi}|\varphi\rangle=\varphi\,|\varphi\rangle,\;\hat{\pi}|\varphi\rangle=i\frac{\partial}{\partial\varphi}|\varphi\rangle,\;\hat{I}=\int
d\varphi\,|\varphi\rangle\langle\varphi|,$ (3.87)

then the matrix elements of $\hat{U}_{J}$ in the coordinate representation
express in terms of the following functional integral

$\displaystyle\langle\varphi_{+}|\hat{U}_{J}(t_{+},t_{-})\,|\varphi_{-}\rangle$
$\displaystyle\quad=\\!\\!\\!\\!\\!\\!\int\limits_{\phi(t_{\pm})=\varphi_{\pm}}\\!\\!\\!\\!\\!\\!{\mathcal{D}}\phi\;\exp\Bigl{\\{}iS[\phi]+i\int_{t_{-}}^{t_{+}}dt\,J^{T}(t)\phi(t)\Bigr{\\}}.$
(3.88)

Since the action (3.1) is quadratic in the field $\varphi$, the latter
integral is Gaussian, so it can be calculated explicitly. We will do this in
the next section with the use of the saddle point method.

###  3.6 Bogoliubov transformations

In the previous subsection we made the choice (3.68) of basis functions which
implies a simple form of the commutation relations (3.81) for the
creation/annihilation operators. It will be useful to study the
transformations, preserving these commutation relations. For this purpose, let
us define a new set of creation/annihilation operators $\hat{b}$,
$\hat{b}^{\dagger}$ as a linear combination of the initial ones,

$\displaystyle\begin{bmatrix}[l]\,\hat{b}\\\
\,\hat{b}^{\dagger}\end{bmatrix}=\begin{bmatrix}[l]&U&V\\\
&V^{*}&U^{*}\end{bmatrix}\begin{bmatrix}[l]\,\hat{a}\\\
\,\hat{a}^{\dagger}\end{bmatrix},$ (3.89)

where $U$, $V$ are referred to as the matrices of Bogoliubov transformations.
Demanding that the commutation relations of the new cration/annihilation
operators coincides with those of the initial ones (3.81), one obtains the
equality

$\begin{bmatrix}[l]&U&V\\\ &V^{*}&U^{*}\end{bmatrix}\begin{bmatrix}[r]0&&I&\\\
-I&&0&\end{bmatrix}\begin{bmatrix}&U^{T}&V^{\dagger}\,\\\
&V^{T}&U^{\dagger}\end{bmatrix}=\begin{bmatrix}[r]0&&I&\\\
-I&&0&\end{bmatrix}$ (3.90)

so that $U$ and $V$ should satisfy

$\displaystyle UU^{\dagger}-VV^{\dagger}=I,\qquad UV^{T}-VU^{T}=0.$ (3.91)

Thus, the field operator $\hat{\phi}(t)$ has two equivalent decompositions

$\displaystyle\hat{\phi}(t)=v(t)\,\hat{a}+v^{*}(t)\,\hat{a}^{\dagger}=\tilde{v}(t)\,\hat{b}+\tilde{v}^{*}(t)\,\hat{b}^{\dagger},$
(3.92)

where the new set of the basis functions $\tilde{v}(t)$, $\tilde{v}^{*}(t)$ is
related to the initial one via the following relation

$\begin{bmatrix}\,v&&v^{*}\,\end{bmatrix}=\begin{bmatrix}\,\tilde{v}&&\tilde{v}^{*}\,\end{bmatrix}\begin{bmatrix}[l]&U&V\\\
&V^{*}&U^{*}\end{bmatrix}$ (3.93)

or, in more explicit form

$v=\tilde{v}\,U+\tilde{v}^{*}V^{*}.$ (3.94)

Equality (3.90) leads to the following formula for the inverse matrix of the
Bogoliubov transformation coefficients

$\displaystyle\begin{bmatrix}[l]&U&V\\\
&V^{*}&U^{*}\end{bmatrix}^{-1}=\begin{bmatrix}[r]U^{\dagger}&-V^{T}\\\
-V^{\dagger}&U^{T}\end{bmatrix}$ (3.95)

so that (3.94) can be inverted as

$\tilde{v}=v\,U^{\dagger}-v^{*}V^{\dagger}.$ (3.96)

Now, let us solve the inverse problem. Namely, suppose we have two sets of the
basis functions $v(t)$, $v^{*}(t)$ and $\tilde{v}(t)$, $\tilde{v}^{*}(t)$,
such that the commutation relation of the corresponding creation/annihilation
operators are of the canonical form (3.81). The question is what are the
Bogoliubov coefficients relating these two sets? To find the explicit form of
the coefficients, let us introduce the following inner product in the space of
solutions of the equation (3.5)

$(\phi_{1},\phi_{2})=i\,\phi_{1}^{\dagger}(W\phi_{2})-i(W\phi_{1})^{\dagger}\phi_{2}.$
(3.97)

This is constant if $\phi_{1}$, $\phi_{2}$ solve (3.5) due to the Wronskian
property (3.22), together with the fact that the operator $F$, defining
equations of motion (and the Wronskian $W$) is real. The inner product (3.97)
is usually referred to as the Klein-Gordon type inner product. The choice
(3.68) of the basis functions implies the following normalization with respect
to this inner product

$(v,v)=-(v^{*},v^{*})=I,\quad(v^{*},v)=0,$ (3.98)

and the same for $\tilde{v}(t)$, $\tilde{v}^{*}(t)$. Projecting the equality
(3.94) onto $\tilde{v}$, and using the property
$(v_{1},v_{2})^{*}=-(v_{1}^{*},v_{2}^{*})$, one obtains the explicit
expressions for the Bogoliubov coefficients

$U=(\tilde{v},v),\qquad V=(\tilde{v},v^{*}).$ (3.99)

If both the old and the new sets of basis functions satisfy the Neumann
conditions with different frequency matrices $\omega$ and $\tilde{\omega}$ at
the initial moment of time

$(iW-\omega)v(t_{-})=0,\quad(iW-\tilde{\omega})\tilde{v}(t_{-})=0,$ (3.100)

one can find the explicit expressions for $U$, $V$ in terms of $\omega$ and
$\tilde{\omega}$. Let us first write down the normalization conditions (3.98)
explicitly

$(v,v)=v^{\dagger}(\omega+\omega^{\dagger})v=I,\;(v^{*},v)=v^{T}(\omega-\omega^{T})v=0,$
(3.101)

where all quantities are evaluated at $t=t_{-}$. The same equations hold for
$\tilde{v}(t)$ and $\tilde{v}^{*}(t)$. The second equation implies that the
matrices $\omega$ and $\tilde{\omega}$ are symmetric (in view of invertibility
of the matrix $v(t)$ at a generic moment of time), whereas the first equation
allows one to fix the initial value of the basis functions as

$v(t_{-})=\frac{1}{\sqrt{2\omega_{\mathrm{re}}}},\quad\tilde{v}(t_{-})=\frac{1}{\sqrt{2\tilde{\omega}_{\mathrm{re}}}},$
(3.102)

where $\omega_{\mathrm{re}}$ and $\tilde{\omega}_{\mathrm{re}}$ are the real
parts of $\omega$ and $\tilde{\omega}$, respectively. Using (3.99) with the
inner product defined in (3.97) one finds the following expressions for
Bogoliubov coefficients relating two sets of Neumann basis functions with
different frequency matrices,

$\displaystyle
U=\frac{1}{\sqrt{2\tilde{\omega}_{\mathrm{re}}}}(\omega+\tilde{\omega}^{\dagger})\frac{1}{\sqrt{2\omega_{\mathrm{re}}}},$
(3.103a) $\displaystyle
V=\frac{1}{\sqrt{2\tilde{\omega}_{\mathrm{re}}}}(\tilde{\omega}^{\dagger}-\omega^{\dagger})\frac{1}{\sqrt{2\omega_{\mathrm{re}}}}.$
(3.103b)

###  3.7 Fock space and the coherent states

Once the basis functions $v(t)$, $v^{*}(t)$ are chosen, we can define the Fock
space, associated to the corresponding creation/annihilation operators.
Namely, introducing the vacuum state $|0\rangle$ as

$\hat{a}|0\rangle=0,$ (3.104)

one defines the Fock space as a linear space spanned by

$|A_{1},A_{2},\ldots A_{n}\rangle\equiv\hat{a}^{\dagger A_{1}}\hat{a}^{\dagger
A_{2}}\ldots\hat{a}^{\dagger A_{n}}|0\rangle.$ (3.105)

Next, let us obtain the coordinate representation of the Fock states. For this
purpose, we rewrite (3.79) explicitly for $t=t_{-}$ as

$\displaystyle\begin{bmatrix}[l]&\hat{a}\\\
&\hat{a}^{\dagger}\end{bmatrix}=\frac{1}{\sqrt{2\omega_{\mathrm{re}}}}\begin{bmatrix}[r]&\omega^{*}&iI&\\\
&\omega^{\phantom{*}}&-iI&\end{bmatrix}\begin{bmatrix}&\hat{\phi}&\\\
&\hat{\pi}&\end{bmatrix},$ (3.106)

where $\omega$ is given in terms of the positive frequency basis function

$\displaystyle\omega=\bigl{(}iWv\bigr{)}v^{-1}\bigr{|}_{t=t_{-}},$ (3.107)

and rewrite the definition (3.104) of the vacuum state in the coordinate
representation (3.87) as

$\frac{1}{\sqrt{2\omega_{\mathrm{re}}}}\left(\frac{\partial}{\partial\varphi}+\omega^{*}\varphi\right)\langle\varphi|0\rangle=0.$
(3.108)

Therefore, up to $\pi$-dependent normalization the wavefunction of a vacuum
reads

$\langle\varphi|0\rangle=\bigl{(}\det\omega_{\mathrm{re}}\bigr{)}^{\frac{1}{4}}\,\exp\biggl{\\{}-\frac{1}{2}\varphi^{T}\omega^{*}\varphi\biggr{\\}}.$
(3.109)

Coordinate representation of the excited states can be found by using the
definition (3.105) and the expression for $\hat{a}^{\dagger}$ in the
coordinate representation.

Similarly, one can define the coherent states $|\alpha\rangle$ as eigenstates
of the annihilation operator

$\hat{a}|\alpha\rangle=\alpha|\alpha\rangle.$ (3.110)

Projecting the definition on the coordinate representation basis vector
$|\varphi\rangle$, one obtains the equation

$\left(\frac{\partial}{\partial\varphi}+\omega^{*}\varphi\right)\langle\varphi|\alpha\rangle=\sqrt{2\omega_{\mathrm{re}}}\langle\varphi|\alpha\rangle,$
(3.111)

whose integration gives the (unnormalized) solution

$\langle\varphi|\alpha\rangle=\exp\biggl{\\{}-\frac{1}{2}\varphi^{T}\omega^{*}\varphi+\alpha^{T}\\!\sqrt{2\omega_{\mathrm{re}}}\,\varphi-\frac{1}{2}\alpha^{T}\alpha\biggr{\\}}.$
(3.112)

For this normalization we have the following expression for the Fock states in
terms of coherent state

$|A_{1},A_{2},\ldots
A_{n}\rangle=\left.\frac{\partial^{n}}{\partial\alpha^{A_{1}}\,\partial\alpha^{A_{2}}\ldots\partial\alpha^{A_{n}}}|\alpha\rangle\right|_{\alpha=0}.$
(3.113)

Coherent states allows to perform a partition of unity as

$\hat{I}=\int
d\alpha^{*}\,d\alpha\,e^{-\alpha^{\dagger}\alpha}|\alpha\rangle\langle\alpha|.$
(3.114)

##  4 Generating functional in the path integral formalism

We begin our derivation of the in-in Green’s function generating functional
for the theory, defined in the previous section, by the physical motivation
and the definition of an arbitrary Gaussian initial state. After that, we
derive the corresponding two-component Green’s functions. As we will observe,
there is an ambiguity in the definition of these Green’s function,
parameterized by a matrix, defining initial conditions for the modes, employed
in the mode expansion of the field operators. There is no any a priory
preferred choice, fixing this ambiguity. However, being motivated by the
simple harmonic oscillator in a thermal state, we make a choice of the modes
such that the resulting Green’s function has the form and the properties, very
close to that of the Green’s functions for the equilibrium system in a thermal
state. Further, we introduce the notion of the quasi-thermal state, which is a
very particular case of the Gaussian state, in which the properties of the
Green’s functions become even more closer to those of the thermal ones, in
particular, satisfying the Kubo-Martin-Schwinger (KMS) condition.

###  4.1 Gaussian states

Our goal is to find the explicit and useful form of the generating functional

$Z[J_{1},J_{2}]=\operatorname{tr}\left[\hat{U}_{J_{1}}(T,0)\,\hat{\rho}\,\,\hat{U}^{\dagger}_{-J_{2}}(T,0)\right]$
(4.1)

of in-in correlation functions

$\displaystyle\operatorname{tr}\left[\hat{\rho}\,\bar{\mathrm{T}}\bigl{(}\hat{\phi}(t^{\prime}_{1})\ldots\phi(t^{\prime}_{m})\bigr{)}\,\mathrm{T}\bigl{(}\hat{\phi}(t_{1})\ldots\phi(t_{n})\bigr{)}\right]$
$\displaystyle=\frac{i^{m-n}}{Z}\left.\frac{\delta^{n+m}Z[J_{1},J_{2}]}{\delta
J_{1}(t_{1})\ldots\delta J_{1}(t_{n})\,\delta
J_{2}(t^{\prime}_{1})\ldots\delta
J_{2}(t^{\prime}_{m})}\right|_{J_{1}=J_{2}=0}.$ (4.2)

where $\hat{U}_{J}$ are the evolution operators subject to equation (3.85)
with different sources $J_{1}$ and $-J_{2}$, whereas $\mathrm{T}$ and
$\bar{\mathrm{T}}$ denote chronological and anti-chronological ordering,
respectively. The relation between (4.1) and the correlation functions (4.2)
obviously follows from (3.86). The basic elements are the two-point
correlation functions, namely

$\displaystyle iG_{\mathrm{T}}(t,t^{\prime})$
$\displaystyle\equiv\operatorname{tr}\left[\hat{\rho}\,\mathrm{T}\bigl{(}\hat{\phi}(t)\hat{\phi}(t^{\prime})\bigr{)}\right],$
(4.3a) $\displaystyle iG_{\bar{\mathrm{T}}}(t,t^{\prime})$
$\displaystyle\equiv\operatorname{tr}\left[\hat{\rho}\,\bar{\mathrm{T}}\bigl{(}\hat{\phi}(t)\hat{\phi}(t^{\prime})\bigr{)}\right],$
(4.3b) $\displaystyle iG_{<}(t,t^{\prime})$
$\displaystyle\equiv\operatorname{tr}\left[\hat{\rho}\,\hat{\phi}(t)\hat{\phi}(t^{\prime})\right],$
(4.3c)

where $G_{\mathrm{T}}$, $G_{\bar{\mathrm{T}}}$, and $G_{<}$ are Feynman, anti-
Feynman and Wightman Green’s functions, respectively.

The density matrix $\hat{\rho}$ is assumed to be the Hermitian positive-
definite operator of unit trace. Inserting the partition of unity in the
coordinate representation to the definition (4.1) of the generating functional
three times, and using the path integral representation (3.5) of the evolution
operator, one obtains the following expression for the generating functional

$\displaystyle Z[J_{1},J_{2}]=\int
d\varphi_{+}\,d\varphi_{-}\;\rho(\varphi_{+},\varphi_{-})$
$\displaystyle\;\;\times\int{\mathcal{D}}\phi_{1}\,{\mathcal{D}}\phi_{2}\;\exp\Biggl{\\{}iS[\phi_{1}]-iS[\phi_{2}]$
$\displaystyle\;\;+i\int_{0}^{T}dt\,\left.\bigl{(}J_{1}^{T}\phi_{1}+J_{2}^{T}\phi_{2}\bigr{)}\Biggr{\\}}\,\right|_{\begin{subarray}{c}\phi_{1}(T)=\phi_{2}(T),\\\
\phi_{1}(0)=\varphi_{+},\,\phi_{2}(0)=\varphi_{-}\end{subarray}}\\!\\!\\!\\!,$
(4.4)

where the integration over $\phi_{1,2}(t)$ runs with the indicated boundary
conditions and we introduce the notation for the coordinate representation of
the density matrix
$\rho(\varphi_{+},\varphi_{-})=\langle\varphi_{+}|\mathinner{\hat{\rho}}|\varphi_{-}\rangle$.

Now, we restrict ourselves with the Gaussian density matrices, i.e. those
whose coordinate representation has the form of the Gaussian exponent

$\rho({\boldsymbol{\varphi}})=\frac{1}{Z}\exp\left\\{-\frac{1}{2}{\boldsymbol{\varphi}}^{T}{\boldsymbol{\varOmega}}\,{\boldsymbol{\varphi}}+{\boldsymbol{j}}^{T}\\!{\boldsymbol{\varphi}}\right\\},\;{\boldsymbol{\varphi}}=\begin{bmatrix}\,\varphi_{+}\\\
\,\varphi_{-}\end{bmatrix},$ (4.5)

where the matrix ${\boldsymbol{\varOmega}}$, and the vector ${\boldsymbol{j}}$
play the role of the parameters of $\hat{\rho}$, and normalization constant
$1/Z$ is independent of ${\boldsymbol{\varphi}}$. The Hermitian property of
the density matrix,
$\langle\varphi_{+}|\mathinner{\hat{\rho}}|\varphi_{-}\rangle=\langle\varphi_{-}|\mathinner{\hat{\rho}}|\varphi_{+}\rangle^{*}$,
which in the coordinate representation reads

$\rho(\varphi_{+},\varphi_{-})=\rho^{*}(\varphi_{-},\varphi_{+}),$ (4.6)

implies the following conditions on ${\boldsymbol{\varOmega}}$ and
${\boldsymbol{j}}$

${\boldsymbol{X}}\,{\boldsymbol{\varOmega}}\,{\boldsymbol{X}}={\boldsymbol{\varOmega}}^{*},\quad{\boldsymbol{X}}\,{\boldsymbol{j}}={\boldsymbol{j}}^{*},\quad{\boldsymbol{X}}\equiv\begin{bmatrix}\;\,0&&I\;\\\
\;I&&0\,\end{bmatrix},$ (4.7)

or, in a more explicit block-matrix form

${\boldsymbol{j}}=\begin{bmatrix}[l]\;\,j\\\
\;\,j^{*}\end{bmatrix},\;{\boldsymbol{\varOmega}}=\begin{bmatrix}[l]\;\,R\,&S\\\
\;\,S^{*}&R^{*}\end{bmatrix},\;R=R^{T},\;S=S^{\dagger}.$ (4.8)

Normalizability of $\hat{\rho}$ implies that the real part of the sum $R+S$ is
positive-definite. The case in which the matrix $S$ is non-vanishing
corresponds to the mixed states, i.e. such that
$\hat{\rho}^{2}\neq\hat{\rho}$. The role of the linear term in the exponential
in (4.5) is two-fold. Firstly, ${\boldsymbol{j}}$ defines non-vanishing mean
value of the field operator. Secondly, it can also be used to introduce non-
linearities to the density matrix, namely, by differentiating it with respect
to ${\boldsymbol{j}}$. The typical example of the (pure) Gaussian state is the
vacuum state (3.104), i.e. $\hat{\rho}=|0\rangle\langle 0|$, associated with
some choice of the annihilation operator, for which $R=\omega^{*}$, $S=0$, and
${\boldsymbol{j}}$ = 0. Another example of the pure Gaussian state is the
coherent state (3.110) whose density matrix reads
$\hat{\rho}=|\alpha\rangle\langle\alpha|$, and $R=\omega^{*}$, $S=0$ again,
but $j=\sqrt{2\omega_{\mathrm{re}}}\alpha^{*}$.

###  4.2 In-in boundary value problem

Substituting the general Gaussian density matrix to (4.1), one obtains

$Z[J_{1},J_{2}]=\\!\\!\\!\\!\\!\\!\int\limits_{\begin{subarray}{c}\phi_{1}(T)=\phi_{2}(T)\\\
{\boldsymbol{\phi}}(0)={\boldsymbol{\varphi}}\end{subarray}}\\!\\!\\!\\!\\!\\!{\mathcal{D}}\phi_{1}\,{\mathcal{D}}\phi_{2}\,\exp\Biggl{\\{}iS[\phi_{1}]-iS[\phi_{2}]\\\
\,\,\,\,+i\int_{0}^{T}\\!dt\bigl{(}J_{1}^{T}\phi_{1}+J_{2}^{T}\phi_{2}\bigr{)}\\!-\\!\frac{1}{2}{\boldsymbol{\varphi}}^{T}{\boldsymbol{\varOmega}}\,{\boldsymbol{\varphi}}\\!+\\!{\boldsymbol{j}}^{T}\\!{\boldsymbol{\varphi}}\Biggr{\\}},$
(4.9)

where we put the boundary points $\varphi_{\pm}$, $\varphi^{\prime}$,
appearing in (4.1) to the functional integration measure, and also omit the
constant normalization factor of the density matrix.

We will compute the integral (4.9) representing the generating functional,
with the use of saddle point method. The latter turns out to be exact since
the integral has the Gaussian form. First of all, we introduce the notations
for the block-matrix operators acting on columns of fields and sources (2.7)
introduced in Section 2,

$\displaystyle{\boldsymbol{F}}=\begin{bmatrix}\;F&0\\\
0&-F\,\end{bmatrix},\quad{\boldsymbol{W}}=\begin{bmatrix}\;W&0\\\
0&-W\,\end{bmatrix},$ (4.10)

so that the sum of the actions for $\phi_{1}$ and $\phi_{2}$ in (4.9) can be
rewritten in the joint form

$\displaystyle{\boldsymbol{S}}[{\boldsymbol{\phi}}]$
$\displaystyle=\frac{1}{2}\int_{0}^{T}dt\,{\boldsymbol{\phi}}^{T}\overset{\leftrightarrow}{{\boldsymbol{F}}}{\boldsymbol{\phi}}$
$\displaystyle=\frac{1}{2}\int_{0}^{T}dt\,{\boldsymbol{\phi}}^{T}{\boldsymbol{F}}{\boldsymbol{\phi}}+\frac{1}{2}{\boldsymbol{\phi}}^{T}{\boldsymbol{W}}{\boldsymbol{\phi}}\,\biggr{|}_{0}^{T}.$
(4.11)

This allows us to treat the underlying equations of motion, Green’s functions,
etc. in exactly the same way as the original theory with the action (3.1),
except that now the field content is doubled. In terms of the new notations
the expression for the generating functional is given by Eqs.(2.8)–(2.10) of
Section 2.

The saddle point equation obtained by varying the exponential of this double-
field action (2.8) with respect to all fields including the boundary values at
$t=0$ and $t=T$ reads

$\displaystyle\delta\left\\{i{\boldsymbol{S}}[{\boldsymbol{\phi}}]+i\int_{0}^{T}dt\,{\boldsymbol{J}}^{T}{\boldsymbol{\phi}}-\frac{1}{2}{\boldsymbol{\varphi}}^{T}{\boldsymbol{\varOmega}}\,{\boldsymbol{\varphi}}+{\boldsymbol{j}}^{T}\\!{\boldsymbol{\varphi}}\right\\}$
$\displaystyle\qquad\quad=\int_{0}^{T}dt\,\delta{\boldsymbol{\phi}}^{T}({\boldsymbol{F}}{\boldsymbol{\phi}}+{\boldsymbol{J}})+i\,\delta{\boldsymbol{\phi}}^{T}\,{\boldsymbol{W}}{\boldsymbol{\phi}}\bigr{|}_{t=T}$
$\displaystyle\qquad\quad-\delta{\boldsymbol{\varphi}}^{T}\Bigl{[}(i{\boldsymbol{W}}+{\boldsymbol{\varOmega}}){\boldsymbol{\phi}}\bigr{|}_{t=0}-{\boldsymbol{j}}\,\Bigr{]}=0.$
(4.12)

Independent variation of the fields $\delta{\boldsymbol{\phi}}(t)$ in the
interior of the time interval gives equations of motion

${\boldsymbol{F}}{\boldsymbol{\phi}}(t)+{\boldsymbol{J}}(t)=0,$ (4.13)

whereas the variation of the boundary values $\delta{\boldsymbol{\phi}}(T)$
and $\delta{\boldsymbol{\phi}}(0)=\delta{\boldsymbol{\varphi}}$ supply these
equations with the boundary conditions. They read as the following matrix
relations

$\displaystyle(i{\boldsymbol{W}}+{\boldsymbol{\varOmega}})\,{\boldsymbol{\phi}}\bigr{|}_{t=0}={\boldsymbol{j}},$
(4.14)
$\displaystyle\begin{bmatrix}\,I&&I\,\end{bmatrix}\,{\boldsymbol{W}}{\boldsymbol{\phi}}\bigr{|}_{t=T}=0,\;\begin{bmatrix}\,I&-I\,\end{bmatrix}{\boldsymbol{\phi}}\bigr{|}_{t=T}=0,$
(4.15)

where we took into account that in view of $\phi_{1}(T)=\phi_{2}(T)$ the
variation
$\delta{\boldsymbol{\phi}}^{T}(T)=\delta\phi_{1}^{T}(T)\begin{bmatrix}\,I&I\,\end{bmatrix}$,
and the boundary conditions at $t=T$ reduce to the equality of both the fields
and their time derivatives of both $\phi_{1}$ and $\phi_{2}$.

To solve the boundary value problem above, we first find the Green’s function
subject to the homogeneous version of the above boundary conditions, i.e.
those of vanishing ${\boldsymbol{j}}$

$\displaystyle{\boldsymbol{F}}{\boldsymbol{G}}(t,t^{\prime})=\boldsymbol{I}\,\delta(t-t^{\prime}),$
(4.16)
$\displaystyle(i{\boldsymbol{W}}+{\boldsymbol{\varOmega}}){\boldsymbol{G}}(t,t^{\prime})\bigr{|}_{t=0}=0,$
(4.17) $\displaystyle\begin{aligned}
&\begin{bmatrix}\,I&&I\,\end{bmatrix}{\boldsymbol{W}}{\boldsymbol{G}}(t,t^{\prime})\bigr{|}_{t=T}=0,\\\
&\begin{bmatrix}\,I&-I\,\end{bmatrix}{\boldsymbol{G}}(t,t^{\prime})\bigr{|}_{t=T}=0.\end{aligned}$
(4.18)

We can construct the Green’s function ${\boldsymbol{G}}$ solving the problem
above, out of the basis functions ${\boldsymbol{v}}_{\pm}$. These basis
function should solve the homogeneous equation, and satisfy the same boundary
conditions as those of the Green’s function,

$\displaystyle{\boldsymbol{F}}{\boldsymbol{v}}_{\pm}(t)=0,\quad(i{\boldsymbol{W}}+{\boldsymbol{\varOmega}})\,{\boldsymbol{v}}_{-}(t)\bigr{|}_{t=0}=0,$
(4.19) $\displaystyle\begin{aligned}
&\begin{bmatrix}\,I&&I\,\end{bmatrix}{\boldsymbol{W}}{\boldsymbol{v}}_{+}(t)\bigr{|}_{t=T}=0,\\\
&\begin{bmatrix}\,I&-I\,\end{bmatrix}{\boldsymbol{v}}_{+}(t)\bigr{|}_{t=T}=0.\end{aligned}$
(4.20)

Applying the generic Green’s function expression (3.3) to the case of the
doubled field content, we obtain the Green’s function ${\boldsymbol{G}}$ in
terms of these basis functions

$\displaystyle{\boldsymbol{G}}(t,t^{\prime})$
$\displaystyle=-{\boldsymbol{v}}_{+}(t)\,{\boldsymbol{\Delta}}_{-+}^{-1}\,{\boldsymbol{v}}^{T}_{-}(t^{\prime})\,\theta(t-t^{\prime})$
$\displaystyle\quad+{\boldsymbol{v}}_{-}(t)\,{\boldsymbol{\Delta}}_{+-}^{-1}\,{\boldsymbol{v}}^{T}_{+}(t^{\prime})\,\theta(t^{\prime}-t)$
(4.21) $\displaystyle{\boldsymbol{\Delta}}_{-+}$
$\displaystyle={\boldsymbol{v}}_{-}^{T}\,{\boldsymbol{W}}{\boldsymbol{v}}_{+}-({\boldsymbol{W}}{\boldsymbol{v}}_{-})^{T}\,{\boldsymbol{v}}_{+}=-{\boldsymbol{\Delta}}_{+-}^{T}.$
(4.22)

###  4.3 Neumann type basis functions and Green’s function representation

However, we do not have the explicit form of basis functions
${\boldsymbol{v}}_{\pm}$. We will construct ${\boldsymbol{v}}_{\pm}$ with the
help of another set basis functions ${\boldsymbol{v}}$, ${\boldsymbol{v}}^{*}$
subject to much simpler boundary conditions

$\displaystyle{\boldsymbol{F}}{\boldsymbol{v}}(t)=0,\quad(i{\boldsymbol{W}}-{\boldsymbol{\omega}}){\boldsymbol{v}}(t)\bigr{|}_{t=0}=0,$
(4.23)
$\displaystyle{\boldsymbol{\omega}}=\begin{bmatrix}\;\,\omega\;&0\phantom{{}^{*}}\\\
\,\;0\;&\omega^{*}\end{bmatrix}.$ (4.24)

Since ${\boldsymbol{W}}$ and ${\boldsymbol{\omega}}$ are block-diagonal, the
basis functions ${\boldsymbol{v}}$, ${\boldsymbol{v}}^{*}$ can be chosen
block-diagonal too, namely

${\boldsymbol{v}}=\begin{bmatrix}[l]&v&&0\\\
&0&&v^{*}\end{bmatrix},\qquad{\boldsymbol{v}}^{*}=\begin{bmatrix}[l]&v^{*}&0&\\\
&0&v&\end{bmatrix}.$ (4.25)

With a real operator $F$ the blocks of these matrices solve the equations
$Fv(t)=0$ and $Fv^{*}(t)=0$, subject to the complex conjugated boundary
conditions

$\displaystyle(iW-\omega)v(t)\bigr{|}_{t=0}=0,\quad(iW+\omega^{*})v^{*}(t)\bigr{|}_{t=0}=0.$
(4.26)

Thus, $v$ and $v^{*}$ are simply the basis functions for the single field
$\phi_{+}$ or $\phi_{-}$ subject to the Neumann boundary conditions introduced
above. We assume that $\omega$ is the symmetric matrix with a positive-
definite real part.

The answer for the basis function ${\boldsymbol{v}}_{+}$ in terms of $v$ and
$v^{*}$ can be easily constructed as

$\displaystyle{\boldsymbol{v}}_{+}={\boldsymbol{v}}+{\boldsymbol{v}}^{*}{\boldsymbol{X}}=\begin{bmatrix}\;v&&v^{*}\\\
\;v&&v^{*}\end{bmatrix},\;{\boldsymbol{X}}=\begin{bmatrix}\,0&&I\,\\\
\,I&&0\,\end{bmatrix},$ (4.27)

while the calculation of ${\boldsymbol{v}}_{-}$ requires more efforts. We will
obtain the answer for ${\boldsymbol{v}}_{-}$ with the use of Bogoliubov
coefficients relating two sets of different Neumann basis functions (3.103) by
treating ${\boldsymbol{v}}_{-}$ as the negative frequency basis function
complex conjugated to its positive frequency counterpart
${\boldsymbol{v}}_{-}^{*}$ satisfying at $t=0$ the boundary condition
$(i{\boldsymbol{W}}-{\boldsymbol{\varOmega}}^{*})\,{\boldsymbol{v}}_{-}^{*}|_{t=0}=0$.
Thus, in accordance with (3.103) the answer for ${\boldsymbol{v}}_{-}$ reads

${\boldsymbol{v}}_{-}={\boldsymbol{v}}^{*}\boldsymbol{U}^{T}-{\boldsymbol{v}}\boldsymbol{V}^{T}.$
(4.28)

where $\boldsymbol{U}$, $\boldsymbol{V}$ are the corresponding Bogoliubov
coefficients

$\displaystyle\boldsymbol{U}=\frac{1}{\sqrt{2{\boldsymbol{\varOmega}}_{\mathrm{re}}}}({\boldsymbol{\varOmega}}+{\boldsymbol{\omega}})\frac{1}{\sqrt{2{\boldsymbol{\omega}}_{\mathrm{re}}}},$
(4.29a)
$\displaystyle\boldsymbol{V}=\frac{1}{\sqrt{2{\boldsymbol{\varOmega}}_{\mathrm{re}}}}({\boldsymbol{\varOmega}}-{\boldsymbol{\omega}}^{*})\frac{1}{\sqrt{2{\boldsymbol{\omega}}_{\mathrm{re}}}}.$
(4.29b)

Here we assume the normalization
${\boldsymbol{v}}(0)=1/\sqrt{2{\boldsymbol{\omega}}_{\mathrm{re}}}$,
${\boldsymbol{v}}_{-}(0)=1/\sqrt{2{\boldsymbol{\varOmega}}_{\mathrm{re}}}$,
and denote the real parts of ${\boldsymbol{\omega}}$ and
${\boldsymbol{\varOmega}}$ respectively as
${\boldsymbol{\omega}}_{\mathrm{re}}$ and
${\boldsymbol{\varOmega}}_{\mathrm{re}}$.

Finally, let us consider the details of the Green’s function
${\boldsymbol{G}}(t,t^{\prime})$ defined by (4.21) for a particular form of
${\boldsymbol{v}}_{\pm}$ we have just built. Its matrix
${\boldsymbol{\Delta}}_{-+}$ given by (4.22) reads

$\displaystyle
i{\boldsymbol{\Delta}}_{-+}=\frac{1}{\sqrt{2{\boldsymbol{\varOmega}}_{\mathrm{re}}}}\Bigl{[}({\boldsymbol{I}}-{\boldsymbol{X}})\,{\boldsymbol{\omega}}+{\boldsymbol{\varOmega}}\,({\boldsymbol{I}}+{\boldsymbol{X}})\Bigr{]}\frac{1}{\sqrt{2{{\boldsymbol{\omega}}_{\mathrm{re}}}}}.$
(4.30)

Next, let us consider separately the first term in (4.21). After the
calculation presented in Appendix B one obtains for it the following form

${\boldsymbol{v}}_{+}(t)\,(i{\boldsymbol{\Delta}}_{-+})^{-1}\,{\boldsymbol{v}}^{T}_{-}(t^{\prime})={\boldsymbol{v}}_{+}(t)\,{\boldsymbol{v}}^{\dagger}(t^{\prime})+{\boldsymbol{v}}_{+}(t)\,{\boldsymbol{\nu}}\,{\boldsymbol{v}}_{+}^{T}(t^{\prime}),$
(4.31)

where we introduce the following symmetric matrix

${\boldsymbol{\nu}}=\Bigl{[}{\boldsymbol{I}}\\!+\\!{\boldsymbol{X}}\\!-\\!\sqrt{2{\boldsymbol{\omega}}_{\mathrm{re}}}\,{\boldsymbol{X}}\,({\boldsymbol{\omega}}\\!+\\!{\boldsymbol{\varOmega}})^{-1}{\boldsymbol{X}}\sqrt{2{\boldsymbol{\omega}}_{\mathrm{re}}}\Bigr{]}^{-1}\\!\\!-{\boldsymbol{X}}.$
(4.32)

Recalling that the second term of the expression (4.21) can be obtained from
the first one by the simultaneous transposition and exchange of time
arguments, and observing that the second term in (4.31) is symmetric under
this transformation, we find that two theta functions sum up to identity, so
that the final expression for the Green’s function reads

$i{\boldsymbol{G}}(t,t^{\prime})=i{\boldsymbol{G}}_{0}(t,t^{\prime})+{\boldsymbol{v}}_{+}(t)\,{\boldsymbol{\nu}}\,{\boldsymbol{v}}_{+}^{T}(t^{\prime}),$
(4.33)

where ${\boldsymbol{G}}_{0}$ is defined as

$\displaystyle i{\boldsymbol{G}}_{0}(t,t^{\prime})$
$\displaystyle={\boldsymbol{v}}_{+}(t)\,{\boldsymbol{v}}^{\dagger}(t^{\prime})\,\theta(t-t^{\prime})+{\boldsymbol{v}}^{*}(t)\,{\boldsymbol{v}}_{+}^{T}(t^{\prime})\,\theta(t^{\prime}-t)$
$\displaystyle={\boldsymbol{v}}(t)\,{\boldsymbol{v}}^{\dagger}(t^{\prime})\,\theta(t-t^{\prime})+{\boldsymbol{v}}^{*}(t)\,{\boldsymbol{v}}^{T}(t^{\prime})\,\theta(t^{\prime}-t)$
$\displaystyle\quad+{\boldsymbol{v}}^{*}(t)\,{\boldsymbol{X}}\,{\boldsymbol{v}}^{\dagger}(t^{\prime}),$
(4.34)

and interpreted as the Green function, corresponding to the vacuum state,
having the density matrix $\hat{\rho}_{0}=|0\rangle\langle 0|$, associated
with the basis functions $v(t)$, $v^{*}(t)$. Indeed, from (3.109) one observes
that the matrix ${\boldsymbol{\varOmega}}$, defining the vacuum density matrix
$\hat{\rho}_{0}$, coincides with ${\boldsymbol{\omega}}^{*}$, i.e.
${\boldsymbol{\varOmega}}={\boldsymbol{\omega}}^{*}$. In this case,
${\boldsymbol{\nu}}$ vanishes due to its definition (4.32), so from (4.33) we
find that ${\boldsymbol{G}}={\boldsymbol{G}}_{0}$.

Substituting the generating functional obtained to (4.2), we observe that for
vanishing ${\boldsymbol{j}}$ the block-matrix components of ${\boldsymbol{G}}$
are composed of the Feynman, anti-Feynman and Wightman Green’s functions
(4.3), namely

${\boldsymbol{G}}(t,t^{\prime})=\begin{bmatrix}[l]G_{\mathrm{T}}(t,t^{\prime})&G_{<}(t,t^{\prime})\\\
G_{>}(t,t^{\prime})&G_{\bar{\mathrm{T}}}(t,t^{\prime})\end{bmatrix},$ (4.35)

where $G_{>}(t,t^{\prime})\equiv G_{<}^{T}(t^{\prime},t)$, and the explicit
form of the block components can be read off form (4.33).

Now we have to find the solution ${\boldsymbol{\phi}}(t)$ of the boundary
value problem (4.13)–(4.15) in order to substitute it into the exponential of
(2.8). The only inhomogeneous boundary conditions in this problem are the
Neumann conditions (4.14), so that the solution is given by the double-field
version of (3.58) with the substitutions $j_{+}\mapsto 0$ (remember that there
is no $j_{+}$ at the point $t=T$) and $j_{-}\mapsto-{\boldsymbol{j}}$. Thus it
reads

${\boldsymbol{\phi}}(t)=i{\boldsymbol{G}}(t,0)\,{\boldsymbol{j}}-\int_{0}^{T}dt^{\prime}\,{\boldsymbol{G}}(t,t^{\prime}){\boldsymbol{J}}(t^{\prime}).$
(4.36)

Substituting it to the exponential of (2.8) then gives Eq. (2.11) advocated in
Section 2.

$\displaystyle
Z[{\boldsymbol{J}}]=\text{const}\times\exp\Biggl{\\{}-\frac{i}{2}\int
dt\,dt^{\prime}\,{\boldsymbol{J}}^{T}\\!(t){\boldsymbol{G}}(t,t^{\prime}){\boldsymbol{J}}(t)$
$\displaystyle\qquad\quad-\int
dt\,{\boldsymbol{J}}^{T}\\!(t){\boldsymbol{G}}(t,0)\,{\boldsymbol{j}}+\frac{i}{2}{\boldsymbol{j}}^{T}{\boldsymbol{G}}(0,0)\,{\boldsymbol{j}}\Biggr{\\}},$
(4.37)

where all time integrations run from $t=0$ to $t=T$. Here the restriction of
${\boldsymbol{G}}(t,t^{\prime})$ to ${\boldsymbol{G}}(t,0)$ does not lead to
essential simplification whereas ${\boldsymbol{G}}(0,0)$ has, as shown in
Appendix C, the following explicit and simple form in terms of the parameters
of the density matrix

$\displaystyle
i{\boldsymbol{G}}(0,0)=\frac{{\boldsymbol{I}}+{\boldsymbol{X}}}{2{\boldsymbol{\varOmega}}_{\mathrm{re}}},$
(4.38)

where the “ratio” of matrices ${\boldsymbol{I}}+{\boldsymbol{X}}$ and
${\boldsymbol{\varOmega}}_{\mathrm{re}}$ is unambiguous because these matrices
are commuting in view of the special form of ${\boldsymbol{\varOmega}}$
subject to the relation
${\boldsymbol{X}}{\boldsymbol{\varOmega}}{\boldsymbol{X}}={\boldsymbol{\varOmega}}^{*}$.

###  4.4 Keldysh rotation

For further convenience it is useful to perform the change of the basis in the
doubled field space $\phi_{+}$, $\phi_{-}$ and introduce the so-called
classical and quantum fields $\phi_{c}$ and $\phi_{q}$ [45, 46],

${\boldsymbol{\phi}}_{K}(t)=\begin{bmatrix}\;\phi_{c}(t)\,\\\
\phi_{q}(t)\end{bmatrix}=\boldsymbol{C}{\boldsymbol{\phi}}(t),\;\boldsymbol{C}\equiv\begin{bmatrix}[c]\;\frac{1}{2}I&&\,\frac{1}{2}I\;\\\
I&&-I\end{bmatrix}.$ (4.39)

This transformation is called Keldysh rotation. In the new basis, the Green’s
function ${\boldsymbol{G}}$ takes the form

${\boldsymbol{G}}_{K}(t,t^{\prime})\\!=\\!\boldsymbol{C}{\boldsymbol{G}}(t,t^{\prime})\,\boldsymbol{C}^{T}\\!\\!=\\!\begin{bmatrix}\,G_{K}(t,t^{\prime})&G_{R}(t,t^{\prime})\,\\\
\,G_{A}(t,t^{\prime})&0\end{bmatrix}\\!.$ (4.40)

Here $G_{R}$ and $G_{A}$ are the retarded and advanced Green’s functions,
respectively, having the following operator form

$\displaystyle G_{R}(t,t^{\prime})$
$\displaystyle=-i\operatorname{tr}\Bigl{(}\hat{\rho}\,\bigl{[}\hat{\phi}(t),\hat{\phi}(t^{\prime})\bigr{]}\Bigr{)}\,\theta(t-t^{\prime})$
$\displaystyle=-i\bigl{[}v(t)v^{\dagger}(t^{\prime})-v^{*}(t)v^{T}(t^{\prime})\bigr{]}\,\theta(t-t^{\prime}),$
(4.41) $\displaystyle G_{A}(t,t^{\prime})$
$\displaystyle=G_{R}^{T}(t^{\prime},t).$ (4.42)

They are consistent with the classical definition (3.38), in particular,
because of independence of the commutator average of the state $\hat{\rho}$.
The block $G_{K}$ is called Keldysh Green’s function and contains the
information about the state. In view of operator averages (4.3) it expresses
as the mean value of the anti-commutator of fields and, due to (4.40),
explicitly reads in terms basis functions as

$\displaystyle iG_{K}$
$\displaystyle(t,t^{\prime})=\tfrac{1}{2}\operatorname{tr}\Bigl{(}\hat{\rho}\,\bigl{\\{}\hat{\phi}(t),\hat{\phi}(t^{\prime})\bigr{\\}}\Bigr{)}$
$\displaystyle=\begin{bmatrix}\;v(t)&&v^{*}(t)\,\end{bmatrix}\bigl{(}{\boldsymbol{\nu}}+\tfrac{1}{2}{\boldsymbol{X}}\bigr{)}\begin{bmatrix}\;v^{T}(t^{\prime})\,\\\
v^{\dagger}(t^{\prime})\end{bmatrix}.$ (4.43)

###  4.5 Special choice of basis functions and particle interpretation

Thus far, the matrix $\omega$, which defines the Neumann boundary conditions
for the basis functions $v$, $v^{*}$, is not fixed except the requirement of
symmetry under transposition and positive definiteness of its real part. In
this section we make a convenient choice of $\omega$ which leads to the
expressions for the Green’s functions admitting particle interpretation with
well-defined notion of average occupation number.

For this purpose, it is useful to rewrite the Keldysh Green’s function in
terms of non-anomalous and anomalous particle averages

$\displaystyle\nu=\operatorname{tr}\bigl{[}\hat{\rho}\,\hat{a}^{\dagger}\hat{a}\bigr{]},\qquad\kappa=\operatorname{tr}\bigl{[}\hat{\rho}\,\hat{a}\,\hat{a}\bigr{]},$
(4.44)

so that from (4.4) $G_{K}$ becomes

$\displaystyle i\,G_{K}(t,t^{\prime})$
$\displaystyle=\begin{bmatrix}\;v(t)&&v^{*}(t)\;\end{bmatrix}$
$\displaystyle\times\begin{bmatrix}\;\kappa&\nu^{*}+\frac{1}{2}I\;\\\
\;\nu+\frac{1}{2}I&\kappa^{*}\;\end{bmatrix}\begin{bmatrix}\;v^{T}(t^{\prime})\;\\\
v^{\dagger}(t^{\prime})\end{bmatrix}.$ (4.45)

Note that the matrix $\kappa$ is symmetric, whereas $\nu$ is Hermitian.
Comparing with (4.4) we find the connection between particle averages and the
matrix ${\boldsymbol{\nu}}$

${\boldsymbol{\nu}}=\begin{bmatrix}\;\kappa&&\nu^{*}\\\
\;\nu&&\kappa^{*}\end{bmatrix}.$ (4.46)

Thus, we see that block-diagonal components of ${\boldsymbol{\nu}}$ are
responsible for anomalous averages. To ascribe the particle interpretation to
the creation/annihilation operators, we will try to choose the matrix
$\omega$, defining the corresponding basis functions $v(t)$ and $v^{*}(t)$ so
that the diagonal blocks of ${\boldsymbol{\nu}}$, defining the anomalous
averages $\kappa$, vanish. Moreover, this choice will simplify the expressions
for the Green’s functions, since they contain the terms, containing $\kappa$.
For example, with a nonzero $\kappa$ the Wightman function reads

$G_{>}(t,t^{\prime})=v(t)\,(\nu^{*}+I)\,v^{\dagger}(t^{\prime})+v^{*}(t)\,\nu\,v^{T}(t^{\prime})\\\
+v(t)\,\kappa\,v^{T}(t^{\prime})+v^{*}(t)\,\kappa^{*}\,v^{\dagger}(t^{\prime}).$
(4.47)

To make the matrix ${\boldsymbol{\nu}}$ block off-diagonal consider the
expression (4.32) and note that the only block diagonal contribution is
contained in the identity matrix ${\boldsymbol{I}}$ and, possibly, in the term
involving $({\boldsymbol{\omega}}+{\boldsymbol{\varOmega}})^{-1}$. Thus, we
want to choose $\omega$ such that block-diagonal contribution of the latter
exactly cancels those of ${\boldsymbol{I}}$. Using the block matrix inversion
formula666The useful form of the block matrix inversion formula is
$\displaystyle\begin{bmatrix}\;A\;&B\;\\\ \;C\;&D\;\end{bmatrix}^{-1}=$
$\displaystyle\begin{bmatrix}\;(A-BD^{-1}C)^{-1}&0\\\
0&(D-CA^{-1}B)^{-1}\end{bmatrix}$
$\displaystyle\times\begin{bmatrix}I&-BD^{-1}\\\ -CA^{-1}&I\end{bmatrix}$ with
$A=R+\omega$, $B=S$, $C=S^{*}$, $D=R^{*}+\omega^{*}$. we have the condition of
the vanishing block-diagonal part of ${\boldsymbol{\nu}}$,

$R+\omega-S(R^{*}+\omega^{*})^{-1}S^{*}-2\omega_{\mathrm{re}}=0.$ (4.48)

We will focus on the case in which $R$ and $S$ are real. The formalism
described below can be easily extended to the complex $R$, but it seems that
there is no straightforward extension to general complex (Hermitian) $S$.
Introducing the dimensionless quantities

$r=\omega^{-1/2}\,R\,\omega^{-1/2},\qquad s=\omega^{-1/2}\,S\,\omega^{-1/2},$
(4.49)

the equation (4.48) can be rewritten as $r+I-s(r+I)^{-1}s=2$ and further
simplified by introducing the new variable
$\tilde{s}=(r+I)^{-1/2}s\,(r+I)^{-1/2}$ and solving for $\tilde{s}$, so that
it takes the following form

$r^{2}=s^{2}+I.$ (4.50)

This is implicit equation on $\omega$, due to the above definition of $r$ and
$s$. Its explicit form reads

$R\omega^{-1}R=S\omega^{-1}S+\omega,$ (4.51)

which can be solved in the form advocated in Section 2

$\omega=R^{1/2}\sqrt{I-\sigma^{2}}R^{1/2},\quad\sigma\equiv
R^{-1/2}SR^{-1/2}.$ (4.52)

Note that the assumption of positive definiteness of $\omega$ implies that
$I-\sigma^{2}=(I-\sigma)(I+\sigma)$ is positive definite. Recalling that
$R+S=R^{1/2}(I+\sigma)R^{1/2}$ should be positive definite for normalizability
of the density matrix, it is easy to see that
$I-\sigma=R^{-1/2}(R-S)R^{-1/2}$, or equivalently $R-S$ should be positive
definite too. Then, the substitution of the obtained expression for $\omega$
to (4.32) gives the desired block-diagonal matrix form of ${\boldsymbol{\nu}}$
advocated in Section 2

$\displaystyle{\boldsymbol{\nu}}=\begin{bmatrix}\;0&&\nu\;\\\
\;\nu&&0\;\end{bmatrix},\quad\nu\equiv\frac{1}{2}\varkappa\left(\sqrt{\frac{I-\sigma}{I+\sigma}}-I\right)\varkappa^{T},$
(4.53)
$\displaystyle\varkappa\equiv\big{[}\omega^{1/2}R^{-1}\omega^{1/2}\big{]}^{1/2}\omega^{-1/2}R^{1/2}=\bigl{(}\varkappa^{T}\bigr{)}^{-1}$
(4.54)

where the matrix $\varkappa$ introduced above is orthogonal. Therefore, as a
consequence of positive definiteness of $I+\sigma$ and $I-\sigma$, the matrix
$\nu$ is necessarily real. As shown in Appendix D, for the density matrix to
be positive definite, the matrix $\sigma$ should be negative definite, so
$\nu$ is positive definite.

Substituting it to (4.33), one immediately obtains simple expressions for the
Green’s functions. In particular, for Wightmann and Feynman functions one has

$\displaystyle iG_{\mathrm{T}}(t,t^{\prime})=$ $\displaystyle
v(t)\,v^{\dagger}(t^{\prime})\,\theta(t-t^{\prime})+v^{*}(t)\,v^{T}(t^{\prime})\,\theta(t^{\prime}-t)$
$\displaystyle+\,v(t)\,\nu\,v^{\dagger}(t^{\prime})+v^{*}(t)\,\nu\,v^{T}(t^{\prime}),$
(4.55) $\displaystyle iG_{>}(t,t^{\prime})=$ $\displaystyle
v(t)\,\bigl{(}\nu+I\bigr{)}\,v^{\dagger}(t^{\prime})+v^{*}(t)\,\nu\,v^{T}(t^{\prime}),$
(4.56)

while the others can be expressed through them in a straightforward way.

It will be useful to express ${\boldsymbol{\varOmega}}$ in terms of $\nu$.
Disentangling ${\boldsymbol{\varOmega}}$ from (4.32), and then using the
explicit form (4.53) of ${\boldsymbol{\nu}}$, corresponding to the special
choice of Neumann basis functions with (4.52), we obtain the following
expression

${\boldsymbol{\varOmega}}=\omega^{1/2}\begin{bmatrix}\;\frac{2\nu^{2}+2\nu+I}{2\nu+I}&\;-\frac{2\nu(\nu+I)}{2\nu+I}&\\\
\;-\frac{2\nu(\nu+I)}{2\nu+I}&\;\frac{2\nu^{2}+2\nu+I}{2\nu+I}\end{bmatrix}\omega^{1/2}.$
(4.57)

###  4.6 Euclidean density matrix state

Now, let us focus on the particular Gaussian state, which is obtained from the
Euclidean path integral, namely

$\displaystyle\rho_{E}(\varphi_{+},\varphi_{-};J_{E}]=\frac{1}{Z}\\!\\!\\!$
$\displaystyle\int\limits_{\phi(\tau_{\pm})=\varphi_{\pm}}D\phi\,\exp\biggl{\\{}-S_{E}[\phi]$
$\displaystyle-\int_{0}^{\beta}d\tau\,J_{E}(\tau)\phi(\tau)\biggr{\\}},$
(4.58)

Here $S_{E}$ is the quadratic action of Euclidean field theory within time