Initial value formulation of a quantum damped harmonic oscillator

The in-in formalism and its influence functional generalization are widely used to describe the out-of-equilibrium dynamics of unitary and open quantum systems, respectively. In this paper, we build on these techniques to develop an effective theory of a quantum damped harmonic oscillator and use it to study initial state-dependence, decoherence, and thermalization. We first consider a Gaussian initial state and quadratic influence functional and obtain general equations for the Green's functions of the oscillator. We solve the equations in the specific case of time-local dissipation and use the resulting Green's functions to obtain the purity and unequal-time two-point correlations of the oscillator. We find that the dynamics must include a non-vanishing noise term to yield physical results for the purity and that the oscillator decoheres in time such that the late-time density operator is thermal. We show that the frequency spectrum or unequal-time correlations can, however, distinguish between the damped oscillator and an isolated oscillator in thermal equilibrium, and obtain a generalized fluctuation-dissipation relation for the damped oscillator. We briefly consider time-nonlocal dissipation as well, to show that the fluctuation-dissipation relation is satisfied for a specific choice of dissipation kernels. Lastly, we develop a double in-out path integral approach to go beyond Gaussian initial states and show that our equal-time results for time-local dissipation are in fact non-perturbative in the initial state.


I. INTRODUCTION
Green's functions are extremely useful for calculating correlation functions of quantum systems, both in-and out-of-equilibrium.In an equilibrium quantum field theory (QFT), for example, we typically assume that the field is in the ground state in the infinite past and future and that any interaction is turned on and off adiabatically.Specifying both initial and final conditions picks out the Feynman Green's function as the primary Green's function, in terms of which we can obtain any time-ordered correlation function of the field.In an outof-equilibrium QFT, on the other hand, we are typically interested in finite-time correlations, with the field initialized in any state at a finite initial time.Specifying initial conditions now picks out the retarded Green's function as the primary Green's function, in terms of which we can obtain any field correlations.
The Green's functions of a given quantum system are most readily obtained through the path integral approach.The standard in-out path integral, for example, allows us to obtain the generating functional and hence correlation functions in an equilibrium QFT.Its out-ofequilibrium generalization or the in-in path integral [1][2][3][4][5] instead allows us to obtain the generating functional and hence correlation functions in an out-of-equilibrium QFT.While the original in-in path integral only describes unitary dynamics, i.e., the dynamics of a closed quantum system that is initialized in any state and evolves with a time-independent or time-dependent Hamiltonian, it can be further generalized to describe the non-Hamiltonian dynamics of an open quantum system by introducing an influence functional [6,7].The influence functional is the path integral analog of the quantum master equation [8,9] and simplifies the calculation of certain quantities in open quantum systems, such as Green's functions and unequal-time correlations [10].
The influence functional method has been used in a variety of problems since it was first developed.Amongst exactly solvable ones are the quantum damped harmonic oscillator (DHO) and linear Brownian motion, for which a standard reference is the textbook [11].An incomplete list of other problems for which it has been used is the study of quantum Brownian motion in different environments [12,13], quantum transport in interacting nanojunctions [14], decoherence in interacting QFTs [15,16] and inflation [17], entanglement in primordial correlations [18,19], coarse-graining in interacting QFTs [20,21], and open holographic QFTs [22].Given the versatility of the method, in this paper, we revisit one of the simplest out-of-equilibrium quantum systems described by an influence functional -a quantum DHOwith a goal of developing an effective field theory-inspired approach to the problem.We are thus interested in constraining parameters that appear in the influence functional on physical grounds, without knowledge of the full microscopic model.
We first initialize the oscillator in a Gaussian state and evolve it with a general but quadratic influence functional while remaining agnostic to the source of dissipation.We then specialize to time-local dissipation, find exact solutions for the Green's functions in this case, and show that terms in the influence functional are constrained, due in large part to fluctuation-dissipation relations of environment degrees of freedom.We find, in particular, that the influence functional must contain a nonvanishing noise term, which leads, when set to zero, to a non-physical late-time purity of the oscillator.Using purity as an indicator of decoherence, we next show that the DHO decoheres in time and settles into a thermal state at a temperature defined by the dissipation parameters.We show that the frequency spectrum or unequal-time correlations can, however, distinguish between the DHO and an isolated oscillator in thermal equilibrium.We use unequal-time correlations in the late-time limit to obtain a generalized fluctuation-dissipation relation and show that it reduces to the usual relation only in the hightemperature regime.We further consider time-nonlocal dissipation with a specific choice of dissipation kernels, where the fluctuation-dissipation relation is satisfied at any temperature.Finally, we develop a double in-out path integral approach that allows us to obtain equaltime correlations for any initial state and show that our equal-time results for time-local dissipation are in fact non-perturbative in the initial state.
The paper is organized as follows.In Sec.II, we review the in-in formalism in the context of a harmonic oscillator, deriving the corresponding generating functional in an Appendix for completeness, and discuss the initial state, which we choose to be Gaussian, and the influence functional, which we choose to be quadratic.In Sec.III, we write the generating functional in terms of Green's functions, use them to obtain n-point correlations, again relegating details to an Appendix, find the initial conditions they must satisfy, and obtain their equations of motion.In Sec.IV, we restrict to time-local dissipation and obtain exact solutions for the Green's functions in this case, highlighting the contribution from the noise term.We use these solutions to understand how the oscillator's purity evolves in time and whether it satisfies the fluctuation-dissipation relation at late times.In Sec.V, we consider the fluctuation-dissipation relation for time-nonlocal dissipation with a specific choice of dissipation kernels.In Sec.VI, we develop a double in-out path integral approach that allows us to go beyond Gaussian initial states and use it to show that our results for equal-time correlations and purity in the case of timelocal dissipation hold for any initial state.We end with a discussion in Sec.VII.

II. IN-IN GENERATING FUNCTIONAL
In this section, we obtain the generating functional of a quantum harmonic oscillator that is initialized in a Gaussian initial state and undergoes non-unitary/dissipative evolution described by a quadratic dissipative action; also see Ref. [23] for a recent review, Ref. [24] for an earlier review, and Ref. [25] for related recent work.We denote the initial state of the oscillator at time t 0 by the density operator ρ(t 0 ) and first ignore dissipation, so that the oscillator evolves unitarily with the Hamiltonian Ĥ = 1 2m P 2 + 1 2 mω 2 X2 , m being the mass of the oscillator, ω its frequency [26], X the position operator, and P the momentum operator.Say we evolve ρ(t 0 ) to a final time t f , which we take to be later than any times of interest, in the presence of a source J + (t) on the forward branch (from t 0 to t f ) and J − (t) on the backward branch (from t f to t 0 ).The in-in generating functional is then defined as Z[J + , J − ] = Tr [ρ(t f )] J + ,J − and finite-time correlation functions of the Heisenberg picture operator X(t) can be obtained by taking functional derivatives of Z[J + , J − ] with respect to the two sources J ± (t), which are typically set to zero at the end of the calculation.
It is convenient to write Z[J + , J − ] in the path integral representation as that allows us to express it in terms of Green's functions.Let us thus define eigenkets and eigenvalues of the Schrödinger picture operator XS , which we denote with |x⟩ and x, so that XS |x⟩ = x|x⟩.As shown in Appendix A, one finds that where ) in position basis with the functions x ± (t) evaluated at t 0 , S[x ± ] is the action of the oscillator, the shorthand t stands for t f t0 dt, the δ-function at the end imposes the boundary condition at the turn-around point, and we have set ℏ to unity.The action S[x ± ] is given by 1  2 t ( ẋ± ) 2 − ω 2 (x ± ) 2 , where we have set the mass m to unity for simplicity and without loss of generality and the dot denotes a derivative with time.Note that Z[J, J] = 1 since in this case the forward and backward evolution exactly cancels out or, in other words, Z[J, J] is simply Tr [ρ(t f )] for evolution in the presence of a source, which is normalized to unity.
Let us now introduce dissipation in the dynamics.Dissipation is described by an influence functional that leads to time-nonlocal terms on each of the forward and backward branches of evolution and additionally ties together terms on the two branches.We thus rewrite the generating functional of Eq. ( 1) as where we have added the influence functional S IF [x + , x − ] to the action.We must still have that Z diss.[J, J] = 1 since Tr [ρ(t f )] for evolution in the presence of a source is normalized to unity for both unitary and non-unitary evolution.
In order to explicitly calculate the generating functional, we need an ansatz for the initial state and influence functional, which we discuss in the two subsections below.

