Nonperturbative renormalization of quantum thermodynamics from weak to strong couplings

By solving the exact master equation of open quantum systems, we formulate the quantum thermodynamics from weak to strong couplings. The open quantum systems exchange matters, energies and information with their reservoirs through quantum particles tunnelings that are described by the generalized Fano-Anderson Hamiltonians. We find that the exact solution of the reduced density matrix of these systems approaches a Gibbs-type state in the steady-state limit for both the weak and strong system-reservoir coupling strengths. When the couplings become strong, thermodynamic quantities of the system must be renormalized. The renormalization effects are obtained nonperturbatively after exactly traced over all reservoir states through the coherent state path integrals. The renormalized system Hamiltonian is characterized by the renormalized system energy levels and interactions, corresponding to the quantum work done by the system. The renormalized temperature is introduced to characterize the entropy production counting the heat transfer between the system and the reservoir. We further find that only with the renormalized system Hamiltonian and other renormalized thermodynamic quantities, can the exact steady state of the system be expressed as the standard Gibbs state. Consequently, the corresponding exact steady-state particle occupations in the renormalized system energy levels obey the Bose-Einstein and the Fermi-Dirac distributions for bosonic and fermionic systems, respectively. Thus, the conventional statistical mechanics and thermodynamics are thereby rigorously deduced from quantum dynamical evolution.


I. INTRODUCTION
Understanding the physical process of thermalization within the framework of quantum mechanical principle has been a long-standing problem.Thermodynamics and statistical mechanics are built with the hypothesis of equilibrium [1,2], that is, over a sufficiently long time, a macroscopic system which is very weakly coupled with a thermal reservoir can always reach thermal equilibrium, and its equilibrium statistical distribution does not depend on the initial state of the system.Over a century and a half, investigating the foundation of statistical mechanics and thermodynamics has been focused on two basic questions [3]: (i) how does macroscopic irreversibility emerge from microscopic reversibility? and (ii) how does the system relax to thermal equilibrium with its environment from an arbitrary initial state?Rigorously solving these problems from the dynamical evolution of quantum systems, namely, finding the underlie of disorder and fluctuations from the deterministic dynamical evolution, has been a big challenge in physics [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15][16][17][18][19].Obviously, the foundation of thermodynamics and statistical mechanics and the answers to these questions rely on a deep understanding of the dynamics of systems interacting with their environments, i.e., the nonequilibrium evolution of * Correspondence author: wzhang@mail.ncku.edu.twopen quantum systems.In 1980's, Caldeira and Leggett investigated the problem of thermalization from the study of the quantum Brownian motion, a Brownian particle coupled to a thermal reservoir made by a continuous distribution of harmonic oscillators [5].They used the Feynman-Vernon influence functional approach [4] to explore the dynamics of quantum Brownian motion, and found the equilibrium thermal state approximately [5].Later, Zurek studied extensively this nontrivial problem from the quantumto-classical transition point of view.Zurek revealed the fact that thermalization is realized through decoherence dynamics as a consequence of entanglement between the system and the reservoir [7].Thermalization in these investigations is demonstrated for quantum Brownian motion for initial Gaussian wave packets at high temperature limit [5,7].However, the thermalization with arbitrary initial state of the system at arbitrary initial temperatures of one or multiple reservoirs for arbitrary system-reservoir coupling strengths have not been obtained.
In the last decade, we have derived the exact master equation of open quantum systems [55][56][57][58][59][60][61][62] by extending the Feynman-Vernon influence functional theory into the coherent state representation [63].The open quantum systems we have studied are a large class of nano-scale open quantum systems that exchange matters, energies and information with their reservoirs through the particle tunneling processes.We also solved the exact master equation of these systems with arbitrary initial states at arbitrary initial reservoir temperatures.Thus, a rather general picture of thermalization processes has been obtained [10,12,18].In this paper, we shall explore the thermodynamic laws and statistical mechanics principles from the dynamical evolution of open quantum systems for both the weak and strong coupling strengths, based on the exact solution of the exact master equation we obtained.
In fact, the difficulty for building the strong coupling quantum thermodynamics is twofold [38][39][40][41][42][43][44][45][46][47][48][49][50]: (i) How to systematically determine the internal energy from the system Hamiltonian which must be modified by the strong coupling with its reservoirs?(ii) How to correctly account the entropy production when the system evolves from nonequilibrium state to the steady state?We find that the nature of solving the above difficulty is the renormalization of both the system Hamiltonian and the system density matrix during the nonequilibrium evolution through the system-reservoir couplings.The system-reservoir couplings also result in the dissipation and fluctuation dynamics in open quantum systems, which are indeed renormalization effects of the systemreservoir interactions.The renormalization effects can be obtained nonperturbatively after exactly traced over all reservoir states.They are manifested in the dynamical evolution of the reduced density matrix with dissipation and fluctuation, and accompanied by the renormalized system Hamiltonian.We develop such a nonperturbative renormalization theory of quantum thermodynamics from weak to strong couplings in this paper.
The rest of the paper is organized as follows.In Sec.II, we begin with the simple open quantum system of a nanophotonic system coupled with a thermal reservoir.The renormalized system Hamiltonian is obtained in the derivation of the exact master equation for the reduced density matrix.The exact solution of the reduced density matrix is also obtained analytically from the exact master equation.Its steady state approaches to a Gibbs state so that quantum thermodynamics emerges naturally.However, we find that the exact solution of the particle occupation in the system at strong coupling does not agree to the Bose-Einstein distribution with the initial reservoir temperature.This indicates that the corresponding equilibrium temperature must also be renormalized when the reduced density matrix is influenced by the system-reservoir interaction through the dissipation and fluctuation dynamics of the system.By introducing the renormalized temperature as the derivative of the renormalized system energy with respect to the von Neumann entropy in terms of the reduced density matrix, we overcome the inconsistency.Thus, the self-consistent renormalized quantum statistics and renormalized quantum thermodynamics are formulated for both the weak and strong coupling strengths.
In Sec.III, we extend such study to more general open quantum systems coupled to multi-reservoirs through particle exchange (tunneling) processes described by generalized Fano-Anderson Hamiltonians.These open systems are typical nano-scale systems that have been studied for various quantum transport in mesoscopic physics.
Here both systems and reservoirs are made of many bosons or many fermions, not limiting to the prototypical open system of a harmonic oscillator coupling to a oscillator reservoir introduced originally by Feynman [4] and by Caldeira and Leggett [5,64].From the exact master equation and its exact solution in the steady state for such class of open quantum systems [55][56][57], we develop the renormalization theory of quantum thermodynamics for both the weak and strong coupling strengths in general.We further take an electronic junction system (a single electronic channel coupled two reservoirs with different initial temperatures and chemical potentials) as a specific nontrivial application.It is a nontrivial example because other approaches proposed for strong coupling quantum thermodynamics in the last few years keep the reservoir temperature unchanged [38][39][40][41][42][43][44][45][46][47][48][49] so that these approaches become invalid for multiple reservoirs when the total system (the system plus all reservoirs) reaches a final equilibrium state.We demonstrate the consistency of the Fermi-Dirac statistics with our renormalized quantum thermodynamic in this nontrivial application.
In Sec.IV, we discuss further the generalization of this nonperturbative renormalization theory for quantum thermodynamics to more complicated interacting open quantum systems.We take the non-relativistic quantum electrodynamics (QED) derived from the fundamental quantum field theory as an example, and considered electrons as the open system and all photonic modes (electromagnetic field) as the reservoir.The system-reservoir interaction is the fundamental electron-photon interaction.We perform the nonperturbative renormalization by integrating out exactly the infinite number of electromagnetic field degrees of freedom.We obtain the re-duced density matrix in terms of only the system degrees of freedom in the same way as we derived the exact master equation for the generalized Fano-Anderson Hamiltonian in Sec.III.The resulting renormalization theory are given by the reduced density matrix for electrons and the nonperturbative renormalized electron Hamiltonian, which can be systemically computed in terms of twoelectron propagating Green functions in principle.Thus, we show that although our renormalized quantum thermodynamics theory is formulated from the exact solvable open quantum systems, it can apply to arbitrary open quantum systems even though the final exact analytical solution is hardly found.In fact, similar situation also exists for the equilibrium statistical mechanics, namely, one cannot solve exactly all equilibrium physical systems, in particular the strongly correlated systems such as the Hubbard model and the general quantum Heisenberg spin model [65] even though the reservoir effect can be ignored there.Therefore, approximations and numerical methods remained to be developed further for the study of renormalized nonequilibrium dynamics within the framework we developed in this paper.A conclusion is given in Sec.V.In Appendices, we provide the necessary analytical derivations of the solutions used in the paper.

