An operator derivation of the Feynman-Vernon theory, with applications to the generating function of bath energy changes and to anharmonic baths

We present a derivation of the Feynman-Vernon approach to open quantum systems in the language of super-operators. We show that this gives a new and more direct derivation of the generating function of energy changes in a bath, or baths. This generating function is given by a Feynman-Vernon-like influence functional, with only time shifts in some of the kernels. We further show that the approach can be extended to anharmonic baths by an expansion in cumulants. Every non-zero cumulant of certain environment correlation functions thus gives a kernel in a higher-order term in the Feynman-Vernon action.


I. INTRODUCTION
When a quantum system interacts with an environment (open quantum system or OQS), the state of the system is influenced by the environment in a fundamental manner. Decoherence, for example, is a consequence of quantum entanglement between the system and the environment [1]. The understanding of effects on the system induced by the environments is essential to the quantum information technology [2] and quantum thermodynamics [3]. A direct inclusion of the environment in a first principle investigation is usually not practically feasible. On the other hand, the structure and details of a large environment can only partially be reflected in the dynamics of the system. In the super-operator approach going back to Nakajima and Zwanzig [4,5] one starts from equations of motion (von Neumann-Liouville equations) of the total density matrix of the system and the environment, and projects that to an effective dynamics for the system density matrix [6][7][8]. In the alternative approach of Feynman and Vernon the influence of an environment on the system is expressed in terms of a influence functional [9]. The development of the reduced density matrix of the system is then given by a double path integral, where the influence functional couples the two paths. If used exactly, both approaches mathematically specify completely positive dynamic maps describing the evolution of the system. Once these are obtained, the system dynamics can be investigated without the knowledge of environment dynamics [7,8].
Even after one has found a closed form expression of a dynamical map or the influence functional, calculating the time-evolution of the system using them is technically challenging. Tractable methods using various approximations have been developed. Quantum master equations (Lindblad equations) [10,11] based on the Born-Markovian approximation are the most popular, and can be derived in both approaches [7,12,13]. Non-Markovian methods are also developed [14][15][16]. For the "spin-boson" problem of one two-state system (a qubit) interacting with a bosonic bath the noninteracting blip approximation (NIBA) was developed [17], and shown to be equivalent to relaxation after a polaron transform [18,19]. Many exact numerical algorithms have also been developed, such as the hierarchical equation of motion (HEOM) [15,20], the quasi-adiabatic propagator path integral (QuAPI) [21], and the Stochastic Liouvillian algorithm [22].
The thermodynamics of an OQS describes how a quantum system exchanges energy, particles and other quantities with one or several reservoirs. While in Quantum Markov dynamics the energy interchanged with a reservoir, which we call heat, can be expressed in terms of Lindblad operators acting on the system density matrix [23], in general that is not so. Nevertheless, it has recently been shown by several groups that the generating function of heat can be computed by the path integral technique in a Feynman-Vernon-like approach [24][25][26][27]. In this formulation appears a new influence functional depending on the generating function parameters. However, the terms in this new influence functional are in fact the same as in Feynman-Vernon, with only a time shift in the argument in some of the kernels [28]. We will in this work derive the same result in the super-operator formalism where it emerges in a straight-forward manner without cancellations in intermediate steps of the calculation.
We will also show that the super-operator approach can be extended beyond ideal Bose gas. Cumulants of specific bath correlation functions (to be discussed below) then enter as kernels in higher-order order terms in the Feynman-Vernon action. For instance, a non-zero third-order bath correlation function gives the kernel of a third-order term in the Feynman-Vernon action. This result appears more difficult to obtain in the path integral formulation.
Several of the results in this paper can be found in the literature, but not in the context of current concerns in quantum thermodynamics, and/or not in a form easily assimilated in different communities. The super-operator expression for evolution operator of the reduced density matrix of the system (Eq. (9) below) is for instance given (in the Schrödinger picture) as Eq. (3.508) on page 187 in the monograph of Breuer and Petruccione [7]. Two recent contributions that use a similar plus/minus (left/right) representation of the super-operator as we do are [29] and [30]; the latter paper also extends the analysis beyond harmonic baths. Time shifts in kernels describing a statistics of heat appear in the theory of heat transport through a Josephson junction developed in [31], though in a particular setting, and for a partially classical model.
We present the equivalence of super-operator and Feynman-Vernon approach in a coherent manner through the generating function of heat as an example, and the extension to anharmonic baths. The paper is organized as follows. Section II contains a brief summary of the path integral and super-operator approaches to the development of the system density matrix, and Section III extends the discussion to generating function of heat. Section IV contains a systematic and general derivation of the super-operator approach, and Section V the extension to anharmonic baths. Section VI sums up and discussed the results. Appendix A defines the super-operator time-ordering used in the main body of the paper, and Appendix B contains the details of the pair correlation functions in harmonic baths. Appendix C shows the third-and fourth-order kernels in high-order influence functional that result from non-zero third-order and fourth-order cumulants in the bath.

