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

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.

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 nonequilibrium system with the initial state in the form of arXiv:2309.03687v1[hep-th] 7 Sep 2023 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].
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 Ψ [φ(x)] of a generic free fields model in flat spacetime can be built by the path integral over the field histories ϕ(τ, x) on a half-space interpolating between a given 3-dimensional configuration φ(x) on the boundary plane of τ = 0 and the vanishing value of these fields at the Euclidean time τ → −∞.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 |ψ⟩ along with the set of their weights w ψ participating in the construction of the density matrix ρ = ψ w ψ |ψ⟩⟨ψ|.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 |ψ⟩.These are the Wheeler-DeWitt equations Ĥµ |ψ⟩ = 0 which are quantum Dirac constraints in gravity theory selecting the physical states [8], µ 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 the factor Z being a partition function which provides the normalization tr ρ = 1 [9].Important feature of this formal projector is that a detailed construction of the delta function of noncommuting operators Ĥµ (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 [9,10] 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 [11].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 τ ).Then the two-point kernel ρ E (φ + , φ − ) = ⟨φ + | ρE |φ − ⟩ of the density matrix in question is defined by where S E [ϕ] is the Euclidean action of the field perturbations ϕ(τ ) on top of the given background, defined on the period of the Euclidean time, τ − ≤ τ ≤ τ + , the functional integration runs over field histories interpolating between their values φ ± -the arguments of the density matrix kernel.Z is the partition function given by the path integral over the periodic histories with the period providing the normalization tr ρE = 1.Hermiticity of this density matrix, which in view of its reality reduces to its symmetry ρ E (φ + , φ − ) = ρ E (φ − , φ + ), implies that the background solution is a bounce that has a reflection symmetry with respect to the middle turning point at τ++τ− 2 , and the turning points τ ± 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τ , with the quantum Hamiltonian ĤS (τ ) 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 τ -ordered exponent Because of the properties of the turning points (zero derivatives of the background field) the Euclidean background can be smoothly matched at τ ± 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) dt Ĥ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 Û (T, 0)ρ E Û † (T, 0) is depicted on Fig. 1.It shows the Euclidean spacetime instanton with the topology R 1 × S 3 , R 1 = [τ − , τ + ], bounded at the turning points τ ± by two minimal surfaces Σ ± with a vanishing extrinsic curvature.This instanton represents the density matrix ρE and connects  the Lorenzian spacetime branches.These branches correspond to the unitary and anti-unitary evolution from Σ ± in some finite interval of the Lorentzian time 0 ≤ t ≤ T . 1 The pictorial representation of the cosmological partition function Z in view of cancellation of unitary evolution factors, tr Û (T, 0)ρ E Û † (T, 0) = tr ρE = 1, contains only the Euclidean part of Fig. 1.It is represented by the closed cosmological instanton with the identified surfaces Σ + = Σ − and their 3-dimensional field configurations φ + = φ − (following from the identification of the arguments in tr ρE = dφ ρ E (φ, φ)).The origin of this instanton having a donut topology S 1 × 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 Ψ (φ + ) and its complex conjugated counterpart 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 1 Of course, the second Lorentzian branch could have been attached to the middle turning point 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.based on the instanton with S 4 -topology of Fig. 4. The latter originates by glueing together two 4-dimensional hemispheres (discs D 4 ± ) along their common equatorial boundary.
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 [12][13][14].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.
The notion of the Euclidean quantum gravity density matrix was pioneered in [15].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 [16].Interpreted as a microcanonical density matrix of spatially closed cosmology [9] 2 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 [16], provides the origin of the Higgs-type or R 2 -type inflationary scenario [17] with subplanckian Hubble scale [18] 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 +∞ [16,19].On the other hand, this is the hill-top nature of inflation starting from the maximum of inflaton potential rather than from its minimum [20].Finally, this model suggests that quantum origin of the Universe is the subplanckian phenomenon subject to semiclassical 1/Nperturbation theory in the number of numerous higherspin conformal fields [21].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 [22] 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 [23].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 [24][25][26] and the pioneering calculation of non-gaussianities in cosmological perturbation spectra [27] these works include the calculation of cosmological correlation functions [28,29], the results on cosmological singularity avoidance due to nonlocal effects [30], equivalence of the Euclidean and in-in formalisms in de Sitter QFT [31,32] and even the analysis of initial conditions within Schwinger-Keldysh formalism [33].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 [34], its applications to four-point correlators in inflationary cosmology [35] and numerous conformal field theory and holography ramifications of Schwinger-Keldysh technique (see, for example [34,36] and references therein).However, the success of these works essentially relies on working with the model of spatially flat Universe -extension to spathe 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.
tially 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 [37], testifies at more that 99% confidence level in favor of the closed Universe preferring a positive curvature with K = +1 [38,39].Remarkably, the model of microcanonical initial conditions in early quantum cosmology of [9,16] 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 inin 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 Its evolution in time gives rise to Matsubara technique of thermal Green's functions [40] and thermofield dynamics [41] which satisfies nontrivial analyticity properties in the complex plane of time including periodicity in the direction of imaginary axis -Kubo-Martin-Schwinger (KMS) condition [42,43].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 inin formalism.Therefore we restrict ourselves with the systems having a quadratic action on top of the nonstationary 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.