II. A SIMPLE EXAMPLE FOR STRONG COUPLING QUANTUM THERMODYNAMICS
For simplicity, we begin with a single-mode bosonic open system (such as a microwave cavity in quantum optics or a vibrational phononic mode in solid-state and biological systems) coupled to a thermal reservoir through the energy exchange interaction.The total Hamiltonian of the system, the reservoir and the coupling between them is considered to be described by the Fano Hamiltonian [66]: where a † and b † k (a and b k ) are the creation (annihilation) operators of the bosonic modes in the system and in the reservoir with energy quanta ω s and ω k , respectively.They obey the standard bosonic commutation relations: The parameter V k is the coupling amplitude between the system and the reservoir and can be experimentally tuned to strong coupling [67,68].In fact, all parameters in the Hamiltonian, including the couplings between the system and the reservoir can be time-dependently controlled with the modern nano and quantum technologies.The universality of Fano resonance also makes this simple system useful in nuclear, atomic, molecular and optical physics, as well as condensed matter systems [69].
A. The exact master equation of the system and its exact nonequilibrium solution To study the thermalization of open quantum systems, the reservoir can be initially set in a thermal state where β 0 = 1/k B T 0 and T 0 is the initial temperature of the reservoir, Z E = Tr E [e −β0H E ] is its partition function.
The system can be initially in arbitrary state ρ S (t 0 ) so that the initial total density matrix of the system plus the reservoir is a direct product state [4,5], After the initial time t 0 , both the system and the reservoir evolve into an entangled nonequilibrium state ρ tot (t) which obeys the Liouville-von Neumann equation in quantum mechanics [70], Because the system and the reservoir together form a closed system, the Liouville-von Neumann equation is the same as the Schrödinger equation of quantum mechanics for the evolution of pure quantum states.But the Liouville-von Neumann equation is more general because it is also valid for mixed states.Quantum states of the system are completely determined by the reduced density matrix ρ S (t).It is defined by the partial trace over all the reservoir states: The equation of motion for ρ S (t), which is called the master equation, determines the quantum evolution of the system at later time t (> t 0 ).In the literature, one usually derives the master equation using various approximations, such as the memory-less dynamical maps, the Born-Markovian approximation, and secular approximation, etc. [71][72][73][74][75][76].But these methods are invalid for strong coupling open quantum systems with strong non-Markovian dynamics.In the past decade, we have developed a very different approach to rigorously derive the exact master equation for a large class of open quantum systems [55][56][57][58][59][60][61][62].Explicitly, we have derived the exact master equation for Eq. ( 1) by exactly tracing over all the reservoir states from the solution of the Liouvillevon Neumann equation [57,[77][78][79].The trace over all the reservoir states is a nonperturbative renormalization to the reduced density matrix of the system and to the system Hamiltonian simultaneously.We complete this partial trace by integrating out exactly all the reservoir degrees of freedom through the coherent state path integrals [57,63].The resulting exact master equation for the reduced density matrix accompanied with the renormalized system Hamiltonian is given by In this exact master equation, the first term describes a unitary evolution of the reduced density matrix with the renormalized Hamiltonian This renormalized Hamiltonian contains all the energy corrections to the system arisen from the system-reservoir interaction through the nonequilibrium evolution.The second and the third terms describe the non-unitary evolution of the reduced density matrix, which contain all non-Markovian dissipation and fluctuation dynamics induced by the back-reactions between the system and the reservoir through the system-reservoir interaction (see Eqs. ( 7) and ( 8) given later).Physically, the second and the third terms in the above master equation also characterize the emergence of disorder and fluctuations induced by the system-reservoir interaction.This is because if the system is initially in a pure quantum state, it contains zero disorder at beginning (its initial entropy is zero).The index r denotes renormalized physical quantities hereafter.The energy renormalization, the dissipation and fluctuation dynamics described in the exact master equation Eq. ( 6) are characterized by the non-Markovian renormalized frequency ω r s (t, t 0 ), the non-Markovian dissipation coefficient γ(t, t 0 ) and the non-Markovian fluctuation coefficient γ(t, t 0 ), respectively.All these non-Markovian coefficients are nonperturbatively and exactly determined by the following relations [57,[77][78][79] Here u(t, t 0 ) and v(τ, t) are the two non-equilibrium Green functions obeying the integro-differential Dyson equations, The non-Markovianity is manifested by the above timeconvolution equation of motion for these non-equilibrium Green functions.The integral kernels in the above con-volution equations are given by which characterize the time correlations between the system and the reservoir through the system-reservoir interaction.The frequency-dependent function is called as the spectral density, which fully encapsulates the fundamental dissipation (relaxation) and fluctuation (noise or dephasing) effects induced by the systemreservoir interaction.Finally, the initial temperature dependent function is the initial particle distribution in the reservoir.
An arbitrary initial state of the system can be expressed as where ) is a pure state, otherwise it is a mixed state.The exact solution of the exact master equation Eq. (6a) has be found [12,18], where As a self-consistent check of the above solution, we calculate the average particle number in the system from the above solution, n(t) ≡ Tr S [a † aρ S (t)], also using the Heisenberg equation of motion directly, n(t) ≡ Tr S + E [a † (t)a(t)ρ tot (t 0 )].Both calculations give the same result [56][57][58]61]: where u(t, t 0 ) and v(t, t) are determined by Eq. ( 8).Based on the above exact formalism, for a given spectral density J(ω), if no localized bound state exists [10,80], the general solution of Eq. (8a) is where [ω−ωs−∆(ω)] 2 +π 2 J 2 (ω) shows the system spectrum broadening due to the coupling to the reservoir, and ∆(ω) = P dω J(ω ) ω−ω is the principal-value integral of the self-energy correction to the system, Σ(ω) = dω J(ω ) ω−ω = ∆(ω) − iπJ(ω).In fact, ∆(ω) gives the system frequency (or energy) shift.As a result, in the steady state limit, we have Then the exact solution of the particle distribution and the reduced density matrix are reduced to Equation ( 17) is the exact steady-state solution of the system coupled to a thermal reservoir for all coupling strengths for the open system Eq.( 1).All the influences of the reservoir on the system through the systemreservoir interaction have been taken into account in this solution.Remarkably, the above results show that the exact solution of the steady state is independent of the initial state of the system and is determined by the particle distribution, as a consequence of thermalization [12,18].Note that the above exact master equation formalism remains the same for initial states involved initial correlations between the system and the reservoir, with the only modification of the correlation function v(τ, t), as we have shown in Refs.[58,62,81].This exact master equation formalism has also been extended to open quantum systems including external deriving fields [57,68].

B. Renormalization of quantum thermodynamics
Now we can study quantum thermodynamics for all the coupling strengths from the above exact solution.First, the master equation Eq. (6) shows that the Hamiltonian of the system must be renormalized from H S to H r S given by the energy (or frequency) shift from ω s to ω r s (t) during the nonequilibrium dynamical evolution.This is a nonperturbative renormalization effect of the system-reservoir coupling on the system.The renormalized frequency ω r s (t) and its steady-state value ω r s = ω r s (t → ∞) can be exactly calculated from Eqs. (7a) and (8a).Here we take the Ohmic spectral density J(ω) = ηω exp(−ω/ω c ) [82] in the practical calculation.The result is presented in Fig. 1(a) and (b).It shows that different system-reservoir coupling strengths η will cause different renormalized system energies, resulting in different cavity frequency shifts, see Fig. 1(a).In Fig. 1(b), we plot the steady-state values of the renormalized cavity frequency as a function of the system-reservoir coupling strength η/η c , where η c = ω s /ω c is a critical coupling strength for the Ohmic spectral density [10,12].When η > η c , the system-reservoir coupling would generate a localized mode (localized bound state) such that the cavity system cannot approach to the equilibrium with the reservoir, as we will discuss later [12,18].n(ω r s , T 0 ) The renormalized system energy (cavity frequency shift) ω r s (t) for three different system-reservoir coupling strengths, η = 0.01ηc, 0.5ηc, 0.9ηc.It is calculated from Eqs. (7a) and (8a) for the Ohmic spectral density J(ω) = ηω exp(−ω/ωc), where the cutoff frequency ωc = 10ωs is taken, and ηc = ωs/ωc is a critical coupling for the Ohmic spectral density [10,12].(b) The steady-state renormalized frequency shift ω r s = ω r s (t → ∞) as a function of the systemreservoir coupling strength η/ηc.(c) The steady-state particle distribution nexact(t → ∞) of Eq. (17a) (the blue-dashed line), the Bose-Einstein distribution without the energy (frequency) renormalization n(ωs, T0) (the black-dot line) and with the energy renormalization n(ω r s , T0) (the green-dasheddot line), respectively.The system is initially set in a pure Fock state |n0 with n0 = 5, and the reservoir initial temperature T0 = 10 ωs.
In Fig. 1(c), we plot the exact solution n exact (t → ∞) of Eq. (17a) as a function of the coupling strength η/η c (the blue-dashed line).We compare the result with the Bose-Einstein distribution without the energy (frequency) renormalization, n(ω s , T 0 ) = 1/[e ωs/k B T0 − 1] (see the black-dot line), also compare to Bose-Einstein distribution with the energy renormalization, n(ω r s , T 0 ) = 1/[e ω r s /k B T0 − 1] (see the green-dashed-dot line).As one can see, the exact solution n exact (t → ∞) derivates significantly from the Bose-Einstein distribution without the energy renormalization, i.e., n(ω s , T 0 ), as η increases.This derivation shows how the system-reservoir coupling strength changes the intrinsic thermal property of the system.On the other hand, the Bose-Einstein distribution with the renormalized energy, given by n(ω r s , T 0 ), changes with the changes of η, similar to the exact solution n exact (t → ∞).But there is still a quantitative disagreement between the exact solution n exact (t → ∞) and the Bose-Einstein distribution n(ω r s , T 0 ) with the renormalized cavity photon energy ω r s .To understand further the origin of the above difference, let us recall that the exact solution ρ exact S (t → ∞) of Eq. ( 17c) is indeed a Gibbs-type state.This indicates that the exact particle distribution n exact (t → ∞) should obey a Bose-Einstein distribution for all coupling strengths.To find such distribution that agrees with the exact solution Eq. (17a), one possibility is to renormalize the temperature, because no other thermal quantity can be modified in the Gibbs state for this photonic system.In the literature, it is commonly believed that the reservoir is large enough so that its temperature should keep invariant [38,41].However, the initial decoupled states Eq. (3) of the system plus the reservoir is not an equilibrium state of the total system.After the initial time t 0 , both the system and reservoir evolve into a correlated (entangled) nonequilibrium state ρ tot (t).When the system and the reservoir reach the equilibrium state, there must be a fundamental way to show whether the new equilibrium state is still characterized by the initial reservoir temperature.
As a self-consistent check, let us denote the final steady-state equilibrium temperature as T f .Then, according to the equilibrium statistical mechanics, the steady state of the total system (the system plus the reservoir) should be where , and H tot is the total Hamiltonian of the system plus the reservoir, including the coupling interaction between them, i.e., Eq. (1).Taking a trace over the reservoir states from the above steady state of the total density matrix, we have rigorously proven [62] that (also see the detailed derivation given in Appendix A) This result is the same as the solution of Eq. ( 17).The latter is the steady state of the exact time-dependent solution Eq. ( 13) solved from the exact master equation Eq. (6a) for arbitrary coupling.This shows that the equilibrium state Eq. ( 18) which is originally proposed in statistical mechanics is indeed valid for both the weak and strong coupling between the system and the reservoir.Furthermore, the exact particle distribution can be obtained from the dynamical evolution of exact master equation or from the Heisenberg equation of motion directly, as shown by Eq. ( 15).Thus, we have This gives a further self-consist justification to the above conclusion.The result presented in Fig. 1(c) shows that n(ω r s , T 0 ) = n exact (t → ∞) except for the weak coupling.This indicates that in general, T f = T 0 , namely the final equilibrium temperature of the total system cannot be the same as the initial equilibrium temperature of the reservoir when the total system reaches the new equilibrium state, except for the very weak coupling strength.Now the question is how to determine this steady-state equilibrium temperature T f when the system and the reservoir finally reach the equilibrium state.According to the axiomatic description of thermodynamics [83], the equilibrium temperature of a system is defined as the change of its internal energy with respect to the change of its thermal entropy.This temperature definition in thermodynamics does not assume a weak coupling between the system and the reservoir because no statistical mechanics is used in this definition.It is the fundamental definition of the temperature for arbitrary two coupled thermodynamic systems when they reach the equilibrium each other, from which the zeroth law of thermodynamics is derived [83].Now, the average energy of the system at arbitrary time, i.e., the nonequilibrium internal energy of the system, is given by the renormalized Hamiltonian Eq. (6b) with the exact solution of the reduced density matrix ρ S (t) of Eq. ( 13): Also, we define the von Neumann entropy with the exact reduced density matrix of Eq. ( 13) as the nonequilibrium thermodynamic entropy of the system [50,70,83]: where k B is the Boltzmann constant.Note that this entropy is defined for the exact reduced density matrix obtained after traced over exactly all the reservoir states.
It also encapsulates all the renormalization effects of the system-reservoir interaction to the system state distributions.Thus, we introduce a renormalized nonequilibrium thermodynamic temperature [84] which is defined as This is a direct generalization of the concept of the equilibrium temperature to nonequilibrium states in open quantum systems.When the system and its reservoir reach the equilibrium steady state, no matter the systemreservoir coupling is strong or weak, we can fundamentally obtain the final equilibrium temperature T f ≡ T r = T r (t → ∞) from the dynamical evolution of the open quantum system.