II. THEORY OF OPEN QUANTUM SYSTEMS: A BRIEF SUMMARY
Consider a system S in contact with an environment B. Their Hamiltonians are denoted as H s and H b , respectively and they are coupled through a interaction Hamiltonian V sb . The whole system is assumed to be initially in a product state ρ(t i ) = ρ s (t i ) ⊗ ρ b (t i ) and evolves by a unitary transformation: where ρ is the density operator of the whole system and U (t; t i ) = e −iH(t−ti) is a usual time-evolution operator with the total Hamiltonian H = H s + H b + V sb . Units of time are chosen such that = 1. The state of the system is defined as a reduced density ρ s = tr b ρ(t) where tr b traces out the degrees of freedom of the environments. The theory of OQS seeks a completely positive operation M(t; t i ) (dynamical map, or quantum map) defined by [11] ρ Alternatively, a quantum map can be considered as the given, and then one of the many possible couplings to an environment giving rise to same map after tracing out the environment is called an environmental representation of that map [32].
In either case, once the map is obtained, we can evaluate any quantity associated with the system. For example the transition probability from an initial pure state of the system |i to a final pure state of the system |f is given by In the path integral approach the dynamical map (2) is written in coordinate basis as where Q andQ are coordinates of the system at the initial time t i and final time t f . Feynman and Vernon wrote M in a path integral form where S[Q] is the classical action of a system trajectory without interactions with the environment, and Q(t) and Q(t) are forward and time-reversed trajectories. The effects of the environment are fully included in the influence functional F [Q, Q ; t f , t i ] [9,33]. General properties of influence functionals were discussed in [9], and we will return to those below. Exact expressions can be obtained when the functional integrals over the environment variables can be done in closed form. In practice this means that integrals have to be Gaussian, which is the case when the environment is a Bose gas with Hamiltonian H b = k p 2 k /2m k + m k ω 2 k q 2 k /2 initially in Gibbs states ρ b (t i ) = e −βHb /Z b , and the coupling takes a bilinear form V sb = X s ⊗ ( x c k q k ) where X s is an operator of the system, and c k is the coupling strength. Such an environment is called a harmonic bath at inverse temperature β. We will later discuss the case when the environment consists of two or more harmonic baths, each with its own temperature. As already discussed in [9], the coupling constants c k may depend on time. While formally that only leads to a redefinition of the system operator X s , in practice it is convenient to assume that c k vanishes at the very beginning and at the very end of the process [34]; as otherwise boundary terms lead to antinomies, well-known in the literature since some time [35,36].
Under the above assumptions the influence functional can be written in exponential form where Φ is known as the influence action or Feynman-Vernon action, and κ r (τ ) and κ i (τ ) are known as dissipation and noise kernel. These functions are In the super-operator approach one instead directly evaluates Eq. (1). Rewriting Eq. (1) using the Liouville operator L(t)• = −i[V sb (t), •] in the interaction picture, the map (in the interaction picture) can be written as where ← − T is time-ordering super-operator which chronologically orders the super-operators (see Appendix A.) The symbol • in above and in the following is the "slot" on which the super-operator acts, and represents any operator, including density operator. Using Wick's theorem, we find the map The two super-operators together can be expressed with commutator and anti-commutator as ( The kernels κ r and κ i in (9) and (6) are the same as those in Eq. (6), and will be shown to be the equilibrium pair coorelation functions of the ideal Bose gas.
While the two methods use different mathematical objects, one with paths Q(t) andQ(t), and the other with super-operators X + s (t) and X − s (t), Eqs. (6) and (9) clearly show similarity. They are essentially the same if two quantities are replaced as . In both approaches, the degrees of freedom of the environment are all integrated out analytically. However, the path integral for the system coordinate still remains in Eq. (5). For the operator approach, we must apply the super-operator ← − T e iΦ(t f ,ti) to the initial density. Neither operation is trivial but certain approximations or numerical methods have been developed [15] Extension of these methods to a system interacting with multiple environments is straight-forward if the environments do not interact between themselves. Indeed, General property of influence functionals 2 of Feynman and Vernon states that "If a number of [environments] act on [the system] and if F k is the influence of the k'th [environment] alone, then the total influence of all [the environments] is given by the product of the individual influences" [9]. In the super-operator approach the same statement follows from the observation that if the Lioville operator is a sum, say •], and if the environment operators in V sh (t) and V sc (t) commute and act on parts of the environment that start in a product state (different baths), then the time ordering of environment operators in (8) can be done separately.