Schwinger-Keldysh technique for models with density matrix state
We consider a generic system with the action S[ϕ] quadratic in dynamical variables ϕ = ϕ I (t), the index I including both the discrete tensor labels and in fieldtheoretical context also the spatial coordinates, Here A = A T ≡ A IJ , B ≡ B IJ and C = C T ≡ C IJ are the matrices acting in the vector space of ϕ J , the superscript T denoting the transposition, ϕ being a column and ϕ 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 nonstationary physical setup.This action will be considered as a quadratic part of the full nonlinear action in field perturbations ϕ 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 ρ.This is the functional of two sources Here the trace is taken over the Hilbert space of the canonically quantized field φ and Û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)ϕ(t) ≡ −J I (t)ϕ I (t) with the source J T (t) = J I (t).In the (2.3) We will consider the class of density matrices whose kernel in the coordinate representation ⟨φ + | ρ |φ − ⟩ = ρ(φ + , φ − ) has the following Gaussian form -exponentiated quadratic and linear forms in φ ± , where we assembled φ ± into the two-component column multiplets (denoted by boldfaced letters) φ, did the same with the coefficients j of the linear form and introduced the 2 × 2 block-matrix Ω acting in the space of such two-component multiplets The blocks of this matrix R = R IJ , S = S IJ and their complex and Hermitian conjugated versions, S † ≡ S T * , should satisfy these transposition and conjugation properties in order to provide Hermiticity of the density matrix.The same concerns the "sources" j ± in the definition of j, j + = j * − ≡ j.Transposition operation above applies also to two-component objects, so that φ T = [φ T + φ T − ].Below we will denote 2 × 2 block matrices and relevant 2block component columns and rows by boldfaced letters.
Such a choice of the density matrix is motivated by the fact that for a block-diagonal Ω it reduces to a pure quasi-vacuum state, its "source" j allows one to induce nonzero mean value of the quantum field and by the differentiation with respect to 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 Ω 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 ϕ 1 (t) and ϕ 2 (t) as well as for their sources, In terms of these notations the generating functional reads where the total action is obviously with the actions S[ϕ 1,2 ] given by (2.1) in the integration range from t = 0 to t = T and the total integration measure over ϕ and φ (2.10) Here dφ and Dϕ denote respectively the integration measures over variables at a given moment of time and the integration measures Dϕ = t dϕ(t) over time histories subject to indicated boundary conditions.Calculation of this Gaussian path integral leads to the expression 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 pointthe point of stationarity of the action with respect to variations of both the histories ϕ(t) and their boundary data φ 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 ϕ(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 G(t, t ′ ) subject to homogeneous version of these boundary conditions.The Green's function has a block-matrix form typical of Schwinger-Keldysh inin formalism composed of the Feynman G T (t, t ′ ), anti-Feynman GT(t, t ′ ) and off-diagonal Wightman Green's functions blocks (4.3), which are related to one another by the equalities ) and satisfy respectively inhomogeneous and homogeneous wave equations with the operator F -the Hessian of the action (2.1),