FIG. 2. (a)-(c)
The nonequilibrium dynamical evolution of the internal energy, the entropy and the corresponding renormalized temperature, given respectively by Eqs. ( 20)-( 22) for different coupling strengths η/ηc = 0.3, 0.5, 0.8 (correspond to the blue-dot-dashed, green-dashed, red-dot lines, respectively).(d) The steady-state values of the renormalized frequency ω r s and renormalized temperature T r as a function of the coupling strength η.(e) The steady-state particle distribution as a function of coupling strength η.It shows that the exact solution nexact(t → ∞) of Eq. (17a) (the blue-dashed line) and the Bose-Einstein distribution n(ω r s , T r ) with the renormalized frequency and the renormalized temperature (the red-dot line) agree perfectly to each other.The greendashed-dot line is n(ω r s , T0) without the renormalized temperature, which cannot describe the exact solution solved from the exact master equation.Other parameters are taken as the same as that in Fig. 1.
From the exact solution of the reduced density matrix ρ exact S (t) of Eq. ( 13), we calculate the time-dependence of the internal energy U S (t), the entropy S S (t) and the dynamical renormalized temperature T r (t) for different coupling strengths.The corresponding results are presented in Fig. 2(a)-(c), respectively.It shows explicitly how the nonequilibrium internal energy, entropy and renormalized temperature evolve differently for different system-reservoir coupling strengths.Their steady-state values also approach different points for different coupling strengths.The different steady-state internal energies and entropies associated with different systemreservoir coupling strengths result in different steady-state temperatures.This indicates that the reservoir cannot remain unchanged from the initial reservoir temperature.This new feature has not been discovered or noticed in all previous investigations of strong-coupling quantum thermodynamics [38][39][40][41][42][43][44][45][46][47][48][49].
In Fig. 2(d), we plot the steady-state renormalized temperature, T r = T r (t → ∞), as a function of the coupling astrength η/η c .Using this renormalized temperature, we further plot the Bose-Einstein distribution with both the renormalized energy and the renormalized temperature: n(ω r s , T r ) = 1/[e ω r s /k B T r − 1], see the red-dot line in Fig. 2(e).Remarkably, it perfectly reproduces the exact solution of Eq. (17a), i.e., In other words, in the steady state, the exact solution of the steady-state particle occupation solved from the exact dynamics of the open quantum system obeys the standard Bose-Einstein distribution only for the renormalized Hamiltonian Eq. (6b) with the renormalized temperature Eq. ( 22).This provides a very strong proof that strong coupling quantum thermodynamics must be renormalized for both the system Hamiltonian and the temperature.Furthermore, in terms of the renormalized Hamiltonian Eq. (6b) and the renormalized temperature Eq. ( 22), the steady state Eq.(17c) can be expressed as the standard Gibbs state, where Z r S = Tr S [e −β r H r S ] is the renormalized partition function, and β r = 1/k B T r is the inverse renormalized temperature in the steady state.This is a direct proof of how the statistical mechanics, as a consequence of disorder or randomness in the nature, emerges from the exact dynamical evolution of quantum mechanics.
Moreover, one can check that in the very weak coupling regime η η c , ∆(ω) → 0 and D(ω) → δ(ω − ω s ) so that the steady state solution of Eq. [17a] is directly reduced to n(ω s , T 0 ) [12,18], and This reproduces the expected solution of the statistical mechanics in the weak coupling regime.Figures 1 and 2 also show that ω r s → ω s and T r → T 0 at very weak coupling.Thus, the equilibrium hypothesis of thermodynamics and statistical mechanics is deduced rigorously from the dynamics of quantum systems.This solves the long-standing problem of how thermodynamics and statistical mechanics emergence from quantum dynamical evolution [3].
On the other hand, η c = ω s /ω c is a critical coupling strength for Ohmic spectral density that when η > η c , the system exists a dissipationless localized bound state (localized mode) at frequency ω b = ω s + ∆(ω b ) with J(ω b ) = 0 [10].Once such a localized mode exists, the spectral function D(ω) of the system in Eq. ( 16) is mod-ified as where is the localized bound state wavefunction.Then, the asymptotic value of the Green function u(t → ∞, t 0 ) never vanishes.As a result, the steady state of the reduced density matrix Eq. ( 13) cannot be reduced to Eq. ( 17).It always depends on the initial state distribution ρ mm of Eq. ( 12).In other words, the system cannot be thermalized with the reservoir [12,18,84].The nonequilibrium dynamical evolution of the internal energy U S (t), the entropy S S (t) and the renormalized nonequilibrium temperature T r (t) for different system initial states |n0 = |5 , |10 and |15 (corresponding to the blue dasheddot line, the grees dashed line and the red solid line, respectively).The left, the midden and the right panels correspond to (a) the weak coupling (η ηc), (b) the strong coupling η < ηc), and (c) the ultra-strong coupling η > ηc).The initial bath temperature T0 = 30 ωs.Other parameters are taken the same as that in Fig. 1.
In Fig. 3, we plot the nonequilibrium dynamical evolution of the internal energy, the entropy production and the renormalized temperature with different initial states for the very weak coupling (η = 0.01η c η c ), the strong coupling (η = 0.5η c < η c ) and the ultra-strong coupling (η = 1.2η c > η c ) cases.The results show that when η > η c , different initial states of the system lead to different steady states.That is, the equilibrium hypothesis of the classical thermodynamics and statistical mechanics is broken down at ultra-strong coupling.Also note that after consider the renormalization of the system Hamiltonian, the dynamics of the internal energy and the renormalized temperature are significantly changed, in particular in the strong coupling regime, in the comparison with our previous study [84] where no renormalization of the system Hamiltonian is taken into account.On the other hand, regarding the system and the reservoir as a many-body system, the existence of the localized bound state in the regime η > η c corresponds to a realization of the many-body localization [13].When the coupling strength η crosses the critical value η c , the transition from thermalization to many-body localization occurs [13].Our exact solution provides the foundation of this transition between thermodynamics and many-body localization.