III. GENERATING FUNCTION OF HEAT
Once the dynamical map is found, we know the state of the system precisely. However, the information on the state of environments is completely buried in the map. If we want to investigate any quantity associated with the environments or correlation between the system and environments, the knowledge of the system density alone is not enough. In order to make our story concrete, we consider a system interacting with a hot and a cold bath. Their Hamiltonians are denoted as H s , H h , and H c , respectively, and the interaction Hamiltonians between the system and the baths are V sh and V sc . The initial state of the whole system is assumed to be a product state and the baths are at thermal equilibrium where Z is a partition function. E (n) and |E (n) are eigenvalue and the corresponding eigenket of H . The system is initially in an arbitrary state ρ s (t i ) = i µ i |i i| where µ i and |i are eigenvalues and eigenkets of the density. Now we want know the change in the energy of the cold bath, ∆E c , over time period t f − t i . The probability distribution of ∆E c may be written as where µ i and |i are the eigenvalue and eigenket of ρ s (t i ) and the final states |f can be any basis set. The generating function of heat is the Fourier transform of parameter ν of the probability distribution P (∆E c ; t f ) with respect to variable ∆E c . One finds is an operator in the Hilbert space of the system. Direct comparison of Eqs (2) and (14) shows that Γ(ν, t) is quite similar to ρ s (t). In fact, when ν = 0, they coincide. The only difference is that one of the time evolution operators in Eq. (14) is rotated by e iνHc . The resemblance suggests that the generating function can be computed with the methods developed for OQS. Following the procedure discussed in the previous section, we first write Γ with a map as where We generally assume that the system interacts separately with the hot and the cold bath, and thus the Liouville operator is split to two parts, L x = −i[V sx (t), ·], x = c, h. If the system parts of the two operators (in the interaction picture) always commute the expression further factorizes into the product of the traces over the two baths separately. The energy change associated with a particular transition from |i to |f can be expressed like the transition probability (3): In order to use the path integral approach, we express Eq. (15) in the coordinate representation in the same way as Eq. (4), The map has been derived using the path integral [28] for When ν = 0, Eq. (16)  Based on the correspondence between the path integral method and the super operator method, we expect that the map defined in Eq. (16) is given by M(ν; t f , t i ) = e iΦh(t f ,ti)+iΦc(ν;t f ,ti) with Φ h (t f ; t i ) remains the same as Eq. (9). In the next section, we derive Eq. (15), and show that it appears more directly in the super-operator approach.