14)
The block-matrix Green's function G(t, t ′ ), as is usually done in boundary value problems, can be built in terms of the full set of basis functions v ± 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 ϕ 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 v ± (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, where W is what we call the Wronskian operator which participates in the Wronskian relation for the operator F , which is valid for arbitrary two complex fields ϕ 1,2 (t), (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 ϕ 1,2 = 0, (2.18) We will call the boundary conditions (2.15) and associated with them Green's functions the Neumann ones 3 .
Important point of the definition (2.15) is that the frequency matrix ω (remember that in the generic setup this is a matrix ω IJ acting in the vector space of ϕ 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 â † and â in terms of which the Heisenberg operator φ(t) is decomposed as the sum of positive-negative basis functions v(t) and v * (t), φ(t) = v(t) â + v * (t) â † .Then there exist non-anomalous and anomalous particle averages with respect to the density matrix, and the requirement of vanishing anomalous average κ = 0 allows one to assign the average ν the interpretation of the set of occupation numbers associated with ρ.This requirement serves as the equation for the frequency matrix ω which, as it is shown in Section 4, can be explicitly solved for a special case of the real matrix Ω.This solution reads and gives the expression for the occupation number matrix in terms of the single symmetric matrix σ after the orthogonal rotation by the orthogonal matrix κ, As shown in Appendix D, the existence of this particle interpretation with a positive definite matrix ν fully matches with conditions of normalizability, boundedness and positivity of the density matrix incorporating positive definiteness of matrices I ± σ and negative definiteness of σ.
With the normalization of these distinguished basis functions v to unity 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 Here the terms of the type v(t) v † (t ′ ) should be understood as the matrix products ) (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 ϕ I components, whereas the subscript indices A enumerate the basis functions in their full linear independent set).Correspondingly, , 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 = p, A = d 3 p/(2π) 3/2 and ν AB = ν p,p ′ = (exp(βω p ) − 1) −1 δ(p − p ′ ) represents expected occupation number for Bose-Einstein statistics at inverse temperature β (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 κ), and its eigenvalues would serve as occupation numbers in the generic nonequilibrium state with the initial density matrix.

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 where integration runs over histories ϕ(τ ) in Euclidean time τ on the segment [τ − , τ + ], interpolating between the arguments of the density matrix φ ± .In what follows we will assume that τ − = 0 and τ + = β.The Euclidean action, supplied by the Euclidean source J E (τ ) 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 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, (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, β].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 τ++τ− 2 = β 2 , and the other one corresponds to the (identified) points of nucleation τ ± 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 β = τ + − τ − 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 ( The resulting density matrix looks like the expression (2.4) amended by the quadratic form in the Euclidean source, with the special expressions for the matrix Ω E and the source j E .The matrix Ω E reads 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 ϕ T (τ ) , etc.The column j E is given by the following integral (2.32) Using (2.11) in the generating functional with the Euclidean density matrix (2.30) 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 inin formalism (2.34)

Reflection symmetry and analyticity properties
The obtained expression for Z[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 where the Lorentzian and Euclidean histories are also related by the same continuation rule ϕ(t)| t=−iτ = ϕ E (τ ).This, in particular, implies that the coefficients of the operators F E and F are related by so that F | t=−iτ = −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 τ 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 τ = 0 where all real field variables have zero τ -derivatives and can be smoothly continued to the imaginary axis of τ = it where they become again real functions of real t.This leads to the above continuation rule with real With this analytic continuation rule, the expression (2.34) for Z[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 nega-tive frequency basis functions.This expression reads where the z-integration runs respectively over t or τ in the domain depending on which of these Lorentzian or Euclidean time variables is in the argument of the following block matrix Green's function G(z, z ′ ) and the corresponding source (2.38) Here the Euclidean and Lorentzian-Euclidean blocks of the total Green's function express in terms of the relevant Euclidean and Lorentzian-Euclidean Wightman functions 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 ± (τ ) basis functions.All these functions are subject to Neumann boundary conditions (2.15) and For ω fixed by the above condition of particle interpretation, leading to the expressions (2.20)-(2.22), the Euclidean basis functions u ± have a remarkable property.They satisfy at opposite ends of the Euclidean segment τ ∓ the same boundary conditions4 If one smoothly continues the operator F E beyond the segment τ − ≤ τ ≤ τ + , then it becomes periodic with the period β (which is possible because τ ± 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 ± of this operator become quasi-periodic -u ± (τ + β) expresses as a linear combination of the same basis functions u ± (τ ) (no mixing between u − and u + occurs in their monodromy matrix).As shown in Section 4, with the normalization u ± (0) = 1/ √ 2ω this quasi-periodicity property reads in terms of the occupation number matrix (2.21) (2.47) 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 τ = 0 (2.48) Therefore, similarly to quasi-periodicity the basis functions u ± (±τ ) are also related by the analogue of the anti-diagonal monodromy matrix L, u The above relations introduce the analytic structure which allows one to express all basis and Green's functions on the Euclidean-Lorentzian domain C in terms of one analytic function V (z) of the complexified time variable z = t − iτ .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 ≡ F (t, d/dt) = −F E | τ =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 + (τ ) also turns out to be the analytic continuation of the latter, (2.50) Therefore, the operators F and −F E as well as the full set of their basis functions v(t) and u ± (τ ) can be represented respectively as the boundary values at the real and imaginary axes of the complex z-plane of the complex operator F C and the solution V (z) of its homogeneous wave equation, ) The function V (z) gives rise to basis functions as 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 + (τ ) from the imaginary axes z = −iτ , and in view of reality of u + (τ ) for real τ it has the property 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

55)
Due to inverse matrix factors of positive and negative basis functions here the Lorentzian Wightman functions G > (t, t ′ ) (given by the expression (2.25)) and which is nothing but Kubo-Martin-Schwinger condition [42,43].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.The tubular Riemann surface of complex time z = t − iτ whose main sheet is compactified in τ to the circle of curcumference β 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 [31,32] might be extended to a generic reflection-symmetry background underlying our definition of the Euclidean density matrix.

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.

Condensed notations
The elements of the field space will be denoted as ϕ I (t), where the index I is, in fact, a multi-index, and contains both the dependence on spatial coordinates denoted as x and discrete spin-tensor labels i, I = (x, i).Thus, we can equivalently write the fields in the form, emphasizing its dependence on the spatial coordinates ϕ I (t) = ϕ i (t, 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 ± , 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 ) where A ij (t, x, x ′ ) is the kernel of the operator.Thus, the first term in (3.1) has the following explicit structure 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 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 ϕ have the form 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 allows one to rewrite the quadratic action (3.1) in even more condensed form Here the Wronskian operator W is defined by (2.16) and the origin of the boundary term at t ± is the result of integration by parts, which is also associated with the Wronskian relation (2.17).

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 π canonically conjugated to the field ϕ where L is the Lagrangian of the action (3.1).The corresponding Hamiltonian has the form (3.9) Together with the Poisson bracket {ϕ I , π J } = δ I J it defines the dynamics of the system.The Hamiltonian equations of motion read Transition to the Lagrangian formalism by expressing π in terms of ϕ and φ 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 ± 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 and can be rewritten in shortened notations as Here α ±A constitute a set of constants, specifying particular initial conditions.Using (3.8), we find the corresponding solution for the momentum so, the evolution of phase space variables can be rewritten in the joint form as (3.14) Now, we can equip the space of initial conditions, consisting of α ± , with the Poisson bracket structure inherited from the Poisson brackets of ϕ and π.Substituting (3.14) into the left hand side of we have in condensed notations where I denotes the identity matrix.The identity above fixes the pairwise Poisson brackets of α ± .Let us denote the right hand side of this equality, playing the role of the Poisson bivector in the Darboux coordinates, as Introducing also the matrix D as inverse to the matrix of the pairwise Poisson brackets Thus, one can express the inverse of M(t) in terms of its transpose, namely (3.21) Thus, for ϕ 1,2 -solutions of (3.5) l.h.s.vanishes, so we have It is easy to see, that each element of (3.18) has the form (3.22) as above, where the role of solutions ϕ 1 , ϕ 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 D is skew-symmetric, since P T = −P.In terms of the block elements of D this means that 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 D is real skew-symmetric, so there is a timeindependent linear transformation S, bringing it to the canonical form, i.e. S T DS = P. Without the loss of generality one set D = P by default 5 .However, for the reasons which will become clear soon (see equation (4.21) below), we will assume that D has the following more general form where In terms of the basis functions, the vanishing of the diagonal blocks of D implies that v + , v − are chosen such that This can always be done by an appropriate transformation of the basis functions, possibly mixing v + and v − .Consequently, the pairwise Poisson brackets of α + and α − take the form ) As noted above, one can go further, and set ∆ +− = −∆ T −+ = I.Now, let us modify the Hamiltonian by introducing time-dependent sources J ϕ , J π for the field and its conjugate momentum The modified equations of motion can be written as where the subscript J of ϕ, π 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 α + , α − in its definition time-dependent Then, we substitute the result to the modified e.o.m. and obtain where we exploit the fact that M(t) satisfies the system (3.10).Using the equality (3.20) for the inverse of the matrix M(t) and integrating the equation on α + (t) and α − (t) we obtain ) where α + 0 and α − 0 are integration constants.Substitution back to (3.30) gives the solution to the equations (3.29) where the initial conditions ϕ 0 (t), π 0 (t) are related to constants of integration by and represent the solution to homogeneous equation, i.e. for vanishing sources J ϕ and J π .Now, let us focus on the case of vanishing momentum source and also redefine the field source for the convenience The corresponding e.o.m. in the Lagrange form reads From (3.33), one obtains the explicit form of the solution for ϕ(t), which is where G R (t, t ′ ) is called the retarded Green's function and expressed through the top-left block of the matrix (3.38)The fact that ∆ ++ = ∆ −− = 0 is crucial in obtaining this simple expression for G R .From (3.37) we find that G R satisfies the equation and is uniquely determined by the condition 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 The fact that solution (3.37) is expressed through the retarded Green's function means that ϕ(t) is subject to the following initial (rather than boundary) value problem (3.42)

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 It is symmetric under the simultaneous transposition and exchange of the time arguments, i.e.G T S (t, t ′ ) = G S (t ′ , t).Unlike the retarded Green's function it is defined nonuniquely 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 one can use the Dirichlet Green's function subject to the same boudary conditions so that the solution reads Similarly, in solving Neumann boundary problem one defines the corresponding Neumann Green's function demanding ) and obtains the solution as (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 ϕ 2 → ϕ(t ′ ), ϕ 1 → G(t ′ , t), where ϕ(t ′ ) solves (3.36) and G(t ′ , t) is some Green's function, solving F G(t ′ , t) = δ(t − t ′ ).Next, integrating both sides of the equality over t ′ from t − to t + , we obtain Now, suppose we are to solve (3.36) supplemented by inhomogeneous boundary conditions (in contrast to homogeneous ones (3.44)) for some constants φ + , φ − .Substituting these conditions to (3.50) together with Dirichlet Green's function G → G D , satisfying (3.45), we observe that the third line vanishes, so we get where we introduce the notation for the two-component row as the transposition of the newly introduced column W denotes the Wronskian operator (2.16) acting from the right on the second argument of G D (t, t ′ ) at the total boundary of the time domain at t ± (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, The quantity w(t) introduced above has the following important property.Namely, evaluating both sides of (3.52) at the boundary points t = t ± , and using (3.45) we observe that Similarly, one can consider inhomogeneous Neumann boundary conditions with some boundary sources j + and j − .Substitution of this condition and Neumann Green's function G → G N , satisfying (3.48), to (3.50) gives the solution to (3.36) with the boundary conditions above Here g T N (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 + ), (3.59)

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 φ ± .Substituting the solution of this problem (3.52) into (3.47)one obtains a linear equation on φ ± , which can be solved as where the matrices ω and Ω read ) Substituting these φ ± back into (3.52)gives which implies, after comparing with (3.49), the following expression for the Neumann Green's functrion Here we use the notations (3.53)-(3.54)introduced above.Substituting t = t ± to the both sides of the equality and using (3.56), we get the equality that allows us to express the Dirichlet Green's function from (3.64) via the Neumann one as ), (3.66)where we use the notation (3.59) for the row 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 g N at both ends of the time segment (double bar denoting the restriction of both arguments to t ± ) (3.67)

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 ± (t) real, and such that the matrix D defined by (3.24) has a canonical form, D = P. Together with the reality of ϕ(t) this implies also the reality of the corresponding integration constants α ± .Next, we combine these basis functions and integration constants into the following complex conjugated pairs The equation (3.20) takes the form 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, 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. ϕ(t), π(t) → φ(t), π(t), whereas the Poisson bracket is replaced by the commutator times the factor i, so that we obtain the equal-time quantum commutation relations [ φ(t), π(t)] = i Î, where Î 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 Here Ĥ(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 ϕ and its conjugate momentum π is literally the same at classical and quantum levels π(t) = W φ(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) Similarly, our quantization procedure implies that the integration constants α, α * are in one-to-one correspondence to the creation/annihilation operators â, â † .According to (3.72) the operators φ(t), π(t) are decomposed in the creation/annihilation operators as that can be inverted similar to (3.74) as The fact that â and â † are indeed Hermitian conjugate to each other immediately follows from the Hermicity of φ(t).Indeed, comparing φ(t) to it's conjugate 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 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, dt ĤS (t) , (3.82) where ĤS (t) is the Hamiltonian in the Schroedinger picture, so that its time dependence is only due to timedependent coefficients A, B, and C. The operators φ, π in the Schroedinger picture are identified with the Heisenberg ones, evaluated at the initial time φ = φ(t − ), π = π(t − ). (3.83) In the presence of the source, H → H −J T ϕ, the solution (3.78) to the Heisenberg equation generalizes to Since the action (3.1) is quadratic in the field φ, 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.

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 b, b † as a linear combination of the initial ones, 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 Thus, the field operator φ(t) has two equivalent decompositions where the new set of the basis functions ṽ(t), ṽ * (t) is related to the initial one via the following relation or, in more explicit form v = ṽ U + ṽ * V * .
(3.94) Equality (3.90) leads to the following formula for the inverse matrix of the Bogoliubov transformation coefficients so that (3.94) can be inverted as Now, let us solve the inverse problem.Namely, suppose we have two sets of the basis functions v(t), v * (t) and ṽ(t), ṽ * (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) This is constant if ϕ 1 , ϕ 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 and the same for ṽ(t), ṽ * (t).Projecting the equality (3.94) onto ṽ, and using the property , one obtains the explicit expressions for the Bogoliubov coefficients U = (ṽ, v), V = (ṽ, v * ). (3.99) If both the old and the new sets of basis functions satisfy the Neumann conditions with different frequency matrices ω and ω at the initial moment of time one can find the explicit expressions for U , V in terms of ω and ω.Let us first write down the normalization conditions (3.98) explicitly ) where all quantities are evaluated at t = t − .The same equations hold for ṽ(t) and ṽ * (t).The second equation implies that the matrices ω and ω 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 where ω re and ωre are the real parts of ω and ω, 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, where ω is given in terms of the positive frequency basis function whose integration gives the (unnormalized) solution For this normalization we have the following expression for the Fock states in terms of coherent state (3.114)

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.

