# Nonperturbative renormalization of quantum thermodynamics
from weak to strong couplings

Wei-Ming Huang  Department of Physics and Center for Quantum Information
Science, National Cheng Kung University, Tainan 70101, Taiwan  Wei-Min Zhang
Correspondence author<EMAIL_ADDRESS>Department of Physics and
Center for Quantum Information Science, National Cheng Kung University, Tainan
70101, Taiwan  Physics Division, National Center for Theoretical Sciences,
Taipei 10617, Taiwan

(This revised manuscript was compiled on Mar. 25, 2022)

###### Abstract

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 the systems in arbitrary initial states as well as 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. In the very weak system-
reservoir coupling limit, the renormalized system Hamiltonian and the
renormalized temperature are reduced to the original bare Hamiltonian of the
system and the initial temperature of the reservoir. Thus, the conventional
statistical mechanics and thermodynamics are thereby rigorously deduced from
quantum dynamical evolution. In the last, this nonperturbative renormalization
method is also extended to general interacting open quantum systems.

Quantum thermodynamics $|$ Quantum statistical mechanics $|$ Open quantum
systems $|$ Exact master equation $|$ Nonequilibrium dynamics

##  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
Landau1969 ; Kubo1991 , 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 Huang1987 : (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 Huang1987 ; Landau1969 ; Kubo1991 ; Feynman1963 ; Leggett1983a ;
Srednicki94 ; Zurek2003 ; Gemmer2004 ; Jarzynski2011 ; Zhang2012 ; Kosloff2013
; Nandkishore15 ; Xiong2015 ; MillenNJP2016 ; EspositoNJP17 ; Binder2018 ;
Deffner2019 ; Xiong2020 ; Talkner2020 . 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 open 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
Leggett1983a . They used the Feynman-Vernon influence functional approach
Feynman1963  to explore the dynamics of quantum Brownian motion, and found the
equilibrium thermal state approximately Leggett1983a . Later, Zurek studied
extensively this nontrivial problem from the quantum-to-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 Zurek2003 . Thermalization in these investigations is
demonstrated for quantum Brownian motion for initial Gaussian wave packets at
high temperature limit Leggett1983a ; Zurek2003 . 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.

On the other hand, in the last two decades, experimental investigations on
nano-scale quantum heat engines have attracted tremendous attentions on the
realization of thermalization and the formulation of quantum thermodynamics
Allahverdyan2000 ; Scully2003 ; Esposito2010a ; Scully2011 ; Trotzky2012 ;
Jezouin2013 ; KZhang2014 ; Bergenfeldt2014 ; Langen2015 ; Pekola2015 ;
David2015 ; Kaufman2016 ; Anders2016 ; Ronagel2016 ; Ochoa15 . Besides
searching new thermal phenomena arisen from quantum coherence and quantum
entanglement, an interesting question also appeared naturally is what happen
when microscopic systems couple strongly with reservoirs. Since then, many
effects have been devoted on the problems of how thermodynamic laws emerge
from quantum dynamics and how these laws may be changed when the system-
reservoir couplings become strong Campisi2009 ; Subas2012 ; Esposito2015 ;
Seifert2016 ; Carrega2016 ; Ochoa2016 ; Jarzynski2017 ; Marcantoni2017 ;
Bruch2018 ; Perarnau2018 ; Hsiang2018 ; Anders2018 ; Strasberg2019 ;
Newman2020 ; Ali2020a ; Rivas2020 . In particular, how to properly define the
thermodynamic work and heat in the quantum mechanical framework becomes an
important issue when the system and reservoirs strongly coupled together
Esposito2010b ; Binder2015 ; Esposito2015 ; Alipour2016 ; Strasberg2017 ;
Perarnau2018 . Due to various assumptions and approximations one inevitably
taken in addressing the open quantum system dynamics, no consensus has been
reached in building quantum thermodynamics at strong coupling.

In the last decade, we have derived the exact master equation of open quantum
systems Tu2008 ; Jin2010 ; Lei2012 ; Yang2015 ; Yang2017 ; Zhang2018 ; Lai2018
; Huang2020  by extending the Feynman-Vernon influence functional theory into
the coherent state representation Zhang1990 . 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
Zhang2012 ; Xiong2015 ; Xiong2020 . 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 Seifert2016 ; Carrega2016 ; Ochoa2016 ;
Jarzynski2017 ; Marcantoni2017 ; Bruch2018 ; Perarnau2018 ; Hsiang2018 ;
Anders2018 ; Strasberg2019 ; Newman2020 ; Ali2020a ; Rivas2020 : (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 system-reservoir 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
Feynman1963  and by Caldeira and Leggett Leggett1983a ; Leggett1983b . From
the exact master equation and its exact solution in the steady state for such
class of open quantum systems Tu2008 ; Jin2010 ; Lei2012 , 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 Seifert2016 ; Carrega2016 ; Ochoa2016 ; Jarzynski2017 ;
Marcantoni2017 ; Bruch2018 ; Perarnau2018 ; Hsiang2018 ; Anders2018 ;
Strasberg2019 ; Newman2020 ; Rivas2020  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 reduced 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 two-electron
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 Nagaosa2010  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 Fano1961 :

$\displaystyle H_{\rm tot}\\!=$ $\displaystyle
H_{{}_{\\!S}}\\!+\\!H_{{}_{\\!E}}\\!+\\!H_{{{}_{\\!S}}{{}_{\\!E}}}$
$\displaystyle=$
$\displaystyle\hbar\omega_{s}a^{\dagger}a\\!+\\!{\sum}_{k}\hbar\omega_{k}b^{\dagger}_{k}b_{k}\\!+\\!{\sum}_{k}\hbar(V_{k}a^{\dagger}b_{k}\\!+\\!V_{k}^{*}b^{\dagger}_{k}a),$
(1)

where $a^{\dagger}$ and $b^{\dagger}_{k}$ ($a$ and $b_{k}$) are the creation
(annihilation) operators of the bosonic modes in the system and in the
reservoir with energy quanta $\hbar\omega_{s}$ and $\hbar\omega_{k}$,
respectively. They obey the standard bosonic commutation relations:
$[a,a^{\dagger}]=1$ and
$[b_{k},b^{\dagger}_{k^{\prime}}]=\delta_{kk^{\prime}}$, etc. The parameter
$V_{k}$ is the coupling amplitude between the system and the reservoir and can
be experimentally tuned to strong coupling Putz14 ; Chiang21 . 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 Miroshnichenko2010 .

###  II.1 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

$\displaystyle\rho_{{}_{\\!E}}(t_{0})\\!=\\!e^{-\beta_{0}H_{{}_{\\!E}}}/Z_{{}_{\\!E}},$
(2)

where $\beta_{0}\\!=\\!1/k_{B}T_{0}$ and $T_{0}$ is the initial temperature of
the reservoir, $Z_{{}_{\\!E}}\\!=\\!{\rm
Tr}_{{}_{\\!E}}[e^{-\beta_{0}H_{{}_{\\!E}}}]$ is its partition function. The
system can be initially in arbitrary state $\rho_{{}_{\\!S}}(t_{0})$ so that
the initial total density matrix of the system plus the reservoir is a direct
product state Feynman1963 ; Leggett1983a ,

$\displaystyle\rho_{\rm
tot}(t_{0})=\rho_{{}_{\\!S}}(t_{0})\otimes\frac{\\!e^{-\beta_{0}H_{{}_{\\!E}}}}{Z_{{}_{\\!E}}}.$
(3)

After the initial time $t_{0}$, both the system and the reservoir evolve into
an entangled nonequilibrium state $\rho_{\rm tot}(t)$ which obeys the
Liouville-von Neumann equation in quantum mechanics Neumann55 ,

$\displaystyle\frac{d}{dt}\rho_{\rm tot}(t)=\frac{1}{i\hbar}[H_{\rm
tot},\rho_{\rm tot}(t)].$ (4)

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 $\rho_{{}_{\\!S}}(t)$. It is defined by the partial trace over all the
reservoir states:

$\displaystyle\rho_{{}_{\\!S}}(t)\\!=\\!{\rm Tr}_{{}_{\\!E}}[\rho_{\rm
tot}(t)].$ (5)

The equation of motion for $\rho_{{}_{\\!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. Lindblad1976 ;
GKS1976 ; Kolodynski2018 ; Breuer2008 ; Paavola2009 ; Rajesh2015 . 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 Tu2008 ; Jin2010 ; Lei2012 ; Yang2015 ; Yang2017 ;
Lai2018 ; Zhang2018 ; Huang2020 . Explicitly, we have derived the exact master
equation for Eq. (1) by exactly tracing over all the reservoir states from the
solution of the Liouville-von Neumann equation Wu2010 ; Xiong2010 ; Lei2011 ;
Lei2012 . 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 Lei2012 ; Zhang1990 . The resulting exact master equation for the
reduced density matrix accompanied with the renormalized system Hamiltonian is
given by

$\displaystyle\frac{d}{dt}\rho_{{}_{\\!S}}$
$\displaystyle(t)\\!=\frac{1}{i\hbar}\big{[}H^{r}_{{}_{\\!S}}(t),\rho_{{}_{\\!S}}(t)\big{]}$
$\displaystyle\\!+\\!\gamma(t,t_{0})\big{\\{}2a\rho_{{}_{\\!S}}(t)a^{\dagger}\\!-\\!a^{\dagger}a\rho_{{}_{\\!S}}(t)\\!-\\!\rho_{{}_{\\!S}}(t)a^{\dagger}a\big{\\}}$
$\displaystyle\\!+\\!\widetilde{\gamma}(t,t_{0})\big{\\{}a^{\dagger}\rho_{{}_{\\!S}}(t)a\\!+\\!a\rho_{{}_{\\!S}}(t)a^{\dagger}\\!-\\!a^{\dagger}a\rho_{{}_{\\!S}}(t)\\!-\\!\rho_{{}_{\\!S}}(t)aa^{\dagger}\big{\\}}.$
(6a) In this exact master equation, the first term describes a unitary
evolution of the reduced density matrix with the renormalized Hamiltonian
$\displaystyle H^{r}_{{}_{\\!S}}(t)=\hbar\omega^{r}_{s}(t,t_{0})a^{\dagger}a.$
(6b)

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 $\omega^{r}_{s}(t,t_{0})$, the non-Markovian
dissipation coefficient $\gamma(t,t_{0})$ and the non-Markovian fluctuation
coefficient $\widetilde{\gamma}(t,t_{0})$, respectively. All these non-
Markovian coefficients are nonperturbatively and exactly determined by the
following relations Wu2010 ; Xiong2010 ; Lei2011 ; Lei2012

$\displaystyle\omega^{r}_{s}(t,t_{0})=-{\rm Im}[\dot{u}(t,t_{0})/u(t,t_{0})],$
(7a) $\displaystyle\gamma(t,t_{0})=-{\rm Re}[\dot{u}(t,t_{0})/u(t,t_{0})],$
(7b) $\displaystyle\tilde{\gamma}(t,t_{0})=\dot{v}(t,t)-2v(t,t){\rm
Re}[\dot{u}(t,t_{0})/u(t,t_{0})].$ (7c)

Here $u(t,t_{0})$ and $v(\tau,t)$ are the two non-equilibrium Green functions
obeying the integro-differential Dyson equations,

$\displaystyle\frac{d}{dt}u(t,t_{0})\\!+\\!i\omega_{s}u(t,t_{0})\\!+\\!\\!\int^{t}_{t_{0}}\\!\\!d\tau
g(t,\tau)u(\tau,t_{0})=0,$ (8a) $\displaystyle
v(\tau,t)=\\!\int^{\tau}_{t_{0}}\\!\\!d\tau_{1}\\!\int^{t}_{t_{0}}\\!\\!d\tau_{2}u(\tau,\tau_{1})\widetilde{g}(\tau_{1},\tau_{2})u^{*}(t,\tau_{2}).$
(8b)

The non-Markovianity is manifested by the above time-convolution equation of
motion for these non-equilibrium Green functions. The integral kernels in the
above convolution equations are given by

$\displaystyle g(t,\tau)=\\!\\!\int^{\infty}_{0}\\!\\!\\!d\omega
J(\omega)e^{-i\omega(t-\tau)},$ (9a)
$\displaystyle\widetilde{g}(\tau_{1},\tau_{2})=\\!\\!\int^{\infty}_{0}\\!\\!\\!d\omega
J(\omega)\overline{n}(\omega,T_{0})e^{-i\omega(\tau_{1}-\tau_{2})},$ (9b)

which characterize the time correlations between the system and the reservoir
through the system-reservoir interaction. The frequency-dependent function

$\displaystyle J(\omega)\equiv\sum_{k}|V_{k}|^{2}\delta(\omega-\omega_{k})$
(10)

is called as the spectral density, which fully encapsulates the fundamental
dissipation (relaxation) and fluctuation (noise or dephasing) effects induced
by the system-reservoir interaction. Finally, the initial temperature
dependent function

$\displaystyle\overline{n}(\omega_{k},T_{0})=\\!{\rm
Tr}_{{}_{\\!E}}[b^{\dagger}_{k}b_{k}\rho_{{}_{\\!E}}(t_{0})]=1/[e^{\hbar\omega_{k}/k_{B}T_{0}}-1]$
(11)

is the initial particle distribution in the reservoir.

An arbitrary initial state of the system can be expressed as

$\displaystyle\rho_{{}_{\\!S}}(t_{0})\\!=\\!\sum^{\infty}_{m,m^{\prime}=0}\rho_{mm^{\prime}}|m\rangle\langle
m^{\prime}|,$ (12)

where $|m\rangle=\frac{1}{\sqrt{m!}}(a^{\dagger})^{m}|0\rangle$ is the bosonic
Fock state. If $\rho_{mm^{\prime}}=c_{m}c^{*}_{m^{\prime}}$, then
$\rho_{{}_{\\!S}}(t_{0})$ is a pure state, otherwise it is a mixed state. The
exact solution of the exact master equation Eq. (6a) has be found Xiong2015 ;
Xiong2020 ,

$\displaystyle\rho^{\rm exact}_{{}_{\\!S}}(t)=$
$\displaystyle\\!\\!\\!\\!\sum^{\infty}_{m,m^{\prime}=0}\\!\\!\\!\rho_{mm^{\prime}}\\!\\!\\!\\!\\!\\!\\!\sum_{k=0}^{\rm
min\\{m,m^{\prime}\\}}\\!\\!\\!\\!\\!\\!\\!d_{k}(t)A^{\dagger}_{mk}(t)\widetilde{\rho}[v(t,t)]A_{m^{\prime}k}(t),$
(13)

where

$\displaystyle\widetilde{\rho}[v(t,t)]=\sum_{n=0}^{\infty}\frac{[v(t,t)]^{n}}{[1+v(t,t)]^{n+1}}|n\rangle\langle
n|,$ (14a) $\displaystyle
A^{\dagger}_{mk}(t)=\frac{\sqrt{m!}}{(m-k)!\sqrt{k!}}\Big{[}\frac{u(t,t_{0})}{1+v(t,t)}a^{\dagger}\Big{]}^{m-k},$
(14b) $\displaystyle
d_{k}(t)=\\!\Big{[}1\\!-\\!\frac{|u(t,t_{0})|^{2}}{1+v(t,t)}\Big{]}^{k}.$
(14c)

As a self-consistent check of the above solution, we calculate the average
particle number in the system from the above solution,
$\overline{n}(t)\equiv{\rm Tr}_{{}_{\\!S}}[a^{\dagger}a\rho_{{}_{\\!S}}(t)]$,
also using the Heisenberg equation of motion directly,
$\overline{n}(t)\equiv{\rm
Tr}_{{{}_{\\!S}}+{{}_{\\!E}}}[a^{\dagger}(t)a(t)\rho_{\rm tot}(t_{0})]$. Both
calculations give the same result Jin2010 ; Lei2012 ; Yang2015 ; Zhang2018 :

$\displaystyle\overline{n}(t)$ $\displaystyle={\rm
Tr}_{{}_{\\!S}}[a^{\dagger}a\rho_{{}_{\\!S}}(t)]={\rm
Tr}_{{{}_{\\!S}}+{{}_{\\!E}}}[a^{\dagger}(t)a(t)\rho_{\rm tot}(t_{0})]$
$\displaystyle=u^{*}(t,t_{0})\overline{n}(t_{0})u(t.t_{0})+v(t,t),$ (15)

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(\omega)$,
if no localized bound state exists Zhang2012 ; Expl , the general solution of
Eq. (8a) is

$\displaystyle u(t,t_{0})=\int_{0}^{\infty}\\!\\!\\!\\!d\omega
D(\omega)e^{i\omega(t-t_{0})}\xrightarrow{t\rightarrow\infty}0$ (16)

where
$D(\omega)\\!=\\!\frac{J(\omega)}{[\omega-\omega_{s}-\Delta(\omega)]^{2}+\pi^{2}J^{2}(\omega)}$
shows the system spectrum broadening due to the coupling to the reservoir, and
$\Delta(\omega)={\cal P}\big{[}\\!\int
d\omega^{\prime}\frac{J(\omega^{\prime})}{\omega-\omega^{\prime}}\big{]}$ is
the principal-value integral of the self-energy correction to the system,
$\Sigma(\omega)\\!=\\!\int
d\omega^{\prime}\frac{J(\omega^{\prime})}{\omega-\omega^{\prime}}\\!=\\!\Delta(\omega)-i\pi
J(\omega)$. In fact, $\Delta(\omega)$ gives the system frequency (or energy)
shift. As a result, in the steady state limit, we have
$A^{\dagger}_{mk}(t)\xrightarrow{t\rightarrow\infty}\delta_{mk}$ and
$d_{k}(t)\xrightarrow{t\rightarrow\infty}1$. Then the exact solution of the
particle distribution and the reduced density matrix are reduced to

$\displaystyle\overline{n}_{\rm exact}\\!(t\\!\rightarrow\\!\infty)$
$\displaystyle=\lim_{t\rightarrow\infty}v(t,t)=\\!\\!\int^{\infty}_{0}\\!\\!\\!\\!d\omega
D(\omega)\overline{n}(\omega,T_{0}),$ (17a) $\displaystyle\rho^{\rm
exact}_{{}_{\\!S}}\\!(t\\!\rightarrow\\!\infty)$
$\displaystyle=\lim_{t\rightarrow\infty}\sum_{n=0}^{\infty}\frac{[v(t,t)]^{n}}{[1+v(t,t)]^{n+1}}|n\rangle\langle
n|$ (17b)
$\displaystyle=\lim_{t\rightarrow\infty}\frac{\exp\big{\\{}\ln\\!\big{[}\\!\frac{v(t,t)}{1+v(t,t)}\big{]}a^{\dagger}a\big{\\}}}{1+v(t,t)}.$
(17c)

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 system-reservoir
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 Xiong2015 ; Xiong2020 .

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(\tau,t)$,
as we have shown in Refs. Tan2011 ; Huang2020 ; Yang2015 . This exact master
equation formalism has also been extended to open quantum systems including
external deriving fields Lei2012 ; Chiang21 .

###  II.2 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
$\hbar\omega_{s}$ to $\hbar\omega^{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
$\omega^{r}_{s}(t)$ and its steady-state value
$\omega^{r}_{s}=\omega^{r}_{s}(t\\!\rightarrow\\!\infty)$ can be exactly
calculated from Eqs. (7a) and (8a). Here we take the Ohmic spectral density
$J(\omega)\\!=\\!\eta\omega\exp(-\omega/\omega_{c})$ Leggett1987  in the
practical calculation. The result is presented in Fig. 1(a) and (b). It shows
that different system-reservoir coupling strengths $\eta$ 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 $\eta/\eta_{c}$, where $\eta_{c}=\omega_{s}/\omega_{c}$ is a critical
coupling strength for the Ohmic spectral density Zhang2012 ; Xiong2015 . When
$\eta>\eta_{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 Xiong2015 ; Xiong2020
.

Figure 1: (a) The renormalized system energy (cavity frequency shift)
$\hbar\omega^{r}_{s}(t)$ for three different system-reservoir coupling
strengths, $\eta=0.01\eta_{c},0.5\eta_{c},0.9\eta_{c}$. It is calculated from
Eqs. (7a) and (8a) for the Ohmic spectral density
$J(\omega)\\!=\\!\eta\omega\exp(-\omega/\omega_{c})$, where the cutoff
frequency $\omega_{c}=10\omega_{s}$ is taken, and
$\eta_{c}\\!=\\!\omega_{s}/\omega_{c}$ is a critical coupling for the Ohmic
spectral density Zhang2012 ; Xiong2015 . (b) The steady-state renormalized
frequency shift $\omega^{r}_{s}=\omega^{r}_{s}(t\\!\rightarrow\\!\infty)$ as a
function of the system-reservoir coupling strength $\eta/\eta_{c}$. (c) The
steady-state particle distribution $\overline{n}_{\rm
exact}(t\\!\rightarrow\\!\infty)$ of Eq. (17a) (the blue-dashed line), the
Bose-Einstein distribution without the energy (frequency) renormalization
$\overline{n}(\omega_{s},T_{0})$ (the black-dot line) and with the energy
renormalization $\overline{n}(\omega^{r}_{s},T_{0})$ (the green-dashed-dot
line), respectively. The system is initially set in a pure Fock state
$|n_{0}\rangle$ with $n_{0}=5$, and the reservoir initial temperature
$T_{0}=10\hbar\omega_{s}$.

In Fig. 1(c), we plot the exact solution $\overline{n}_{\rm
exact}(t\\!\rightarrow\\!\infty)$ of Eq. (17a) as a function of the coupling
strength $\eta/\eta_{c}$ (the blue-dashed line). We compare the result with
the Bose-Einstein distribution without the energy (frequency) renormalization,
$\overline{n}(\omega_{s},T_{0})=1/[e^{\hbar\omega_{s}/k_{B}T_{0}}-1]$ (see the
black-dot line), also compare to Bose-Einstein distribution with the energy
renormalization,
$\overline{n}(\omega^{r}_{s},T_{0})=1/[e^{\hbar\omega^{r}_{s}/k_{B}T_{0}}-1]$
(see the green-dashed-dot line). As one can see, the exact solution
$\overline{n}_{\rm exact}(t\\!\rightarrow\\!\infty)$ derivates significantly
from the Bose-Einstein distribution without the energy renormalization, i.e.,
$\overline{n}(\omega_{s},T_{0})$, as $\eta$ 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 $\overline{n}(\omega^{r}_{s},T_{0})$,
changes with the changes of $\eta$, similar to the exact solution
$\overline{n}_{\rm exact}(t\\!\rightarrow\\!\infty)$. But there is still a
quantitative disagreement between the exact solution $\overline{n}_{\rm
exact}(t\\!\rightarrow\\!\infty)$ and the Bose-Einstein distribution
$\overline{n}(\omega^{r}_{s},T_{0})$ with the renormalized cavity photon
energy $\hbar\omega^{r}_{s}$.

To understand further the origin of the above difference, let us recall that
the exact solution $\rho^{\rm exact}_{{}_{\\!S}}\\!(t\\!\rightarrow\\!\infty)$
of Eq. (17c) is indeed a Gibbs-type state. This indicates that the exact
particle distribution $\overline{n}_{\rm exact}(t\\!\rightarrow\\!\infty)$
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 Seifert2016 ; Jarzynski2017 . 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 $\rho_{\rm 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

$\displaystyle\rho_{\rm tot}(t\rightarrow\infty)=\frac{1}{Z_{\rm
tot}}e^{-\beta_{f}H_{\rm tot}},$ (18)

where $\beta_{f}=1/k_{B}T_{f}$, and $H_{\rm 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 Huang2020
that (also see the detailed derivation given in Appendix A)

$\displaystyle\rho_{{}_{\\!S}}(t\rightarrow\infty)\\!$ $\displaystyle=\\!{\rm
Tr}_{{}_{\\!E}}\Big{[}\frac{e^{-\beta^{r}H_{\rm tot}}}{Z_{\rm tot}}\Big{]}$
$\displaystyle=\\!\frac{\exp\big{\\{}\ln\\!\big{[}\\!\frac{\overline{n}(t\\!\rightarrow\\!\infty)}{1+\overline{n}(t\\!\rightarrow\\!\infty)}\big{]}a^{\dagger}a\big{\\}}}{1+\overline{n}(t\\!\rightarrow\\!\infty)}.$
(19)

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
$\overline{n}(t\\!\rightarrow\\!\infty)={\rm
Tr}_{{{}_{\\!S}}+{{}_{\\!E}}}[a^{\dagger}a\rho_{\rm
tot}(t\\!\rightarrow\\!\infty)]={\rm
Tr}_{{}_{\\!S}}[a^{\dagger}a\rho_{{}_{\\!S}}(t\\!\rightarrow\\!\infty)]=\overline{n}_{\rm
exact}(t\\!\rightarrow\\!\infty)$. This gives a further self-consist
justification to the above conclusion.

The result presented in Fig. 1(c) shows that
$\overline{n}(\omega^{r}_{s},T_{0})\neq\overline{n}_{\rm
exact}(t\\!\rightarrow\\!\infty)$ except for the weak coupling. This indicates
that in general, $T_{f}\neq 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 Callen1985 , 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 Callen1985 .

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
$\rho_{{}_{\\!S}}(t)$ of Eq. (13):

$\displaystyle U_{{}_{\\!S}}(t)\equiv{\rm
Tr}_{{}_{\\!S}}[H^{r}_{{}_{\\!S}}\\!(t)\rho_{{}_{\\!S}}(t)].$ (20)

Also, we define the von Neumann entropy with the exact reduced density matrix
of Eq. (13) as the nonequilibrium thermodynamic entropy of the system Ali2020a
; Neumann55 ; Callen1985 :

$\displaystyle S_{{}_{\\!S}}(t)=\\!-k_{B}{\rm
Tr}_{{}_{\\!S}}[\rho_{{}_{\\!S}}(t)\ln\rho_{{}_{\\!S}}(t)],$ (21)

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 Ali2020
which is defined as

$\displaystyle T^{r}\\!(t)\equiv\frac{\partial U_{{}_{\\!S}}(t)}{\partial
S_{{}_{\\!S}}(t)}\bigg{|}_{\omega^{r}_{s}}\\!\\!={\rm
Tr}_{{}_{\\!S}}\bigg{[}H^{r}_{{}_{\\!S}}\\!(t)\frac{d\rho_{{}_{\\!S}}(t)}{dS_{{}_{\\!S}}(t)}\bigg{]}.$
(22)

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 system-reservoir
coupling is strong or weak, we can fundamentally obtain the final equilibrium
temperature $T_{f}\equiv T^{r}\\!=\\!T^{r}(t\\!\rightarrow\\!\infty)$ from the
dynamical evolution of the open quantum system.

Figure 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
$\eta/\eta_{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 $\omega^{r}_{s}$ and renormalized temperature $T^{r}$ as a function
of the coupling strength $\eta$. (e) The steady-state particle distribution as
a function of coupling strength $\eta$. It shows that the exact solution
$\overline{n}_{\rm exact}(t\\!\rightarrow\\!\infty)$ of Eq. (17a) (the blue-
dashed line) and the Bose-Einstein distribution
$\overline{n}(\omega^{r}_{s},T^{r})$ with the renormalized frequency and the
renormalized temperature (the red-dot line) agree perfectly to each other. The
green-dashed-dot line is $\overline{n}(\omega^{r}_{s},T_{0})$ 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 $\rho^{\rm
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 system-reservoir
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
Seifert2016 ; Carrega2016 ; Ochoa2016 ; Jarzynski2017 ; Marcantoni2017 ;
Bruch2018 ; Perarnau2018 ; Hsiang2018 ; Anders2018 ; Strasberg2019 ;
Newman2020 ; Rivas2020 .

In Fig. 2(d), we plot the steady-state renormalized temperature,
$T^{r}\\!=\\!T^{r}(t\\!\rightarrow\\!\infty)$, as a function of the coupling
astrength $\eta/\eta_{c}$. Using this renormalized temperature, we further
plot the Bose-Einstein distribution with both the renormalized energy and the
renormalized temperature:
$\overline{n}(\omega^{r}_{s},T^{r})=1/[e^{\hbar\omega^{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.,

$\displaystyle\overline{n}_{\rm
exact}(t\\!\rightarrow\\!\infty)=\overline{n}(\omega^{r}_{s},T^{r})=\frac{1}{e^{\hbar\omega^{r}_{s}/k_{B}T^{r}}-1}.$
(23)

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,

$\displaystyle\rho^{\rm exact}_{{}_{\\!S}}\\!(t\\!\rightarrow\\!\infty)\\!$
$\displaystyle=\\!\\!\sum_{n=0}^{\infty}\\!\frac{[\overline{n}(\omega^{r}_{s},T^{r})]^{n}}{[1+\overline{n}(\omega^{r}_{s},T^{r})]^{n+1}}|n\rangle\langle
n|$ $\displaystyle=\\!\frac{1}{Z^{r}_{{}_{\\!S}}}e^{-\beta^{r}\\!H^{r}_{s}},$
(24)

where $Z^{r}_{{}_{\\!S}}\\!=\\!{\rm
Tr}_{{}_{\\!S}}[e^{-\beta^{r}H^{r}_{{}_{\\!S}}}]$ is the renormalized
partition function, and $\beta^{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
$\eta\\!\ll\\!\eta_{c}$, $\Delta(\omega)\\!\rightarrow\\!0$ and
$D(\omega)\\!\rightarrow\\!\delta(\omega-\omega_{s})$ so that the steady state
solution of Eq. [17a] is directly reduced to $\overline{n}(\omega_{s},T_{0})$
Xiong2015 ; Xiong2020 , and

$\displaystyle\rho^{\rm
exact}_{{}_{\\!S}}\\!(t\\!\rightarrow\\!\infty)\\!\xrightarrow{\eta\ll\eta_{c}}$
$\displaystyle\sum_{n=0}^{\infty}\\!\frac{[\overline{n}(\omega_{s},T_{0})]^{n}}{[1\\!+\\!\overline{n}(\omega_{s},T_{0})]^{n+1}}|n\rangle\langle
n|$ $\displaystyle=\\!\frac{1}{Z_{{}_{\\!S}}}e^{-\beta_{0}H_{{}_{\\!S}}}.$
(25)

This reproduces the expected solution of the statistical mechanics in the weak
coupling regime. Figures 1 and 2 also show that
$\hbar\omega^{r}_{s}\\!\rightarrow\\!\hbar\omega_{s}$ and
$T^{r}\\!\rightarrow\\!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 Huang1987 .

On the other hand, $\eta_{c}=\omega_{s}/\omega_{c}$ is a critical coupling
strength for Ohmic spectral density that when $\eta>\eta_{c}$, the system
exists a dissipationless localized bound state (localized mode) at frequency
$\omega_{b}=\omega_{s}+\Delta(\omega_{b})$ with $J(\omega_{b})=0$ Zhang2012 .
Once such a localized mode exists, the spectral function $D(\omega)$ of the
system in Eq. (16) is modified as

$\displaystyle
D(\omega)\\!=Z(\omega_{b})\delta(\omega\\!-\\!\omega_{b})+\\!\frac{J(\omega)}{[\omega\\!-\\!\omega_{s}\\!-\\!\Delta(\omega)]^{2}\\!+\\!\pi^{2}J^{2}(\omega)},$
(26)

where
$Z(\omega_{b})=[1-\partial\Sigma(\omega)/\partial\omega]^{-1}\big{|}_{\omega=\omega_{b}}$
is the localized bound state wavefunction. Then, the asymptotic value of the
Green function $u(t\\!\rightarrow\\!\infty,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
$\rho_{mm^{\prime}}$ of Eq. (12). In other words, the system cannot be
thermalized with the reservoir Xiong2015 ; Xiong2020 ; Ali2020 .

Figure 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
$|n_{0}\rangle=|5\rangle,|10\rangle$ and $|15\rangle$ (corresponding to the
blue dashed-dot 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 ($\eta\ll\eta_{c}$), (b) the strong coupling $\eta<\eta_{c}$),
and (c) the ultra-strong coupling $\eta>\eta_{c}$). The initial bath
temperature $T_{0}=30\hbar\omega_{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 ($\eta=0.01\eta_{c}\ll\eta_{c}$),
the strong coupling ($\eta=0.5\eta_{c}<\eta_{c}$) and the ultra-strong
coupling ($\eta=1.2\eta_{c}>\eta_{c}$) cases. The results show that when
$\eta>\eta_{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 Ali2020  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 $\eta>\eta_{c}$
corresponds to a realization of the many-body localization Nandkishore15 .
When the coupling strength $\eta$ crosses the critical value $\eta_{c}$, the
transition from thermalization to many-body localization occurs Nandkishore15
. Our exact solution provides the foundation of this transition between
thermodynamics and many-body localization.

###  II.3 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 $\hbar\omega^{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 Zemansky1997 . The other part is the change of the
density state $\rho_{{}_{\\!S}}(t)$ which corresponds to quantum heat
associated with the entropy production. Consequently,

$\displaystyle\frac{dU_{{}_{\\!S}}(t)}{dt}$ $\displaystyle=\\!{\rm
Tr}_{{}_{\\!S}}\Big{[}\rho_{{}_{\\!S}}(t)\frac{dH^{r}_{{}_{\\!S}}\\!(t)}{dt}\Big{]}+{\rm
Tr}_{{}_{\\!S}}\Big{[}H^{r}_{{{}_{\\!S}}}\\!(t)\frac{d\rho_{{}_{\\!S}}(t)}{dt}\Big{]}$
$\displaystyle=\frac{dW_{s}(t)}{dt}+\frac{dQ_{s}(t)}{dt}.$ (27)

This is the first law of nonequilibrium quantum thermodynamics. Thus, the
quantum work and quantum heat can be naturally determined by

$\displaystyle\frac{dW_{{}_{\\!S}}(t)}{dt}$ $\displaystyle=\\!{\rm
Tr}_{{}_{\\!S}}\Big{[}\rho_{{}_{\\!S}}(t)\frac{dH^{r}_{{}_{\\!S}}(t)}{dt}\Big{]}=\overline{n}(t)\frac{d(\hbar\omega^{r}_{s}(t))}{dt},$
(28a) $\displaystyle\frac{dQ_{{}_{\\!S}}(t)}{dt}$ $\displaystyle={\rm
Tr}_{{}_{\\!S}}\Big{[}H^{r}_{{}_{\\!S}}(t)\frac{d\rho_{{}_{\\!S}}(t)}{dt}\Big{]}=T^{r}\\!(t)\frac{dS_{{}_{\\!S}}(t)}{dt},$
(28b)

where $\overline{n}(t)={\rm Tr}_{{}_{\\!S}}[a^{\dagger}a\rho_{{}_{\\!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 Esposito2010b ; Binder2015 ;
Esposito2015 ; Alipour2016 ; Strasberg2017 ; Perarnau2018 ; Rivas2020 . 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

$\displaystyle H_{{}_{\\!S}}^{r}(t)$
$\displaystyle=\hbar\omega^{r}_{s}(t,t_{0})a^{\dagger}a=\hbar\omega_{s}a^{\dagger}a+\delta\omega_{{}_{\\!S}}(t,t_{0})a^{\dagger}a$
$\displaystyle=H_{{}_{\\!S}}+\delta H_{{}_{\\!S}}(t).$ (29)

From Eq. (7a) and Eq. (8a), we have

$\displaystyle\omega^{r}_{s}(t,t_{0})$
$\displaystyle=\omega_{{}_{\\!S}}\\!+\delta\omega_{{}_{\\!S}}(t,t_{0})$
$\displaystyle=\omega_{{}_{\\!S}}\\!+\frac{1}{2}{\rm
Im}\bigg{[}\\!\int^{t}_{t_{0}}\\!\\!d\tau
g(t,\tau)u(\tau,t_{0})/u(t,t_{0})\bigg{]}.$ (30)

Here $H_{{}_{\\!S}}=\hbar\omega_{s}a^{\dagger}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(\tau,t_{0})$ approximately with the free-
particle (zero-th order) Green function
$u_{0}(\tau,t_{0})=e^{-i\omega_{{}_{\\!S}}(\tau-t_{0})}$ (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 Xiong2010 . But once we have the exact master
equation with the exact solution, such approximated master equation is no
longer needed.

Figure 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 $\Sigma(\omega)\\!=\\!\int
d\omega^{\prime}\frac{J(\omega^{\prime})}{\omega-\omega^{\prime}}$ is the
self-energy arisen from the coupling between the system and the reservoir, and
$J(\omega)\equiv\sum_{k}|V_{k}|^{2}\delta(\omega-\omega_{k})$. The
renormalized system energy $\hbar\omega^{r}_{s}$ and the dissipation
coefficient $\gamma$ 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

$\displaystyle
dQ_{{}_{\\!S}}(t)/dt=\hbar\omega^{r}_{s}\big{[}-2\gamma(t,t_{0})n(t)+\widetilde{\gamma}(t,t_{0})\big{]},$
(31)

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)$ Ali2020a ;
Callen1985 :

$\displaystyle F_{{}_{\\!S}}(t)$
$\displaystyle=U_{{}_{\\!S}}(t)-T^{r}\\!(t)S_{{}_{\\!S}}(t)$
$\displaystyle\stackrel{{\scriptstyle
t\rightarrow\infty}}{{\longrightarrow}}-(1/\beta^{r})\ln Z^{r}_{{}_{\\!S}}.$
(32)

From the above solution, we have

$\displaystyle
dF_{{}_{\\!S}}(t)=dW_{{}_{\\!S}}(t)-S_{{}_{\\!S}}(t)dT^{r}\\!(t),$ (33)

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 Callen1985 . 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),

$\displaystyle C=$
$\displaystyle\frac{dQ_{{}_{\\!S}}}{dT^{r}}=T^{r}\frac{dS_{{}_{\\!S}}}{dT^{r}}=\frac{\partial
U_{{}_{\\!S}}}{\partial T^{r}}\bigg{|}_{\omega^{r}_{s}},$ (34)

where the third thermodynamic law is justified from the specific heat at
arbitrary coupling strength: $C\sim(T^{r})^{3}$ as $T^{r}\\!\rightarrow\\!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).

Figure 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 $\hbar\omega_{s}^{2}$), 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 $\eta/\eta_{c}=0.3,0.5,0.8$, respectively. Other parameters
are taken the same as that in Fig. 1.

##  III The more general formulation of quantum thermodynamics for all
couplings

###  III.1 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
Seifert2016 ; Carrega2016 ; Ochoa2016 ; Jarzynski2017 ; Marcantoni2017 ;
Bruch2018 ; Perarnau2018 ; Hsiang2018 ; Anders2018 ; Strasberg2019 ;
Newman2020 ; Ali2020a ; Rivas2020 , 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
$H_{S}=\sum_{i}\varepsilon_{i}a^{\dagger}_{i}a_{i}$. As a specific example,
consider the system be an individual system and the reservoir be a many-body
system. The system Hamiltonian can be generally expressed as

$\displaystyle H_{{}_{\\!S}}=\frac{{\bf P}^{2}}{2m}+V({\bf
Q})=\sum_{i}\varepsilon_{i}|\psi_{i}\rangle\langle\psi_{i}|=\sum_{i}\varepsilon_{i}a^{\dagger}_{i}a_{i}.$
(35)

In Eq. (35), the second equality is the spectral decomposition of the system
Hamiltonian: $H_{{}_{\\!S}}|\psi_{i}\rangle=\varepsilon_{i}|\psi_{i}\rangle$,
and the last equality uses the second quantization language:
$|\psi_{i}\rangle=a^{\dagger}_{i}|0\rangle$ and $a_{i}|0\rangle=0$, and
$|0\rangle$ is the vacuum state. The particle creation and annihilation
operators $a^{\dagger}_{i}$ and $a_{i}$ obey the standard bosonic commutation
and fermionic anticommutation relations:
$[a_{i},a^{\dagger}_{j}]_{\mp}=a_{j}a^{\dagger}_{j}\mp
a^{\dagger}_{j}a_{i}=\delta_{ij}$ when the system being boson and fermion
systems, respectively.

Similarly, the Hamiltonian of a reservoir can also be written as
$H_{{}_{\\!E}}=\sum_{k}\epsilon_{k}b^{\dagger}_{k}b_{k}$, where $\epsilon_{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
Thouless1972 . Then the reservoir Hamiltonian can be expressed approximated as

$\displaystyle H_{{}_{\\!E}}$
$\displaystyle\simeq\\!\sum_{j}\\!\Big{[}\frac{\bm{p}^{2}_{j}}{2m_{j}}\\!+\\!U({\bm{q}}_{j})\\!+\\!\overline{\sum_{j^{\prime}}\\!V({\bm{q}}_{j},{\bm{q}}_{j^{\prime}})}\,\Big{]}$
$\displaystyle=\sum_{k}\epsilon_{k}|\psi_{k}\rangle\langle\psi_{k}|=\sum_{k}\epsilon_{k}b^{\dagger}_{k}b_{k},$
(36)

where $\overline{\sum_{j^{\prime}}V({\bm{q}}_{j},{\bm{q}}_{j^{\prime}})}$
represents the effective mean-field potential of many-body interactions, and
$\big{[}\frac{\bm{p}^{2}_{j}}{2m_{j}}+U({\bm{q}}_{j})+\overline{\sum_{j^{\prime}}V({\bm{q}}_{j},{\bm{q}}_{j^{\prime}})}\,\big{]}|\psi_{k}\rangle=\epsilon_{k}|\psi_{k}\rangle$
gives the quasiparticle continuous spectrum of the reservoir. The reservoir
particle creation and annihilation operators $b^{\dagger}_{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,

$\displaystyle
H_{{{}_{\\!S}}{{}_{\\!E}}}=\sum_{ik}\big{(}V_{ik}a^{\dagger}_{i}b_{k}+V^{*}_{ik}b^{\dagger}_{k}a_{i}\big{)},$
(37)

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 Hang1996 ; Miroshnichenko2010 ; Jin2010 ; Lei2012 . The
coupling strengths $V_{ik}$ are proportional to the quasiparticle wavefunction
overlaps between the system and reservoirs and therefore are tunable through
nanotechnology manipulations Hang1996 ; Miroshnichenko2010  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

$\displaystyle H_{\rm tot}(t)=$
$\displaystyle\,H_{{}_{\\!S}}(t)+\sum_{\alpha}H^{\alpha}_{{}_{\\!E}}(t)+\sum_{\alpha}H^{\alpha}_{{{}_{\\!S}}{{}_{\\!E}}}(t)$
$\displaystyle=$
$\displaystyle\sum_{i}\varepsilon_{i}(t)a^{\dagger}_{i}a_{i}\\!+\\!\\!\sum_{\alpha
k}\epsilon_{\alpha k}(t)b^{\dagger}_{\alpha k}b_{\alpha k}$
$\displaystyle+\\!\sum_{\alpha ik}\\!\big{[}V_{\alpha
ik}(t)a^{\dagger}_{i}b_{\alpha k}\\!+\\!V^{*}_{\alpha
ik}(t)b^{\dagger}_{\alpha k}a_{i}\big{]},$ (38)

which describes the system one concerned couples to multiple reservoirs. This
is a generalization of the Fano-Anderson Hamiltonian we introduced Zhang2012 ;
Zhang2018 . The index $\alpha$ 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 Feynman1963  as well as by Caldeira and Leggett Leggett1983b  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:

$\displaystyle\rho_{{}_{\\!S}}(t)\\!=\\!{\rm Tr}_{{}_{\\!E}}\big{[}{\cal
U}(t,t_{0})\big{(}\rho_{{}_{\\!S}}(t_{0}){\prod}_{\alpha}\\!\\!\\!\otimes\rho^{\alpha}_{{}_{\\!E}}(t_{0})\big{)}{\cal
U}^{\dagger}(t,t_{0})\big{]},$ (39)

where ${\cal U}(t,t_{0})={\cal
T}_{\rightarrow}\exp\big{\\{}\\!-\\!\frac{i}{\hbar}\\!\int^{t}_{t_{0}}\\!H_{\rm
tot}(t^{\prime})dt^{\prime}\big{\\}}$ is the time evolution of the total
system, and ${\cal T_{\rightarrow}}$ is the time-ordering operator. Here the
system is initially in an arbitrary state $\rho_{{}_{\\!S}}(t_{0})$. All
reservoirs can be initially in their own equalibrium thermal states,
$\rho^{\alpha}_{{}_{\\!E}}(t_{0})=e^{-\beta_{\alpha
0}(H^{\alpha}_{{}_{\\!E}}-\mu_{\alpha 0}\hat{N}^{\alpha})}/Z_{\alpha}$ which
can have different initial temperature $\beta_{\alpha 0}=1/k_{B}T_{\alpha 0}$
and different chemical potentials $\mu_{\alpha 0}$ for different reservoir
$\alpha$. Here $\hat{N}^{\alpha}$ is the total particle number operator of
reservoir $\alpha$. After trace out all the environmental states through the
coherent state path integrals Zhang1990 , the resulting exact master equation
of the system is indeed a generalization of Eq. (6a) to multi-level open
systems Tu2008 ; Jin2010 ; Lei2012 ; Yang2015 ; Yang2017 ; Zhang2018 ;
Huang2020 ,

$\displaystyle\frac{d}{dt}{\rho}_{{}_{\\!S}}(t)=$
$\displaystyle\frac{1}{i\hbar}\big{[}H^{r}_{{}_{\\!S}}(t),\rho_{{}_{\\!S}}(t)\big{]}\\!+\\!\\!{\sum}_{ij}\\!\big{\\{}\gamma_{ij}\left(t,t_{0}\right)\\!\\!\big{[}2a_{j}\rho_{{}_{\\!S}}(t)a_{i}^{\dagger}$
$\displaystyle-\\!a_{i}^{\dagger}a_{j}\rho_{{}_{\\!S}}(t)\\!-\\!\rho_{{}_{\\!S}}(t)a_{i}^{\dagger}a_{j}\big{]}\\!+\\!\widetilde{\gamma}_{ij}(t,t_{0})\big{[}a_{i}^{\dagger}\rho_{{}_{\\!S}}(t)a_{j}$
$\displaystyle~{}~{}~{}{\pm}a_{j}\rho_{{}_{\\!S}}(t)a_{i}^{\dagger}\mp
a_{i}^{\dagger}a_{j}\rho_{{}_{\\!S}}(t)\\!-\\!\rho_{{}_{\\!S}}(t)a_{j}a_{i}^{\dagger}\big{]}\big{\\}}.$
(40)

where the upper and lower signs of $\pm$ 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,

$\displaystyle
H^{r}_{{}_{\\!S}}(t)={\sum}_{ij}\varepsilon^{r}_{s,ij}(t,t_{0})a_{i}^{\dagger}a_{j}$
(41)

and the dissipation and fluctuations coefficients $\gamma_{ij}(t,t_{0})$ and
$\widetilde{\gamma}_{ij}(t,t_{0})$ in Eq. (40). These time-dependent
coefficients are determined nonperturbatively and exactly by the following
relations,

$\displaystyle\varepsilon^{r}_{ij}(t,t_{0})=\\!-\hbar{\rm
Im}\big{[}\dot{\bm{u}}(t,t_{0})\bm{u}^{-1}(t,t_{0})\big{]}_{ij},$ (42a)
$\displaystyle\gamma_{ij}(t,t_{0})=\\!-{\rm
Re}\big{[}\dot{\bm{u}}(t,t_{0})\bm{u}^{-1}(t,t_{0})\big{]}_{ij},$ (42b)
$\displaystyle\widetilde{\gamma}_{ij}(t,t_{0})=\dot{\bm{v}}_{ij}(t,t)\\!-\\!\big{[}\dot{\bm{u}}(t,t_{0})\bm{u}^{-1}(t,t_{0})\bm{v}(t,t)\\!+\\!\text{h.c.}\big{]}_{ij}.$
(42c)

where ${\bm{u}}(t,t_{0})$ and ${\bm{v}}(t,t)$ are $N\times 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})\equiv\langle[a_{i}(t),a^{\dagger}_{j}(t_{0})]_{\pm}$ which
obeys the equation of motion Zhang2012 ; Tu2008 ; Jin2010 ; Lei2012 ,

$\displaystyle\frac{d}{dt}\bm{u}(t,t_{0})-$
$\displaystyle\frac{1}{i\hbar}{\bm{\varepsilon}}(t)\bm{u}(t,t_{0})+\\!\\!\int_{t_{0}}^{t}\\!\\!\\!\\!dt^{\prime}\bm{g}(t,t^{\prime})\bm{u}(t^{\prime},t_{0})=0.$
(43a) The nonequilibrium correlation Green function ${\bm{v}}(t,t)$ obeys the
nonequilibrium fluctuation-dissipation relation Zhang2012 ,
$\displaystyle{\bm{v}}(\tau,t)\\!=\\!\\!\\!\int^{\tau}_{t_{0}}\\!\\!dt_{1}\\!\int^{t}_{t_{0}}\\!\\!dt_{2}\bm{u}(\tau,t_{1})\,\widetilde{\bm{g}}(t_{1},t_{2})\,\bm{u}^{\dagger}(t,t_{2}).$
(43b)

The integral memory kernels $\bm{g}(t,t^{\prime})$ and
$\widetilde{\bm{g}}(t_{1},t_{2})$ are the system-reservoir time correlations
and are given by

$\displaystyle g_{ij}(t,t^{\prime})=\\!\sum_{\alpha
k}\frac{1}{\hbar^{2}}V_{\alpha ik}(t^{\prime})V^{*}_{\alpha
jk}(t)\exp\\!\bigg{\\{}\\!\\!-\\!\frac{i}{\hbar}\\!\\!\int_{t}^{t^{\prime}}\\!\\!\\!\\!d\tau\epsilon_{\alpha
k}(\tau)\bigg{\\}},$ (44a)
$\displaystyle\widetilde{g}_{ij}(t_{1},t_{2})=\\!\sum_{\alpha
k}\frac{1}{\hbar^{2}}V_{\alpha ik}(t_{2})V^{*}_{\alpha
jk}(t_{1})\big{\langle}b^{\dagger}_{\alpha k}(t_{0})b_{\alpha
k}(t_{0})\big{\rangle}_{E}$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\times\\!\exp\\!\bigg{\\{}\\!\\!-\\!\frac{i}{\hbar}\\!\\!\int_{t_{2}}^{t_{1}}\\!\\!\\!\\!d\tau\epsilon_{\alpha
k}(\tau)\bigg{\\}}.$ (44b)

Here the initial reservoir correlation function,

$\displaystyle\big{\langle}b^{\dagger}_{\alpha k}(t_{0})b_{\alpha
k}(t_{0})\big{\rangle}_{E}$ $\displaystyle=f(\epsilon_{\alpha k},T_{\alpha
0},\mu_{\alpha 0})$ $\displaystyle=\frac{1}{[e^{(\epsilon-\mu_{\alpha
0})/k_{B}T_{\alpha 0}}{\mp}1]},$ (45)

determines the initial particle distribution of the bosons or fermions in the
initial thermal reservoir $\alpha$ with the chemical potential $\mu_{\alpha
0}$ and the temperature $T_{\alpha 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
$g_{ij}(t,t^{\prime})\\!=\\!\\!\\!\int\\!\\!d\epsilon
J_{ij}(\epsilon)e^{-i\epsilon(t-t^{\prime})}$,
$\widetilde{g}_{ij}(t_{1},t_{2})\\!=\\!\\!\\!\int\\!\\!d\epsilon
J_{ij}(\epsilon)f(\epsilon,T_{\alpha},\mu_{\alpha})e^{-i\epsilon(t_{1}-t_{2})}$,
where

$\displaystyle J_{ij}(\epsilon)=\frac{1}{\hbar^{2}}{\sum}_{\alpha k}V_{\alpha
ik}V^{*}_{\alpha jk}\delta(\epsilon-\epsilon_{\alpha
k})={\sum}_{\alpha}J_{\alpha,ij}(\epsilon),$ (46)

and $J_{\alpha ij}(\epsilon)$ is the spectral density matrix of reservoir
$\alpha$.

###  III.2 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 Xiong2020  and its exact steady state is
(see a detailed derivation in Appendix B)

$\displaystyle\rho^{\rm
exact}_{{}_{\\!S}}(t\\!\rightarrow\\!\infty)=\frac{\exp\Big{\\{}\sum_{ij}\Big{(}\\!\ln\frac{\overline{\bm{n}}}{I\pm\overline{\bm{n}}}\Big{)}_{ij}a^{\dagger}_{i}a_{j}\Big{\\}}}{[\det(I\pm\overline{\bm{n}})]^{\pm
1}}$ (47)

which is a generalized Gibbs-type state. Here
$\overline{n}_{ij}=\lim_{t\rightarrow\infty}n_{ij}(t)$ is the one-particle
density matrix defined as Jin2010 ; Yang1962

$\displaystyle n_{ij}(t)\equiv{\rm
Tr}_{{}_{\\!S}}[a^{\dagger}_{i}a_{j}\rho_{{}_{\\!S}}(t)]=\rho^{(1)}_{ij}(t).$
(48)

The solution Eq. (47) remains the same for initial system-reservoir correlated
states with a modification of $\widetilde{\bm{g}}(t_{1},t_{2})$ in Eq. (43b)
to include the initial correlations between the system and reservoirs Yang2015
; Huang2020 . Thus, the nonequilibrium internal energy, entropy and particle
number can be defined by

$\displaystyle U_{{}_{\\!S}}(t)\\!\equiv\\!{\rm
Tr}_{{}_{\\!S}}[H^{r}_{{}_{\\!S}}\\!(t)\rho_{{}_{\\!S}}(t)]=\sum_{ij}\varepsilon^{r}_{ij}(t)n_{ij}(t),$
(49a) $\displaystyle S_{{}_{\\!S}}(t)\equiv\\!-\\!k_{B}{\rm
Tr}_{{}_{\\!S}}[\rho_{{}_{\\!S}}(t)\ln\rho_{{}_{\\!S}}(t)],$ (49b)
$\displaystyle N_{{}_{\\!S}}(t)\equiv\\!{\rm
Tr}_{{}_{\\!S}}[a^{\dagger}_{i}a_{i}\rho_{{}_{\\!S}}(t)]\\!=\\!\
\sum_{i}n_{ii}(t).$ (49c)

They are related to each other and may form the fundamental equation for
quantum thermodynamics Callen1985 ; Ali2020a :
$U_{{}_{\\!S}}(t)=U_{{}_{\\!S}}(\varepsilon^{r}_{s}(t),S_{{}_{\\!S}}(t),N_{{}_{\\!S}}(t))$.
Here energy levels play a similar role as the volume Zemansky1997 . Thus,

$\displaystyle
dU_{{}_{\\!S}}(t)\\!=\\!dW_{{}_{\\!S}}(t)+T^{r}(t)dS_{{}_{\\!S}}(t)+\mu^{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,

$\displaystyle\frac{dW_{{}_{\\!S}}(t)}{dt}=\\!{\rm
Tr}_{{}_{\\!S}}\Big{[}\rho_{{}_{\\!S}}(t)\frac{dH^{r}_{{}_{\\!S}}(t)}{dt}\Big{]}\\!=\\!{\sum}_{ij}n_{ij}(t)\frac{d\varepsilon^{r}_{s,ij}(t)}{dt}.$
(51)

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)),

$\displaystyle dQ_{{}_{\\!S}}(t)+dW^{c}_{{}_{\\!S}}(t)$
$\displaystyle\\!=\\!{\sum}_{ij}\\!\varepsilon^{r}_{s,ij}(t)dn_{ij}(t)$
$\displaystyle=\\!T^{r}\\!(t)dS_{{}_{\\!S}}(t)\\!+\\!\mu^{r}\\!(t)dN_{{}_{\\!S}}(t).$
(52)

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 $\hbar\omega$ so that
$\mu^{r}(t)\\!=\\!0$. From the above formulation, we can define the
renormalized temperature and renormalized chemical potential by

$\displaystyle T^{r}\\!(t)=\\!\frac{\partial U_{{}_{\\!S}}(t)}{\partial
S_{{}_{\\!S}}(t)}\bigg{|}_{\varepsilon^{r}_{s}(t),N_{{}_{\\!S}}(t)},~{}\mu^{r}(t)=\\!\frac{\partial
U_{{}_{\\!S}}(t)}{\partial
N_{{}_{\\!S}}(t)}\bigg{|}_{\varepsilon^{r}_{s}(t),S_{{}_{\\!S}}(t)}.$ (53)

As a result, Eq. (47) can be also written as the standard Gibbs state,

$\displaystyle\rho^{\rm
exact}_{{}_{\\!S}}(t\\!\rightarrow\\!\infty)=\frac{1}{Z^{r}}\exp\big{\\{}\\!-\\!\beta^{r}(H^{r}_{{}_{\\!S}}\\!-\\!\mu^{r}\hat{N})\big{\\}},$
(54)

which is given in terms of the renormalized Hamiltonian
$H^{r}_{{}_{\\!S}}(t)$, the renormalized temperature $T^{r}(t)$ and the
renormalized chemical potential $\mu^{r}(t)$ at steady state, and
$\hat{N}=\sum_{i}a^{\dagger}_{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.

###  III.3 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
Hang1996 ; Tu2008 ; Jin2010 , see Fig. 6(a). The total Hamiltonian is

$\displaystyle H_{\rm tot}\\!=$
$\displaystyle\sum_{\sigma}\varepsilon_{\sigma}a^{\dagger}_{\sigma}a_{\sigma}\\!+\\!\sum_{\alpha,k,\sigma}\epsilon_{\alpha
k}c^{\dagger}_{\alpha k\sigma}c_{\alpha k\sigma}$
$\displaystyle+\\!\sum_{\alpha,k,\sigma}(V_{\alpha
k}a^{\dagger}_{\sigma}c_{\alpha k\sigma}\\!+\\!V^{*}_{\alpha
k}c^{\dagger}_{\alpha k\sigma}a_{\sigma}).$ (55)

The index $\sigma=\uparrow,\downarrow$ labels electron spin, $\alpha=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
$\mu_{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 $\mu^{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 Seifert2016 ;
Carrega2016 ; Ochoa2016 ; Jarzynski2017 ; Marcantoni2017 ; Bruch2018 ;
Perarnau2018 ; Hsiang2018 ; Anders2018 ; Strasberg2019 ; Newman2020 ;
Rivas2020  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\rangle,|\uparrow\rangle,|\downarrow\rangle,|d\rangle$ (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,

$\displaystyle\rho(t)=\begin{pmatrix}\rho_{00}(t)&0&0&0\\\
0&\rho_{\uparrow\uparrow}(t)&\rho_{\uparrow\downarrow}(t)&0\\\
0&\rho_{\downarrow\uparrow}(t)&\rho_{\downarrow\downarrow}(t)&0\\\
0&0&0&\rho_{dd}(t)\end{pmatrix}.$ (56)

If the dot is initially empty, the $4\times 4$ reduced density matrix has been
solved exactly from the exact master equation Tu2011 ; Yang2018 :

$\displaystyle\rho_{00}(t)\\!=\\!\det[I\\!-\\!\bm{v}(t,t)],~{}~{}\rho_{dd}\\!=\\!\det[\bm{v}(t)],$
$\displaystyle\rho_{\uparrow\uparrow}(t)\\!=\\!v_{\uparrow\uparrow}(t)\\!-\\!\rho_{33}(t),~{}~{}\rho_{\downarrow\downarrow}(t)\\!=\\!v_{\downarrow\downarrow}(t)\\!-\\!\rho_{33}(t),$
$\displaystyle\rho_{\uparrow\downarrow}(t)\\!=\\!v_{\uparrow\downarrow}(t)\\!=\\!\rho^{*}_{\downarrow\uparrow}(t).$
(57)

Here the $2\times 2$ matrix Green function $\bm{v}(t)\equiv\bm{v}(t,t)$ is
determined by the Green function $\bm{u}(t,t^{\prime})$. We take reservoir
spectra as a Lorentzian form, then the spectral densities
$J_{\alpha}(\epsilon)$ can be expressed as Meir1993 ; Tu2008 :

$\displaystyle
J_{\alpha,ij}(\epsilon)=\frac{\Gamma_{\alpha}d^{2}}{\epsilon^{2}+d^{2}}\delta_{ij}~{}~{}(i,j=\uparrow,\downarrow),$
(58)

where $\Gamma_{\alpha}$ is the tunneling rate (the coupling strength) between
the quantum dot and the lead $\alpha$. For simplicity, we also ignore the
spin-flip tunneling. The exact solution of the reduced density matrix is
rather simple,

$\displaystyle\rho(t)=\det[1-\bm{v}(t)]\exp\Big{\\{}\bm{a}^{\dagger}\ln\frac{\bm{v}(t)}{1-\bm{v}(t)}\bm{a}\Big{\\}}.$
(59)

Here $\bm{a}^{\dagger}=(a^{\dagger}_{\uparrow},a^{\dagger}_{\downarrow})$, and

$\displaystyle
v_{ii}(t)=\\!\\!\int^{t}_{t_{0}}\\!\\!\\!dt_{1}\\!\\!\int^{t}_{t_{0}}\\!\\!\\!\\!dt_{2}\,u_{ii}(t,t_{1})\,\widetilde{g}(t_{1},t_{2})\,u^{*}_{ii}(t,t_{2}),$
(60a) $\displaystyle v_{\uparrow\downarrow}=0,$ (60b)
$\displaystyle\frac{d}{dt}u_{ii}(t,t_{0})+i\varepsilon_{i}u_{ii}(t,t_{0})+\\!\\!\int_{t_{0}}^{t}\\!\\!\\!dt^{\prime}\\!\\!\int\\!\\!d\epsilon\frac{\Gamma
d^{2}e^{-\epsilon(t-t^{\prime})}}{\epsilon^{2}+d^{2}}u_{ii}(t,t^{\prime})$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}=0,$
(60c)

for $i=\uparrow,\downarrow$, and $\Gamma=\Gamma_{L}+\Gamma_{R}$. As a result,
the nonequilibrium internal energy, the entropy and the total average particle
number can be found analytically,

$\displaystyle U_{{}_{\\!S}}(t)=$
$\displaystyle\varepsilon^{r}_{\uparrow}(t)v_{\uparrow\uparrow}(t)+\varepsilon^{r}_{\downarrow}(t)v_{\downarrow\downarrow}(t),$
(61a) $\displaystyle S_{{}_{\\!S}}(t)=$ $\displaystyle-
v_{\uparrow\uparrow}(t)\ln
v_{\uparrow\uparrow}(t)-v_{\downarrow\downarrow}(t)\ln
v_{\downarrow\downarrow}(t)$
$\displaystyle-(1-v_{\uparrow\uparrow}(t))\ln(1-v_{\uparrow\uparrow}(t))$
$\displaystyle-(1-v_{\downarrow\downarrow}(t))\ln(1-v_{\downarrow\downarrow}(t)),$
(61b) $\displaystyle N_{{}_{\\!S}}(t)=$ $\displaystyle
v_{\uparrow\uparrow}(t)+v_{\downarrow\downarrow}(t).$ (61c)

From the above solution, the corresponding renormalized energy, renormalized
temperature and renormalized chemical potential can be calculated
straightforwardly.

Figure 6: (a) A schematic plot of the single-electron transistor device.
(b)-(g) The nonequilibrium evolution of the energy levels
$\varepsilon^{r}_{\uparrow,\downarrow}(t)$, the particle occupation in each
level $\overline{n}_{\uparrow,\downarrow}(t)$, the internal energy
$U_{{}_{\\!S}}(t)$, the entropy $S_{{}_{\\!S}}(t)$, the renormalized
temperature $T^{r}(t)$ and chemical potential $\mu^{r}(t)$ at different
coupling strength
$\Gamma=0.2\varepsilon_{\uparrow},0.8\varepsilon_{\uparrow}$, respectively.
(h) The steady-state value of the renormalized energy levels
$\varepsilon^{r}_{\uparrow,\downarrow}$, the renormalized temperature $T^{r}$
and the renormalized chemical potential $\mu^{r}$ as a function of the
coupling strength, and (i) the comparison of the renormalized Fermi-Dirac
distribution $f(\varepsilon^{r}_{\uparrow,\downarrow},T^{r},\mu^{r})$ with the
exact solution of the
$\overline{n}_{\uparrow,\downarrow}(t\\!\rightarrow\\!\infty)$ as a function
of the coupling strength $\Gamma$. Other parameters:
$\varepsilon_{\downarrow}=3\varepsilon_{\uparrow}$,
$k_{B}T_{L,R}=(3,0.1)\varepsilon_{\uparrow}$,
$\mu_{L,R}=(5,2)\varepsilon_{\uparrow}$, and $d=10\varepsilon_{\uparrow}$.

In Fig. 6(b)-(g), we show the nonequilibrium evolution of the renormalized
energy levels $\varepsilon^{r}_{\uparrow,\downarrow}(t)$, the particle
occupations in each levels $\overline{n}_{\uparrow,\downarrow}(t)$, the
internal energy $U_{{}_{\\!S}}(t)$, the entropy $S_{{}_{\\!S}}(t)$, the
renormalized temperature $T^{r}(t)$ and the renormalized chemical potential
$\mu^{r}(t)$ for the coupling strength $\Gamma_{L}=\Gamma_{R}=\Gamma/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
$\varepsilon^{r}_{\uparrow,\downarrow}$, the renormalized temperature $T^{r}$
and the renormalized chemical potential $\mu^{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):
$f(\varepsilon^{r}_{\uparrow,\downarrow},T^{r},\mu^{r})=1/[e^{(\varepsilon^{r}_{\uparrow,\downarrow}-\mu^{r}_{\alpha})/k_{B}T^{r}_{\alpha}}+1]$.
We compare the renormalized Fermi-Dirac distributions with the exact solution
of the occupation numbers
$\overline{n}_{\uparrow,\downarrow}(t\\!\rightarrow\\!\infty)$ 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),

$\displaystyle N_{{}_{\\!S}}(t\rightarrow\infty)$
$\displaystyle=n_{\uparrow}(\mu^{r},T^{r})+n_{\downarrow}(\mu^{r},T^{r})$
$\displaystyle\stackrel{{\scriptstyle\rm
WCL}}{{\rightarrow}}\frac{1}{2}\big{[}n_{\uparrow}(\mu_{L},T_{L})+n_{\uparrow}(\mu_{R},T_{R})\big{]}$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}+\frac{1}{2}\big{[}n_{\downarrow}(\mu_{L},T_{L})+n_{\downarrow}(\mu_{R},T_{R})\big{]}.$
(62)

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
$\mu_{L},T_{L}$ and $\mu_{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
$\mu^{r}=(\mu_{L}+\mu_{R})/2$. From Fig. 7(b), we see that the renormalized
chemical potential in the very weak coupling limit is
$\mu^{r}=(\mu_{L}+\mu_{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
$\mu^{r}=(\mu_{L}+\mu_{R})/2\neq\mu_{L}\neq\mu_{R}$, Eq. (62) shows that
$T_{L}\neq T^{r}\neq T_{R}$ in the very weak-coupling limit, even through
$T_{L}=T_{R}$. Figure 7(c) shows further the case $\mu_{L}=\mu_{R}$ and
$T_{L}\neq T_{R}$. We have $\mu^{r}=\mu_{L}=\mu_{R}$, and from Eq. (62), we
find that $T^{r}\neq(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}\simeq(T_{L}+T_{R})/2$ in the very weak coupling limit, even if
$\mu_{L}\neq\mu_{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 weak-coupling
limit for multi-reservoirs. This justifies further the consistency of our
renormalized theory for quantum thermodynamics at any coupling.

Figure 7: The steady-state renormalized temperature and renormalized chemical
potential of the single-electron transistor as a function of the system-
reservoir coupling strength $\Gamma$ 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 system-reservoir 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 $\bm{u}(t,t_{0})$ of Eq.
(43a). It can be applied to arbitrary system interacting with arbitrary
environment. In Eq. (43a), $\bm{g}(t,t^{\prime})$ is the self-energy
correlation that can be easily generalized to any interacting system using the
nonequilibrium Green function technique in many-body systems Kadanoff1962 ;
Keldysh1965 . Meanwhile, the renormalized temperature $T^{r}(t)$ (the
renormalized chemical potential $\mu^{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

$\displaystyle\rho_{{}_{\\!S}}(t)={\rm Tr}_{{}_{\\!E}}[{\cal
U}(t,t_{0})\rho_{\rm tot}(t_{0}){\cal U}^{\dagger}(t,t_{0})].$ (63)

Here, ${\cal 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 Breuer2008 , or using the coherent state path integrals to
nonperturbatively trace over all the environmental states as we did Tu2008 ;
Jin2010 ; Lei2012 ; Yang2015 ; Zhang2018 ; Huang2020 ; Zhang1990 ; Zhang2022 .
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 Zhang2012 ; Zhang2018

$\displaystyle
H_{{{}_{\\!S}}{{}_{\\!E}}}=\sum_{ijk}\big{[}V_{ij}(k)c^{\dagger}_{i}c_{j}b^{\dagger}_{k}+V^{*}_{ij}(k)c^{\dagger}_{j}c_{i}b_{k}\big{]}.$
(64)

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^{\dagger}_{i},c_{i}$
($b^{\dagger}_{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 Peskin1995 ,

$\displaystyle{\cal L}_{\rm
QED}=\overline{\psi}(i\gamma^{\mu}\partial_{\mu}-m)\psi-\frac{1}{4}F^{\mu\nu}F_{\mu\nu}-e\overline{\psi}\gamma^{\mu}A_{\mu}\psi,$
(65)

where $\psi(x)$ is the fermionic field for electrons, $A_{\mu}(x)$ is the
covariant 4-vector of the electromagnetic (EM) field, $\gamma^{\mu}$ is the
Dirac matrix, and
$F_{\mu\nu}(x)=\partial_{\mu}A_{\nu}(x)-\partial_{\nu}A_{\mu}(x)$ 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
electron-photon interaction in the non-relativistic limit (by ignoring
positrons and choosing a proper gauge). Thus, the non-relativistic QED
Hamiltonian is given by Coulomb_gauge

$\displaystyle\\!\\!\\!H_{\rm QED}\\!=$ $\displaystyle H_{\rm
electron}\\!+\\!H_{\rm photon}\\!+\\!H_{\rm e-p}$ $\displaystyle=$
$\displaystyle\\!\sum_{\bm{p}}\\!\varepsilon_{\bm{p}}c^{\dagger}_{\bm{p}}c_{\bm{p}}\\!+\\!\\!\\!\\!\sum_{\bm{p},\bm{p}^{\prime},\bm{q}}\\!\\!\\!U(\bm{q})c^{\dagger}_{\bm{p}+\bm{q}}c^{\dagger}_{\bm{p}^{\prime}\\!-\\!\bm{q}}c_{\bm{p}^{\prime}}c_{\bm{p}}\\!+\\!\\!\sum_{\bm{k}}\\!\hbar\omega_{\bm{k}}b^{\dagger}_{\bm{k}}b_{\bm{k}}$
$\displaystyle+\\!\hbar\\!\sum_{\bm{p}\bm{k}}\\!\big{[}V({\bm{k}})c^{\dagger}_{\bm{p}}c_{\bm{p}-\bm{k}}b_{\bm{k}}\\!+\\!V^{*}({\bm{k}})c^{\dagger}_{\bm{p}-\bm{k}}c_{\bm{p}}b^{\dagger}_{\bm{k}}\big{]},$
(66)

where $c^{\dagger}_{\bm{p}},c_{\bm{p}}$ and $b^{\dagger}_{\bm{k}},b_{\bm{k}}$
are creation and annihilation operators of electrons and photons with momentum
$\bm{p}$ and $\bm{k}$. The summations over $\bm{p}$ and $\bm{k}$ should be
replaced by the continuous integrals
$\int\\!\\!\frac{d^{3}\bm{p}}{(2\pi)^{3}}$ and
$\int\\!\\!\frac{d^{3}\bm{k}}{(2\pi)^{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 non-relativistic 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 electron-photon interaction by using the
coherent state path integrals Zhang1990  to nonperturbatively trace over all
the environmental states. To do so, we may express the exact reduced density
$\rho_{{}_{\\!S}}(t)$ of Eq. (63) as
$\rho_{{}_{\\!S}}(\bm{\xi}^{\dagger}_{f},\bm{\xi}^{\prime}_{f},t)=\langle\bm{\xi}_{f}|\rho_{{}_{\\!S}}(t)|\bm{\xi}^{\prime}_{f}\rangle$
which is generally given by

$\displaystyle\rho_{{}_{\\!S}}(\bm{\xi}^{\dagger}_{f},\bm{\xi}^{\prime}_{f},t)=\\!\\!\int\\!\\!d\mu$
$\displaystyle(\bm{\xi}_{0})d\mu(\bm{\xi}^{\prime}_{0})\rho_{{}_{\\!S}}(\bm{\xi}_{0},\bm{\xi}^{\prime}_{0},t_{0})$
$\displaystyle~{}~{}~{}\times\\!\mathcal{J}_{\rm
QED}\\!(\bm{\xi}^{\dagger}_{f},\bm{\xi}^{\prime}_{f},t;\bm{\xi}_{0},\bm{\xi}^{\prime{\dagger}}_{0},t_{0}).$
(67)

Here we have utilized the unnormalized fermion coherent states
$|\bm{\xi}\rangle\equiv\exp(\sum_{\bm{p}}\xi_{\bm{p}}c^{{\dagger}}_{\bm{p}})|0\rangle$.
The integral measure
$d\mu(\bm{\xi})=\prod_{\bm{p}}d\xi^{*}_{\bm{p}}d\xi_{\bm{p}}e^{-|\xi_{\bm{p}}|^{2}}$
is defined the Haar measure in Grassmannian space. The vectors
$\bm{\xi}\equiv(\xi_{\bm{p}_{1}},\xi_{\bm{p}_{2}},...)$ is a one-column matrix
and $\xi^{*}_{\bm{p}_{i}}$ is a Grassmannian variable. The propagating
function
$\mathcal{J}\\!(\bm{\xi}^{\dagger}_{f},\bm{\xi}^{\prime}_{f},t;\bm{\xi}_{0},{\bm{\xi}^{\prime}}_{0}^{\dagger},t_{0})$
in Eq. (67), which describes the nonequilibrium evolution of the states of the
electron system from the initial state $\rho_{{}_{\\!S}}(t_{0})$ to the state
at any later time $\rho_{{}_{\\!S}}(t)$, can be obtained analytically after
completing exactly the path integrals over all the photon modes. The result is

$\displaystyle\mathcal{J}(\bm{\xi}^{\dagger}_{t},\bm{\xi}^{\prime}_{t},t;\,$
$\displaystyle\bm{\xi}_{0},{\bm{\xi}^{\prime}_{0}}^{{\dagger}},t_{0})=\int\mathcal{D}[\bm{\xi};\bm{\xi}^{\prime}]\exp\Big{\\{}\frac{i}{\hbar}\big{(}S_{s}[\bm{\xi}^{\dagger},\bm{\xi}]$
$\displaystyle-S^{*}_{s}[\bm{\xi}^{\prime{\dagger}},\bm{\xi}^{\prime}]+S^{\rm
qed}_{\rm
IF}[\bm{\xi}^{\dagger},\bm{\xi};\bm{\xi}^{\prime{\dagger}},\bm{\xi}^{\prime}]\big{)}\Big{\\}},$
(68)

where
$\mathcal{D}[\bm{\xi};\bm{\xi}^{\prime}]=\prod_{\bm{p},t<\tau<t_{0}}d\xi^{*}_{\bm{p}}(\tau)d\xi_{\bm{p}}(\tau)$,
$S_{e}[\bm{\xi}^{\dagger},\bm{\xi}]$ is the bare electron action of the
original free-electron Hamiltonian plus the electron-electron Coulomb
interaction in QED. In the fermion coherent state representation, it is given
by

$\displaystyle S_{e}[\bm{\xi}^{\dagger},\bm{\xi}]=$
$\displaystyle-{i\hbar\over
2}\big{[}\bm{\xi}^{\dagger}_{t}\bm{\xi}(t)+\bm{\xi}^{\dagger}(t_{0})\bm{\xi}_{t_{0}}\big{]}$
$\displaystyle+\\!\\!\int_{t_{0}}^{t}\\!\\!d\tau\Big{\\{}{i\hbar\over
2}\big{[}\dot{\bm{\xi}}^{\dagger}(\tau)\bm{\xi}(\tau)-\bm{\xi}^{\dagger}(\tau)\dot{\bm{\xi}}(\tau)\big{]}$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}-{\cal
H}(\bm{\xi}^{\dagger}(\tau),\bm{\xi}(\tau))\Big{\\}},$ (69)

where ${\cal
H}(\bm{\xi}^{\dagger}(\tau),\bm{\xi}(\tau))=\sum_{\bm{p}}\varepsilon_{\bm{p}}\xi^{*}_{\bm{p}}(\tau)\xi_{\bm{p}}(\tau)+\sum_{\bm{p},\bm{p}^{\prime},\bm{q}}U(\bm{q})\xi^{*}_{\bm{p}+\bm{q}}(\tau)\xi^{*}_{\bm{p}^{\prime}\\!-\\!\bm{q}}(\tau)\xi_{\bm{p}^{\prime}}(\tau)\xi_{\bm{p}}(\tau)$.
The additional action $S^{\rm qed}_{\rm
IF}[\bm{\xi}^{\dagger},\bm{\xi};\bm{\xi}^{\prime{\dagger}},\bm{\xi}^{\prime}]$
in Eq. (68) is an electron action correction arisen from electron-photon
interaction after exactly integrated out all the photon modes Zhang2022 :

$\displaystyle S^{\rm qed}_{\rm
IF}[\bm{\xi}^{\dagger},\bm{\xi};\bm{\xi}^{\prime{\dagger}},\bm{\xi}^{\prime}]=\\!\\!\int_{t_{0}}^{t}\\!\\!\\!$
$\displaystyle
d\tau\Bigg{\\{}i\hbar\sum_{\bm{p}\bm{p}^{\prime}\bm{k}}\bigg{[}\int_{t_{0}}^{\tau}\\!\\!\\!d\tau^{\prime}\\!\Big{\\{}\sigma^{+}_{\bm{p},\bm{k}}(\tau)G_{\bm{k}}(\tau,\tau^{\prime})\sigma^{-}_{\bm{p}^{\prime},\bm{k}}(\tau^{\prime})+\sigma^{\prime-}_{\bm{p}^{\prime},\bm{k}}(\tau)G^{*}_{\bm{k}}(\tau,\tau^{\prime})\sigma^{\prime+}_{\bm{p},\bm{k}}(\tau^{\prime})\\!\Big{\\}}$
$\displaystyle-\\!\\!\\!\int_{t_{0}}^{t}\\!\\!\\!d\tau^{\prime}\Big{\\{}\sigma^{\prime+}_{\bm{p},\bm{k}}(\tau)G_{\bm{k}}(\tau,\tau^{\prime})\sigma^{-}_{\bm{p}^{\prime},\bm{k}}(\tau^{\prime})\\!-\\!\big{[}\sigma^{+}_{\bm{p},\bm{k}}(\tau)\\!-\\!\sigma^{\prime+}_{\bm{p},\bm{k}}(\tau)\big{]}\widetilde{G}_{\bm{k}}(\tau,\tau^{\prime})\big{[}\sigma^{-}_{\bm{p}^{\prime},\bm{k}}(\tau^{\prime})\\!-\\!\sigma^{\prime-}_{\bm{p}^{\prime},\bm{k}}(\tau^{\prime})\big{]}\\!\Big{\\}}\bigg{]}\\!\Bigg{\\}}$
(70)

This is a generalization of the Feynman-Vernon influence functional
Feynman1963  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
$\sigma^{+}_{\bm{p},\bm{k}}(\tau)\equiv\xi^{*}_{\bm{p}}(\tau)\xi_{\bm{p}-\bm{k}}(\tau)$
and
$\sigma^{-}_{\bm{p},\bm{k}}(\tau)\equiv\xi^{*}_{\bm{p}-\bm{k}}(\tau)\xi_{\bm{p}}(\tau)=(\sigma^{+}_{\bm{p},\bm{k}}(\tau))^{\dagger}$,
which correspond to the spin-like variables of the exciton operators
$a^{\dagger}_{\bm{p}}a_{\bm{p}-\bm{k}}$ and
$a^{\dagger}_{\bm{p}-\bm{k}}a_{\bm{p}}$, respectively. The non-local time
correlations in Eq. (70) are given by

$\displaystyle G_{\bm{k}}(\tau,\tau^{\prime})$
$\displaystyle=|V(\bm{k})|^{2}e^{-i\omega_{\bm{k}}(\tau-\tau^{\prime})},$
(71a) $\displaystyle\widetilde{G}_{\bm{k}}(\tau,\tau^{\prime})$
$\displaystyle=|V(\bm{k})|^{2}\overline{n}(\omega_{\bm{k}},T_{0})e^{-i\omega_{\bm{k}}(\tau-\tau^{\prime})}$
(71b)

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. (89) and (91)
in Appendix B. The main difference is that the bare system action Eq. (90) and
the influence functional action Eq. (91) for the generalized 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. (94)
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. (91) 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. (91) shares the same property. Furthermore, it is not
difficult to show that the last two terms in both Eqs. (70) and (91) 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
(91) can combine with the forward and backward bare system actions in Eq. (69)
and (90), 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. (90) plus the first term in the
influence functional action Eq. (91), i.e.,

$\displaystyle{\cal H}^{r}$
$\displaystyle[\bm{\alpha}^{\dagger},\bm{\alpha}]={\cal
H}[\bm{\alpha}^{\dagger},\bm{\alpha}]+\delta{\cal
H}[\bm{\alpha}^{\dagger},\bm{\alpha}]$
$\displaystyle=\\!\sum_{i}\\!\varepsilon_{i}(\tau)\alpha^{*}_{i}(\tau)\alpha_{i}(\tau)\\!-\\!i\hbar\\!\sum_{ij}\\!\\!\int_{t_{0}}^{\tau}\\!\\!\\!\\!d\tau^{\prime}\\!\alpha^{*}_{i}(\tau)g_{ij}(\tau,\tau^{\prime})\alpha_{j}(\tau^{\prime}).$
(72)

Note that the evolution of $\alpha_{j}(\tau)$ along the forward path is
determined by Jin2010 ; Lei2012 ; Yang2015 ; Zhang2018

$\displaystyle\alpha_{j}(\tau)=\bm{u}_{jj^{\prime}}(\tau,t_{0})\alpha_{j}(t_{0})+f_{j}(\tau)$
(73)

where $\bm{u}_{jj^{\prime}}(\tau,t_{0})$ is the propagating Green function
which obeys the integro-differential Dyson equation Eq. (43a). While,
$f_{j}(\tau^{\prime})$ is the noise source associated with the correlation
Green function $\bm{v}(\tau,t)$ of Eq. (43b) that has no contribution to the
system Hamiltonian renormalization Yang2015 ; Zhang2018 . Thus, we can only
take the part of the evolution $\alpha_{j}(\tau^{\prime})$ that has the
contribution to system Hamiltonian renormalization, i.e.

$\displaystyle\alpha_{j}(\tau^{\prime})$
$\displaystyle\propto\bm{u}_{jj^{\prime}}(\tau^{\prime},t_{0})\alpha_{j}(t_{0})$
$\displaystyle\propto[\bm{u}(\tau^{\prime},t_{0})\bm{u}^{-1}(\tau,t_{0})]_{jj^{\prime}}\alpha_{j}(\tau).$
(74)

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

$\displaystyle{\cal H}^{r}[$
$\displaystyle\bm{\alpha}^{\dagger},\bm{\alpha}]=\\!\sum_{ij}\alpha^{*}_{i}(\tau)\varepsilon^{r}_{ij}(\tau,t_{0})\alpha_{ij}(\tau)$
(75)

and

$\displaystyle\varepsilon^{r}_{ij}(\tau,t_{0})$
$\displaystyle=\varepsilon_{i}(\tau)\delta_{ij}\\!+\\!\hbar{\rm
Im}\\!\\!\int_{t_{0}}^{\tau}\\!\\!\\!\\!\\!d\tau^{\prime}[\bm{g}(\tau,\tau^{\prime})\bm{u}(\tau^{\prime},t_{0})\bm{u}^{-1}(\tau,t_{0})]_{ij}$
$\displaystyle=-\hbar{\rm
Im}\big{[}\dot{\bm{u}}(\tau,t_{0})\bm{u}^{-1}(\tau,t_{0})\big{]}_{ij}$ (76)

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. (92) 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

$\displaystyle{\cal H}^{r}($
$\displaystyle\bm{\xi}^{\dagger}(\tau),\bm{\xi}(\tau))=\sum_{\bm{p}}\varepsilon_{\bm{p}}\xi^{*}_{\bm{p}}(\tau)\xi_{\bm{p}}(\tau)$
$\displaystyle+\\!\\!\sum_{\bm{p},\bm{p}^{\prime},\bm{q}}\\!\\!U(\bm{q})\xi^{*}_{\bm{p}+\bm{q}}(\tau)\xi^{*}_{\bm{p}^{\prime}\\!-\\!\bm{q}}(\tau)\xi_{\bm{p}^{\prime}}(\tau)\xi_{\bm{p}}(\tau)$
$\displaystyle\\!-\\!i\hbar\\!\sum_{\bm{p}\bm{p}^{\prime}\bm{k}}\\!\int_{t_{0}}^{\tau}\\!\\!\\!d\tau^{\prime}\\!\xi^{*}_{\bm{p}}(\tau)\xi_{\bm{p}-\bm{k}}(\tau)G_{\bm{k}}(\tau,\tau^{\prime})\xi^{*}_{\bm{p}^{\prime}-\bm{k}}(\tau^{\prime})\xi_{\bm{p}^{\prime}}(\tau^{\prime}).$
(77)

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
$\bm{W}_{\bm{p},\bm{p}^{\prime},\bm{k}}(t,t_{0})\equiv\langle[c^{\dagger}_{\bm{p}-\bm{k}}(t)c_{\bm{p}}(t),c^{\dagger}_{\bm{p}^{\prime}}(t_{0})c_{\bm{p}^{\prime}-\bm{k}}(t_{0})]\rangle$.
Consequently, we have

$\displaystyle\xi^{*}_{\bm{p}-\bm{k}}(t)\xi_{\bm{p}}(t)\propto\sum_{\bm{p}^{\prime\prime}}\bm{W}_{\bm{p},\bm{p}^{\prime},\bm{k}}(t,t_{0})\xi^{*}_{\bm{p}^{\prime}-\bm{k}}(t_{0})\xi_{\bm{p}^{\prime}}(t_{0}).$
(78)

Then the renormalized electron Hamiltonian can be obtained as

$\displaystyle H^{r}_{\rm electron}(t,t_{0})$
$\displaystyle=\sum_{\bm{p}}\varepsilon^{r}_{\bm{p}}(t,t_{0})c^{\dagger}_{\bm{p}}c_{\bm{p}}$
$\displaystyle+\\!\\!\sum_{\bm{p},\bm{p}^{\prime},\bm{q}}\\!\\!U^{r}_{\bm{p}^{\prime}}(\bm{q},t,t_{0})c^{\dagger}_{\bm{p}+\bm{q}}c^{\dagger}_{\bm{p}^{\prime}\\!-\\!\bm{q}}c_{\bm{p}^{\prime}}c_{\bm{p}},$
(79a) where
$\displaystyle\varepsilon^{r}_{\bm{p}}(t,t_{0})=\varepsilon_{\bm{p}}+\sum_{\bm{p}^{\prime}}\delta
U_{\bm{p}^{\prime}}(\bm{q},t,t_{0}),$ (79b) $\displaystyle
U^{r}_{\bm{p}^{\prime}}(\bm{q},t,t_{0})=U(\bm{q})+\delta
U_{\bm{p}^{\prime}}(\bm{q},t,t_{0}),$ (79c) and $\displaystyle\delta
U_{\bm{p}^{\prime}}(\bm{q},t,t_{0})=\hbar{\rm
Im}\\!\sum_{\bm{q}^{\prime}\bm{q}^{\prime\prime}}\\!\int_{t_{0}}^{t}\\!\\!\\!d\tau$
$\displaystyle
G_{\bm{q}}(t,\tau)\bm{W}_{\bm{q}^{\prime},\bm{q}^{\prime\prime},\bm{q}}(\tau,t_{0})$
$\displaystyle\times\bm{W}^{-1}_{\bm{q}^{\prime\prime},\bm{p}^{\prime},\bm{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 $\bm{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 $\bm{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
$\bm{u}(t,t_{0})$ has been solved analytically in our previous work Zhang2012
.

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 non-relativistic 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
Thouless1972  and in quantum field theory Peskin1995 . 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 strongly-correlated 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 Nagaosa2010 .
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 Zhang1992
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 Zhang1990 . 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 system 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 frequency 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.

###### Acknowledgements.

This work is supported by Ministry of Science and Technology of Taiwan,
Republic of China under Contract No. MOST-108-2112-M-006-009-MY3. WMZ proposed
the idea and formulated the theory, WMH performed the numerical calculations,
WMH and WMZ discussed and analyzed the results, WMZ interpreted the physics
and wrote the manuscript.

##  Appendix A The derivation of Eq. (19) from Eq. (18)

In this appendix, we shall derive rigorously the reduced density matrix from
the Gibbs state of the total system at the steady state by trace over all
reservoir states.

The total system is initially in an decoupled state between the system and the
reservoir, given by Eq. (3). After a long-time nonequilibrium evolution, the
total system approaches to the steady equilibrium state which is the Gibbs
state with a final equilibrium temperature denoted by
$\beta_{f}=1/k_{B}T_{f}$, i.e. Eq. (18):

$\displaystyle\rho_{\rm tot}(t\\!\rightarrow\\!\infty)=\frac{1}{Z_{\rm
tot}}e^{-\beta_{f}H_{\rm tot}},$ (80)

where $H_{\rm tot}$ is the total Hamiltonian of Eq. (1), and $Z_{\rm tot}={\rm
Tr}_{{{}_{\\!S}}+{{}_{\\!E}}}[e^{-\beta_{f}H_{\rm tot}}]$ is the corresponding
partition function of the total system. This is also a direct consequence of
statistical mechanics, namely when the total system is in equilibrium, its
state is given by the Bolzertmann distribution, i.e., Eq. (80). Because
$H_{\rm tot}$ is a bilinear operator in terms of the bosonic creation and
annihilation operators, Eq. (80) becomes a Gaussian function in the coherent-
state representation. Explicitly, Eq. (80) can be expressed as Huang2020 ,

$\displaystyle\langle z,\bm{z}|\rho_{\rm
tot}|z^{\prime},\bm{z}^{\prime}\rangle\\!=\\!\frac{1}{Z_{\rm
tot}}\exp\bigg{\\{}\\!\\!\begin{pmatrix}z^{\dagger}&\\!\bm{z}^{\dagger}\end{pmatrix}\\!\\!\begin{pmatrix}\Omega_{{{}_{\\!S}}{{}_{\\!S}}}&\bm{\Omega}_{{{}_{\\!S}}{{}_{\\!E}}}\\\
\bm{\Omega}_{{{}_{\\!E}}{{}_{\\!S}}}&\bm{\Omega}_{{{}_{\\!E}}{{}_{\\!E}}}\end{pmatrix}\\!\\!\begin{pmatrix}z^{\prime}\\\
\bm{z}^{\prime}\end{pmatrix}\\!\\!\bigg{\\}},$ (81)

where

$\displaystyle\begin{pmatrix}\Omega_{{{}_{\\!S}}{{}_{\\!S}}}&\bm{\Omega}_{{{}_{\\!S}}{{}_{\\!E}}}\\\
\bm{\Omega}_{{{}_{\\!E}}{{}_{\\!S}}}&\bm{\Omega}_{{{}_{\\!E}}{{}_{\\!E}}}\end{pmatrix}=\exp\\!\bigg{[}\\!-\\!\beta_{f}\begin{pmatrix}\hbar\omega_{s}&\hbar\bm{V}\\\
\hbar\bm{V}^{\dagger}&\hbar\bm{\omega}\end{pmatrix}\\!\bigg{]}.$ (82)

Here, we have used the combined bosonic coherent state of the system plus the
reservoir
$|z,\bm{z}\rangle=\exp(za^{\dagger}+\sum_{k}z_{k}b^{\dagger}_{k})|0\rangle$.
We also used boldface to denote matrices and vectors. As an explicit example,
the vector $\bm{z}\equiv(z_{k_{1}},z_{k_{2}},...)$ and $z_{k_{i}}$ is a
complex variable for reservoir boson mode $k_{i}$.

Taking the trace over all the reservoir modes, we obtain the reduced density
matrix in the coherent state representation,

$\displaystyle\langle
z|\rho_{{}_{\\!S}}(t\\!\rightarrow\\!\infty)|z^{\prime}\rangle=\\!\\!\int\\!\\!d\mu(\bm{z})\langle
z,\bm{z}|\rho_{\rm tot}(t\\!\rightarrow\\!\infty)|z^{\prime},\bm{z}\rangle,$
(83)

where $d\mu(\bm{z})=\prod_{k}\frac{dz^{*}_{k}dz_{k}}{2\pi i}e^{-|z_{k}|^{2}}$.
Use the Gaussian integral to complete the integration of Eq. (83), we have

$\displaystyle\langle
z|\rho_{{}_{\\!S}}(t\\!\rightarrow\\!\infty)|z^{\prime}\rangle=\frac{1}{Z_{{}_{\\!S}}}e^{z^{*}\Omega_{{}_{\\!S}}z^{\prime}}$
(84)

where $Z_{{}_{\\!S}}=[1-\Omega_{{}_{\\!S}}]^{-1}$ and
$\Omega_{{}_{\\!S}}=\Omega_{{{}_{\\!S}}{{}_{\\!S}}}+\bm{\Omega}_{{{}_{\\!S}}{{}_{\\!E}}}[\bm{I}-\bm{\Omega}_{{{}_{\\!E}}{{}_{\\!E}}}]^{-1}\bm{\Omega}_{{{}_{\\!E}}{{}_{\\!S}}}$.
On the other hand, the average particle number in the steady state is
$\overline{n}={\rm
tr}_{{}_{\\!S}}[a^{\dagger}a\rho_{{}_{\\!S}}(t\\!\rightarrow\\!\infty)]=\\!\\!\int\\!\\!d\mu(\bm{z})\langle
z|a^{\dagger}a\rho_{{}_{\\!S}}(t\\!\rightarrow\\!\infty)|z\rangle$. The result
is

$\displaystyle\overline{n}(t\\!\rightarrow\\!\infty)$
$\displaystyle=\Omega_{{}_{\\!S}}[1-\Omega_{{}_{\\!S}}]^{-1}$ (85)

Thus, the reduced density matrix becomes

$\displaystyle\langle
z|\rho_{{}_{\\!S}}(t\\!\rightarrow\\!\infty)|z^{\prime}\rangle=\frac{1}{Z_{{}_{\\!S}}}e^{z^{*}\frac{\bar{n}(t\\!\rightarrow\\!\infty)}{1+\bar{n}(t\\!\rightarrow\\!\infty)}z^{\prime}}.$
(86)

Using the fact that $\langle
z|e^{fa^{\dagger}a}|z^{\prime}\rangle=e^{z^{*}e^{f}z^{\prime}}$, the above
reduced density matrix can be written as a operator,

$\displaystyle\rho_{{}_{\\!S}}(t\\!\rightarrow\\!\infty)=\frac{\exp\\{\ln[\frac{\bar{n}(t\\!\rightarrow\\!\infty)}{1+\bar{n}(t\\!\rightarrow\\!\infty)}]a^{\dagger}a\\}}{1+\bar{n}(t\\!\rightarrow\\!\infty)},$
(87)

which is the solution of Eq. (19). This solution is obtained directly from Eq.
(80). As we have found, the exact solution of
$\bar{n}(t\\!\rightarrow\\!\infty)$ does not agree with the Bose-Einstein
distribution if one takes the steady-state temperature as the initial
reservoir temperature $T_{0}$ at strong coupling, even though the energy
correction arisen from the strong coupling with the reservoir is properly
included. Therefore, $T_{f}\neq T_{0}$ and temperature renormalization is
necessary for strong coupling quantum thermodynamics. Note that a similar
proof claimed in Subas2012  without providing any details should be incorrect.

##  Appendix B Exact solution and the steady state of open quantum systems
with Eq. (38)

In this Appendix, we present the general solution of open quantum systems with
Eq. (38). Without loss of the generality, the initial state of reservoir
$\alpha$ can be assumed to be in a thermal state at temperature $T_{\alpha
0}$, and the system can be an arbitrary state,

$\displaystyle\rho_{\rm
tot}(t_{0})=\rho_{{}_{\\!S}}(t_{0})\otimes\rho_{{}_{\\!E}}(t_{0}),~{}~{}~{}\rho_{{}_{\\!E}}(t_{0})=\prod_{\alpha}\frac{1}{Z_{\alpha}}e^{-\beta_{\alpha}H^{\alpha}_{{{}_{\\!E}}}}.$

Here $Z_{\alpha}=\prod_{k}(1\mp e^{-\beta_{\alpha 0}(\varepsilon_{\alpha
k}-\mu_{\alpha 0})})^{\mp}$ is the partition function of reservoir $\alpha$.
Different reservoir $\alpha$ could have differential initial chemical
potential $\mu_{\alpha 0}$ and different initial temperature $\beta_{\alpha
0}=1/k_{B}T_{\alpha 0}$. The up and down signs $\mp$ 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
$\rho_{{}_{\\!S}}(t_{0})$ of the system, the reduced density matrix at later
time $t$ is defined by $\rho_{{}_{\\!S}}(t)={\rm Tr}_{{}_{\\!E}}[\rho_{\rm
tot}(t)]$, which can be solved from Eq. (39) in general. To find the exact
solution, we take the coherent state representation Zhang1990  again. Then,
the reduced density matrix $\rho_{{}_{\\!S}}\left(t\right)$ of Eq. (39) can be
expressed as
$\rho_{{}_{\\!S}}\\!(\bm{\alpha}_{t},\bm{\alpha}_{t}^{\prime},t)=\langle\bm{\alpha}_{t}|\rho_{{}_{\\!S}}(t)|\bm{\alpha}^{\prime}_{t}\rangle$.
Here we have used the unnormalized coherent state defined as:
$|\bm{\alpha}\rangle=\exp(\sum_{i}\alpha_{i}a^{\dagger}_{i})|0\rangle,~{}d\mu\left(\bm{\alpha}\right)=\prod_{i}g_{i}d\alpha^{*}_{i}d\alpha_{i}e^{-\left|\alpha_{i}\right|^{2}}$,
the vector $\bm{\alpha}\equiv(\alpha_{1},\alpha_{2},...)$ is one-column
matrix, and $\alpha_{i}$ are complex variables for bosons and Grassmannian
variables for fermions with $g_{i}=1/2\pi i$ and $1$ in the Haar measure,
respectively. Thus, the reduced density matrix in the coherent state
representation can be expressed as

$\displaystyle\rho_{{}_{\\!S}}\\!(\bm{\alpha}_{t},\bm{\alpha}_{t}^{\prime},t)=\\!\\!\int\\!\\!d\mu\left(\bm{\alpha}_{0}\right)d\mu\left(\bm{\alpha}_{0}^{\prime}\right)\rho_{{}_{\\!S}}\left(\bm{\alpha}_{0},\bm{\alpha}_{0}^{\prime},t_{0}\right)$
$\displaystyle\times\mathcal{J}\\!\\!\left(\bm{\alpha}_{t},\bm{\alpha}_{t}^{\prime},t;\bm{\alpha}_{0},\bm{\alpha}_{0}^{\prime},t_{0}\right).$
(88)

The propagating function
$\mathcal{J}\\!\\!\left(\bm{\alpha}^{\dagger}_{t},\bm{\alpha}_{t}^{\prime},t;\bm{\alpha}_{0},{\bm{\alpha}^{\prime}_{0}}^{{\dagger}},t_{0}\right)$
in Eq. (88) can be obtained analytically after integrated exactly over all the
environmental degrees of freedom using the coherent state path integrals. The
result is Tu2008 ; Jin2010 ; Lei2012

$\displaystyle\mathcal{J}(\bm{\alpha}^{\dagger}_{t},\bm{\alpha}^{\prime}_{t},t;$
$\displaystyle\bm{\alpha}_{0},{\bm{\alpha}^{\prime}_{0}}^{{\dagger}},t_{0})=\int\\!\mathcal{D}[\bm{\alpha};\bm{\alpha}^{\prime}]\exp\Big{\\{}\frac{i}{\hbar}\big{(}S_{s}[\bm{\alpha}^{\dagger},\bm{\alpha}]$
$\displaystyle-S^{*}_{s}[\bm{\alpha}^{\prime{\dagger}},\bm{\alpha}^{\prime}]+S_{\rm
IF}[\bm{\alpha}^{\dagger},\bm{\alpha};\bm{\alpha}^{\prime{\dagger}},\bm{\alpha}^{\prime}]\big{)}\\!\Big{\\}},$
(89)

where $\mathcal{D}[\bm{\alpha};\bm{\alpha}^{\prime}]$ is the path integral
measure over the parameter space of the system coherent states
$|\bm{\alpha}\rangle$. The bare system action in the coherent state basis is
given by

$\displaystyle S_{s}[\bm{\alpha}^{\dagger},\bm{\alpha}]=$
$\displaystyle-{i\hbar\over
2}\big{[}\bm{\alpha}^{\dagger}_{t}\bm{\alpha}(t)+\bm{\alpha}^{\dagger}(t_{0})\bm{\alpha}_{t_{0}}\big{]}$
$\displaystyle+\\!\\!\int_{t_{0}}^{t}\\!\\!d\tau\Big{\\{}{i\hbar\over
2}\big{[}\dot{\bm{\alpha}}^{\dagger}(\tau)\bm{\alpha}(\tau)-\bm{\alpha}^{\dagger}(\tau)\dot{\bm{\alpha}}(\tau)\big{]}$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}-{\cal
H}(\bm{\alpha}^{\dagger}(\tau),\bm{\alpha}(\tau))\Big{\\}}$ (90)

with the Hamiltonian function ${\cal
H}(\bm{\alpha}^{\dagger}(\tau),\bm{\alpha}(\tau))=\sum_{i}\varepsilon_{i}(\tau)\alpha^{*}_{i}(\tau)\alpha_{i}(\tau)$.
The actions $S_{s}[\bm{\alpha}^{\dagger},\bm{\alpha}]$ and
$S^{*}_{s}[\bm{\alpha}^{\prime{\dagger}},\bm{\alpha}^{\prime}]$ in Eq. (89
describe the forward and backward evolution of the system. The influence
functional action $S_{\rm
IF}[\bm{\alpha}^{\dagger},\bm{\alpha};\bm{\alpha}^{\prime{\dagger}},\bm{\alpha}^{\prime}]$
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 Feynman1963 . 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 Tu2008 ;
Jin2010 ; Lei2012 . The resulting influence functional action for the open
quantum system Eq. (38) is

$\displaystyle S_{\rm
IF}[\bm{\alpha}^{\dagger},\bm{\alpha};\bm{\alpha}^{\prime{\dagger}},\bm{\alpha}^{\prime}]=\\!\\!\int_{t_{0}}^{t}\\!\\!\\!d\tau\bigg{\\{}i\hbar\sum_{ij}$
$\displaystyle\Big{[}\int_{t_{0}}^{\tau}\\!\\!\\!d\tau^{\prime}\\!\Big{\\{}\alpha^{*}_{i}(\tau)g_{ij}(\tau,\tau^{\prime})\alpha_{j}(\tau^{\prime})+\alpha^{\prime*}_{i}(\tau)g^{*}_{ij}(\tau,\tau^{\prime})\alpha^{\prime}_{j}(\tau^{\prime})$
$\displaystyle\mp\\!\\!\int_{t_{0}}^{t}\\!\\!\\!d\tau^{\prime}\\!\Big{\\{}\alpha^{\prime*}_{i}(\tau)g_{ij}(\tau,\tau^{\prime})\alpha_{j}(\tau^{\prime})\\!-\\!\big{[}\alpha^{*}_{i}(\tau)\mp\alpha^{\prime*}_{i}(\tau)\big{]}\widetilde{g}_{ij}(\tau,\tau^{\prime})\big{[}\alpha_{j}(\tau^{\prime})\mp\alpha^{\prime}_{j}(\tau^{\prime})\big{]}\\!\Big{\\}}\Big{]}\\!\bigg{\\}},$
(91)

where the up and down signs of $\mp$ correspond respectively to the system
being bosonic and fermionic. The system-reservoir correlation functions
$g_{ij}(\tau,\tau^{\prime})$ and $\widetilde{g}_{ij}(\tau,\tau^{\prime})$ 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
${\alpha^{*}_{i},\alpha_{i}}$ so that the path integrals of Eq. (89) can be
solved exactly. After a tedious calculation, we obtain Tu2008 ; Jin2010 ;
Lei2012

$\displaystyle\mathcal{J}(\bm{\alpha}^{\dagger}_{t},$
$\displaystyle\bm{\alpha}^{\prime}_{t},t;\bm{\alpha}_{0},{\bm{\alpha}^{\prime}_{0}}^{{\dagger}},t_{0})=\big{(}{\rm
det}[\bm{w}(t)]\big{)}^{\pm 1}$
$\displaystyle\times\exp\Big{\\{}\bm{\alpha}_{t}^{{\dagger}}\bm{J}_{1}(t,t_{0})\bm{\alpha}_{0}+{\bm{\alpha}^{\prime}}_{0}^{{\dagger}}\bm{J}^{\dagger}_{1}(t,t_{0})\bm{\alpha}^{\prime}_{t}$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\pm{\bm{\alpha}_{0}^{\prime}}^{\dagger}\bm{J}_{3}(t,t_{0})\bm{\alpha}_{0}\pm\bm{\alpha}_{t}^{{\dagger}}\bm{J}_{2}(t)\bm{\alpha}^{\prime}_{t}\Big{\\}},$
(92)

where $\bm{w}\left(t\right)=[\bm{1}\pm\bm{v}\left(t,t\right)]^{-1}$,
$\bm{J}_{1}(t,t_{0})=\bm{w}(t)\bm{u}(t,t_{0})$,
$\bm{J}_{2}(t)=\bm{1}-\bm{w}(t)$ and
$\bm{J}_{3}(t,t_{0})=\bm{1}-\bm{u}^{\dagger}(t,t_{0})\bm{w}(t)\bm{u}(t,t_{0})$,
and the up and down signs of $\pm$ correspond to the system being bosonic or
fermionic systems again. The functions $\bm{u}(t,t_{0})$ and $\bm{v}(t,t)$ are
the nonequilibrium dissipative particle-propagating Green’s function and
fluctuated particle-correlation Green’s function, respectively. They are
determined by Eq. (43). Equations (88) and (92) give the exact solution of the
reduced density matrix in coherent state representation for the open quantum
systems by Eq. (38). From Eqs. (88)-(92), we have derived rigorously the exact
master equation Eq. (40) Tu2008 ; Jin2010 ; Lei2012 .