IV. SUPER-OPERATOR APPROACH FOR THE GENERATING FUNCTION
We derive Eq. (15) by evaluating the super-operator expression of map For simplicity we have assumed that the system parts of the interaction Hamiltonians commute, which allow us to focus on the trace over the cold bath.
First we rewrite H b with creation and annihilation operators, a † k and a k : and the interaction Hamiltonian where c k = c k / √ 2m k ω k is coupling strength. The system part of the coupling X s is arbitrary. Using the interaction picture, the Liouville super-operator is defined by where X s (t) = e iHs(t−ti) X s e −iHs(t−ti) and e iHb(t−ti) Y −iHb(t−ti) b and for mathematical convenience, we introduced the following super-operators Expanding the exponential function in Eq. (IV) ×C d1,··· ,dn (ν; t 1 , · · · , t n ) where multi-time correlation functions of the environment are defined as For the cross correlation C +− (t1, t2), only t1 on the chronological branch shifts and thus the time difference also shifts. The situation for C −+ (t1, t2) is similar to C +− (t1, t2) except that t2 shifts by ν instead of t1. For t1 < t2, the direction of τ is reversed for C ++ (t1, t2) and C −− (t1, t2). However, the sign of t1 −t2 does not affect the situation of C +− (t1, t2) and C −+ (t1, t2).
where · · · ti indicates expectation value tr c {· · · ρ c (t i )}. Since H c is quadratic in a and a † , all odd order correlation functions vanish. For the even order terms, we apply the Wick's theorem for operators e iνHc Y d1 c (t)e −iνHc C d1,···d2n (t 1 , · · · , t 2n ) = all possible pairing all pairs where all pairing indicates the sum of all possible combinations of pairs. The map is now expressed with the pair correlation functions as The four pair correlation functions C ++ , C +− , C −+ , and C −− can be expressed with the standard pair correlation function κ(τ ) as shown in Fig. 1. (See Appendix B.) The correlation functions between the two times on the same branch are where τ = t 1 − t 2 . The cross correlation functions are Substituting these correlation functions into (12) we obtain Eq. (15). Note that only the difference between the map for ρ derived and discussed in Section II and the map for the generating function is the cross correlations.

V. ANHARMONIC BATHS AND CLUSTER EXPANSIONS
A second advantage of the super-operator formulation is in the derivation of corrections to the Feynman-Vernon theory. The starting point is then the dynamical map (14) with the multi-time correlation functions of the environment (13), but without assuming Wick's theorem. The outcome will be that multi-time cumulants of the environment (discussed below) translate into kernels of higher-than-quadratic contributions to the Feynman-Vernon action.
For ordinary operator correlation functions, successive orders of cumulants are defined inductively as Owing to the time-ordering super-operator ← − T and indexes d j , C d1,··· ,dn (t 1 , · · · , t n ) defined in Eq. (13) behaves like an ordinary multi-time correlation function and the relations (9) hold. Hence, C d1,···d N (t 1 , · · · , t N ) = all possible groupings groups of one time 10) where N can be even or odd. The first order cumulant (G d 1 ) can be set to zero by a shift. The first non-trivial cumulant is then where we have retained the super-operator notation on the right-hand side. For a bath that satisfies Wick's theorem, this cumulant and all others beyond G d1,d2 2 vanish. The second step is to count the number of groupings in (10) with n 1 groups of one element, n 2 groups of two elements (pairs), n 3 groups of three elements, etc. There are N ! n 1 !n 2 !(2!) n2 n 3 !(3!) n3 · · · such groupings. The correlation functions appear inside the time integral and index sums in (14) and the indices and time variables can therefore be renamed in any way. Each grouping of the same type (same n 1 , n 2 , . . .) hence contributes the same, and the quantum map can be summed in an analogous way to Section IV.
Introducing Q(t) as the coordinate representation of X + s (t) and −Q(t) the coordinate representation of X − s (t) one can show that the contribution to the Feynman-Vernon action from n number of X and m number Y is where the last term is the cumulant of the operator correlation function with the times ordered as required in the super-operator cumulant. One can further sum all contributions of the same order and express them in terms of timeordered sums ζ + (t) = Q(t) +Q(t) and differences ζ − (t) = Q(t) −Q(t). The most important general result one can find this way is for the largest time, the dependence in only through the difference ζ − (t) as also follows from Feynman and Vernon's General property of influence functionals 5 [9]. Ultimately this is a consequence of the super-operator correlation function C d1,··· ,dn (t 1 , · · · , t n ) being independent of the symbol connected to the largest time. Otherwise the general expressions are somewhat unwieldy, and we will here only quote the result to third order where A, B, C and D are combinations of third order bath correlation functions given in Appendix C.