C. Quantum work and quantum heat
Quantum mechanics does not introduce the concepts of work and heat because it deals with closed systems.For open quantum systems, the exchanges of matters, energies and information between the system and the reservoir cause the energy change of the system in the nonequilibrium evolution.This results in the work and chemical work (associated with chemical potential) done on the system or by the system, and the heat flowed into or out of the system.But usually the exchanges of matters, energies and information are correlated and interfered each other.This makes difficulties to define clearly the concepts of work, heat and chemical work in quantum thermodynamics.For the photons and phonons described by Eq. ( 1), no matters exchange between the system and the reservoir so that no chemical work involves (chemical potential is zero here).Thus, the energy change of the system only involves with work and heat.After integrated out exactly the reservoir degrees of freedom, the reduced density matrix Eq. ( 13) and the associated renormalized system Hamiltonian Eq. (6b) can be used to properly define thermodynamic work and heat within the quantum mechanics framework.The chemical work will be considered when we study fermion systems, as we will discuss in the next section.
As it is shown from Eq. ( 20), the nonequilibrium change of the internal energy in time contains two parts.One is the change (i.e. the renormalization) of the system Hamiltonian H r S (t) (through the renormalization of the energy level ω r S (t)) which corresponds to the quantum work done on the system.Note that in quantum mechanics, the concept of volume in a physical system is mainly characterized by energy levels through energy quantization.Thus, the change of volume is naturally replaced by the change of energy levels, which results in a proper definition of work in quantum mechanics [85].The other part is the change of the density state ρ S (t) which corresponds to quantum heat associated with the entropy production.Consequently, This is the first law of nonequilibrium quantum thermodynamics.Thus, the quantum work and quantum heat can be naturally determined by where n(t) = Tr S [a † aρ S (t)] is given by Eq. ( 15).The second equalities in the above equations have used explicitly the renormalized system Hamiltonian given after Eq. (6a) and the definition of the renormalized temperature Eq. ( 22), respectively.In the literature, there are various definitions about work and heat for strong coupling quantum thermodynamics but no consensus has been reached.The main concern is how to correctly include the system-reservoir coupling energy into the internal energy of the system [37,44,48,[51][52][53][54].The difficulty comes from the fact that most of open quantum systems cannot be solved exactly so that it is not clear how to properly separate the contributions of the system-reservoir coupling interaction into the system and into the reservoir, respectively.However, this difficulty can be overcome in our exact master equation formalism.Because we obtain the renormalized system Hamiltonian accompanied with the reduced density matrix after integrated out exactly all the reservoir degrees of freedom, namely exactly solved the partial trace over all the reservoir states.Thus, the renormalized system Hamiltonian contains all possible contributions of the system-reservoir coupling interaction to the system energy.
Explicitly, let us rewrite the renormalized system Hamiltonian Eq. (6b) as From Eq. (7a) and Eq.(8a), we have Here H S = ω s a † a is the bare Hamiltonian of the system.The second term in Eq. ( 30) contains all order contributions of the system-reservoir coupling interaction to the system energy, as shown in Fig. 4. Figure 4(a) is a diagrammatic plot of the bare Hamniltonians of the system, the reservoir and the interaction between them, respectively.Figure 4(b) is the diagrammatic expansion (up to infinite orders) of the retarded Green function of Eq. (8a), from which all order renormalization effects to the system energy change (the system frequency shift) are reproduced.This diagrammatic expansion up to the infinite orders illustrate the nonperturbative renormalized energy arisen from the system-reservoir interaction in our exact master equation theory.
On the other hand, it is interesting to see that if we replace the full solution of the Green function u(τ, t 0 ) approximately with the free-particle (zero-th order) Green function u 0 (τ, t 0 ) = e −iω S (τ −t0) (also for u(t, t 0 )) in Eq. ( 30), the result is just the second-order renormalized energy correction.By applying this same approximation to the dissipation and fluctuation coefficients in Eq. ( 7), it is straightforward to obtain the time-dependent decay rate and noise in the Born-Markovian master equation, as we have shown in our previous work [78].But once we have the exact master equation with the exact solution, such approximated master equation is no longer needed.The diagrammatic Dyson expansion of Eq. (8a) in the energy domain, where Σ(ω) = dω J(ω ) ω−ω is the self-energy arisen from the coupling between the system and the reservoir, and The renormalized system energy ω r s and the dissipation coefficient γ of Eqs.(7a)-(7b) are determined nonperturbatively from Eq. (8a) with a Laplace transformation, which contains all order contributions upto the infinite orders from the system-reservoir coupling Hamiltonian, as shown in this diagrammatic expansion.
In Fig. 5 (a)-(b), we plot the nonequilibrium evolution of dW S (t)/dt and dQ S (t)/dt for different coupling strengths.The negative values of dW S (t)/dt show quantum work done by the system during the quantum mechanical time evolution, and more work is done by the system for the stronger system-reservoir coupling.While, dQ S (t)/dt is negative and then becomes positive in time, which shows that quantum heat flows into the reservoir in the beginning and then flows back to the system in later time.This corresponds to the system dissipate energy very quickly into the reservoir in the very beginning, and then the thermal fluctuations arisen from the reservoir makes the heat flowing back slowly into the system.This heat flowing process can indeed be explained clearly from the exact master equation Eq. (6a) combined with Eq. (28b).It directly results in where the first term is the contribution from dissipation and the second is the contribution of fluctuations in our exact master equation.That is, the heating flow in open quantum systems is a combination effect of dissipation and fluctuation dynamics, which makes the system and the reservoir approach eventually to the equilibrium.This is also a renormalization effect.Furthermore, the quantum Helmholtz free energy of the system is defined by a Legendre transformation from the internal energy U S (t) [50,83]: From the above solution, we have which naturally leads to the consistency that the quantum thermodynamic work done on the system can be identified with the change of the Helmholtz free energy of the system in isothermal processes [83].Moreover, the specific heat calculated from the internal energy and from the Gibbs state with the renormalized Hamiltonian and temperature are also identical, as shown in Fig. 5(c), where the third thermodynamic law is justified from the specific heat at arbitrary coupling strength: C ∼ (T r ) 3 as T r → 0. Thus, a consistent formalism of quantum thermodynamics from the weak coupling to the strong coupling is obtained from a simple open quantum system of Eq. (1).

III. THE MORE GENERAL FORMULATION OF QUANTUM THERMODYNAMICS FOR ALL COUPLINGS A. Multi-level open quantum system couple to multiple reservoirs
The results from the exact solution of the single-mode bosonic open system in the last section show that different from the previous investigations [38][39][40][41][42][43][44][45][46][47][48][49][50], only by introducing the renormalized temperature and incorporating with the renormalized system Hamiltonian, can we obtain the consistent quantum thermodynamics for all coupling strengths.Now we extend this quantum thermodynamics formulation to the more general situation: a multi-level system couples to multiple reservoirs (including both bosonic and fermionic systems) through the particle exchange (tunneling) processes.
In a quasiparticle picture, the Hamiltonian of a microscopic system in the energy eigenbasis can be written as  34) (red lines) and from the partition function given in the Gibbs state Eq. ( 24) for different initial temperatures.The dashed-dot, dot and dashed lines correspond to the different coupling strengths η/ηc = 0.3, 0.5, 0.8, respectively.Other parameters are taken the same as that in Fig. 1.
As a specific example, consider the system be an individual system and the reservoir be a manybody system.The system Hamiltonian can be generally expressed as In Eq. ( 35), the second equality is the spectral decomposition of the system Hamiltonian: H S |ψ i = ε i |ψ i , and the last equality uses the second quantization language: |ψ i = a † i |0 and a i |0 = 0, and |0 is the vacuum state.The particle creation and annihilation operators a † i and a i obey the standard bosonic commutation and fermionic anticommutation relations: [a i , a † j ] ∓ = a j a † j ∓ a † j a i = δ ij when the system being boson and fermion systems, respectively.
Similarly, the Hamiltonian of a reservoir can also be written as where k is usually a continuous spectrum and could have band structure for structured reservoir.For a many-body reservoir in which the particle-particle interaction is not strong enough, the single quasiparticle picture works [86].Then the reservoir Hamiltonian can be expressed approximated as where j V (q j , q j ) represents the effective mean-field potential of many-body interactions, and gives the quasiparticle continuous spectrum of the reservoir.The reservoir particle creation and annihilation operators b † k and b k also obey the standard bosonic commutation or fermionic anticommutation relations.In fact, the system can also be either a simple system or such a many-body system.
To dynamically address statistical mechanics and thermodynamics from quantum mechanical principle, the fundamental system-reservoir interactions are required to contain at least the basic physical processes of energy exchanges, matter exchanges and information exchanges between the system and reservoirs.The simplest realization for such a minimum requirement is the quantum tunneling Hamiltonian, which is also the basic Hamiltonian in the study of quantum transport in mesoscopic physics as well as in nuclear, atomic and condensed matter physics for various phenomena [56,57,69,87].The coupling strengths V ik are proportional to the quasiparticle wavefunction overlaps between the system and reservoirs and therefore are tunable through nanotechnology manipulations [69,87] so that they can be weak or strong coupling.More discussions about fundamental system-reservoir interactions will be given in the next section.Thus, a basic Hamiltonian with the minimum requirement for solving the foundation of quantum thermodynamics and statistical mechanics can be modeled as which describes the system one concerned couples to multiple reservoirs.This is a generalization of the Fano-Anderson Hamiltonian we introduced [10,61].The index α stands for different reservoirs.All parameters in the Hamiltonian can be time-dependently controlled with the current nano and quantum technologies.This is an exact solvable Hamiltonian that involves explicit exchanges of energies, matter and information between the system and reservoirs.It allows us to rigorously solve quantum statistics and thermodynamics from the dynamical evolution of quantum systems.Also note that the above open quantum systems are different from the one proposed by Feynman and Vernon [4] as well as by Caldeira and Leggett [64] in the previous investigations of dissipative quantum dynamics in the sense that their environment is made only by harmonic oscillators and the system-environment coupling is limited to the weak coupling.
We have derived the exact master equation of the open systems with Eq. ( 38) for the reduced density matrix of the system.The formal solution of the total density matrix of the Liouville-von Neumann equation ( 4) can be expressed as: where U(t, t 0 ) = T → exp − i t t0 H tot (t )dt is the time evolution of the total system, and T → is the time-ordering operator.Here the system is initially in an arbitrary state ρ S (t 0 ).All reservoirs can be initially in their own equalibrium thermal states, Here N α is the total particle number operator of reservoir α.After trace out all the environmental states through the coherent state path integrals [63], the resulting exact master equation of the system is indeed a generalization of Eq. (6a) to multi-level open systems [55-59, 61, 62], where the upper and lower signs of ± correspond respectively to the bosonic and fermionic systems.
In the above exact master equation, all the renormalization effects arisen from the system-reservoir interactions have been taken into account when all the environmental degrees of freedoms are integrated out nonperturbatively and exactly in finding the reduced density matrix.These renormalization effects are manifested by the renormalized system Hamiltonian, and the dissipation and fluctuations coefficients γ ij (t, t 0 ) and γ ij (t, t 0 ) in Eq. ( 40).These time-dependent coefficients are determined nonperturbatively and exactly by the following relations, where u(t, t 0 ) and v(t, t) are N ×N nonequilibrium Green function matrices and N is the total number of energy levels in the system.The nonequilibrium retarded Green functions u ij (t, t 0 ) ≡ [a i (t), a † j (t 0 )] ± which obeys the equation of motion [10,[55][56][57], The nonequilibrium correlation Green function v(t, t) obeys the nonequilibrium fluctuation-dissipation relation [10], The integral memory kernels g(t, t ) and g(t 1 , t 2 ) are the system-reservoir time correlations and are given by Here the initial reservoir correlation function, determines the initial particle distribution of the bosons or fermions in the initial thermal reservoir α with the chemical potential µ α0 and the temperature T α0 at initial time t 0 .In the case the energy spectra of the reservoirs and the system-reservoir couplings are time-independent, the memory kernels are simply reduced to , where and J αij ( ) is the spectral density matrix of reservoir α.
B. The theory of quantum thermodynamics from the weak to the strong couplings Again, if there exist no many-body localized bound states, the exact solution of Eq. ( 40) has recently been solved [18] and its exact steady state is (see a detailed derivation in Appendix B) which is a generalized Gibbs-type state.Here n ij = lim t→∞ n ij (t) is the one-particle density matrix defined as [56,88] The solution Eq. ( 47) remains the same for initial systemreservoir correlated states with a modification of g(t 1 , t 2 ) in Eq. (43b) to include the initial correlations between the system and reservoirs [58,62].Thus, the nonequilibrium internal energy, entropy and particle number can be defined by They are related to each other and may form the fundamental equation for quantum thermodynamics [50,83]: U S (t) = U S (ε r s (t), S S (t), N S (t)).Here energy levels play a similar role as the volume [85].Thus, dU S (t) = dW S (t) + T r (t)dS S (t) + µ r (t)dN S (t), (50) as the first law of nonequilibrium quantum thermodynamics.
Explicitly, the quantum work dW S (t) done on the system is arisen from the changes of energy levels, The quantum heat dQ S (t) (also including the chemical work dW c S (t)) comes from the changes of particle distributions and transitions (the one-particle density matrix, see Eq. ( 48)), It shows that dn ij (t) characterizes both the state information exchanges (entropy production) and the matter exchanges (chemical process for massive particles) between the systems and the reservoir.For photon or phonon systems, particle number is the number of energy quanta ω so that µ r (t) = 0. From the above formulation, we can define the renormalized temperature and renormalized chemical potential by As a result, Eq. ( 47) can be also written as the standard Gibbs state, which is given in terms of the renormalized Hamiltonian H r S (t), the renormalized temperature T r (t) and the renormalized chemical potential µ r (t) at steady state, and N = i a † i a i is the particle number operator of the system.Because the exact solution of the steady state is a Gibbs state, thermodynamic laws are all preserved at steady state.This completes our nonperturbative renormalization theory of quantum thermodynamics for all the coupling strengths.
C. An application to a nanoelectronic system with two reservoirs As a practical and nontrivial application: we consider a nanoelectronic system, the single electron transistor made of a quantum dot coupled to a source and a drain.Here the two leads which are treated as two reservoirs [55,56,87], see Fig. 6(a).The total Hamiltonian is The index σ =↑, ↓ labels electron spin, α = L, R labels the left and right leads.The two leads are setup initially in thermal states with different initial temperatures T L,R and chemical potentials µ L.R .This is a prototype with nontrivial feature in the sense that two reservoirs initially have different temperatures and different chemical potentials so that when the system reaches the steady state, there exists only one final temperature and one final chemical potential.That is, one has to introduce the renormalized temperature T r and the renormalized chemical potential µ r to characterize this final equilibrium state when the system and two reservoirs reach equilibrium.While, other approaches proposed for strong coupling quantum thermodynamics in the last few years [38][39][40][41][42][43][44][45][46][47][48][49] keep the reservoir temperature unchanged and therefore must be invalid for such simple but nontrivial open quantum systems.
To explicitly solve the renormalized thermodynamics of the above system, let |0 , | ↑ , | ↓ , |d (the empty state, the spin up and down states and the double occupied state, respectively) be the basis of the 4-dim dot Hilbert space of this quantum dot system.Then the reduced density matrix has the form, If the dot is initially empty, the 4 × 4 reduced density matrix has been solved exactly from the exact master equation [89,90]: Here the 2 × 2 matrix Green function v(t) ≡ v(t, t) is determined by the Green function u(t, t ).We take reservoir spectra as a Lorentzian form, then the spectral densities J α ( ) can be expressed as [55,91]: where Γ α is the tunneling rate (the coupling strength) between the quantum dot and the lead α.For simplicity, we also ignore the spin-flip tunneling.The exact solution of the reduced density matrix is rather simple, Here a † = (a † ↑ , a † ↓ ), and for i =↑, ↓, and Γ = Γ L +Γ R .As a result, the nonequilibrium internal energy, the entropy and the total average particle number can be found analytically, From the above solution, the corresponding renormalized energy, renormalized temperature and renormalized chemical potential can be calculated straightforwardly.
FIG. 6.(a) A schematic plot of the single-electron transistor device.(b)-(g) The nonequilibrium evolution of the energy levels ε r ↑,↓ (t), the particle occupation in each level n ↑,↓ (t), the internal energy U S (t), the entropy S S (t), the renormalized temperature T r (t) and chemical potential µ r (t) at different coupling strength Γ = 0.2ε ↑ , 0.8ε ↑ , respectively.(h) The steady-state value of the renormalized energy levels ε r ↑,↓ , the renormalized temperature T r and the renormalized chemical potential µ r as a function of the coupling strength, and (i) the comparison of the renormalized Fermi-Dirac distribution f (ε r ↑,↓ , T r , µ r ) with the exact solution of the n ↑,↓ (t → ∞) as a function of the coupling strength Γ.Other parameters: In Fig. 6(b)-(g), we show the nonequilibrium evolution of the renormalized energy levels ε r ↑,↓ (t), the particle occupations in each levels n ↑,↓ (t), the internal energy U S (t), the entropy S S (t), the renormalized temperature T r (t) and the renormalized chemical potential µ r (t) for the coupling strength Γ L = Γ R = Γ/2 for different coupling strengths.It shows that in such a nano-scale device, all physical quantities are quickly approach to the steady state.Then, in Fig. 6(h), we plot the steady-state values of the renormalized energy levels ε r ↑,↓ , the renormalized temperature T r and the renormalized chemical potential µ r as a function of the coupling strength, respectively.These renormalized thermodynamical quantities change as the change of the coupling strength.Finally, in Fig. 6(i), we present the corresponding renormalized Fermi-Dirac distributions (the Fermi-Dirac distribution with renormalized energy, the renormalized temperature and the renormalized chemical potential): We compare the renormalized Fermi-Dirac distributions with the exact solution of the occupation numbers n ↑,↓ (t → ∞) which are solved from the exact master equation.The results shows that they completely agree with each other.This provides the proof to the consistency of the renormalized strong coupling quantum thermodynamics for fermionic systems.
To have a clearer physical picture about the renormalized temperature and renormalized chemical potential when the system coupled to two reservoirs, we take various different setups of the initial temperatures and initial chemical potentials of the two reservoirs in Fig. 7. From these results, we can see how the renormalized temperature and the renormalized chemical potential changes for the different setups, even in the weak-coupling regime.To understand these results, we first compare the exact solution with its weak-coupling limit.Since we also take the same spectral density for two reservoirs, we find that in the very weak-coupling limit (WCL), The first equality is the exact solution from Eqs. (61c) and the second equality is obtained with the help (60) in the very weak-coupling limit, where µ L , T L and µ R , T R are the initial chemical potential and temperatures of the left and right reservoirs, respectively.Figure 7(a) shows the results for the two reservoirs that have the same initial temperature and the same initial chemical potential.
Because the two reservoirs are set to have the same spectral density, the two reservoirs are equivalent to one single reservoir in this case.Thus, the renormalized temperature and the chemical potential approach to the initial temperature and the initial chemical potential in the very weak coupling limit, as shown in Fig. 7(a), also as we expected.Figure 7(b) shows the results for the two reservoirs sharing the same initial temperature but having different initial chemical potentials.Naively, one may think that the renormalized temperature in the very weak coupling limit should be the same as the same initial temperature of the two reservoirs and the renormalized chemical potential should be µ r = (µ L +µ R )/2.From Fig. 7(b), we see that the renormalized chemical potential in the very weak coupling limit is µ r = (µ L + µ R )/2, as we expected from energy conservation law.However, the renormalized temperature is a bit larger than the initial temperature.This result can actually be understood from Eq. ( 62).Because µ r = (µ L + µ R )/2 = µ L = µ R , Eq. ( 62) shows that T L = T r = T R in the very weak-coupling limit, even through T L = T R .Figure 7(c) shows further the case µ L = µ R and T L = T R .We have µ r = µ L = µ R , and from Eq. ( 62), we find that T r = (T L + T R )/2 in the very weak coupling limit, as shown in Fig. 7(c).Figure 7(d) shows the high temperature limit in which the chemical potentials play a little role.Thus, we have T r (T L + T R )/2 in the very weak coupling limit, even if µ L = µ R .This is shown in Fig. 7(d).These results demonstrate that only at very high temperature, the renormalized temperature T r = (T L + T R )/2.In other words, in the quantum regime, the renormalized temperature we introduced is necessary even at very weakcoupling limit for multi-reservoirs.This justifies further the consistency of our renormalized theory for quantum thermodynamics at any coupling.7. The steady-state renormalized temperature and renormalized chemical potential of the single-electron transistor as a function of the system-reservoir coupling strength Γ changing from weak to strong for different setups of the initial temperatures and initial chemical potentials of the two leads (reservoirs): (a) Two reservoirs have the same initial temperature and chemical potential; (b) The initial temperatures of the two reservoirs are the same but their initial chemical potentials are different; (c) The initial chemical potentials of the two reservoirs are the same but their initial temperatures are different; and (d) Two reservoirs at high temperature limit.

IV. EXTENSION TO MORE ARBITRARY SYSTEM-RESERVOIR COUPLINGS
The renormalized quantum thermodynamics for arbitrary coupling strength presented in Sec.III, given by Eq. ( 49) to Eq. ( 54), is formulated from the exact master equation ( 40) based on the system-reservoir coupling of Eq. (38).However, this formulation can be directly extended to general open quantum systems with systemreservoir couplings not limiting to the form of Eq. (38).This is because the renormalized Hamiltonian H r S (t) is determined by the nonequilibrium Green function u(t, t 0 ) of Eq. (43a).It can be applied to arbitrary system interacting with arbitrary environment.In Eq. (43a), g(t, t ) is the self-energy correlation that can be easily generalized to any interacting system using the nonequilibrium Green function technique in many-body systems [92,93].Meanwhile, the renormalized temperature T r (t) (the renormalized chemical potential µ r (t)) of Eq. ( 53) are determined by the changes of internal energy of the system with respect to the changes of the von Neumann entropy (the average particle number of the system).These nonequilibrium thermodynamic quantities are well defined by Eq. ( 49).They rely neither on the exact master equation of Eq. ( 40) nor the system-environment coupling in the Hamiltonian of Eq. (38).They are all determined by the reduced density matrix which can be solved from the Liouville-von Neumann equation ( 4), whose formal solution can be expressed as Here, U(t, t 0 ) is the quantum evolution operator, the same as the one given after Eq. ( 39) but the total Hamiltonian can be extended to arbitrary system interacting with arbitrary reservoir.In most of cases, taking the trace over the environmental states is the most difficult problem in open quantum systems.Practically, one can use the perturbation expansion method to calculate the trace over the environmental states order by order approximately [73], or using the coherent state path integrals to nonperturbatively trace over all the environmental states as we did [55-58, 61-63, 94].Here we focus on the nonperturbation procedure.All the renormalization effects of the system-reservoir interactions on the system can be obtained from this procedure.To be specific, let us consider a general fermionic system coupled to a general bosonic reservoir.Notice that Eq. ( 38) describes the exchanges of energies and particles between the system and the reservoir only for both the system and the reservoir that are made of the same type of quasiparticles, either bosons or fermions.When the system is a fermionic system and the reservoir is a bosonic system, the system-reservoir coupling Hamiltonian generally has the following interaction form [10,61] This system-reservoir interaction describes the energy exchange between the system and the reservoir through the transition of a fermion (e.g. an electron) between two states by emitting a boson (a quanta of energy, such as a photon or a phonon) into the reservoir or absorbing a boson from the reservoir.The creation and annihilation operators c † i , c i (b † k , b k ) obey the standard fertmionic anticommutation (bosonic commutation) relationships.In fact, Eq. ( 64) is the general form of the non-relativistic electron-photon interaction that can be derived from the fundamental field theory of quantum electrodynamics (QED).
Explicitly, the QED Lagrangian determines the fundamental electron-photon interaction as follows [95], where ψ(x) is the fermionic field for electrons, A µ (x) is the covariant 4-vector of the electromagnetic (EM) field, γ µ is the Dirac matrix, and is the EM field strength tensor.The first two terms of Eq. ( 65) leads to the free electron and photon Hamiltonians.The last term gives the fundamental electronphoton interaction in the non-relativistic limit (by ignoring positrons and choosing a proper gauge).Thus, the non-relativistic QED Hamiltonian is given by [96] where c † p , c p and b † k , b k are creation and annihilation operators of electrons and photons with momentum p and k.The summations over p and k should be replaced by the continuous integrals d 3 p (2π) 3 and d 3 k (2π) 3 .Also, without loss of generality, we have omitted the indices of electron spin and photon polarization.The first term in the second equality in Eq. ( 66) is the free electron Hamiltonian.The second term is the electron-electron Coulomb interaction arisen from the choice of Coulomb gauge.The third term is the EM field Hamiltonian, and the last two terms are the electron-photon interaction.Note that the electron-phonon interaction in solid-state physics has the same form.Equation (66) can describe most of nonrelativistic physics in the current physics research, unless one is also interested in the phenomena in the smaller scale of nuclear arisen from the weak and strong interactions or the larger scale of universe from gravity.
In the following, we shall find all the nonperturbative renormalization effects on electrons from the electronphoton interaction by using the coherent state path integrals [63] to nonperturbatively trace over all the environmental states.To do so, we may express the exact reduced density ρ S (t) of Eq. ( 63) as Here we have utilized the unnormalized fermion coherent states |ξ ≡ exp( p ξ p c † p )|0 .The integral measure dµ(ξ) = p dξ * p dξ p e −|ξp| 2 is defined the Haar measure in Grassmannian space.The vectors ξ ≡ (ξ p1 , ξ p2 , ...) is a one-column matrix and ξ * pi is a Grassmannian variable.The propagating function J(ξ † f , ξ f , t; ξ 0 , ξ † 0 , t 0 ) in Eq. ( 67), which describes the nonequilibrium evolution of the states of the electron system from the initial state ρ S (t 0 ) to the state at any later time ρ S (t), can be obtained analytically after completing exactly the path integrals over all the photon modes.The result is where D[ξ; ξ ] = p,t<τ <t0 dξ * p (τ )dξ p (τ ), S e [ξ † , ξ] is the bare electron action of the original free-electron Hamil-tonian plus the electron-electron Coulomb interaction in QED.In the fermion coherent state representation, it is given by where H(ξ † (τ ), ξ(τ )) = p ε p ξ * p (τ )ξ p (τ ) + p,p ,q U (q)ξ * p+q (τ )ξ * p −q (τ )ξ p (τ )ξ p (τ ).The additional action S qed IF [ξ † , ξ; ξ † , ξ ] in Eq. ( 68) is an electron action correction arisen from electron-photon interaction after exactly integrated out all the photon modes [94]: This is a generalization of the Feynman-Vernon influence functional [4] to electron-photon interacting systems so that we may also call the action of Eq. ( 70) as the influence functional action.For simplicity, here we have introduced the composite-particle variables σ + p,k (τ ) ≡ ξ * p (τ )ξ p−k (τ ) and σ − p,k (τ ) ≡ ξ * p−k (τ )ξ p (τ ) = (σ + p,k (τ )) † , which correspond to the spin-like variables of the exciton operators a † p a p−k and a † p−k a p , respectively.The nonlocal time correlations in Eq. ( 70) are given by which depict the time-correlations between electrons and photons.The four terms in Eq. ( 70) come from the contributions of the electron-photon interactions to the electron forward propagation, the electron backward propagation, and to the electrons mixed from the forward with backward propagations at the end point time t and at the initial time t 0 , respectively, through the path integrals over all the photon modes.
The above results show that the propagating function Eq. ( 68) of the reduced density matrix for electrons in non-relativistic QED and the corresponding influence functional action Eq. ( 70) have the same form as that for our generalized Fano-Anderson Hamiltonian Eq. ( 38), as shown by Eqs.(B2) and (B4) in Appendix B. The main difference is that the bare system action Eq.(B3) and the influence functional action Eq.(B4) for the general-ized Fano-Anderson Hamiltonian are quadratic with respect to the integrated variables in the path integrals of the propagating function for the reduced density matrix.They can be solved exactly and the resulting propagating function is given by Eq. (B7) in Appendix B.Here the bare electron action Eq. ( 69) and the influence functional action Eq. ( 70) for QED Hamiltonian are highly nonlinear so that the path integrals of the propagating function Eq. ( 68) cannot be carried out exactly.However, the similarity between Eq. ( 70) and Eq.(B4) allows us to find the nonperturbative renormalized electron Hamiltonian in non-rtelativistic QED.
Note that the influence functional action Eq. ( 70) is a complex function.Its real part contains all the corrections to the electron Hamiltonian in non-rtelativistic QED, which results in the renormalization of both the single electron energy and the electron-electron interaction.The imaginary part contains two decoherence sources.One contributes to the energy dissipation into the environment induced by the electron-photon interaction.The other contributes to the fluctuations arisen from the initial states of the thermal photonic reservoirs through the electron-photon interaction.The influence functional action Eq.(B4) shares the same property.Furthermore, it is not difficult to show that the last two terms in both Eqs.(70) and (B4) are pure imaginary so that they only contribute to the dissipation and fluctuation dynamics of the electron system.The first two terms in both Eqs. ( 70) and (B4) can combine with the forward and backward bare system actions in Eq. ( 69) and (B3), respectively, from which we can systematically find the general nonperturbative renormalized Hamiltonian of the system.
Explicitly, let us first examine the generalized Fano-Anderson systems Eq.(38) in Sec.III.The renormalized system Hamiltonian can also be determined by the bare system Hamiltonian function in Eq. (B3) plus the first term in the influence functional action Eq.(B4), i.e., Note that the evolution of α j (τ ) along the forward path is determined by [56][57][58]61] where u jj (τ, t 0 ) is the propagating Green function which obeys the integro-differential Dyson equation Eq. (43a).While, f j (τ ) is the noise source associated with the correlation Green function v(τ, t) of Eq. (43b) that has no contribution to the system Hamiltonian renormalization [58,61].Thus, we can only take the part of the evolution α j (τ ) that has the contribution to system Hamiltonian renormalization, i.e.
The second line in the above expression also shows how the memory effect is taken into account.Using this result, Eq. ( 72) can be rewritten by and This is the same solution for the renormalized system Hamiltonian given by Eqs. ( 41) and (42a) that we obtained after completely solved the propagating function Eq. (B5) and derived the exact master equation Eq. ( 40).Thus, we can find the renormalized electron Hamiltonian in non-relativistic QED from the electron influence functional action Eq. ( 70) in the same way.The result is In Eq. ( 77), the first two terms are the bare electron Hamiltonian in Eq. ( 69), and the last term comes from the the first term in the electron influence functional action Eq. ( 70), as the renormalization effect arisen from the electron-photon interaction after integrated out all the photonic modes.Moreover, we can similarly introduce the two-electron propagating Green function Then the renormalized electron Hamiltonian can be obtained as where U r p (q, t, t 0 ) = U (q) + δU p (q, t, t 0 ), and δU p (q, t, t 0 ) = Im q q t t0 dτ G q (t, τ )W q ,q ,q (τ, t 0 ) × W −1 q ,p ,q (t, t 0 ).(79d) The correction to the single electron energy in Eq. (79b) comes from the operator normal ordering of the renormalized electron-electron interaction in the last term of Eq. (77).By calculating the two-electron propagating Green function W (t, t 0 ) from the total non-relativistic QED Hamiltonian Eq. ( 66), the renormalized electron Hamiltonian and the electron reduced density matrix can be obtained.From the renormalized electron Hamiltonian Eq. ( 79) and the reduced density matrix given by Eq. ( 67)-( 70), the renormalized quantum thermodynamics, Eq. ( 49) to Eq. ( 54) formulated in Sec.III, can be directly extended to complicated interacting open quantum system.Of course, in practical, the two-electron propagating Green function W (t, t 0 ) is very difficult to be calculated, not like the systems of Eq. ( 38) discussed in Sec.III where the general solution of the single particle Green function u(t, t 0 ) has been solved analytically in our previous work [10].Nevertheless, Eqs. ( 67) to (70) provide a full nonequilibrium electron-electron interaction theory rigorously derived from the non-relativistic QED theory after we integrated out exactly all the photonic modes.It can describe various nonequilibrium physical phenomena in many-body electronic systems based on the nonrelativistic QED theory, where all the renormalization effects arisen from electron-photon interaction has been taken into account in the reduced density matrix.In practical, the reduced density matrix of Eq. ( 67) with Eqs.(68) to (70) are still hard to be solved exactly, because the contributions from all the photon modes has been included exactly and it goes far beyond the perturbation expansion one usually used in many-body systems [86] and in quantum field theory [95].In particular, when the Coulomb interaction dominates the electron-electron interaction, the system become a strongly correlated electronic system.Then further nonperturbative approximations and numerical methods have to be introduced to find properly the renormalized Hamiltonian and the reduced density matrix of the open system for the strong coupling quantum thermodynamics.
In fact, the same problem also exists in the equilibrium physics, namely one cannot solve all the equilibrium physical problems exactly even though the Gibbs state is well defined under the equilibrium hypothesis of statistical mechanics.The typical example is the stronglycorrelated electron systems, such as Hubbard model or quantum Heisenberg spin model, which are the approximation of the above nonequilibrium electron-electron interaction QED theory.But so far one is still unable to solve them exactly [65].Therefore, how to practically solve nonequilibrium quantum thermodynamics for arbitrary system-environment interactions remains to be a challenge problem.The closed time-path Green functions technique with the loop expansion to quantum transport phenomena developed by one of us long time ago [97] could be a possible method for solving nonperturbatively the nonequilibrium quantum thermodynamics of strong interacting many-body systems, and we leave this problem for further investigation.