A. Initial state
We choose a Gaussian initial state for the oscillator, parameterized as ρ[x + , x − , t 0 ] = N exp {iS 0 [x + , x − ]}, N being a normalization constant chosen so that Tr [ρ(t 0 )] = 1 and [27,28] with the functions x ± (t) evaluated at t 0 , as mentioned earlier.Note that B here must be real since ρ(t 0 ) is Hermitian.A, on the other hand, can be complex and we denote its real and imaginary parts with A R and A I .The normalization N is then found to be N = (A I + B)/π with the condition that A I + B > 0. x 0 and p 0 are the initial one-point functions X(t 0 ) and P (t 0 ) , where we have used angular brackets to denote the expectation value in ρ(t 0 ).Lastly, A and B are related to the initial two-point correlators, that we denote as c xx , c xp , and c pp for convenience, The subscript 'c' here denotes connected correlators, for example, X2 c = X2 − X 2 , and {•, •} is the anticommutator.We can also calculate the purity of our initial state, which we denote Pu(t 0 ) = Tr ρ2 (t 0 ) , and relate it to the initial correlators as Pu(t xp .Since purity must be between zero and one, we also have the condition that B ≤ 0, with the initial state being pure for B = 0 and mixed for B < 0.

B. Influence functional
We choose the influence functional to be quadratic and parameterize it as where γ 1 (t, t ′ ) and γ 2 (t, t ′ ), which we refer to as dissipation kernels, are complex functions.The complex conjugates and signs in Eq. ( 7) ensure that the generating functional is real or, equivalently, the density operator is Hermitian at all times.The dissipation kernels must further satisfy γ 1 (t, t ′ ) = γ 1 (t ′ , t) and γ 2 (t, t ′ ) = −γ * 2 (t ′ , t) since the integration measure is symmetric under the interchange of t and t ′ .We find in the next section that only two real functions contribute to dissipation and show later in the paper that they too are not independent of one another.We also note that in a microscopic derivation of the influence functional, we would typically take the oscillator to be coupled to some environment degrees of freedom, with the full system plus environment evolving unitarily.Performing a partial trace over the environment, we would find that the dissipation kernels here are related to correlation functions of the environment and additionally satisfy γ We will not impose this constraint here, except when mapping to a specific microscopic model.

III. GREEN'S FUNCTIONS
In this section, we write the generating functional in terms of Green's functions and obtain the equations of motion and initial conditions that they satisfy.To introduce Green's functions, it is convenient to first express the integrand in Eq. ( 2) as an exponent.Since we have already written the initial state and influence functional as exponents in Secs.II A and II B, we only need to further express the δ-function in Eq. ( 1) as an exponent.Following [29], the δ-function can be written as where ) and the extra factor of two is needed to cancel the factor of half that arises from evaluating the δ-function at a limit of the integral.We will also integrate by parts the kinetic term in the action S[x ± , J ± ], to move both time derivatives to a single x ± .This generates boundary terms at t 0 and t f , of which those at t f cancel out between the plus and minus branches given that x + (t f ) = x − (t f ) and ẋ+ (t f ) = ẋ− (t f ), where the second condition is shown to follow from the first in Appendix A. We also note in Appendix A that, in the presence of dissipation, ẋ+ (t f ) = ẋ− (t f ) only holds for dissipation kernels γ 1 (t, t ′ ) and γ 2 (t, t ′ ) that are not proportional to d 2 δ(t − t ′ )/dt ′2 .The boundary terms at t 0 remain, and we keep track of them below.
Putting everything together, we can now write the generating functional in Eq. ( 2) as The vectors x(t) and J s (t) here are given by with the extra factors of 2 in front of the δ(t ′ − t 0 ) terms again arising from evaluating the δ-function at a limit of the integral.Let us now define a 2 × 2 matrix G(t, t ′ ) that satisfies the following Green's function equation where and I is the 2 × 2 identity matrix.This leads to four equations of motion that couple the functions G ±± (t, t ′ ).Since the ϵ that appears in O(t, t ′ ) is arbitrarily small, all terms proportional to 1/ϵ in Eq. ( 13) must cancel out, giving us the constraints for all t 0 < t ′ < t f .Note that these constraints follow directly from the requirement that x + (t f ) = x − (t f ).We similarly expect the derivatives at t f to also satisfy the same constraints, so that for all t 0 < t ′ < t f , in analogy with the condition that ẋ+ (t f ) = ẋ− (t f ).One way to see this is to first write the solution for the classical field that satisfies the equation in the presence of any source Then realizing that the boundary conditions are carried by the classical field, i.e., ), all four constraints in Eqs. ( 15)-( 18) follow.We note again, however, that this argument only holds for dissipation kernels that are not proportional to d 2 δ(t − t ′ )/dt ′2 .All four constraints in Eqs. ( 15)-( 18) are needed to write the generating functional in standard form.To do so, we first shift x(t) in Eq. (9) to x s (t) = x(t)−i t ′ G(t, t ′ )J s (t ′ ), then integrate by parts to move the derivatives from x ± s (t) to G ±± (t, t ′ ), and lastly make use of the four constraints above to cancel the boundary terms.The generating functional can then be written as where the normalization Z 0 is chosen such that Z diss.[J, J] = 1 and it turns out to be unity in the case that the initial state is a coherent state.Further, since we can freely interchange the t and t ′ integrals in the exponent, the functions G ±± (t, t ′ ) must also satisfy the conditions in addition to the constraints written earlier.
Writing the generating functional in the form of Eq. ( 19) makes it easy to obtain correlation functions as usual, and we discuss this further in the first subsection below.In the next two subsections, we obtain the initial conditions and equations of motion for G ±± (t, t ′ ), and sketch how to solve the resulting equations.