VI. DISCUSSION
In this paper we have compared the path integral and super-operator approaches to the theory of open quantum system (OQS). We have pointed out that both approaches lead to equivalent descriptions of a system interacting with one or several harmonic oscillator baths, but that the routes to the result are qualitatively different. In the super-operator approach the kernels in the description are found to be certain pair correlation functions of the bath (or baths), and the main assumption is Wick's theorem, reducing any correlation function to sums or products of pair correlation functions. In the path integral approach, the result on the hand follow from integrating over the initial and final points of the propagator of an harmonic oscillator (one of the degrees of freedom of the bath) acted upon by a linear drive (a linear interaction with the system), and after a fair amount of cancellation.
We have here shown that same holds for the generating function of heat: both approaches give the same result, but the super-operator approach is more direct. In particular, the fact that the generating function of heat can be expressed with the same kernels as for the system density matrix (Feynman-Vernon theory), with only a time shift in the terms mixing the forward and time-revered paths, follows in a more straight-forward manner in the super-operator approach.
We have also shown that the super-operator approach extends in a natural way to interactions with environments where Wick's theorem does not hold. Cumulants of correlation functions of the environment, which vanish when Wick's theorem holds, hence translate to kernels in higher-order terms in the Feynman-Vernon action. In the text we have discussed that the resulting higher-order theory of the influence functional satisfies the general properties stated by Feynman and Vernon. Considerations of when the higher-order terms are comparable or more important than the Feynman-Vernon terms are left for future work.
We end by summarize the assumptions that go and do not go into the higher-order theory. First, we assume that the system and the environment start out in a product state. Second, we assume that it is possible to write the system-environment interaction as V sb = k X k s ⊗ B k b , where X k s and B k b are operators on respectively the system and the environment, and where all the B k b commute. Third, we assume that the initial state of the environment is a product state compatible with the interaction. By the latter we mean that if the full environment Hilbert space is a product space H b = H 1 b ⊗ H 2 b · · · and the operators B j b act on H j b , then the initial environment density matrix factorizes as as is a unit trace positive Hermitian operator on H j b . One class of models that fullfill the above is when the system interacts with one or several baths which start out independent, and which do not interact between themselves. In the other direction, in each bath the environmental degrees of freedom can be either Bosonic or Fermionic (or both), and the Hamiltonians can be arbitrary. The initial state of each bath does not even have to be in equilibrium. We suspect that such a general-looking result will find applications also outside the current realm of theory of open quantum system.

VII. ACKNOWLEDGEMENT
This work was initiated at the Nordita program "New Directions in Quantum Information" (Stockholm, April 2019). EA thanks Dr Dmitry Golubev for discussions. RK thanks Garrett Higginbotham and Saarth Anjali Chitale for helpful discussion.
where the forward and backward evolution operators are defined by where ← − T and − → T are chronological and anti-chronological time ordering operator. The evolution of a density operator involves both forward and backward evolution operators as Managing the order of operators is a bit complicated due to the presence of two evolutions. There is a simpler expression using the time line shown in Fig. 2. We note that the evolution of the density operator is determined by the Liouville-von Neumann equation where the Liouville super-operator is defined by L(t) = where the time-ordering super-operator ← − T orders super-operators such as L(t) chronologically. It automatically orders regular operators along the time line shown in Fig. 2. As an example, consider t 1 > t 2 , which automatically orders V (t) chronologically if it is on the left of ρ(t i ) and anti-chronologically on the right.
permuted within themselves which can be done in n! × (N − n)! ways. The contribution from n forward paths (Q) and m backward paths (Q) is thus n time-ordered and m reverse time-ordered integrals multiplying the corresponding correlation function, which is (11) in main text.
Eq. (C-2) can be analyzed further by considering in the mixed terms the three ranges of v: less than u; between u and s, and larger than s. Renaming the variables so that times are always ordered s > u > v this gives Collecting terms with the same last entries one sees that this is The third-order terms hence satisfy the general property of the Feynman-Vernon action that if Q(s) =Q(s) for all s greater than τ , then the action does not depend on Q(s) orQ(s) for s > τ . This conclusion also holds more generally: starting from (11) in main text one can first insert u 1 in any of the intervals [t f ; s 1 ], [s 1 ; s 2 ], . . . , [s n ; t i ], then u 2 in the same interval as u 1 or any one further down the list, and so on. Each such insertion can be identified by a sequence Z τ1 , Z τ2 , . . . , Z τ N where each symbol is Q orQ, and the times are ordered τ 1 > τ 2 > · · · > τ N . Consider now two cases that only differ by the symbol Z τ1 . The first case has n symbols Q and N − n symbolsQ (0 < n ≤ N ), and arises from inserting (u 1 , u 2 , . . where the combined amplitudes can be written out as By comparison, the standard Feynman-Vernon action can be written in a similar way as where The contributions from the fourth order cumulants are analogously found to be