Gaussian states
Our goal is to find the explicit and useful form of the generating functional of in-in correlation functions where ÛJ are the evolution operators subject to equation (3.85) with different sources J 1 and −J 2 , whereas T and 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 where G T , GT, and G < are Feynman, anti-Feynman and Wightman Green's functions, respectively.The density matrix ρ 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.88) of the evolution operator, one obtains the following expression for the generating functional where the integration over ϕ 1,2 (t) runs with the indicated boundary conditions and we introduce the notation for the coordinate representation of the density matrix Now, we restrict ourselves with the Gaussian density matrices, i.e. those whose coordinate representation has the form of the Gaussian exponent where the matrix Ω, and the vector j play the role of the parameters of ρ, and normalization constant 1/Z is independent of φ.The Hermitian property of the density matrix, ⟨φ implies the following conditions on Ω and j or, in a more explicit block-matrix form Normalizability of ρ 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 ρ2 ̸ = ρ.The role of the linear term in the exponential in (4.5) is two-fold.Firstly, j defines nonvanishing 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 j.The typical example of the (pure) Gaussian state is the vacuum state (3.104), i.e. ρ = |0⟩⟨0|, associated with some choice of the annihilation operator, for which R = ω * , S = 0, and j = 0. Another example of the pure Gaussian state is the coherent state (3.110) whose density matrix reads ρ = |α⟩⟨α|, and R = ω * , S = 0 again, but j = √ 2ω re α * .

In-in boundary value problem
Substituting the general Gaussian density matrix to (4.4), one obtains where we put the boundary points φ ± , φ ′ , appearing in (4.4) 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, so that the sum of the actions for ϕ 1 and ϕ 2 in (4.9) can be rewritten in the joint form 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 Independent variation of the fields δϕ(t) in the interior of the time interval gives equations of motion whereas the variation of the boundary values δϕ(T ) and δϕ(0) = δφ supply these equations with the boundary conditions.They read as the following matrix relations where we took into account that in view of ϕ 1 (T ) = ϕ 2 (T ) the variation δϕ T (T ) = δϕ T 1 (T ) I I , and the boundary conditions at t = T reduce to the equality of both the fields and their time derivatives of both ϕ 1 and ϕ 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 j We can construct the Green's function G solving the problem above, out of the basis functions v ± .These basis function should solve the homogeneous equation, and satisfy the same boundary conditions as those of the Green's function, Applying the generic Green's function expression (3.43) to the case of the doubled field content, we obtain the Green's function G in terms of these basis functions