A. n-point correlations
As shown in Appendix B for the case of unitary dynamics, one-and two-point correlation functions of X(t) are easily obtained by taking functional derivatives of Z[J + , J − ] with respect to J ± .We can write similar expressions for the case when the oscillator is coupled to an environment, with the full system plus environment evolving unitarily, and then trace out the environment.Correlation functions in the presence of dissipation can, therefore, be obtained by taking functional derivatives of Z diss.[J + , J − ] instead.Generalizing the calculation in Appendix B to the dissipative case and additionally beyond one-and two-point correlations gives T X(t n+1 ) . . .X(t n+m ) T X(t 1 ) . . .
in the presence of a source J(t), where the first n operators are time-ordered, denoted by T , and the next m are anti-time-ordered, denoted by T .We now want to relate n-point correlations to the Green's functions by making use of Eq. (19).Let us first consider the one-point function X(t) .Setting n = 1 and m = 0 in Eq. ( 23) and using Eq. ( 19), we find that We expect ⟨ X(t)⟩ to match the solution to the classical equation of motion.Let us next consider the two-point correlations T X(t) X(t ′ ) and X(t ′ ) X(t) .These yield the functions G ++ (t, t ′ ) and G +− (t, t ′ ) plus a product of one-point expectation values and, therefore, We can similarly show that T X(t) , where T denotes antitime-ordering.Since these are also equal to the Hermitian conjugates of the above expressions, we further have that We can now write the functions G ±± (t, t ′ ) in a more convenient form.Let us denote G +− (t, t ′ ) as Then we can write and G −− (t, t ′ ) as the complex conjugate of Eq. (27).Note that the functions G ±± (t, t ′ ) written as above satisfy all constraints and conditions written earlier in this section.By using Eqs.( 27) and (28) and similar expressions for G −− (t, t ′ ) and G −+ (t, t ′ ), we can verify that they additionally satisfy which is needed to ensure that Z diss.[J, J] = 1 and will also turn out to be useful later in this section.
Let us lastly consider the n-point correlation T X(t 1 ) . . .X(t n ) obtained by setting m = 0 in Eq. (23).Using again Eq. ( 19) for Z diss.[J + , J − ] and expanding out the J + (t i ) derivatives on the right hand side gives us a product of all G ++ (t i , t j ) plus a sum of disconnected correlators for even n and only a sum of disconnected correlators for odd n.Taking the disconnected correlators to the left hand side then turns it into a time-ordered product of X(t) − X(t) at times t 1 , . . ., t n .The timeordered product of X(t)− X(t) , therefore, obeys Wick's theorem.We can also obtain this result in a simpler way, and it is thus worth discussing, by writing Eq. ( 19) in a slightly different form.On expanding out the shifted sources on the right hand side of Eq. ( 19) and rearranging the resulting expression, we can write it as where ⟨ X(t)⟩, given in Eq. ( 24), is calculated in the presence of a source J(t) and J c (t) = J + (t) − J(t) −J − (t) + J(t) .Taking derivatives of Eq. ( 30) with respect to J + (t i ) and setting J ± (t) = J(t) now produces the time-ordered product of X(t)− X(t) at times t 1 , . . ., t n on the left hand side and a product of G ++ (t i , t j ) for even n and zero for odd n on the right hand side, directly giving us Wick's theorem.

B. Initial conditions
In the previous subsection, we showed that connected two-point correlations of X(t) can be identified with the functions G ±± (t, t ′ ).We next want to argue that the Heisenberg picture operator P (t) is given by P (t) = Ẋ(t), so that we can additionally relate the connected twopoint correlations that involve P (t) with time derivatives of G ±± (t, t ′ ).If we substitute the form of G ++ (t, t ′ ) given in Eq. ( 27) into its equation of motion given by Eq. ( 13) and equate the δ-function on both sides, we find that it yields the Wronskian condition , the Wronskian condition further implies that X(t), Ẋ(t) = i which, together with the commutation relation X(t), P (t) = i, suggests that indeed P (t) = Ẋ(t).We will restrict to those influence functionals that preserve P (t) = Ẋ(t) in this paper, but note that this is not guaranteed for more general influence functionals, specifically those that arise from derivative system-environment interactions.
With P (t) = Ẋ(t), we can directly transcribe the initial conditions given in Eqs. ( 4), (5), and ( 6) to conditions on G ±± (t, t ′ ), in particular on G < (t, t ′ ), From the Wronskian condition, or the commutator of X(t) and P (t) at the initial time, we additionally have that As shown later, the three initial conditions in Eqs. ( 31), (32), and ( 33) are sufficient to solve the equations of motion for the symmetrized two-point correlation that we introduce in the next subsection.
It is also worth noting that the initial conditions written here are consistent with the δ(t ′ − t 0 ) terms in the equations of motion in Eq. ( 13) and we show this next for completeness.Consider specifically the equations of motion of G ++ (t, t ′ ) and G +− (t, t ′ ) for t 0 ≤ t, t ′ < t f , which are obtained by plugging the 11 and 12 components of O(t, t ′ ), given in Eqs. ( 11) and ( 12), into Eq.( 13).On dropping the terms containing a δ-function at t f since they cancel out and further dropping the dissipative terms assuming they are of a form that does not affect the initial conditions, the G ++ (t, t ′ ) and Let us first set t ′ = t 0 + ϵ, for some small parameter ϵ, in Eq. ( 35) and integrate the equation over t from t 0 to t 0 + 2ϵ, taking care of factors of 1/2 that arise from evaluating the δ-function at a limit of the integral.Using the forms of G ±± (t, t ′ ) given in Eqs. ( 27), (28), and the following text in the resulting equation and further taking the limit ϵ → 0 gives Let us next integrate Eq. ( 36) over t from t 0 to t 0 + ϵ.
Using again the forms of G ±± (t, t ′ ) given in Eqs. ( 27), (28), and the following text in the resulting equation and further taking the limit ϵ → 0 gives for any t ′ .For the specific choice of t ′ = t 0 and using G < (t, t ′ ) = G > * (t, t ′ ), adding Eqs. ( 37) and (38) gives Now equating first the imaginary parts and then the real parts on both sides of this equation, and using the definitions of c xx and c xp from Eqs. ( 4) and ( 5), yields the initial conditions given in Eqs.(31) and (32).Similarly, setting t ′ = t 0 and subtracting Eq. ( 37) from Eq. (38) gives the commutator or Wronskian condition at the initial time given in Eq. (34).Lastly, let us differentiate Eq. ( 38) with respect to t ′ and set t ′ = t 0 + ϵ.Using the previous initial conditions to simplify the resulting equation, taking the limit ϵ → 0, and using the definition of c pp from Eq. ( 6) then gives the initial condition in Eq. (33).

C. Equations of motion
We next consider the equations of motion for G ±± (t, t ′ ) that are obtained as noted in the previous subsection.We drop the terms containing a δ-function at t 0 since these simply impose initial conditions on G ±± (t, t ′ ) and those containing a δ-function at t f since they cancel out as noted earlier.Then G ++ (t, t ′ ) and G +− (t, t ′ ) satisfy the following equations for t 0 < t, t ′ < t f , while G −− (t, t ′ ) and G −+ (t, t ′ ) satisfy the complex conjugates of the above equations.
To solve these equations, it is simplest to first decouple them by rotating to a new basis, where ξ(t) is a vector containing ξ ± (t).The differential operators and Green's functions in the two bases are then related by T , as found using Eqs.( 9) and ( 13), respectively, where G ξ (t, t ′ ) is a 2 × 2 matrix containing the functions G ξ,±± (t, t ′ ).It is instructive to write these functions explicitly, where the last identity follows from Eq. ( 29).G ξ,++ (t, t ′ ) is thus the symmetrized correlation and G ξ,+− (t, t ′ ) is the retarded Green's function of the theory.
We can now obtain the equations of motion for the functions G ξ,±± (t, t ′ ) in a similar way as we did for the functions G ±± (t, t ′ ).Since G ξ,−− (t, t ′ ) vanishes identically, its equation of motion yields the constraint where the subscripts I indicate imaginary parts as before.We take this to imply that γ 1I (t, t ′ ) = γ 2I (t, t ′ ), although one can envision a more general class of functions that satisfy this constraint.The equations for G ξ,++ (t, t ′ ) and G ξ,+− (t, t ′ ) then simplify to which can be solved for specific dissipation kernels.Note that only two real functions, γ 1R (t, t ′ ) + γ 2R (t, t ′ ) and γ 1I (t, t ′ ), contribute to dissipation.The solution to the second equation above gives us the retarded Green's function, G ξ,+− (t, t ′ ), of the theory.The first equation, on the other hand, yields the symmetrized correlation, G ξ,++ (t, t ′ ), and can be solved by first writing it as the sum of a homogeneous piece, that we denote G ξ,++[h] (t, t ′ ), and a particular solution or noise piece, that we denote G ξ,++[n] (t, t ′ ).While the homogeneous piece is the solution to Eq. ( 48) with zero on the right hand side and satisfies the initial conditions written in the previous subsection, the noise piece is obtained by convolving the source term on the right hand side of Eq. ( 48) with the retarded Green's function and vanishes at the initial time (i.e., at t ′ = t 0 ).We will choose a specific form for the dissipation kernels in the next section that will allow us to explicitly solve for G ξ,+− (t, t ′ ) and G ξ,++ (t, t ′ ) in that case.