V. CONCLUSION AND PERSPECTIVE
In conclusion, we formulate the renormalization theory of quantum thermodynamics and quantum statistical mechanics based on the exact dynamic evolution of quantum mechanics for both weak and strong coupling strengths.For a class of generally solvable open quantum systems described by the generalized Fano-Anderson Hamiltonians, we show that the exact steady state of open quantum systems coupled to reservoirs through the particle exchange processes is a generalized Gibbs state.The renormalized system Hamiltonian and the reduced density matrix are obtained nonperturbatively when we traced over exactly all the reservoir states through the coherent state path integrals [63].Using the renormalized system Hamiltonian and introducing the renormalized temperature, the exact steady state of the reduced density matrix can be expressed as the standard Gibbs state.The corresponding steady-state particle distributions obey the Bose-Einstein and the Fermi-Dirac distributions for bosonic and fermionic systems, respectively.In the very weak coupling limit, the renormalized sys-tem Hamiltonian and the renormalized temperature are reduced to the original bare Hamiltonian of the system and the initial temperature of the reservoir if it couples to a single reservoir.Thus, classical thermodynamics and statistical mechanics emerge naturally from the dynamics of open quantum systems.Thermodynamic laws and statistical mechanics principle are thereby deduced from the dynamical evolution of quantum dynamics.If open quantum systems contain dissipationless localized bound states, thermalization cannot be reached.This is the solution to the long-standing problem in thermodynamics and statistical mechanics that one has been trying to solve from quantum mechanics for a century.
On the other hand, the renormalization theory presented in this work is nonperturbative.The traditional renormalization theory in quantum field theory and in many-body physics are build on perturbation expansions with respect to the interaction Hamiltonian.As we have systematically shown, the system Hamiltonian renormalization and the reduced density matrix are finally expressed in terms of the nonequilibrium Green functions.We have nonperturbatively derived the equation of motion for these nonequilibrium Green functions and obtained the general nonperturbation solution.We can easily reproduce the traditional perturbation renormalization theory by expanding order by order our solution with respect to the system-reservoir interaction.Furthermore, this nonperturbative renormalization theory also corresponds to the one-step renormalization in the framework of Wilson renormalization group framework.The renormalization group is build through subsequent integrations of physical degrees of freedom from the higher energy scale to lower energy scale.For open quantum systems, instead of integrating out the higher energy degrees of freedom, the dynamics is fully determined by nonperturbatively integrating out all the reservoir degrees of freedom at once but including all energy levels from the low energy scale to the high energy scale of the reservoirs.Therefore, the underlying physical picture of our renormalization roots on the different physical basis.If the open quantum system interacts hierarchically with many reservoirs, then hierarchically tracing over all the reservoirs' states would lead to a new renormalization group theory to open quantum systems that count all influences of hierarchical reservoirs on the system.
As a consequence of the renormalized Hamiltonian and renormalized temperature, we find that the system can become colder or hotter, as the coupling increases.For fermion systems, as the coupling increases, the renormalized energy levels can be increased or decreased, depending on the dot energy levels are greater than or less than the center energy of the Lorentz-type spectral density, but the renormalized temperature is always increased (becomes hotter).For boson systems with the Ohmic-type spectral density, both the renormalized frequency and the renormalized temperature always decrease (becomes colder) as the coupling increases, while for a Lorentz-type spectral density, the renormalized fre-quency and temperature will simultaneously decrease or increase, which is quite different from fermion systems.This reveals the very flexible controllability for energy and heat transfers between systems and reservoirs, and provides potential applications in building quantum heat engines in strong coupling quantum thermodynamics.temperature β α0 = 1/k B T α0 .The up and down signs ∓ correspond to the reservoir being bosonic and fermionic systems, respectively.
After the initial time t 0 , both the system and all reservoirs will evolve into a fully non-equilibrium state.For an arbitrary initial state ρ S (t 0 ) of the system, the reduced density matrix at later time t is defined by ρ S (t) = Tr E [ρ tot (t)], which can be solved from Eq. ( 39) in general.To find the exact solution, we take the coherent state representation [63] again.Then, the reduced density matrix ρ S (t) of Eq. ( 39) can be expressed as ρ S (α t , α t , t) = α t |ρ S (t)|α t .Here we have used the unnormalized coherent state defined as: |α = exp( i α i a † i )|0 , dµ (α) = i g i dα * i dα i e −|αi| 2 , the vector α ≡ (α 1 , α 2 , ...) is onecolumn matrix, and α i are complex variables for bosons and Grassmannian variables for fermions with g i = 1/2πi and 1 in the Haar measure, respectively.Thus, the reduced density matrix in the coherent state representation can be expressed as ρ S (α t , α t , t) = dµ (α 0 ) dµ (α 0 ) ρ S (α 0 , α 0 , t 0 ) ×J(α t , α t , t; α 0 , α 0 , t 0 ) .(B1) The propagating function J α † t , α t , t; α 0 , α 0 † , t 0 in Eq. (B1) can be obtained analytically after integrated exactly over all the environmental degrees of freedom using the coherent state path integrals.The result is [55][56][57] J (α † t , α t , t;α 0 , α 0 † , t s [α † , α ] in Eq. (B2 describe the forward and backward evolution of the system.The influence functional action S IF [α † , α; α † , α ] represents all the influences of the reservoirs on the system after integrated out exactly all the environmental degrees of freedom.This procedure is called as the influence functional approach, proposed originally by Feynman and Vernon [4].We extended the Feynman-Vernon's influence functional in terms of the coherent state path integrals so that the influence functional theory can be applied to both bosonic and fermionic environments [55][56][57].The resulting influence functional action for the open quantum system Eq.( 38) is where the up and down signs of ∓ correspond respectively to the system being bosonic and fermionic.The systemreservoir correlation functions g ij (τ, τ ) and g ij (τ, τ ) are given by Eq. (44).As one can see, both the bare system action and the action correction are quadratic with respect to the variables α * i , α i so that the path integrals of Eq. (B2) can be solved exactly.After a tedious calcula-tion, we obtain [55][56][57] J (α † t ,α t , t; α 0 , α 0 † , t 0 ) = det[w(t)]
Based on the general solution of the nonequilibrium Green functions we obtained recently [10], if there is no localized bound states (modes), the dissipative propagating Green function vanishes in the steady-state limit, namely u(t → ∞, t 0 ) = 0. (B6) This solution is valid for arbitrary continuous spectral density matrix of multiple reservoirs J αij (ω) that cover every point of the whole energy frequency domain.Then the coefficients in the propagating function of the reduced density matrix, Eq. (B5), is largely simplified: J 1 (t, t 0 ) = 0, J 2 (t) = 1−w(t) = ±v(t, t)/[1±v(t, t)] and J 3 (t, t 0 ) = 1.Thus, the propagating function is simply reduced to (B10) This shows that the steady-state reduced density matrix is independent of its initial states, as a consequence of thermalization.Equation (B10) directly results in the operator form of the steady-state reduced density matrix where a ≡ (a 1 , a 2 , • • • , a N ) T is a one-column matrix operator.This is the exact steady-state solution of Eq. ( 47), which is also recently derived from the general solution of the reduced density matrix [18].As one can see, the solution of Eq. ( 19) derived alternatively in Appendix A is a special case of the above general solution.We should also point out that the above solution remains the same for the initial coupled system-reservoir state [58,62].

1 (
FIG. 3.The nonequilibrium dynamical evolution of the internal energy U S (t), the entropy S S (t) and the renormalized nonequilibrium temperature T r (t) for different system initial states |n0 = |5 , |10 and |15 (corresponding to the blue dasheddot line, the grees dashed line and the red solid line, respectively).The left, the midden and the right panels correspond to (a) the weak coupling (η ηc), (b) the strong coupling η < ηc), and (c) the ultra-strong coupling η > ηc).The initial bath temperature T0 = 30 ωs.Other parameters are taken the same as that in Fig.1.

FIG. 4 .
FIG. 4. (a)The diagrammatic spectra of the Hamiltonians of the system, the reservoir and the their interaction.(b) The diagrammatic Dyson expansion of Eq. (8a) in the energy domain, where Σ(ω) = dω J(ω ) ω−ω is the self-energy arisen from the coupling between the system and the reservoir, and J(ω) ≡ k |V k | 2 δ(ω − ω k ).The renormalized system energy ω r s and the dissipation coefficient γ of Eqs.(7a)-(7b) are determined nonperturbatively from Eq. (8a) with a Laplace transformation, which contains all order contributions upto the infinite orders from the system-reservoir coupling Hamiltonian, as shown in this diagrammatic expansion.

FIG. 5 .
FIG. 5. (a)-(b)The nonequilinrium evolution of quantum work and quantum heat changings with respect to the time, dW S (t)/dt and dQ S (t)/dt (in the unit of ω 2 s ), for different coupling strengths.(c) The steady-state specific heat as a function of the renormalized temperature calculated from the derivative of the internal energy with respect to the renormalized temperature Eq. (34) (red lines) and from the partition function given in the Gibbs state Eq.(24) for different initial temperatures.The dashed-dot, dot and dashed lines correspond to the different coupling strengths η/ηc = 0.3, 0.5, 0.8, respectively.Other parameters are taken the same as that in Fig.1.

< l a t e x i t s h a 1 _
b a s e 6 4 = " N Y Z t 0 m Z J m I n D 9 b X u u i C z n o H O T Z 4 = " > A A A B 7 H i c b V B N S 8 N A E J 3 U r 1 q / q h 6 9 L B b B U 0 l E 1 G P R i 8 c q T V t oS 9 l s J + 3 S z S b s b o Q S + h u 8 e F D E q z / I m / / G b Z u D t j 4 Y e L w 3 w 8 y 8 I B F c G 9 f 9 d g p r 6 x u b W 8 X t 0 s 7 u 3 v 5 w D K / w 5 k j n x X l 3 P h a t B S e f O Y Y / c D 5 / A M 4 5 j r I = < / l a t e x i t > T R < l a t e x i t s h a 1 _ b a s e 6 4 = " 2 X 0 N W Y 4 f 7 4 b y + H q P U U B E V C I R t t A = " > A A A B 7 n i c b V D L S g N B E O y N r x h f U Y 9 e B o P g K e x K U I 9 B L x 6 j m A c k S 5 i d z C Z D 5 r H M z A p h y U d 4 8 a C y c O f e e M 3 P n 2 p H n x s I w X i a 0 y a n p m d n M 3 P z C 4 t L y S n Z 1 7 SI O E + 6 w i h N 6 I a / Z V s w 8 N 2 A V 4 Q q P 1 S L O L N / 2 W N X u H s t 4 9 Y b x 2 A 2 D c 9 G L W N O 3 r g K 3 4 z q W I O q y 4 S e t f k O w W 9 E / G Q x a 2 Z y R N 9 T Q x 4 G Z g l x h 7 3 0 j 2 D k t l s L s M x p o I 4 S D B D 4 Y A g j C H i z E 9 N V h w k B E X B N 9 4 j g h V 8 U Z B p g n b U J Z j D I s Y r s 0 X 9 G u n r I B 7 a V n r N Q O n e L R z0 m p Y 5 s 0 I e V x w v I 0 X c U T 5 S z Z 3 7 z 7 y l P e r U e r n X r 5 x A p c E / u X b p T 5 X 5 2 s R a C D I 1 W D S z V F i p H V O a l L o l 5 F 3 l z / U p U g h 4 g 4 i d s U 5 4 Q d p R y 9 s 6 4 0 s a p d v q 2 l 4 q 8 q U 7 J y 7 6 S 5 C d 7 k L a n B 5 s 9 2 j o O L / b x 5 k D f L 1 O n y E G p k s I k t 7 F I / D 1 F A E S V U y J v j A Y 9 4 0 s 6 0 n j b U 7 j 5 T t Y n P F e v 4 N r T 7 D 3 H q m K g = < / l a t e x i t > µ L < l a t e x i t s h a 1 _ b a s e 6 4 = " h M R b d I p D f 0 B K F X n S e J Q X E B a n q T I e S W f B / u Y d S 0 9 x t x G t Z u I 1 J J b j j N i / d J P M / + p E L R w 9 H M k a b K r J l 4 y o z k p c I v k q 4 u b q l 6 o 4 O f j E C d y l e E D Y k s r J O 6 t S E 8 r a x d s a M v 4 m M w U r 9 l a S G + F d 3 J I a r P 9 s 5 z S o 7 R f 0 g 4 J e p k 6 X r y F H G t v Y w R 7 1 8 x B F n K C E K n n 7 u M U d 7 p W y c q F c K U m q k v p c s Y l v Q 7 n 5 A N 4 f m I w = < / l a t e x i t > T L < l a t e x i t s h a 1 _ b a s e 6 4 = " J g / 5 l 5 o D S x p E s p T h M 1 p r o 3 n 1 t m w v 3 m 3 r a e 5 2 w X 9 g 8 I r I l b j j Ni / d B + Z / 9 W Z W j R O s W x r k F R T a h l T H S 9 c c t s V c 3 P 3 U 1 W a H F L i D D 6 h u C L M r f K j z 6 7 V Z L Z 2 0 1 t m 4 6 8 2 0 7 B m z 4 v c H G / m l j R g / / s 4 f 4 L d + a q / W F 3 Y p k m v o r N K m M I0 5 m i e S 1 j B B m q o k 3 e M W 9 z h 3 t l y c u f S u e 6 k O l 2 F Z g J f l n P z D u w 3 l u c = < / l a t e x i t > L < l a t e x i t s h a 1 _ b a s e 6 4 = " l GI F y t G e / G Z v N 0 x F l a V D 1 t 4 V J c I = " > A A A C B X i c b Z D L S s N A F I Y n 9 V b r L e p S F 4 N F q J u S S F G X R T c u K 9 g L t K F M p p N 2 6 O T C z I l Y Q j Z u f B U 3 L h Rx 6 z u 4 8 2 2 c p h G 0 9 Y e B j / + c w 5 n z u 5 H g C i z r y y g s L a + s r h X X S x u b W 9 s 7 5 u 5 e S 4 W x 0 ) = D[α; α ] exp i S s [α † , α] − S * s [α † , α ] + S IF [α † , α; α † , α ] , (B2)where D[α; α ] is the path integral measure over the parameter space of the system coherent states |α .The bare system action in the coherent state basis is given byS s [α † , α] = − i 2 α † t α(t) + α † (t 0 )α t0 τ )α(τ ) − α † (τ ) α(τ ) − H(α † (τ ), α(τ )) (B3)with the Hamiltonian functionH(α † (τ ), α(τ )) = i ε i (τ )α * i (τ )α i (τ ).The actions S s [α † , α] and S *