Neumann type basis functions and Green's function representation
However, we do not have the explicit form of basis functions v ± .We will construct v ± with the help of another set basis functions v, v * subject to much simpler boundary conditions Since W and ω are block-diagonal, the basis functions v, v * can be chosen block-diagonal too, namely With a real operator F the blocks of these matrices solve the equations F v(t) = 0 and F v * (t) = 0, subject to the complex conjugated boundary conditions Thus, v and v * are simply the basis functions for the single field ϕ + or ϕ − subject to the Neumann boundary conditions introduced above.We assume that ω is the symmetric matrix with a positive-definite real part.
The answer for the basis function v + in terms of v and v * can be easily constructed as while the calculation of v − requires more efforts.We will obtain the answer for v − with the use of Bogoliubov coefficients relating two sets of different Neumann basis functions (3.103) by treating v − as the negative frequency basis function complex conjugated to its positive frequency counterpart v * − satisfying at t = 0 the boundary condition (iW − Ω * ) v * − | t=0 = 0. Thus, in accordance with (3.103) the answer for v − reads where U , V are the corresponding Bogoliubov coefficients Here we assume the normalization v(0 √ 2Ω re , and denote the real parts of ω and Ω respectively as ω re and Ω re . Finally, let us consider the details of the Green's function G(t, t ′ ) defined by (4.21) for a particular form of v ± we have just built.Its matrix ∆ −+ given by (4.22) reads 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 ) where we introduce the following symmetric matrix 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 where G 0 is defined as and interpreted as the Green function, corresponding to the vacuum state, having the density matrix ρ0 = |0⟩⟨0|, associated with the basis functions v(t), v * (t).Indeed, from (3.109) one observes that the matrix Ω, defining the vacuum density matrix ρ0 , coincides with ω * , i.e.
Ω = ω * .In this case, ν vanishes due to its definition (4.32), so from (4.33) we find that G = G 0 .Substituting the generating functional obtained to (4.2), we observe that for vanishing j the block-matrix components of G are composed of the Feynman, anti-Feynman and Wightman Green's functions (4.3), namely where G > (t, t ′ ) ≡ G T < (t ′ , t), and the explicit form of the block components can be read off form (4.33).Now we have to find the solution ϕ(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 doublefield version of (3.58) with the substitutions j + → 0 (remember that there is no j + at the point t = T ) and j − → −j.Thus it reads Substituting it to the exponential of (2.8) then gives Eq. (2.11) advocated in Section 2.
where all time integrations run from t = 0 to t = T .
Here the restriction of G(t, t ′ ) to G(t, 0) does not lead to essential simplification whereas G(0, 0) has, as shown in Appendix C, the following explicit and simple form in terms of the parameters of the density matrix where the "ratio" of matrices I +X and Ω re is unambiguous because these matrices are commuting in view of the special form of Ω subject to the relation XΩX = Ω * .

Keldysh rotation
For further convenience it is useful to perform the change of the basis in the doubled field space ϕ + , ϕ − and introduce the so-called classical and quantum fields ϕ c and ϕ q [44,45], This transformation is called Keldysh rotation.In the new basis, the Green's function G takes the form Here G R and G A are the retarded and advanced Green's functions, respectively, having the following operator form They are consistent with the classical definition (3.38), in particular, because of independence of the commutator average of the state ρ.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

. Special choice of basis functions and particle interpretation
Thus far, the matrix ω, 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 ω 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 so that from (4.43) Note that the matrix κ is symmetric, whereas ν is Hermitian.Comparing with (4.43) we find the connection between particle averages and the matrix ν Thus, we see that block-diagonal components of ν are responsible for anomalous averages.To ascribe the particle interpretation to the creation/annihilation operators, we will try to choose the matrix ω, defining the corresponding basis functions v(t) and v * (t) so that the diagonal blocks of ν, defining the anomalous averages κ, vanish.Moreover, this choice will simplify the expressions for the Green's functions, since they contain the terms, containing κ.For example, with a nonzero κ the Wightman function reads To make the matrix ν block off-diagonal consider the expression (4.32) and note that the only block diagonal contribution is contained in the identity matrix I and, possibly, in the term involving (ω + Ω) −1 .Thus, we want to choose ω such that block-diagonal contribution of the latter exactly cancels those of I. Using the block matrix inversion formula 6 we have the condition of the vanishing block-diagonal part of ν, (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 the equation (4.48) can be rewritten as r + I − s(r + I) −1 s = 2 and further simplified by introducing the new variable s = (r + I) −1/2 s (r + I) −1/2 and solving for s, so that it takes the following form This is implicit equation on ω, due to the above definition of r and s.Its explicit form reads which can be solved in the form advocated in Section 2 6 The useful form of the block matrix inversion formula is Note that the assumption of positive definiteness of ω implies that I − σ 2 = (I − σ)(I + σ) is positive definite.Recalling that R + S = R 1/2 (I + σ)R 1/2 should be positive definite for normalizability of the density matrix, it is easy to see that I − σ = 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 ω to (4.32) gives the desired block-diagonal matrix form of ν advocated in Section 2 where the matrix κ introduced above is orthogonal.Therefore, as a consequence of positive definiteness of I + σ and I − σ, the matrix ν is necessarily real.As shown in Appendix D, for the density matrix to be positive definite, the matrix σ should be negative definite, so ν 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 while the others can be expressed through them in a straightforward way.
It will be useful to express Ω in terms of ν.Disentangling Ω from (4.32), and then using the explicit form (4.53) of ν, corresponding to the special choice of Neumann basis functions with (4.52), we obtain the following expression Now, let us focus on the particular Gaussian state, which is obtained from the Euclidean path integral, namely Here S E is the quadratic action of Euclidean field theory within time limits τ ± which we will chose to be τ + = β and τ − = 0, , (4.59) (4.60) The partition function Z in the normalization factor is such that tr ρ E = 1 for vanishing source J E = 0. Hermiticity of the density matrix implies the following (sufficient) condition on the coefficient matrices A E , B E , and C E as the functions of τ which are not necessary real.Nevertheless, we restrict ourselves to the real case below.The source J E is included in the path integral to be able to introduce nonlinear terms of the Euclidean action, leading to the non-Gaussinities of the resulting density matrix.We take the path integral (4.58) over the Euclidean fields ϕ by using the saddle point method.The boundary conditions of the integral fix the endpoints ϕ(β) = φ + , ϕ(0) = φ − , so we have the boundary problem with the Dirichlet boundary conditions Using the Dirichlet Green's function G D for vanishing boundary conditions and substituting to (3.52), one expresses the solution of (4.62) as follows where we introduce the notations, similarly to those of the Lorentzian context (3.53)-(3.54),for the row w T E (τ ) obtained by the transposition of the column w E (τ ) in where we disregard the source independent prefactor, all τ -integrations run from 0 to β, whereas the matrix Ω and the source j, introduced in (4.5) take the following particular form Now, one can substitute the density matrix (4.66), defined by the parameters (4.67) to the general expression for the generating functional (4.37).This leads to Note that the kernel of the third integral here is the periodic Euclidean Green's function 70) corresponding to the fact that with the Lorentzian sources switched off the functional Z[0, J E ] represents the Euclidean path integral over periodic fields ϕ(τ ) on the time interval with the identified boundary points τ ± .The expression for this Green's function seemingly dependent via G(0, 0) on Lorentzian objects is in fact independent of them.This property is based on the relation (4.38) and derived in Appendix C.