IV. TIME-LOCAL DISSIPATION
In this section, we specialize to the simple and wellstudied case of time-local dissipation to understand whether the oscillator thermalizes.Time-local dissipation is obtained within the Caldeira-Leggett model, for example, when specializing to an Ohmic spectral density with infinite cutoff and the high-temperature regime [8,30].First, to reduce the integrals on the left hand side of Eqs. ( 48) and ( 49) to a simple damping term, we choose dissipation kernels such that γ 1R (t, t ′ ) + γ 2R (t, t ′ ) = γdδ(t − t ′ )/dt ′ , where γ is a (real) positive constant.Second, to reduce the integral on the right hand side of Eq. ( 48) to a local source term, we choose γ 1I (t, t ′ ) = (α/2)δ(t − t ′ ), where α is a (real) constant.Eqs. ( 48) and ( 49) then simplify to At the moment, we have not imposed any conditions on the constants γ and α, except that γ ≥ 0. We note that in the Caldeira-Leggett model with time-local dissipation, the dissipation kernels are given by γ 1R (t, , where β env is the inverse of the Boltzmann constant times the temperature of the thermal environment.We can, therefore, specialize to this model by leaving γ as is and replacing α with 4γ/β env in our results below. We solve the equations of motion ( 50) and ( 51) in the first subsection below.In the next two subsections, we obtain the purity of the oscillator and check whether the fluctuation-dissipation relation is satisfied in the latetime limit, and show that γ and α can in fact not be independent of one another.

A. Green's functions
Let us first solve Eq. ( 51) for the retarded Green's function of the theory, G ξ,+− (t, t ′ ).This is easily solved in Fourier space [31], so that G ξ,+− (t, t ′ ) is given by the following inverse Fourier transform, The integrand has two simple poles at ω ′ = −iγ ± ω 2 − γ 2 , both on the negative imaginary axis for ω > γ or ω < γ and assuming that γ > 0. Closing the contour from below gives where Γ = ω 2 − γ 2 , and we can check that this expression has the appropriate ω = γ and γ = 0 limits as well.
We next solve Eq. ( 50) for the symmetrized two-point correlation, G ξ,++ (t, t ′ ).Eq. ( 50) is a nonhomogeneous differential equation and, as mentioned earlier, its solution consists of a homogeneous piece G ξ,++[h] (t, t ′ ) and a particular solution or noise piece that we consider in turn.The homogeneous piece is the solution to the equation and, as also mentioned earlier, satisfies the initial conditions imposed on the full solution G ξ,++ (t, t ′ ).Now since G ξ,++ (t, t ′ ) is symmetric under the interchange of t and t ′ , it satisfies the same equation of motion in both time coordinates.Also noting that G ξ,++ (t, t ′ ) is a real function, the solution to Eq. ( 55) must be of the form where a is a complex constant, b is a real constant, and h(t) and h * (t) are solutions to the equation ḧ(t)+2γ ḣ(t)+ ω 2 h(t) = 0.
The constants a and b in Eq. ( 56) can be fixed by making use of the initial conditions in Eqs. ( 31), (32), and (33), which directly translate into initial conditions on G ξ,++ (t, t ′ ) and, therefore, on G ξ,++[h] (t, t ′ ).We thus have three equations in three constants, which yield a R , a I , and b in terms of c xx , c xp , and c pp and additionally h(t 0 ), ḣ(t 0 ), h * (t 0 ), and ḣ * (t 0 ) [32].We finally need a solution for h(t) and h * (t).The solution for h(t) can be written as where G 1 (t) and G 2 (t) are solutions to the same equation as that of h(t) with the initial conditions G 1 (t 0 ) = i, Ġ1 (t 0 ) = 0 and G 2 (t 0 ) = 0, Ġ2 (t 0 ) = −i, and the solution for h * (t) is similarly given by h * (t) = −iG 1 (t)h * (t 0 ) + iG 2 (t) ḣ * (t 0 ).Now plugging in the resulting expressions for the constants a R , a I , and b and the solutions for h(t) and h * (t) back into Eq.(56) gives us the final expression for Lastly, we write down expressions for the functions G 1 (t) and G 2 (t), To summarize, the homogeneous part of the symmetrized two-point correlation is obtained by plugging Eqs. ( 59) and (60) into Eq.( 58), which gives We see that G ξ,++[h] (t, t ′ ) is proportional to e −γ(t+t ′ −2t0) and thus vanishes in the limits that γ(t − t 0 ) ≫ 1 and γ(t ′ − t 0 ) ≫ 1 for any choice of initial conditions c xx , c xp , and c pp .
The noise contribution in Eq. ( 54) is simpler to calculate and, as mentioned earlier, obtained by convolving the source term on the right hand side of Eq. ( 50) with the retarded Green's function, We will choose t ≥ t ′ without loss of generality, in which case we find on using Eq. ( 53) for G ξ,+− (t, t ′ ), Note that this vanishes at the initial time (i.e., at t ′ = t 0 ), as mentioned earlier.We see that G ξ,++[n] (t, t ′ ) also has a piece proportional to e −γ(t+t ′ −2t0) which vanishes in the limits that γ(t − t 0 ) ≫ 1 and γ(t ′ − t 0 ) ≫ 1, but additionally has a piece proportional to e −γ(t−t ′ ) that need not vanish unless γ(t − t ′ ) ≫ 1 as well.This will be important for the fluctuation-dissipation relation that we discuss later in this section.
With the Green's functions in hand, we can also use Eq.(24) to write an expression for the one-point function ⟨ X(t)⟩.Substituting the expressions given after Eq. (10) for the shifted sources into Eq.( 24), we see that ⟨ X(t)⟩ can be expressed as which, in the absence of a source and using Eqs.( 59) and (60) becomes Note that this vanishes in the limit that γ(t − t 0 ) ≫ 1.

B. Purity
As discussed in Sec.II A, the purity of an oscillator in a Gaussian state can be written in terms of its two-point correlators.At the initial time t 0 , we found that Pu(t 0 ) = 0.5/ c xx c pp − c 2 xp .We can similarly write the purity at any time t in terms of two-point correlators at that time using with the correlators in turn obtained from G ξ,++ (t, t ′ ) of the previous subsection, We do not write explicit expressions for the two-point correlators, but they can be obtained using the results from the previous subsection.It is instructive, however, to consider the final expression for the purity, which we find to be In the late-time limit, in particular, the above expression reduces to lim The late-time purity is thus independent of the initial conditions c xx , c xp , and c pp , similar to the late-time correlations, and is in fact constant.Since purity must be between 0 and 1 at all times, we find a constraint on the three parameters: 0 ≤ 2γω/α ≤ 1, showing that the dissipation parameters γ and α are not independent of one another and additionally that α is positive.We also see that α cannot be set to zero as the late-time purity would otherwise diverge.In other words, the α → 0 limit of Eq. ( 70) reads which grows exponentially in time and is, therefore, not physical.The fact that α needs to be non-zero if γ is non-zero has to do with properties of the environment that is causing dissipation of the oscillator in the first place.
We can go even further and reconstruct the late-time density operator.Since our initial state is Gaussian and the dynamics are linear, we expect the density matrix ⟨x + |ρ(t)|x − ⟩ at any time t to be of the same form as the initial density matrix in Sec.II A, except that all correlations -x 0 , p 0 , c xx , c xp , and c pp -must be replaced by those at time t.Upon doing so using the expressions in the previous subsection (in the absence of a source), we find that the late-time density matrix is given by Let us compare this to the density operator for an isolated harmonic oscillator in a thermal state, ρβ (t) = e −β Ĥ , where β is the inverse of the Boltzmann constant times the temperature of the oscillator.In position basis, ρβ (t) is given by Eqs. ( 73) and (74) together suggest that the DHO thermalizes at the temperature β = (2/ω) tanh −1 (2γω/α).We can arrive at the same conclusion by directly comparing the late-time purity of the DHO in Eq. ( 71) with the purity of an oscillator in thermal equilibrium, given by Pu(t) = tanh(βω/2), as well.Note, however, that although the late-time density operator is thermal, this does not necessarily imply that the fluctuationdissipation relation is satisfied, and we discuss this further in the next subsection.
Before closing this subsection, we also plot in Fig. 1 the two-point correlators found from Eqs. (67), (68), and (69), on substituting for G ξ,++ (t, t ′ ) from the previous subsection, and the purity in Eq. ( 70), for three different choices of initial state and one set of dissipation parameters for illustration.In units of the oscillator frequency ω, we first consider a coherent initial state with c xx = 1/(2ω), c xp = 0, and c pp = ω/2, therefore Pu(t 0 ) = 1.We next consider a squeezed initial state with c xx = e −2r /(2ω), c xp = 0, c pp = e 2r ω/2, and r = 1, therefore Pu(t 0 ) = 1 again.And we lastly consider a thermal initial state with c xx = coth(βω/2)/(2ω), c xp = 0, c pp = coth(βω/2)ω/2, and β = 1/(2ω), therefore Pu(t 0 ) ≈ 0.24.We choose the dissipation parameters to be γ = ω/2 and α = 10ω 2 in all three cases and set ωt 0 = 0 for simplicity.We also show the result for an isolated harmonic oscillator at temperature β = (2/ω) tanh −1 (2γω/α) ≈ 0.2/ω for comparison, using Eq.(75) below to find the correlators and the expression from the previous paragraph for the purity.It is evident from the figure that the late-time behavior of the system is independent of the choice of initial state and equal-time results match those of an oscillator in thermal equilibrium.

C. Fluctuation-dissipation relation
The fluctuation-dissipation relation goes beyond demanding that the density operator be thermal as it is a statement about unequal-time correlators of the system.Let us first review the fluctuation-dissipation relation for an isolated harmonic oscillator in a thermal state, ρβ (t) = e −β Ĥ .In our notation (see Sec. III A), the two-point correlation G < β (t, t ′ ) is given by e βω e iω(t−t ′ ) + e −iω(t−t ′ ) 2ω(e βω − 1) , and G > β (t, t ′ ) = G < * β (t, t ′ ).Since Eq. (75) only depends on the difference of times, we will write G < β (t, t ′ ) and G > β (t, t ′ ) as G < β (T ) and G > β (T ) instead, where T = t − t ′ .Denoting the Fourier transforms of G < β (T ) and Let us now consider the Fourier transform of the antisymmetric two-point correlation G > β (T ) − G < β (T ), which is equal to twice the real part of the Fourier transform of G > β (T ) − G < β (T ) θ(T ).We define Similarly, we denote the Fourier transform of the symmetric two-point correlation with S β (ω ′ ), so that and it is given by G> β (ω ′ ) + G< β (ω ′ ) /2.Using the KMS relation, we can show that A β (ω ′ ) and S β (ω ′ ) satisfy known as the fluctuation-dissipation relation.We now want to check whether the quantum DHO that we studied in the previous subsections satisfies the fluctuation-dissipation relation in the late-time limit.Consider first A(ω ′ ) for the DHO, given by the real part of the Fourier transform of G ξ,+− (t, t ′ ) with respect to T = t − t ′ .Reading off the Fourier transform of G ξ,+− (t, t ′ ) from Eq. ( 52), we find that Next, consider S(ω ′ ), given by the Fourier transform of G ξ,++ (t, t ′ ).We noted in Sec.IV A that the homogeneous part of G ξ,++ (t, t ′ ) vanishes in the late-time limit, and therefore we only need to consider its noise part.
The Fourier transform of G ξ,++[n] (t, t ′ ) in the late-time limit can in turn be obtained from Eq. ( 62) by first replacing G ξ,+− (t, t ′′ ) and G ξ,+− (t ′ , t ′′ ) with their integral representation in Eq. ( 52), then setting the limits of the t ′′ integral to be t 0 = −∞ to t f = ∞, and then using ∞ −∞ dt ′′ e i(ω ′ +ω ′′ )t ′′ = 2πδ(ω ′ + ω ′′ ) to collapse one of the two frequency integrals.The resulting expression for S(ω ′ ) is given by and A(ω ′ ) and S(ω ′ ) for the DHO, therefore, satisfy Comparing Eqs. ( 78) and (81), we see that they match for multiple values of ω ′ only in the high-temperature regime where βω ′ ≪ 1 and coth (βω ′ /2) ≈ 2/(βω ′ ), and the temperature to which the oscillator thermalizes is given by β = 4γ/α.This agrees with the temperature found towards the end of the previous subsection, β = (2/ω) tanh −1 (2γω/α), in the high-temperature limit.It also matches the temperature of the environment for the Caldeira-Leggett model in the regime mentioned at the beginning of this section, where we noted that α = 4γ/β env in this model.
Eq. ( 81) can also be interpreted as a generalized fluctuation-dissipation relation for the quantum DHO, which when comparied with Eq. ( 78) gives an effective frequency-dependent temperature, [33][34][35].We see that the effective temperature agrees with the temperature β obtained earlier at the isolated frequency ω ′ = ω; the physical significance of this, however, is unclear to us.We emphasize that the generalized fluctuation-dissipation relation in Eq. ( 81) can be used to distinguish a quantum DHO from an isolated harmonic oscillator in a thermal state.Alternatively, the frequency spectrum or unequal-time measurements of the system can also be used to distinguish between the two cases even in the high-temperature regime.We show the unequal-time correlations explicitly in Fig. 2 for clarity, for the same choice of parameters as in Fig. 1 with a coherent initial state.
It is also worth noting that the late-time limit of G ξ,++ (t, t ′ ), which is simply the second term in Eq. (63), also agrees with the symmetric two-point correlation obtained in [11,[36][37][38][39] by first assuming that the fluctuation-dissipation relation holds for the quantum DHO, then finding the inverse Fourier transform of coth (βω ′ /2) ρ(ω ′ ), and finally taking the hightemperature and additionally weak-coupling (ω > γ) limit of the result.

V. TIME-NONLOCAL DISSIPATION
In this section, we briefly consider time-nonlocal dissipation with a specific choice of dissipation kernels, such that the fluctuation-dissipation relation is satisfied in the late-time limit at any temperature.We choose our dissipation kernels to match those in the Caldeira-Leggett model, specializing again to an Ohmic spectral den-sity with infinite cutoff but without restricting to the high-temperature regime [8,30].In this case, we have which simply evaluates to −(πγ/β 2 env ) cosech 2 [π(t − t ′ )/β env ].Using the techniques of the previous section, we can now solve for the corresponding anti-symmetric and symmetric two-point correlations of the time-nonlocal DHO.
G ξ,+− (t, t ′ ) is, in fact, the same as before, and is given by Eq. ( 53).Its Fourier transform with respect to t − t ′ is thus simply given by Eq. ( 79), which we write again for convenience, The homogeneous piece of G ξ,++ (t, t ′ ) is also unchanged, given by Eq. ( 61), and again vanishes in the late-time limit.The noise piece of G ξ,++ (t, t ′ ), on the other hand, is given by with γ 1I (t, t ′ ) given by Eq. ( 82).Since we are only interested in checking the fluctuation-dissipation relation in this section, we only need to consider the Fourier transform of G ξ,++ [n] (t, t ′ ) in the late-time limit.This can be obtained from Eq. ( 84) by first replacing G ξ,+− (t, t ′′ ), G ξ,+− (t ′ , t ′′′ ), and γ 1I (t, t ′ ) with their integral representations in Eqs. ( 52) and (82), then setting the limits of the t ′′ and t ′′′ integrals to be t 0 = −∞ to t f = ∞, and then using the δ-functions obtained by performing the time integrals to collapse two of the three frequency integrals.The resulting expression for the Fourier transform of G ξ,++ (t, t ′ ) in the late-time limit is given by and A(ω ′ ) and S(ω ′ ) for the time-nonlocal DHO, therefore, satisfy Since this matches Eq. ( 78), we conclude that the fluctuation-dissipation relation is satisfied for a timenonlocal DHO with the specific dissipation kernels chosen here, and it thermalizes to the temperature of the thermal environment.As in the time-local case, the frequency spectrum or unequal-time measurements of the system can, however, be used to distinguish a quantum DHO from an isolated harmonic oscillator in a thermal state.Lastly, we comment on how this relates to the time-local case of the previous section.In the hightemperature limit, we can expand the integrand in Eq. (82) using coth(x) ≈ x −1 + x/3 + O(x 3 ) as x → 0. Since the factors of ω ′ cancel in the zeroth-order piece, the kernel reduces to the time-local one with γ 1I (t, t ′ ) = (2γ/β env )δ(t − t ′ ).The next-to-leading order term is equally interesting since it results in a factor of ω ′2 under the integrand, which gives a −(γβ env /6)d 2 δ(t − t ′ )/dt ′2 contribution to the kernel.Note that such a term by itself will need to be handled carefully since we have assumed in our derivation that the dissipation kernels γ 1 (t, t ′ ) and γ 2 (t, t ′ ) are not proportional to d 2 δ(t − t ′ )/dt ′2 .Our derivation should be valid here, however, since this term arises from expanding the full dissipation kernel in Eq. ( 82).We can see that a −(γβ env /6)d 2 δ(t − t ′ )/dt ′2 contribution would also generate time-local dissipation, but with higher-order time derivatives appearing in the noise term.In the spirit of effective field theory, we can, therefore, replace the time-nonlocal dissipation kernel with a time-local one with successively higher-order time derivatives, in the high-temperature limit.Going back to the frequency domain, we see that the series expansion tracks all the way to Eq. ( 86) and the fluctuation-dissipation relation matches that in Eq. ( 81) with α = 4γ/β env to leading order, with higher-order corrections.

VI. BEYOND GAUSSIAN INITIAL STATES
In the previous sections, we folded in the initial conditions arising from the Gaussian density matrix of Eq. ( 3) with the dynamics of the quantum DHO from the outset, and solved for the relevant 2 × 2 Green's function G(t, t ′ ).We then solved for the n-point correlation functions and purity of the system from G(t, t ′ ) itself.In this section, we instead compute the time evolution operator of the density matrix for the quantum DHO, without first convolving it against the initial density operator.This approach not only allows us to cleanly separate the influence of dynamics from that of the initial conditions, but also results in expressions for equal-time n-point functions and purity that are valid for any choice of initial state.
As a start, we recall that the path integral -without dissipation - for an appropriate Lagrangian L, evolves any initial state |Ψ[t 0 ]⟩ (in the Schrödinger picture) forward in time, through the relation This in turn means that the density matrix of a closed system is related to its initial state via , (89) where the time evolution of ρ[t f ] would be the result of Eq. ( 88) acting upon the ket and the bra separately, namely, with the over-bar denoting complex conjugation.Note that the time evolution operator here factorizes into two separate ones.As already mentioned in Sec.II, however, if ρ describes a subsystem embedded in a larger environment, then this is in general no longer the case.Then instead, where the total action S T consists of three separate terms, That is, there is now a piece of the action, S I , that couples the variables x + and x − in the double path integrals as before, rendering the time evolution operator KK nonfactorizable.Its presence is again due to the 'tracing out' of irrelevant or inaccessible states from the full system's density operator.
In the subsections below, we first impose a general constraint on KK and next solve for it in the case of time-local dissipation.We then use the resulting solution to time-evolve the n-point correlations and purity of the system for any initial state.We emphasize that solving for KK rather than ρ[t f ] in this section allows us to obtain results for any initial state, whether Gaussian or non-Gaussian.

A. Probability conservation
As long as our subsystem can be embedded in a larger closed system, the total probability of finding it in some state has to remain unity at all times.And since KK arose from tracing over the irrelevant states of the larger closed system, tracing over the remaining states of the subsystem then corresponds to performing the trace over the entirety of the former.We must, therefore, have Now since the initial density operator is arbitrary, the time evolution operator must obey This condition for KK that preserves Tr [ρ[t f ]] = 1 can be contrasted with the condition Z diss.[J, J] = 1 discussed after Eq. ( 2).

B. DHO: Setup
Let us compute this time evolutionary operator KK for the quantum DHO, defined via the Lagrangians and This is, of course, the time-local case considered in Sec.IV, where ω, γ, and α here are the same as those in that section.That this is the quantum version of the classical DHO system will be explicitly verified when we obtain both the solution and equation of motion of the position operator's expectation value in Eqs. ( 148) and (149) below.Physically speaking, ω 2 − iα is the square of the oscillation frequency of the simple harmonic oscillator, where ω 2 and −α are, respectively, its real and imaginary parts.The strength of friction is controlled by the magnitude of γ.We have included the g term because it is allowed at the quadratic level; in fact, it turns out to be mandatory.Since ω appears only as a squared quantity, we will take it to be non-negative without loss of generality.We will also see that in order to prevent a runaway solution to the quantum statistical expectation value of the position, the friction γ also needs to non-negative.Moreover, to ensure probability conservation in Eq. (98), we will find that g = −α.Finally, for non-zero γ and ω, α needs to be positive to produce a well-defined density operator in the asymptotic future.To summarize, we will find that ω > 0 , (g = −α) < 0 , and γ > 0 . (101) As before, we will remain agnostic to how the model in Eqs. ( 99) and (100) arises, but simply assume that it is a consistent effective theory with time-independent parameters.For technical convenience, we now perform the change-of-variables [40] under which the total action becomes on using Eqs.( 99) and (100) for L and L I .

C. Evaluation of the DHO KK
The double path integrals occurring in the timeevolution operator KK may be evaluated in a similar manner to their single path integral counterpart K. First, we seek the classical solutions ξ ± [t] = ξ ± c [t] that extremize the total action.Then, we perform a change in path integration variables from ξ ± to ξ ± q by expanding around the classical path, namely, ξ ± = ξ ± c + ξ ± q .We will find that the double path integrals over ξ ± q may be determined by probability conservation.
By varying Eq. ( 103) with respect to ξ ± , we find that the classical solutions ξ ± c must obey where the indices I and J run over ± and the 2 × 2 differential operator reads The x + portion of the KK path integrals in Eq. ( 92) propagates the quantum system from (t 0 , x + 0 ) to (t f , x + f ), whereas the x − portion propagates the same system from (t 0 , x − 0 ) to (t f , x − f ).This motivates us to impose the following boundary conditions on the classical solutions, keeping in mind that The solutions ξ ± c may be expressed in terms of the retarded Green's function, which we express here as a function of a single time variable τ , of the matrix differential operator D I (τ ) J in Eq. ( 104), When α = −g, Eq. ( 108) is intimately related to Eqs. ( 50) and (51).Unlike the G ξ,++ and G ξ,+− there, however, we will solve for G J (R) M with retarded boundary conditions and use it to construct the classical trajectories ξ ± c .The 2×2 retarded Green's function admits the integral representation The integrand here is the matrix inverse of the D in Eq. ( 104), but expressed in frequency ω ′ -space.The contour C skirts all four poles of the matrix integrand on the complex ω ′ plane from below so that when, and only when, τ > 0 the result is non-zero.Upon summing over the residues, we find with the relations A direct calculation further reveals that In terms of G, the homogeneous solution portion of the classical trajectories is That Eq. ( 119) solves its equation of motion in (104) is because of Eq. ( 116), and that it obeys the appropriate boundary conditions in Eq. ( 106) is because of Eqs. ( 117) and ( 118).We will also need the first derivatives of these trajectories at the boundary times.Invoking Eq. ( 118) hands us while ξJ c [t f ] is simply the t-derivative of Eq. ( 119) with t replaced with t f .Note that even though the total action involves the time integral over t ∈ [t 0 , t f ], it may be converted into a difference of boundary terms when it is evaluated on the classical solutions ξ ± c .For, upon integration-by-parts, the action in Eq. ( 103) is transformed into where the integral terms vanish because of Eq. ( 104) and we used the boundary conditions in Eq. ( 106).We now explicitly shift the path integration variables using Since the boundary conditions for the path integration variables, ξ ± [t f ] = x + f ± x − f and ξ ± [t 0 ] = x + 0 ± x − 0 , are already accounted for by the ξ ± c in Eq. ( 106), we must impose which in turn implies that Eq. ( 103) now reads where S T [ξ ± c ] is the total action evaluated solely on ξ ± c and given by Eq. ( 122) while S T [ξ ± q ] is evaluated solely on ξ ± q but with the limits implied by Eq. ( 124).There are no cross terms between ξ ± c and ξ ± q as the action is quadratic and the cross terms are, therefore, themselves necessarily linear in ξ q .Since the action itself is extremized by the solution to Eq. ( 104), these linear-in-ξ q terms must vanish upon integrating-by-parts all the derivatives acting on ξ ± q and employing Eq. (124) to set to zero the associated boundary terms.
With the new form of the total action in Eq. ( 126) and keeping in mind the integration limits in Eq. (125), we can now write the time evolution operator in Eq. (92) as Let us now consider again the trace of ρ[t f ], which involves setting x + f = x − f ≡ x in Eq. ( 127), followed by integrating over all real x.This procedure will, however, yield an S T [ξ ± c ] that is quadratic in x, because of the 122).More explicitly, a direct calculation tells us that these quadratic-in-x terms in S T [ξ ± c ] are If we were to substitute the above into the integral on the left hand side of the probability conservation statement in Eq. ( 98), we would end up with a Gaussian integral and would not obtain the required δ-function on the right hand side.We must, therefore, eliminate these x 2 terms completely, which can be done by setting as indicated by the overall factor of (g + α) in Eq. ( 128).This, in fact, recovers the time-local version of the γ 1I (t, t ′ ) = γ 2I (t, t ′ ) condition that we deduced from Eq. (47).Notice also that we had to perform an explicit calculation of the density operator's time evolution op-erator KK in this section to obtain the condition in Eq. ( 129), whereas it arose directly from the Green's function equation in Sec.III.
With the condition in Eq. ( 129) and the definitions -the T here not to be confused with the T ≡ t−t ′ of Sec.IV -we may now record the following.The homogeneous solution portion of the 2 × 2 Green's function and its inverse, written as a function of a single time variable τ , are where we have identified the Green's function of the classical DHO as and its anti-damped counterpart G −γ [τ ] and G −γ [τ ] given by the same expressions, except with γ replaced by −γ.
The DHO Green's function obeys and the anti-damped one again obeys the same equations, except with γ replaced by −γ.Let us also observe that G + + in Eq. ( 131) is, up to a factor of −i, the same as G ξ,+− in Eq. (53).
The action evaluated on the classical trajectory now reads and setting x + f = x − f ≡ x in Eq. ( 137) gives The probability conservation of Eq. ( 98) applied to the time evolution operator in Eq. ( 127) further yields Employing Eq. ( 138) thus hands us the relation As alluded to earlier, probability conservation allows us to evaluate the remaining double path integrals involving ξ ± q , namely, KK[t f , 0, 0; t 0 , 0, 0].We may, at this point, conclude that the time evolution operator for the DHO system in Eq. ( 127) with g = −α is given by with the classical action specified in Eq. ( 137).Furthermore, armed with this KK for the DHO, we may immediately integrate it against a Gaussian initial density matrix ⟨x + |ρ(t 0 )| x − ⟩ = N exp {iS 0 [x + , x − ]}, with S 0 given in Eq. (3) -i.e., compute Eq. (89) -to discover that the final density matrix takes the same form as the initial one, except that all relevant 'initial' parameters are replaced with the time-dependent final ones.For instance, the initial position and momentum are replaced as and the parameters A and B, related to the initial twopoint correlations c xx , c xp , and c pp by Eqs. ( 4), ( 5), and ( 6), are similarly replaced by their counterparts at t f .

D. DHO 'equal-time' n-point functions
We now turn to calculating the quantum statistical average of a product of n ≥ 1 position operators at time t f .An application of Eqs. ( 89) and (141) gives us To proceed, we consider the generating functional The time-independent J here is similar in spirit to, but not to be confused with, the time-dependent J ± introduced in the previous sections.From Eq. ( 145), the npoint function in Eq. ( 144) may be extracted via or, equivalently, by reading off the coefficient in the Taylor series of Z[J] containing exactly n powers of J. Now, employing Eq. (138) in Eq. ( 145), we would discover that x + 0 −x − 0 = G γ [T ]J is enforced by the δ-functions obtained by integrating over x.The generating functional is then found to be We note, parenthetically, that the x integral in Eq. ( 147) bears a resemblance to the Wigner function in quantum mechanics.Moreover, we have written iJ as a group in the above expressions since n-point correlations are generated by differentiating Z[J] with respect to it.Applying Eq. ( 146) to Eq. ( 147), we first find the onepoint function to not only depend on the initial position and momentum one-point functions, which is, of course, Eq. ( 65) with t replaced with t f , but also to obey the classical DHO equation, since the pair e −γT sin[ΓT ] and e −γT cos[ΓT ] do.Eq. ( 148) is thus the quantum analog of the initial value formulation of the classical DHO differential equation.
Notice that it does not depend on the additional parameter α, the imaginary portion of the DHO's frequencysquared, occurring in our dynamics defined in Eqs. ( 99) and (100) (with g = −α).This suggests that α is tied to the underlying quantum and/or statistical aspect(s) of the DHO system at hand.As in Sec.IV, we choose γ ≥ 0 in order for the e −γT factors in Eq. (148) to not describe runaway trajectories.Also, X[t 0 ] is readily recovered by setting t f = t 0 in Eq. (148), while taking one time derivative before setting t f = t 0 gives us P [t 0 ] , In fact, a direct calculation shows that the general momentum one-point function, is exactly the time derivative of the position one-point function in Eq. (148), i.e.,

P [t
again suggesting that P (t) = Ẋ(t) for the problem we consider.
Turning next to the two-point function, application of Eq. (146) returns This recovers Eq. ( 67) with G ξ,++ (t, t) there given by the sum of Eqs. ( 61) and (63) with t ′ = t, except that the initial two-point correlations here are arbitrary.We also note that the α term is independent of the initial conditions -a feature already present in Eq. (63) -whereas the remaining terms depend linearly on the initial twopoint functions involving both position and momentum operators.In fact, in the asymptotic future, T → ∞, the exp(−2γT ) in G γ [T ] 2 -recall Eq. ( 134) -implies that the entire dependence on the initial two-point functions is damped out and all that remains is the α-term, It is worth reiterating that although our dynamics as encoded within Eq. ( 103) are purely quadratic -i.e., with linear equations of motion -we have not made any assumptions regarding the initial density operator ρ[t 0 ] in Eq. ( 146) and, therefore, regarding the initial n-point functions in Eqs. ( 148) and (154).
E. Purity of the DHO system and the α → 0 limit We now turn to computing the purity of the DHO system without assumptions on the initial state.Using the time evolution in Eq. ( 89) and the form of KK for the DHO in Eq. (141), the purity at time t f is given by Pu where the time evolution operator associated with purity is Let us now take the α → 0 limit of Pu[t f ]'s time evolution operator KKKK.Upon switching to the ξ ± ≡ x + ± x − variables, the exponent is found to be linear in both ξ ± , The integral representation of the Dirac δ-function together with the identity δ[az] = δ[z]/|a|, which holds for integration variable z and constant a, tells us that KKKK is proportional to where we have used that x 158), therefore, simplifies to Plugging Eq. ( 160) back into Eq.( 156) informs us that, in the α = 0 case, the purity at time t f is the initial purity multiplied by an exponential growth factor, ) We conclude that the quantum DHO system cannot have a purely real frequency-squared, i.e., α ̸ = 0 as long as γ > 0. This is because purity is constrained to lie within [0, 1], while the e 2γT factor in Eq. ( 161) indicates that the α = 0 DHO system has an infinite purity in the asymptotic future.This was already noted in Eq. ( 72), but once again we see that it holds for arbitrary initial conditions.

VII. DISCUSSION
The quantum DHO is a prototypical example of an exactly solvable dissipative quantum system.We revisited this simple system in this paper, with a goal of developing an effective field theory-inspired and influence functional-based method to describe its dynamics.We were particularly interested in what constraints one can impose on various terms that can appear in the influence functional.We thus considered a quadratic influence functional without describing the system-environment interaction that led us to it.We first found that the dissipation kernels we introduced are related in such a way that only two real functions are needed to describe the system's dynamics.We next restricted to time-local dissipation and solved for the exact Green's functions of the system in this case.We used the resulting Green's functions to show that late-time correlations and purity are independent of the initial conditions for both, Gaussian and more general initial states, and argued that the two dissipation kernels must in fact be constrained in order to obtain a physically meaningful purity in the late-time limit.We briefly considered time-nonlocal dissipation as well, and showed that the fluctuation-dissipation relation is satisfied for a specific choice of dissipation kernels.
It is worth highlighting that we discussed two different approaches to solve the problem -in Secs.II-V, we focused on the in-in formalism, and in Sec.VI, we focused on a double in-out formalism.We found that each approach has its own advantages.For example, the former allowed the calculation of unequal-time correlations, which are crucial to understand whether the fluctuationdissipation relation holds at late times.The latter, on the other hand, allowed us to go beyond Gaussian initial states and obtain results that are non-perturbative in the initial state.Which approach is "better" to use thus depends on the problem at hand, but it may be interesting to understand how to extend each approach.
We would also like to highlight a few other extensions of the methods developed in this paper that could be useful.First, we restricted our calculation to dissipation kernels γ 1 (t, t ′ ) and γ 2 (t, t ′ ) that are not proportional to d 2 δ(t − t ′ )/dt ′2 .We further specialized first to time-local dissipation with γ 1R (t, t ′ ) + γ 2R (t, t ′ ) = γdδ(t − t ′ )/dt ′ and γ 1I (t, t ′ ) = (α/2)δ(t−t ′ ), with γ and α constant, and then to time-nonlocal dissipation with the specific choice of γ 1I (t, t ′ ) given in Eq. (82).It would be interesting to lift these restrictions to consider more general influence functionals -time-dependent, generally nonlocal, and arising from derivative system-environment interactions.Second, it would be interesting to adapt these methods to describe a dissipative QFT in Minkowski spacetime as well as time-dependent spacetimes such as de Sitter.Such an effective field theory-inspired approach would again complement existing work on open QFTs.And third, it would be interesting to study loop corrections in dissipative QFTs in the presence of interactions and/or nonlinear dissipation, replacing the usual Green's functions in loop corrections with their dissipative analogs.This may, for example, help us to better understand the secular growth that is often found in loop corrections to observables in time-dependent quantum systems and QFTs; see [41][42][43][44] for reviews and [45][46][47][48][49][50][51][52][53][54] for related work and other possible resolutions.
write the initial density matrix element as ρ[x + , x − , t 0 ], with the functions x ± (t) evaluated at t 0 .The generating functional then becomes as given in Eq. ( 1).The δ-function at the turn-around time t f also imposes the constraint ẋ+ (t f ) = ẋ− (t f ) that is needed to cancel the boundary terms at t f , as discussed in Sec.III.This was shown in [55] and we reproduce their argument here for completeness.Following [55], let us first discretize the x + path integral and action as Dx + = lim δt→0 dx + f dx + t f −δt . . .dx + 0 and where the summation runs over all t between t 0 + δt and t f in steps of δt, and similarly discretize the x − path integral and action.Using the δ-function in Eq. (A4), we can also write f as a single integral that we denote R dx f , replacing x ± f everywhere with x f .Let us now specifically consider the x f integral.The terms proportional to ω 2 (and multiplied by x 2 f ) cancel between S[x + ] and S[x − ].We can also set J ± (t f ) = 0 since t f is chosen to be later than any times at which we are interested in calculating correlation functions.Then the x f integral in the discretized version of the generating functional in Eq. (A4) becomes which is simply 2πδ . This can be further simplified to δ ẋ+ f − ẋ− f by Taylor expanding x ± t f −δt = x f − ẋ± f δt.The path integral, therefore, imposes the constraint ẋ+ (t f ) = ẋ− (t f ).
Lastly, we note that the constraint ẋ+ (t f ) = ẋ− (t f ) may not necessarily hold in the presence of dissipation.Specifically, if the dissipation kernels γ 1 (t, t ′ ) and γ 2 (t, t ′ ) that we introduced in Eq. ( 7) are proportional to d 2 δ(t − t ′ )/dt ′2 , then the argument in the previous paragraph needs to be revisited as Eq.(A6) could potentially contain other terms.If, on the other hand, the dissipation kernels are proportional to dδ(t − t ′ )/dt ′ , as considered in the main text, then the constraint ẋ+ (t f ) = ẋ− (t f ) continues to hold.
denoting shifted sources, and the superscript T indicates a transpose.O(t, t ′ ) is a 2 × 2 matrix of differential operators whose 11 and 12 components are given by The other two components of O(t, t ′ ) are related to these by O 21 (t, t ′ ) = −O 12 * (t, t ′ ) and O 22 (t, t ′ ) = −O 11 * (t, t ′ ).