Analytic continuation and KMS condition
The further transformation of the generating functional, which allows one to reveal its new analyticity properties, can be performed due to two assumptions.The first assumption is that the Euclidean action (4. where g N (τ ) is the Euclidean version of the definition (3.59) for the Neumann Green's function.
To proceed further we have to derive several important properties of Euclidean Neumann Green's function which is the part of (4.73), specific to the choice (4.52) of ω.In terms of the Euclidean basis functions it reads as Here u + , u − are the basis functions obeying Neumann boundary conditions and, as usual, where we use the explicit form (4.53) of ν, corresponding to the particular choice of basis functions described in the previous subsection.Equating these two expressions for G N ∥, with due regard to the structure of ∆ N +− in (4.76), we find the two sets of equalities.The first set follows from the diagonal blocks where the basis functions u + , u − are evaluated either at τ = 0 or τ = β.Similarly, from the off-diagonal blocks of (4.78), one gets the formulas, relating the boundary values of the basis functions It is useful to continue the Euclidean equations of motion beyond the interval 0 < τ < β with the period β (which is again possible because τ = 0 and τ = β) are the turning points), Together with (4.61) it also implies

.83)
Once the functions u ± (τ ) satisfy the same homogeneous boundary conditions for both τ = 0 and τ = β (cf.(4.75) and (4.79)), being translated by the period they can only differ by the multiplication with some nonsingular matrices L ± , u ± (τ +β) = u ± (τ )L ± .From (4.81) we obtain their explicit form With the normalization this monodromy simplifies to then G E (τ, τ ′ ) can be expressed as .91) and the Lorenzian Wighmann Green's function (4.56) is an analytic continuation of G > E (τ ).Now, it is time to connect the Euclidean basis functions u ± and the Lorentzian ones v, v * .Specifically, let us show that both sets of functions can be obtained from a single function V (z) of the complex time z = t − iτ , obeying complexified equations of motion (3.5) .93) which reduces to those for v or u + after the substitution z = t or z = −iτ , respectively.Supplementing the latter condition with the normalization one finds i.e. v and u + are analytic continuation of each other.Similarly, complex conjugation of (4.93) and the same assumptions of coefficient functions reality and its reflection symmetry, we find that V * obeys the following boundary condition so that v * and u + can be obtained from V * as

.97)
Thus, assuming that the complexified basis function V (z), z = t − iτ is analytic on 0 ≤ t ≤ T , 0 ≤ τ < β, we have the following transformation law of the basis functions Substituting to (4.56), one has the following condition on Wightmann Green's function which is nothing but KMS condition advocated in Section 4.7.

Harmonic oscillator
In this section we consider harmonic oscillator as the simplest instructive example, which demonstrates the main concepts and quantities, introduced above, together with the convenience of the special choice of the basis functions v, v * .The corresponding action reads where ϕ is one-component field, defining the coordinate of oscillator, and ω 0 is its frequency.We will consider the system in the state, defined by the Euclidean path integral (4.58),where the Euclidean action is an analytic continuation (4.71) of the Lorentzian one Note that for J E = 0 density matrix (4.58) coincides with the thermal density matrix of the inverse temperature β.
The corresponding differential operator defining the Euclidean equation of motion F E ϕ E = 0 and the Wronskian read To exploit the answer (4.66), one should first calculate the Dirichlet Green's function, which can be constructed out of corresponding basis functions u D ± (τ ) satisfying These basis function can be chosen as so that Dirichlet Green's function has the following form Substituting the Green's function obtained to (4.67), one finds the explicit form of the density matrix constituents where we assume ω to be real for the simplicity.
This basis inherits the properties of y ± (τ ) under translation by period, reflection and complex conjugation.In particular u ± (τ + β) = e ∓βε u ± (τ ). (5.23) Comparing with (2.47) one concludes that the parameter ε is connected to ν as ν = 1 e βε − 1 . (5.24) The basis functions u ± (τ ) have significantly different frequency properties depending on whether ε is real or imaginary.Thus, real ε implies where ω is a real number, which coincides with those defined in (2.20), as will be described below.In contrast, imaginary ε leads to the property (u ± (τ )) * = u ± (−τ ), so that the fraction is imaginary 7 , and one can write where the number ω ′ = iω is real.Let us calculate the density matrix (4.58) and examine its properties.To use the answer (4.66), one should first construct the Dirichlet Green's function.The corresponding basis functions u D ± (τ ) obeying u D − (0) = u D + (β) = 0 can be constructed as linear combinations of u ± (τ ).Namely, one defines u D − (τ ) as so that u D − (0) = 0 due to u − (τ ) = u + (−τ ).Due to reflection symmetry of (5.16) one can obtain u The corresponding Wronskian of u D + and u D − reads where we use the relations (5.27)-(5.28) between Dirichlet basis functions and u ± (τ ), and its derivatives at the boundary points (5.30) 7 In deriving this property we use that Substitution of the corresponding Dirichlet Green's function to (4.67) gives where ω is defined in (5.25).Note that for real ε this coincides with (4.57), with (5.24) substituted.For imaginary ε we express it as ε = iq, so Ω has the form where ω ′ is defined in (5.26).Following Appendix D, let us examine the properties of the underlying density matrix, defined by the obtained Ω.For real ε we have R = ω coth βε and S = −ω/ sinh βε, so that R, R + S and R − S have the same sign, and we conclude that the density matrix is bounded, normalizable and positive-definite for ω > 0.
If it is the case, σ ≡ S/R = −1/ cosh βε, so the definition (5.25) is consistent with (2.20), and particle interpretation is allowed.In contrast, for imaginary ε we have R = ω ′ cot βq, S = −ω ′ / sin βq, so R + S and R − S have different signs, so even if the density matrix is normalizable, the particle interpretation is not available.

The case of a pure state: vacuum no-boundary wavefunction
As we have shown above, the Euclidean density matrix prescription in a rather nontrivial way suggests a distinguished choice of the particle interpretation.In context of the pure Hartle-Hawking state this fact is well known and takes place in a much simpler way.Let us briefly discuss this here along with a general demonstration how the transition from a mixed state to the pure one proceeds via the change of spacetime topology of the underlying Euclidean instanton from Fig. 1 to Fig. 3.
The no-boundary state defined by the path integral over the fields on the Euclidean "hemisphere" D 4 + of Fig. 3 (and its reflection dual on D 4 − considered as a factor in the factorizable pure density matrix of Fig. 3) is the vacuum wavefunction (3.109) with the real frequency (3.107), ω = [iW v(t)][v(t)] −1 | t=0 .The relevants positivefrequency basis function v(t), similarly to (2.50), can be regarded as the analytic continuation of a special Euclidean basis function u(τ ), v(t) = u(τ + + it).This basis function is selected by the requirement that it is regular everywhere inside D 4 + , including its pole which we label by τ = 0 [12,48].
To show this one should repeat the calculation of Section 4. where the second equality follows from the analytic continuation rule v(t) = u(τ + + it).Thus, the Hartle-Hawking no-boundary wavefunction of the linearized field modes is the vacuum of particles uniquely defined by a particular choice of positive-frequency basis functions v(t) which in their turn are the analytic continuation of the regular Euclidean basis functions u(τ ), v(t) = u(τ + +it). 9This is a well-known fact [12,48] which in the case of de Sitter cosmology corresponds to the Euclidean de Sitter invariant vacuum [13,14].
It is known that vacuum in-in formalism in equilibrium models can be reached by taking the zero temperature limit β → ∞.It is not quite clear how this limit can be obtained in generic non-equilibrium situations, but it is likely that the transition from mixed Euclidean density matrix to a pure state is always associated with ripping the Euclidean domain into two disjoint manifolds D 4  + and D 4 − depicted in Fig. 3. To show this consider generic situation of the mixed state with the Euclidean density matrix of Fig. 1.This density matrix has a Gaussian form (2.4)-(2.6)with the matrix Ω given by Eq. (2.31) with the Dirichlet Green's function which can be represented in terms of two sets of Dirichlet basis functions u D ± (τ ), u D ± (τ ± ) = 0, (5.35) 8 The point τ = 0 is an internal regular point of a smooth manifold D 4 + , so that this point with τ treated as a radial coordinate turns out to be a regular singularity of the equation F E ϕ(τ ) = 0. Its two linearly independent solutions u ∓ (τ ) have the asymptotic behavior u ∓ ∝ τ µ ∓ with µ − > 0 > µ + , so that only u − (τ ) ≡ u(τ ) is the regular one, while the contribution of the singular u + (τ ) → ∞, τ → 0, should be discarded from the solution ϕ(τ ) [48]. 9 The set u(τ ) is of course defined only up to a linear transformation with some constant matrix L, u(τ ) → u(τ )L, v(t) → v(t)L, but this Bogoliubov transformation does not mix frequencies and therfore does not change particle interpretation.
Now consider the case of a pure state, when the density matrix factorizes into the product of two wavefunctions, or the situation of Ω +− ≡ S = 0.This off-diagonal block of Ω reads as where we used the fact that in view of boundary conditions on u D ± (τ ).Therefore, the requirement of S = 0 implies singularity of u D + (τ − ) which is impossible, because the Green's function G D (τ, τ ′ ) can have a singularity only at the coincidence point of its arguments τ = τ ′ .This means that no Dirichlet Green's function on a smooth connected Euclidean manifold of the topology [τ − , τ + ]×S 3 can generate the density matrix of a pure state.The only remaining option is ripping the bridge between Σ + and Σ − into the union of two disjoint parts D 4 ± by shrinking the middle time slice at τ ≡ τ++τ− 2 to a point.In context of the cosmological model driven by the set of Weyl invariant quantum fields [9,16,22] this option also matches with the interpretation of zero temperature limit β → ∞, because the inverse temperature of the gas of conformal particles in this model is given by the instanton period in units of the conformal time β = 2 τ+ τ dτ /a(τ ) → ∞, which diverges because the cosmological scale factor (the size of the spatial S 3 -section) a(τ ) → 0 at τ → τ .

DISCUSSION AND CONCLUSIONS
Generality of the above formalism allows one to apply it to a wide scope of problems ranging from condensed matter physics to quantum gravity and cosmology.Our goal in future work will be its use in the calculation of the primordial CMB spectrum of cosmological perturbations in the model of microcanonical initial conditions for inflationary cosmology [9,16,20], which was briefly discussed as a motivation for this research.Quasi-thermal nature of this setup was associated in these papers with the fact that the model was based on local Weyl invariant (conformal) matter which, on the one hand, generates the Friedmann background providing the necessary reflection symmetry and, on the other hand, turns out to be effectively in equilibrium, because in the comoving frame it describes a static situation.
Our results show, however, that thermal properties, including particle interpretation with the distinguished positive/negative frequency decomposition, are valid in much more general case.Specifically, the corresponding frequency matrix ω in the initial conditions problem for basis functions (2.15) is shown to be determined by the parameters of Gaussian type density matrix (2.20), and the occupation number matrix ν reads as (2.21)- (2.22).In this setup, the Euclidean density matrix, which incorporates the reflection symmetry property guaranteed by (4.61), plays the role of the particular case.If in addition the Lorentzian action is related to the Euclidean action via the analytic continuation at the turning points of the bounce background (which, of course, respects its reflection symmetry), important analytic properties of correlation functions, including the KMS condition, begin to hold.These are the main results of the paper.They allow one to derive the full set of Lorentzian domain, Euclidean domain and mixed, Lorentzian-Euclidean, Green's functions of the in-in formalism and reveal its rich analytic structure.In particular, the results of Section 4.2 significantly extend those of [49], where the nonequilibrium evolution of Gaussian type density matrices was examined.The discussion of simple application examples of Section 5 shows the relation of the obtained formalism to the stability properties of dynamical systems in Floquet theory and the theory of Bloch functions.These properties, in their turn, are related to the eigenmode properties of the wave operator F E subject to periodic boundary conditions on the bounce instanton within Euclidean time [ τ − , τ + ]-range and deserve further studies.
Prospective nature of rich analytic structure of the Euclidean-Lorentzian in-in formalism consists in the hope that quantum equivalence of purely Euclidean calculations of loop effects with those of the Lorentzian calculations can be extended to generic bounce type backgrounds.This equivalence was proven in [31,32] for the vacuum case of the flat chart of the de Sitter spacetime vs its Euclidean -S 4 instanton.A similar but much simpler equivalence at the one-loop order was observed within covariant curvature expansion in asymptotically flat spacetime for systems with the Poincareinvariant vacuum which is prescribed as the initial condition at asymptotic past infnity [50].This equivalence is realized via a special type of analytic continuation from Euclidean to Lorentzian spacetime, which guarantees unitarity and causality of relevant nonlocal form factors.
Further applications of the in-in formalism in quantum cosmology require its extension to models with local gauge and diffeomorphism invariance (see also [51] for related problem in the context of quantum electrodynamics).What have been built thus far is the formalism in the physical sector of the theory for explicitly disentangled physical degrees of freedom.In cosmological models subject to time parametrization invariance time is hidden among the full set of metric and matter field variables, and disentangling time is a part of the Hamiltonian reduction to the physical sector.This reduction shows that the cosmological background can be devoid of physical degrees of freedom (just like Friedmann equation in FRW-metric models does not involve any physical degree of freedom in the metric sector of the system).This might play a major role in handling a zero mode of the wave operator F E , which necessarily arises on the bounce type background [52] and comprises in the cosmological context one of the aspects of the problem of time in quantum gravity [11].This and the other problems of cosmological applications of the in-in formalism go beyond the scope of this paper and will be the subject of future research.

Figure 1 .
Figure 1.Picture of instanton representing the density matrix.Gray lines depict the Lorentzian Universe nucleating from the instanton at the minimal surfaces Σ− and Σ+.

Figure 2 .
Figure 2. Origin of the partition function instanton from the density matrix instanton by the procedure of gluing the boundaries Σ+ and Σ− -tracing the density matrix.

Figure 3 .
Figure 3. Density matrix of the pure Hartle-Hawking state represented by the union of two no-boundary instantons.

Figure 4 .
Figure 4. Origin of the partition function instanton from the density matrix instanton by the procedure of gluing the boundaries Σ+ and Σ− -tracing the density matrix.
Schroedinger picture (labelled by S) it reads as the chronologically ordered operator T

Figure 5 .
Figure 5. Euclidean-Lorentzian contour C on the Riemann surface of complex time z = t − iτ .Wightman functions are periodic in imaginary (Euclidean) time direction with a period β, 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.
.65) (the last equality implies the symmetry of the Dirichlet Green's function, G T D (τ, τ ′ ) = G D (τ ′ , τ )).Substitution back to (4.58) gives 59) is obtained by analytic continuation of the Lorentzian one (3.1),namely iS[ϕ] t=−iτ = −S E [ϕ] (4.71)This implies the following form of the Euclidean action coefficient functionsA E (τ ) = A(−iτ ), B E (τ ) = −iB(−iτ ), C E (τ ) = −C(−iτ).(4.72)Though this requirement sounds rather restrictive, it can be based on the assumptions discussed in Introduction about the properties of the Euclidean background underlying the quadratic action and sandwiched between the two (identified) turning points at which the analytic match between the Euclidean and Lorentzian branches can be done.Another assumption which we use in what follows is the possibility to make a special choice of the Neumann basis functions, derived above.The first step is to rewrite the second and the third terms in the exponential of the generating functional (4.69) in terms of the Euclidean Neumann Green's function G N (τ, τ ′ ) instead of the Dirichlet one, i.e. (W E+ ω)G N (β, τ ′ ) = (W E − ω * )G N (0, τ ′ ) =0 where ω is the same as in (4.23)-(4.24).This is done using the relations (3.65)-(3.66)(after the replacement ω → −iω associated with the transition to the Euclidean version of Dirichlet and Neumann Green's functions) and the derivation in Appendix C. The result reads as the expression (4.69) with the kernel of the Lorentzian-Euclidean term −G(t, 0) w E (τ ) replaced by G(t, 0) (ω + Ω) g N (τ ) and the new form of the periodic Green's function G E (τ, τ ′ ) in the Euclidean-Euclidean block .76) Note that the boundary conditions on u ± above are exactly the analytic continuation t → −iτ of the boundary conditions (4.26) on v, v * .Now, consider in detail the matrix of boundary values of the Euclidean Neumann Green's function at τ + = β and τ − = 0 77) (double vertical bar denotes here the restruction of the two Green's function arguments to two boundary surfaces thus forming the 2×2 block matrix).Using the Euclidean version of the relation (3.66), we find the alternative form of this matrix view of the reflection symmetry (4.83) of the operator F E the functions u + (τ ) and u − (−τ ) can differ at most by some non-degenerate matrix L, u + (τ ) = u − (−τ ) L. For the normalization (4.85) this implies u + (τ ) = u − (−τ ).(4.87)For the choice (4.85) we have ∆ N +− = −∆ N −+ = I, so that the blocks of the Euclidean and Lorentzian-Euclidean Green's function in (4.73) read

. 9 )
The basis functions, satisfying (4.26) are the linear combinations of e ±iω0t which are the solutions of e.o.m. and read (cf.(3.96) and (3.103)) 6 on D 4 + -the support of the Euclidean action S E (φ) evaluated at the regular solution of equations of motion F E ϕ(τ ) = 0 with the boundary value φ = ϕ(τ + ) at the single boundary Σ + = ∂D4  + .This regular solution is given by the expression proportional to the regular basis function u(τ) of F E on D 4 + , ϕ(τ ) = u(τ )[u(τ + )] −1 φ,(5.33)because the contribution of the complementary basis function dual to the regular u(τ ) should be excluded in view of its singularity at τ = 0. 8 After the substitution into the expression for the action (4.59) its onshell value reduces to the contribution of the single surface term at Σ + , S E (φ) = 1 2 ϕ T (W E ϕ)| Σ+ .As a result S E (φ) = 1 2 φ T ωφ, and the Hartle-Hawking wavefunction Ψ HH (φ) ∝ e −S E (φ) becomes the vacuum state (3.109) with ω = −[W E u(τ + )][u(τ + )] −1 = [iW v(t)][v(t)] −1 t=0 , (5.34) which can be obtained from (2.14) by the replacement (2.27), subject to Dirichlet boundary conditions, t) P.(3.20)Before proceeding further, let us show explicitly that the right hand side(3.19)isindeedindependent of time t.To demonstrate this, we contract l.h.s of the equation(3.5)wherefield ϕ = ϕ 1 , with another field ϕ 2 , and subtract the same quantity, but with F , acting on ϕ 2 (ϕ 1,2 are not necessarily solve e.o.m.).The result can be written as .79) and means that the basis functions u + , u − obey the same Neumann boundary conditions at both boundary values of the Euclidean time (cf.Eq. (4.75)).This also implies the following explicit form of the matrices ∆ .92) This equation reduces to the Lorentzian e.o.m. for z = t and to the Euclidean ones for z = −iτ .Under the assumption that coefficient functions A(t), B(t), and C(t) are real, together with the reflection symmetry (4.83), one can find that V * (z) ≡ (V (z * )) * obeys the same equation.Moreover, the initial conditions (4.26) for v, v * are connected with those (4.75) for u ± via analytic continuation t → −iτ .This motivates us to impose the boundary condition on V as follows