# Numerical study of non-adiabatic quantum thermodynamics of the driven
resonant level model: Non-equilibrium entropy production and higher order
corrections

Kaiyi Tong  School of Science, Westlake University, Hangzhou, Zhejiang 310024,
China  Institute of Natural Sciences, Westlake Institute for Advanced Study,
Hangzhou, Zhejiang 310024, China  Wenjie Dou<EMAIL_ADDRESS>School
of Science, Westlake University, Hangzhou, Zhejiang 310024, China  Institute
of Natural Sciences, Westlake Institute for Advanced Study, Hangzhou, Zhejiang
310024, China

###### Abstract

We present our numerical study on quantum thermodynamics of the resonant level
model subjected to non-equilibrium condition as well as external driving.
Following our previous work on non-equilibrium quantum thermodynamics (Phys.
Rev. B 101, 184304 [2020]), we expand the density operator into a series of
power in the driving speed, where we can determine the non-adiabatic
thermodynamic quantities. Particularly, we calculate the non-equilibrium
entropy production rate as well as higher order non-adiabatic corrections to
the energy and/or population. In the limit of weak system-bath coupling, our
results reduce to the one from the quantum master equation.

##  I Introduction

The study of dynamics and thermodynamics for a quantum system strongly coupled
to a set of baths are of great interests recently,Gemmer  _et al._ (2009);
Kosloff (2013); Esposito  _et al._ (2015a); Vinjanampathy and Anders (2016);
Anders and Esposito (2017); Carrega  _et al._ (2016); Alicki and Kosloff
(2018); Benenti  _et al._ (2017); Campisi  _et al._ (2011, 2009); Liu  _et
al._ (2021) particularly due to its applications in nano technology, quantum
information, and quantum measurement.Poot and van der Zant (2012); Pekola
(2015); Roßnagel  _et al._ (2016); Argun  _et al._ (2017); Klatzow  _et al._
(2019) Quantum thermodynamics address the energy and information flow of a
system consisting of a few atoms or qubits interacting strongly to non-
equilibrium environments.Esposito  _et al._ (2015b); Millen and Xuereb (2016);
Gemmer  _et al._ (2009); Talkner  _et al._ (2007); Brandao  _et al._ (2015);
Jarzynski (2011) Moreover, to make a useful quantum engine, external driving
is usually applied. The finite speed driving can introduce non-adiabatic
effects, including entropy production, friction (or dissipation), and
fluctuation. Solinas  _et al._ (2013); Schmidt  _et al._ (2015); Gogolin and
Eisert (2016); Subaşı  _et al._ (2012); Ness (2017); Strasberg  _et al._
(2017)

The driven resonant-level model has been studied extensively for
thermodynamics in the strongly coupled regimes.Esposito  _et al._ (2015b, a);
Bruch  _et al._ (2016); Haughian  _et al._ (2018) At equilibrium (with one
bath or no current), studies based on different methods (e.g symmetric
splittingBruch  _et al._ (2016); Ochoa  _et al._ (2016), scattering
matrixBruch  _et al._ (2018), and non-equilibrium Green’s functionEsposito
_et al._ (2015a)) arrive at a (somewhat) consistent quantum description of the
thermodynamics. In particular, the first order non-adiabatic corrections to
work, population, and entropy have been identified. Furthermore, the entropy
production is proportional to the frictional work and remains positive at
equilibrium, which is consistent with the second law of thermodynamics. That
being said, out of equilibrium (baths with different temperatures or chemical
potentials), no such formulations have been agreed upon.Semenov and Nitzan
(2020); Bergmann and Galperin (2021); Strasberg and Winter (2021); Bruch  _et
al._ (2018) Particularly the definition and the positivity of entropy
production remain as open questions. Furthermore, there are far less studies
on the higher order non-adiabatic corrections to the thermodynamic quantities,
which is a challenging task.

In the previous studies,Dou  _et al._ (2018, 2020) by expanding the density
operator into a series of power in the driving speed, we formulate a general
description of quantum thermodynamics for a generic model strongly coupled to
one or more baths. We identify the non-adiabatic corrections to thermodynamic
quantities, such as work, population, and entropy. In particular, the non-
adiabatic entropy production rate is given byDou  _et al._ (2020)

$\displaystyle\delta\dot{S}=-k_{B}\sum_{\alpha\nu}\int_{0}^{\infty}\dot{R}_{\alpha}\dot{R}_{\nu}tr(e^{-i\hat{H}t^{\prime}/\hbar}\partial_{\nu}\hat{\rho}_{ss}e^{i\hat{H}t^{\prime}/\hbar}\partial_{\alpha}\ln\hat{\rho}_{ss})dt^{\prime}.$
(1)

Here, $\hat{\rho}_{ss}$ is the steady state density of the total system.
$\hat{H}$ is the total Hamiltonian, which depends on a set of external
parameters $\\{R_{\alpha}\\}$. $\dot{R}_{\alpha}$ is the driving speed. The
above equation can be recasted into a Kubo transformed correlation function
using the Barker-Campell-Hassdrauff formula. Such that we have proven that the
entropy production rate is always positive, which is consistent with the
second law of thermodynamics. At equilibrium, the steady state density reduces
to equilibrium density, such that we can show that the non-adiabatic entropy
is proportional to the friction (or non-adiabatic correction to the work). Out
of equilibrium, however, we have not been able to calculate $\delta\dot{S}$
explicitly due to the difficulty of determining $\hat{\rho}_{ss}$.

In the present manuscript, we calculate this entropy production rate
explicitly for the driven resonant level model, where we discretize the
continuous bath degrees of freedom and diagonalize the total Hamiltonian
numerically. In doing so, we can identify the steady state density and
calculate relevant thermodynamic quantities. We show that the numerical
results recover analytical solutions when available. We further identify the
second order correction to the populations using analytical analysis as well
as the solutions from hierarchical quantum master equation (HQME) Schinabeck
and Thoss (2020); Dou  _et al._ (2020); Jin  _et al._ (2008); Tanimura (2020).

We organize the manuscript as follows. In II, we present our analysis as well
as our numerical methods to calculate thermodynamic quantities for the
resonant level model. In III, we show our numerical results in non-adiabatic
corrections to population, work as well entropy production. We conclude our
work in IV.

##  II Theory

###  II.1 The driven resonant-level model

The driven resonant-level model is probably the simplest yet a heuristic model
for the study of quantum thermodynamics in the strongly coupled regime. The
model consists of a system, two baths, as well as the couplings between them:

$\hat{H}=\hat{H}_{\text{sys}}+\hat{H}_{\text{bath}}+\hat{H}_{\text{int}},$ (2)

The system consists of one level with a time-dependent on-site energy:

$\displaystyle\hat{H}_{\text{sys}}$ $\displaystyle=$
$\displaystyle\epsilon_{d}(t)\hat{d}^{\dagger}\hat{d},$ (3)

The bath consists of a set of non-interacting Fermions:

$\displaystyle\hat{H}_{\text{bath}}$ $\displaystyle=$
$\displaystyle\sum_{k\alpha}\epsilon_{k\alpha}\hat{c}_{k\alpha}^{\dagger}\hat{c}_{k\alpha},$
(4)

Here $\alpha=L$ (or $R$) indicates the left (or right) bath. The couplings
between the system and the baths are bilinear:

$\displaystyle\hat{H}_{\text{int}}$ $\displaystyle=$
$\displaystyle\sum_{k\alpha}V_{k\alpha}(\hat{c}_{k\alpha}^{\dagger}\hat{d}+\hat{d}^{\dagger}\hat{c}_{k\alpha}).$
(5)

To describe the strength of the coupling, we define the hybridization function
$\Gamma_{\alpha}(\epsilon)$

$\Gamma(\epsilon)=\sum_{\alpha}\Gamma_{\alpha}(\epsilon)=\sum_{k\alpha}2\pi|V_{k\alpha}|^{2}\delta(\epsilon-\epsilon_{k\alpha}).$
(6)

Below, we consider the wide-band limit, such that $\Gamma$ is a constant
(independent of energy $\epsilon$).

Figure 1:  An illustrative picture of the driving resonant-level model.

###  II.2 Analytical analysis

We now try to define the quantum thermodynamic quantities for the driven
resonant-level model. The main quantities of interests are dot population,
work, and entropy production:

$\displaystyle N=tr(\hat{\rho}\hat{d}^{\dagger}\hat{d})$ (7) $\displaystyle
W=tr(\hat{\rho}\partial_{t}\hat{H})$ (8) $\displaystyle
S=-k_{B}tr(\hat{\rho}\ln\hat{\rho})$ (9)

Here $\hat{\rho}$ is the density operator. Note that, we have avoided defined
heat and/or energy as they are extensive quantities and requires spectral
treatment.Dou  _et al._ (2018, 2020) Since the Hamiltonian is quadratic, we
can rewrite the above quantities using the one-body density operator
$\hat{\sigma}$:

$\displaystyle N$ $\displaystyle=$ $\displaystyle\langle
d|\hat{\sigma}|d\rangle$ (10) $\displaystyle W$ $\displaystyle=$
$\displaystyle\dot{\epsilon}_{d}\langle
d|\hat{\sigma}|d\rangle=\dot{\epsilon}_{d}N$ (11) $\displaystyle S$
$\displaystyle=$
$\displaystyle-\text{Tr}(\hat{\sigma}\ln\hat{\sigma})-\text{Tr}((1-\hat{\sigma})\ln(1-\hat{\sigma}))$
(12)

Here Tr denotes the one-body trace (trace in orbital representation). Note
that for the resonant level model, work $W$ is proportional to dot population
by a factor of $\dot{\epsilon}_{d}$, where $\dot{\epsilon}_{d}$ is the speed
of driving. Such that we will mainly focus on population and entropy below.

The task is then to calculate the density operator, which can be very
difficult for the time-dependent Hamiltonian. Nevertheless, we have the
equation of motion for the density matrix

$\displaystyle\partial_{t}\hat{\sigma}=\dot{\epsilon}_{d}\partial_{\epsilon_{d}}\hat{\sigma}-\frac{i}{\hbar}[\hat{h},\hat{\sigma}]$
(13)

Here $\hat{h}$ is the one body Hamiltonian. In the limit of slow driving, we
can expand the density operator into a series of power in the driving speed,

$\hat{\sigma}=\hat{\sigma}^{(0)}+\hat{\sigma}^{(1)}+\hat{\sigma}^{(2)}+\cdots,$
(14)

By matching the order in driving speed from both sides of the above equation,
we arrive at a set of equations:

$\displaystyle\partial_{t}\hat{\sigma}^{(0)}$ $\displaystyle=$
$\displaystyle-\frac{i}{\hbar}[\hat{h},\hat{\sigma}^{(0)}],$ (15)
$\displaystyle\partial_{t}\hat{\sigma}^{(n)}$ $\displaystyle=$
$\displaystyle\dot{\epsilon}_{d}\partial_{\epsilon_{d}}\sigma^{(n-1)}-\frac{i}{\hbar}[\hat{h},\sigma^{(n)}],n\geq
1$ (16)

The steady state solution of Eq. 15 gives us the $0$th order density matrix,
$\hat{\sigma}^{(0)}=\hat{\sigma}_{ss}$, such that we can proceed to calculate
the higher order density matrix

$\displaystyle\hat{\sigma}^{(n)}$
$\displaystyle=-\dot{\epsilon}_{d}\int_{0}^{t}e^{-i\hat{h}(t-t^{\prime})/{\hbar}}\partial_{\epsilon_{d}}\hat{\sigma}^{(n-1)}e^{{i\hat{h}(t-t^{\prime})}/{\hbar}}dt^{\prime}$
(17)
$\displaystyle\approx-\dot{\epsilon}_{d}\int_{0}^{\infty}e^{-{i\hat{h}t}/{\hbar}}\partial_{\epsilon_{d}}\hat{\sigma}^{(n-1)}e^{{i\hat{h}t}/{\hbar}}dt.$

In the second line of the above equation, we have used the Markovian
approximation, assuming the timescale of driving is much smaller as compared
with the timescale for the relaxation. The Markovian approximation is
consistent with the slow driving approximation.

####  II.2.1 Zeroth Order Quantities

Since the zeroth order entropy is an extensive quantity, we mainly calculate
the zeroth order dot population here. To do so, we need to diagonalize the
one-body Hamiltonian:

$\hat{h}=\sum_{k\alpha}\epsilon_{k\alpha}|{\psi}_{k\alpha}\rangle\langle{\psi}_{k\alpha}|,$
(18)

The diagonalization can be done analytically through the following
transformation:

$|{\psi}_{k\alpha}\rangle=|{c}_{k\alpha}\rangle+\frac{V_{k\alpha}}{\epsilon_{k\alpha}-\epsilon_{d}+i\Gamma/2}\left(|d\rangle+\sum_{k^{\prime}\alpha^{\prime}}\frac{V_{k^{\prime}\alpha^{\prime}}}{\epsilon_{k\alpha}-\epsilon_{k^{\prime}\alpha^{\prime}}+i\eta)}|{c}_{k^{\prime}\alpha^{\prime}}\rangle\right).$
(19)

With the diagonalized Hamiltonian, we can determine the steady-state single
particle density matrix

$\displaystyle\hat{\sigma}^{(0)}=\sum_{k\alpha}f(\epsilon_{k\alpha}-\mu_{\alpha})|{\psi}_{k\alpha}\rangle\langle{\psi}_{k\alpha}|$
(20)

Here, $f$ is the Fermi function
$f(\epsilon_{k\alpha}-\mu_{\alpha})=\left[1+e^{\beta(\epsilon_{k\alpha}-\mu_{\alpha})}\right]^{-1}$,
and $\mu_{\alpha}$ is the chemical potential for the $\alpha$ lead.

With the analytical results shown above, we can obtain the zeroth order
population:

$\displaystyle
N^{(0)}=\left<d\right|\hat{\sigma}^{(0)}\left|d\right>=\frac{1}{2\pi}\int
d\epsilon A(\epsilon)\bar{f}(\epsilon),$ (21)

Here $A(\epsilon)$ is the spectral function

$\displaystyle
A(\epsilon)=\frac{\Gamma}{(\epsilon-\epsilon_{d})^{2}+(\Gamma/2)^{2}},$ (22)

and $\bar{f}(\epsilon)$ is the weighted Fermi function
$\bar{f}(\epsilon)=\sum_{\alpha}\frac{\Gamma_{\alpha}}{\Gamma}\left[1+e^{\beta(\epsilon-\mu_{\alpha})}\right]^{-1}$.
Details of the deviation are shown in Appendix A.

####  II.2.2 First Order Quantities

Now we turn to the first order correction. With $0$th order
$\hat{\sigma}^{(1)}$, we can proceed to calculate the first order density:

$\hat{\sigma}^{(1)}=-\dot{\epsilon}_{d}\int_{0}^{\infty}e^{-{i\hat{h}t}/{\hbar}}\partial_{\epsilon_{d}}\hat{\sigma}^{(0)}e^{{i\hat{h}t}/{\hbar}}dt.$
(23)

The first order correction to the population is then given by

$N^{(1)}=\left<d\right|\hat{\sigma}^{(1)}\left|d\right>=-\frac{\hbar\dot{\epsilon}_{d}}{4\pi}\int
d\epsilon A^{2}{\partial}_{\epsilon}\bar{f}(\epsilon).$ (24)

Detailed derivation can be found in the appendix of Refs. (35; 40). Note that
the population is linearly proportional to driving speed $\dot{\epsilon}_{d}$.
Such that the first order work depends on $\dot{\epsilon}_{d}$ quadratically,
$W^{(1)}=\gamma^{(1)}\dot{\epsilon}^{2}_{d}$. Here $\gamma^{(1)}$ is the
frictional coefficient:

$\gamma^{(1)}=\frac{N^{(1)}}{\dot{\epsilon}_{d}}=-\frac{\hbar}{4\pi}\int
d\epsilon A^{2}{\partial}_{\epsilon}\bar{f}(\epsilon),$ (25)

which accounts for the dissipative effects.

We could also define a correlation function of the random force to qualify the
fluctuation

$\displaystyle D_{\mu\nu}$ $\displaystyle=$
$\displaystyle\int_{0}^{\infty}\langle
e^{{i\hat{H}t}/{\hbar}}\delta\hat{F}_{\mu}e^{-{i\hat{H}t}/{\hbar}}\delta\hat{F}_{\nu}(0)\rangle_{S}$
(26) $\displaystyle\delta\hat{F}_{\alpha}$ $\displaystyle=$
$\displaystyle{\partial}_{\alpha}\hat{H}-tr({\partial}_{\alpha}\hat{H}\hat{\rho}_{ss})$
(27)

$\langle\cdots\rangle_{S}$ denote the symmetric average of the correlation
functions. Here $\alpha$ and $\mu$ are degrees of freedom for the motion. In
the driven resonant level model (with $\alpha=\mu=\epsilon_{d}$), we can
calculate the fluctuation analytically

$D=\frac{\hbar}{4\pi}\int d\epsilon
A^{2}\bar{f}(\epsilon)(1-\bar{f}(\epsilon)).$ (28)

We then turn to the first order correction to entropy, which can be calculated
through Taylor expansion of Eq. 12:

$S^{(1)}=-k_{B}\text{Tr}\left[\hat{\sigma}^{(1)}\ln\hat{\sigma}^{(0)}-\hat{\sigma}^{(1)}\ln{(1-\hat{\sigma}^{(0)})}\right],$
(29)

The non-adiabatic entropy rate is the derivative of the above equation.

$\displaystyle\dot{S}^{(2)}=\partial_{t}S^{(1)}=\frac{\dot{Q}^{(2)}}{T}+\delta\dot{S},$
(30)

The first term in Eq. 30 is the first order correction to the heat Dou  _et
al._ (2020)

$\dot{Q}^{(2)}=-k_{B}T\dot{\epsilon}_{d}\text{Tr}(\left[{\partial}_{\epsilon_{d}}{\sigma}^{(1)}\ln{{\sigma}^{(0)}}-{\partial}_{\epsilon_{d}}{\sigma}^{(1)}\ln{(1-{\sigma}^{(0)})}\right].$
(31)

The second term $\delta\dot{S}$ in Eq. 30 is the additional entropy production
rate due to the external driving

$\displaystyle\delta\dot{S}=$ $\displaystyle-
k_{B}\dot{\epsilon}_{d}\left[\hat{\sigma}^{(1)}{\partial}_{\epsilon_{d}}\ln\hat{\sigma}^{(0)}-\hat{\sigma}^{(1)}{\partial}_{\epsilon_{d}}\ln{(1-\hat{\sigma}^{(0)})}\right]$
(32)

At equilibrium, we can calculate the entropy production rate explicitly and
show that

$\gamma^{(1)}=\beta D=\frac{\hbar T\delta\dot{S}}{\dot{\epsilon}_{d}}.$ (33)

The first equality in the above equation is the fluctuation-dissipation
theorem. Under non-equilibrium condition, the fluctuation-dissipation theorem
is no longer obeyed and we have not been able to determine the entropy
production rate. In Sec. II.3, we will show how to evaluate entropy production
rate numerically.

####  II.2.3 Second Order Quantities

Evaluation of the second order corrections to thermodynamic quantities is very
tricky, which will require

$\hat{\sigma}^{(2)}=-\dot{\epsilon}_{d}\int_{0}^{\infty}e^{-{i\hat{h}t}/{\hbar}}{\partial}_{\epsilon_{d}}{\sigma}^{(1)}e^{{i\hat{h}t}/{\hbar}}dt.$
(34)

Nevertheless, as shown in Appendix B, we manage to evaluate the second order
correction to the population analytically for the equilibrium case,

$N^{(2)}=\left<d\right|\hat{\sigma}^{(2)}\left|d\right>=\frac{\hbar^{2}\dot{\epsilon}_{d}^{2}}{12\pi}\int
d\epsilon A^{3}{\partial}_{\epsilon}^{2}{f}(\epsilon).$ (35)

Similarly to the first order case, we can define the friction term
$\gamma^{(2)}$ as

$\gamma^{(2)}=\frac{N^{(2)}}{\dot{\epsilon}_{d}^{2}}=\frac{\hbar^{2}}{12\pi}\int
d\epsilon A^{3}{\partial}_{\epsilon}^{2}{f}(\epsilon).$ (36)

Out of equilibrium, we do not have analytical results for the second order
corrections. We could simply replace $f$ by $\bar{f}$ in the above equation as
our trial results. In general, however, the higher order correction can be
evaluated numerically using hierarchical quantum master equation (HQME).

###  II.3 Numerical Methods

The derivations in evaluating the non-adiabatic corrections can be very
lengthy and a compact form of solutions may not available for certain
quantities. In this subsection, we consider a finite number of degrees of
freedom for the baths in a similar fashion to the Ref. (41), such that we can
evaluate non-adiabatic corrections numerically.

To do so, we discrete the baths degrees of freedom, and we build the one-body
Hamiltonian $\hat{h}$ with finite number of levels from the baths. $\hat{h}$
can be written in the matrix form:

$\displaystyle\hat{h}=\begin{pmatrix}\epsilon_{d}&\mathbb{V_{L}}&\mathbb{V_{R}}\\\
\mathbb{V_{L}}^{\dagger}&\hat{\epsilon}_{L}&0\\\
\mathbb{V_{R}}^{\dagger}&0&\hat{\epsilon}_{R}\end{pmatrix}$ (37)

Here, $\mathbb{V_{\alpha}}=\\{V_{k\alpha}\\}$ is a vector representation of
the system-bath coupling, and
$\hat{\epsilon}_{\alpha}=\\{\epsilon_{k\alpha}\\}$ is a matrix representation
of the energies in the baths.

We then proceed to build the steady state density matrix $\hat{\sigma}_{ss}$.
We first define non-interacting density operator

$\displaystyle\hat{\sigma}_{0}=\begin{pmatrix}f(\epsilon_{d}-\frac{1}{2}\mu_{L}-\frac{1}{2}\mu_{R})&0&0\\\
0&f(\hat{\epsilon}_{L}-\mu_{L})&0\\\
0&0&f(\hat{\epsilon}_{R}-\mu_{R})\end{pmatrix}$ (38)

The interacting steady-state density operator is obtained by transforming
$\hat{\sigma}_{0}$ in the basis where $\hat{h}$ is diagonal and zero out all
non-diagonal terms of $U^{\dagger}\hat{\sigma}_{0}U$. Here $U^{\dagger}$
diagonalizes the $\hat{h}$. Such that $\hat{\sigma}_{ss}$ and $\hat{h}$ are
both diagonal in such a basis.

$\displaystyle\hat{h}=\sum_{m}\epsilon_{m}|m\rangle\langle m|$ (39)
$\displaystyle\hat{\sigma}_{ss}=\sum_{m}\sigma_{m}|m\rangle\langle m|$ (40)

Such a scheme is similar to the scattering formulation of the density
operator. Ness (2017)

With the eigenvalues and eigenbasis of $\hat{h}$ and $\hat{\sigma}_{ss}$, the
zeroth order population is easily obtained:

$N^{(0)}=\sum_{m}|\langle d|m\rangle|^{2}\sigma_{m}$ (41)

The first order correction requires $\hat{\sigma}^{(1)}$ (Eq. 23), which can
be evaluated in the eigenbasis of $\hat{h}$:

$\displaystyle\hat{\sigma}^{(1)}=\pi\hbar\dot{\epsilon}_{d}\sum_{mn}|m\rangle\langle
m|\partial_{\epsilon_{d}}\sigma_{ss}|n\rangle\delta(\epsilon_{m}-\epsilon_{n})\langle
n|$ (42)

Using Hellman-Feynman theorem, we can rewrite the above as (see the appendix
in Ref. (35))

$\displaystyle\hat{\sigma}^{(1)}=\pi\hbar\dot{\epsilon}_{d}\sum_{mn}|m\rangle\left<m|d\right>\frac{\sigma_{m}-\sigma_{n}}{\epsilon_{m}-\epsilon_{n}}\delta(\epsilon_{m}-\epsilon_{n})\left<d|n\right>\langle
n|$ (43)

Such that the first order population is

$N^{(1)}=\pi\hbar\dot{\epsilon}_{d}\sum_{mn}|\left<m|d\right>|^{2}\frac{\sigma_{m}-\sigma_{n}}{\epsilon_{m}-\epsilon_{n}}\delta(\epsilon_{m}-\epsilon_{n})|\left<d|n\right>|^{2}$
(44)

In numerical simulations, we replace the delta function by a Gaussian with a
broadening parameter $\eta$:

$\displaystyle\pi\delta(\epsilon_{m}-\epsilon_{n})$ $\displaystyle\rightarrow$
$\displaystyle\sqrt{\frac{\pi}{2t^{2}}}\exp{\left[\frac{(\epsilon_{m}-\epsilon_{n})^{2}}{2\eta^{2}}\right]}$
(45)

We also take care of the denominator $\epsilon_{m}-\epsilon_{n}$ in Eq. 44 as

$\displaystyle\frac{\sigma_{m}-\sigma_{n}}{\epsilon_{m}-\epsilon_{n}}$
$\displaystyle\rightarrow$
$\displaystyle\Re\frac{\sigma_{m}-\sigma_{n}}{\epsilon_{m}-\epsilon_{n}+i\eta}$
(46)

In practice, the broadening parameter $\eta$ is set to be in the order of the
energy spacing to better converge the results.

The correlation function in Eq. 26 can be recast in the eigenbasis (see the
supplementary material in Ref. (40)):

$D=\pi\hbar\dot{\epsilon}_{d}\sum_{mn}|\left<m|d\right>|^{2}\sigma_{m}(1-\sigma_{n})\delta(\epsilon_{m}-\epsilon_{n})|\left<d|n\right>|^{2}$
(47)

which is readily to be calculated. Finally, $\delta\dot{S}$ can be written in
the eigenbasis as

$\delta\dot{S}=\pi\hbar\dot{\epsilon}_{d}\sum_{mn}\langle
m|\hat{\sigma}^{(1)}|n\rangle\langle
n|\partial_{\epsilon_{d}}(\ln\hat{\sigma}_{ss}-\ln(1-\hat{\sigma}_{ss}))|m\rangle$
(48)

Here, $\langle m|\partial_{\epsilon_{d}}\ln\hat{\sigma}_{ss}|n\rangle$ is done
numerically with finite difference. Note that $|m\rangle$ is the adiabatic
basis. To calculate the matrix element properly, we rotate the matrix
$\ln\hat{\sigma}_{ss}$ back to the diabatic basis. After done derivative, we
rotate the matrix back to the adiabatic basis.

##  III Results and discussion

###  III.1 Weak coupling limit: Master Equations

In this subsection, we consider the the weak coupling limit. In such a limit,
we can trace out the bath degrees of freedom and derive an equation of motion
for the dot density $\sigma_{d}$ only Dou  _et al._ (2015); Dou and Subotnik
(2016):

${\partial}_{t}\sigma_{d}=-\dot{\epsilon_{d}}{\partial}_{\epsilon_{d}}\sigma_{d}-\frac{\Gamma}{\hbar}\left[\sigma_{d}-\bar{f}(\epsilon_{d}\right].$
(49)

Again, we can expand the dot density into a series of powers in terms of the
driving speeds,
$\sigma_{d}=\sigma_{d}^{(0)}+\sigma_{d}^{(1)}+\sigma_{d}^{(2)}+\cdots$.
Similar to the procedures in Sec. II.2, we obtain the $n$th order population
as

$\displaystyle N^{(0)}=\bar{f}(\epsilon_{d}),$ (50) $\displaystyle
N^{(1)}=-\frac{\hbar\dot{\epsilon}_{d}}{\Gamma}{\partial}_{\epsilon_{d}}\bar{f}(\epsilon_{d}),$
(51) $\displaystyle
N^{(2)}=\frac{\hbar^{2}\dot{\epsilon}_{d}^{2}}{\Gamma^{2}}{\partial}_{\epsilon_{d}}^{2}\bar{f}(\epsilon_{d}).$
(52)

Similarly, we can also obtain the friction and fluctuation,

$\displaystyle\gamma^{(1)}=-\frac{\hbar}{\Gamma}{\partial}_{\epsilon_{d}}\bar{f}(\epsilon_{d}),$
(53) $\displaystyle D=\frac{\hbar
k_{B}T}{\Gamma}\bar{f}(\epsilon_{d})[1-\bar{f}(\epsilon_{d})]$ (54)

To get the entropy production rate, we define the entropy for the dot as
$S=-\sigma_{d}\ln\sigma_{d}-(1-\sigma_{d})\ln(1-\sigma_{d})$. Similar to the
procedures in II.2, we arrive at the entropy production rate as

$\displaystyle\frac{\hbar
T\delta\dot{S}}{\dot{\epsilon}_{d}}=[\gamma^{(1)}]^{2}/(\beta D).$ (55)

Note that the above equation is only true for the driven resonant level model
in the weak coupling limit. We expect that our numerical results can reduce to
the master equation results in the limit of $\Gamma<kT$. What is missing in
the master equation is mainly the level broadening effects. Hence, we refer to
the master equation results as unbroadened results.

###  III.2 Results

We now present our results in this subsection. In our numerical calculations,
number of the levels is set to be $1000$ for both leads, and level spacing is
set to be uniform. The bandwidth is $100\Gamma$. We have also set
$\hbar\dot{\epsilon}_{d}=0.1\Gamma^{2}$ and the broadening term $\eta$ is
equal to the energy spacing.

####  III.2.1 Zeroth and First Order Results

We first look at the zeroth and first order corrections to the thermodynamic
quantities. In Fig. 2, we plot the dot population $N^{(0)}$ and $N^{(1)}$ as a
function of $\epsilon_{d}$ at equilibrium where $\mu_{L}=\mu_{R}=-2\Gamma$. In
such a case, analytical results are available. Notice the good agreement
between the numerical results and the analytical results for both the zeroth
order and first order populations. In the numerical calculations, the main
source of the error is the discretization of finite bandwidth. This agreement
verify the reliability of the numerical method. The results from master
equation having a sharper feature as compared to the other results, suggesting
the missing of the level broadening effects. Note also that the first order
correction to the population $N^{(1)}$ shows a peak at
$\epsilon_{d}=\mu_{L}=\mu_{R}$.

Figure 2:  $N^{(0)}$ and $N^{(1)}$ as a function of $\epsilon_{d}$ at
equilibrium. $\mu_{L}=\mu_{R}=-2\Gamma$, $\Gamma=k_{B}T$.

In Fig. 3, we plot friction, fluctuation, and entropy production rate as a
function of $\epsilon_{d}$ at equilibrium. In Eq. 33, we have shown that
friction, fluctuation, and entropy production rate are proportional to each
other analytically without non-equilibrium condition. Indeed, our numerical
calculations further verify this statement. Such an agreement will no long
exist under non-equilibrium condition (as shown below).

Figure 3: Friction, fluctuation, and entropy production rate as a function of
$\epsilon_{d}$ at equilibrium. $\mu_{L}=\mu_{R}=-2\Gamma$, $\Gamma=k_{B}T$.

Now, let us look at the non-equilibrium case. In Fig. 4, we plot $N^{(0)}$ and
$N^{(1)}$ as a function of $\epsilon_{d}$ under non-equilibrium condition,
where $\mu_{L}=-\mu_{R}=-2\Gamma$. Again, we see good agreement between the
analytical solution and the numerical calculation. The first order population
$N^{(1)}$ shows two peaks at $e_{d}=\mu_{L}$ and $e_{d}=\mu_{R}$ respectively,
resulting a dip at $\epsilon_{d}=\frac{1}{2}(\mu_{L}+\mu_{R})$. Hence, the
peaks are Fermi resonance in nature. Again, results from quantum master
equation show shaper dips and peaks due to lack of broadening effects.

Figure 4: $N^{(0)}$ and $N^{(1)}$ as a function of $\epsilon_{d}$ under
nonequilibrium condition. $\mu_{L}=-\mu_{R}=-2\Gamma$, $\Gamma=k_{B}T$.

In Fig. 5, we plot friction, fluctuation, and entropy production rate as a
function of $\epsilon_{d}$ under non-equilibrium condition. Now we see that,
unlike the equilibrium case, friction, fluctuation, and entropy production
rate do not agree with each other. In fact, in general all three quantities
can be recast into correlation functions: Dou  _et al._ (2020)

$\displaystyle\delta\dot{S}$ $\displaystyle=$
$\displaystyle{k_{B}}{\bar{\beta}^{2}}\sum_{\alpha\nu}\dot{R}_{\alpha}\dot{R}_{\nu}\int_{0}^{\infty}\langle\delta\hat{\mathcal{F}}_{\alpha}(t)\delta\hat{\mathcal{F}}_{\nu}\rangle_{K}dt$
(56) $\displaystyle\gamma_{\alpha\nu}$ $\displaystyle=$
$\displaystyle\bar{\beta}\int_{0}^{\infty}\langle\delta\hat{{F}}_{\alpha}(t)\delta\hat{\mathcal{F}}_{\nu}\rangle_{K}dt$
(57) $\displaystyle D_{\alpha\nu}$ $\displaystyle=$
$\displaystyle\int_{0}^{\infty}\langle\delta\hat{{F}}_{\alpha}(t)\delta\hat{{F}}_{\nu}\rangle_{S}dt.$
(58)

Here $\langle\cdots\rangle_{K}$ denotes Kubo transformed average of the
correlation functions. $\bar{\beta}$ is the reduced temperature Ness (2017).
$\delta\hat{\mathcal{F}}_{\alpha}$ can be seen as the generalized random
force:

$\displaystyle\delta\hat{\mathcal{F}}_{\alpha}=-\frac{1}{\bar{\beta}}\partial_{\alpha}\ln\hat{\rho}_{ss}=\partial_{\alpha}\hat{H}-\partial_{\alpha}\hat{Y}-Tr(\hat{\rho}_{ss}(\partial_{\alpha}\hat{H}-\partial_{\alpha}\hat{Y}))$
(59)

Here $\hat{Y}$ is an operator that accounts for particle transport. Hershfield
(1993) $\delta\hat{{F}}_{\alpha}$ is the random force defined in Eq. 27.
Hence, friction, fluctuation, and entropy production rate are just correction
functions of the two random force operators. Since $\delta\dot{S}$ and
$D_{\alpha\nu}$ are the correlation functions of the same random force, such
that $\delta\dot{S}$ and $D_{\alpha\nu}$ are positive definite, whereas in
general $\gamma_{\alpha\nu}$ is not positive definite under non-equilibrium
condition. Here, for the resonant level model, we see that
$\gamma_{\alpha\nu}$ remains positive even out of equilibrium. However, For
more complicated model, this is not true. Bode  _et al._ (2012)

Figure 5: Friction, fluctuation, and entropy production rate as a function of
$\epsilon_{d}$ under nonequilibrium condition. $\mu_{L}=-\mu_{R}=-2\Gamma$,
$\Gamma=k_{B}T$.

####  III.2.2 Second Order Results

We now look at the second order correction to the population (or second order
friction). Here, HQME is used to simulate the exact dynamics for the driven
resonant level model, where we can calculate the exact population as a
function of time. We then estimate the second order friction by subtracting
the zeroth and first order correction:

$\gamma^{(2)}=\left(N-N^{(0)}-\dot{\epsilon}_{d}\gamma^{(1)}\right)/\dot{\epsilon}_{d}^{2},$
(60)

In the above equation, we have ignored the third order and higher correction.
In Fig. 6, we plot $\gamma^{(2)}$ as a function of $\epsilon_{d}$ from HQME
with different driving speeds. At equilibrium, HQME results reproduce
analytical results in Eq. 36. Out of equilibrium, we do not have analytical
results for $\gamma^{(2)}$. However, we simply replace $f$ by $\bar{f}$ in Eq.
36 as our trial solution. Such a solution agrees with the HQME results as well
as results from master equations. Indeed, future work must verify this trial
solution for the second order correction out of equilibrium.

Figure 6:  $\gamma^{(2)}$ as a function of $\epsilon_{d}$ under equilibrium
($\mu_{L}=\mu_{R}=0$) and nonequilibrium condition
($\mu_{L}=-\mu_{R}=-2\Gamma$). $\Gamma=k_{B}T$.

##  IV Conclusion

Using the expansion of the density operator in the power of driving speed, we
have identified the non-adiabatic correction to the quantum thermodynamic
quantities in the strongly coupled regimes. With numerical tools, we have
calculated non-adiabatic corrections for the resonant level model. We have
verified our numerical results against analytical results for the zeroth and
first order correction to the population in and out of equilibrium. We then
have calculated the non-adiabatic entropy production out of equilibrium. We
show that, at equilibrium, the friction, random force, and entropy production
agree with each other. Out of equilibrium, friction, random force and entropy
production are all correlation functions between different (generalized)
random forces. Our results agree with the one from the master equation in the
limit of weak couplings. Future work must study the non-adiabatic entropy
production and higher order corrections go beyond simple resonant level model.

###### Acknowledgements.

We thank Jakob Bätge for providing us data from HQME calculations. We also
acknowledge the startup funding from Westlake University.

## Data availability

The data that support the findings of this study are available upon reasonable
request from the authors.

##  Appendix A Evaluation of
$\left<d\right|\delta(\epsilon-\hat{h})\left|d\right>$

In this appendix, we draw certain identities and then evaluate
$\left<d\right|\delta(\epsilon-\hat{h})\left|d\right>$. To start with, we have
the following identity:

$\delta(\epsilon-\hat{h})(\epsilon-\hat{h})=(\epsilon-\hat{h})\delta(\epsilon-\hat{h})=0.$
(61)

Taking ${\partial}_{\epsilon_{d}}$ on both sides, we arrive at

$-{\partial}_{\epsilon_{d}}\hat{h}\delta(\epsilon-\hat{h})+(\epsilon-\hat{h}){\partial}_{\epsilon_{d}}\delta(\epsilon-\hat{h})=-\delta(\epsilon-\hat{h}){\partial}_{\epsilon_{d}}\hat{h}+{\partial}_{\epsilon_{d}}\delta(\epsilon-\hat{h})(\epsilon-\hat{h})=0.$
(62)

Hence,

${\partial}_{\epsilon_{d}}\hat{h}=(\epsilon-\hat{h})^{-1}{\partial}_{\epsilon_{d}}\hat{h}\delta(\epsilon-\hat{h})=\delta(\epsilon-\hat{h}){\partial}_{\epsilon_{d}}\hat{h}(\epsilon-\hat{h})^{-1}.$
(63)

Similarly, we can show that

${\partial}_{\epsilon}\delta(\epsilon-\hat{h})=-(\epsilon-\hat{h})^{-1}\delta(\epsilon-\hat{h})=-\delta(\epsilon-\hat{h})(\epsilon-\hat{h})^{-1}.$
(64)

The above equations will be useful for higher order corrections.

In the diagonal basis, we have

$\delta(\epsilon-\hat{h})=\sum_{k\alpha}\left|\psi_{k\alpha}\right>\delta(\epsilon-\epsilon_{k\alpha})\left<\psi_{k\alpha}\right|,$
(65)

such that

$\displaystyle\left<d\right|\delta(\epsilon-\hat{h})\left|d\right>=\sum_{k\alpha}\left|\left<d|\psi_{k\alpha}\right>\right|^{2}\delta(\epsilon-\epsilon_{k\alpha})$
$\displaystyle=\sum_{k\alpha}\frac{|V_{k\alpha}|^{2}\delta(\epsilon-\epsilon_{k\alpha})}{(\epsilon-\epsilon_{d})^{2}+(\Gamma/2)^{2}}=\frac{1}{2\pi}\sum_{\alpha}\frac{\Gamma_{\alpha}}{(\epsilon-\epsilon_{d})^{2}+(\Gamma/2)^{2}}=\frac{A(\epsilon)}{2\pi}.$
(66)

##  Appendix B Deviation of Eq. 35

At equilibrium, the one-body steady state density is

$\hat{\sigma}_{ss}=f(\hat{h})=\int
f(\epsilon)\delta(\epsilon-\hat{h})d\epsilon,$ (67)

Use the identity in Appendix A, we can show that

$\hat{\sigma}^{(1)}=-\pi\hbar\dot{\epsilon_{d}}\int
d\epsilon\delta(\epsilon-\hat{h}){\partial}_{\epsilon_{d}}\hat{h}\delta(\epsilon-\hat{h}){\partial}_{\epsilon}f(\epsilon),$
(68)

Taking the derivative of $\hat{\sigma}^{(1)}$ with respect to $\epsilon_{d}$,
we have

$\displaystyle{\partial}_{\epsilon_{d}}\hat{\sigma}^{(1)}$ $\displaystyle=$
$\displaystyle-\pi\hbar\dot{\epsilon_{d}}\int
d\epsilon{\partial}_{\epsilon_{d}}\delta(\epsilon-\hat{h}){\partial}_{\epsilon_{d}}\hat{h}\delta(\epsilon-\hat{h}){\partial}_{\epsilon}f(\epsilon)$
(69) $\displaystyle-\pi\hbar\dot{\epsilon_{d}}\int
d\epsilon\delta(\epsilon-\hat{h}){\partial}_{\epsilon_{d}}\hat{h}{\partial}_{\epsilon_{d}}\delta(\epsilon-\hat{h}){\partial}_{\epsilon}f(\epsilon)$
$\displaystyle=$ $\displaystyle-\pi\hbar\dot{\epsilon_{d}}\int
d\epsilon(\epsilon-\hat{h})^{-1}{\partial}_{\epsilon_{d}}\hat{h}\delta(\epsilon-\hat{h}){\partial}_{\epsilon_{d}}\hat{h}\delta(\epsilon-\hat{h}){\partial}_{\epsilon}f(\epsilon)$
$\displaystyle-\pi\hbar\dot{\epsilon_{d}}\int
d\epsilon\delta(\epsilon-\hat{h}){\partial}_{\epsilon_{d}}\hat{h}\delta(\epsilon-\hat{h}){\partial}_{\epsilon_{d}}\hat{h}(\epsilon-\hat{h})^{-1}{\partial}_{\epsilon}f(\epsilon)$

We have used the identities shown in Appendix A. To evaluate
$\hat{\sigma}^{(2)}$, we first calculate the following term,

$\displaystyle\int_{0}^{\infty}e^{i\hat{h}t/\hbar}\int
d\epsilon(\epsilon-\hat{h})^{-1}{\partial}_{\epsilon_{d}}\hat{h}\delta(\epsilon-\hat{h}){\partial}_{\epsilon_{d}}\hat{h}\delta(\epsilon-\hat{h}){\partial}_{\epsilon}f(\epsilon)e^{-i\hat{h}t/\hbar}$
(70) $\displaystyle=$
$\displaystyle\sum_{mnp}\int_{0}^{\infty}dte^{i(\epsilon_{m}-\epsilon_{n})t/\hbar}\int
d\epsilon\left|m\right>(\epsilon-\epsilon_{m})^{-1}\left<m\right|{\partial}_{\epsilon_{d}}\hat{h}\left|p\right>\delta(\epsilon-\epsilon_{p})\left<p\right|{\partial}_{\epsilon_{d}}\hat{h}\left|n\right>\delta(\epsilon-\epsilon_{n}){\partial}_{\epsilon}f(\epsilon)\left<n\right|$
$\displaystyle=$
$\displaystyle\pi\hbar\sum_{mnp}\delta(\epsilon_{n}-\epsilon_{m})\left|m\right>(\epsilon_{n}-\epsilon_{m})^{-1}\left<m\right|{\partial}_{\epsilon_{d}}\hat{h}\left|p\right>\delta(\epsilon_{n}-\epsilon_{p})\left<p\right|{\partial}_{\epsilon_{d}}\hat{h}\left|n\right>{\partial}_{\epsilon}f(\epsilon_{n})\left<n\right|$
$\displaystyle=$ $\displaystyle\pi\hbar\sum_{mnp}\int
d\epsilon\left|m\right>(\epsilon-\epsilon_{m})^{-1}\delta(\epsilon-\epsilon_{m})\left<m\right|{\partial}_{\epsilon_{d}}\hat{h}\left|p\right>\delta(\epsilon-\epsilon_{p})\left<p\right|{\partial}_{\epsilon_{d}}\hat{h}\left|n\right>\delta(\epsilon-\epsilon_{n}){\partial}_{\epsilon}f(\epsilon)\left<n\right|$
$\displaystyle=$ $\displaystyle-\pi\hbar\int
d\epsilon{\partial}_{\epsilon}\delta(\epsilon-\hat{h}){\partial}_{\epsilon_{d}}\hat{h}\delta(\epsilon-\hat{h}){\partial}_{\epsilon_{d}}\hat{h}\delta(\epsilon-\hat{h}){\partial}_{\epsilon}f(\epsilon).$

Similarly,

$\displaystyle\int_{0}^{\infty}e^{i\hat{h}t/\hbar}\int
d\epsilon\delta(\epsilon-\hat{h}){\partial}_{\epsilon_{d}}\hat{h}\delta(\epsilon-\hat{h}){\partial}_{\epsilon_{d}}\hat{h}(\epsilon-\hat{h})^{-1}{\partial}_{\epsilon}f(\epsilon)e^{-i\hat{h}t/\hbar}$
(71) $\displaystyle=$ $\displaystyle-\pi\hbar\int
d\epsilon\delta(\epsilon-\hat{h}){\partial}_{\epsilon_{d}}\hat{h}\delta(\epsilon-\hat{h}){\partial}_{\epsilon_{d}}\hat{h}{\partial}_{\epsilon}\delta(\epsilon-\hat{h}){\partial}_{\epsilon}f(\epsilon).$

Therefore,

$\displaystyle\hat{\sigma}^{(2)}=$
$\displaystyle-\dot{\epsilon_{d}}\int_{0}^{\infty}e^{i\hat{h}t/\hbar}{\partial}_{\epsilon_{d}}\hat{\sigma}^{(1)}e^{-i\hat{h}t/\hbar}$
(72) $\displaystyle=$
$\displaystyle-\pi^{2}\hbar^{2}\dot{\epsilon_{d}}^{2}\int
d\epsilon{\partial}_{\epsilon}\delta(\epsilon-\hat{h}){\partial}_{\epsilon_{d}}\hat{h}\delta(\epsilon-\hat{h}){\partial}_{\epsilon_{d}}\hat{h}\delta(\epsilon-\hat{h}){\partial}_{\epsilon}f(\epsilon)$
$\displaystyle-\pi^{2}\hbar^{2}\dot{\epsilon_{d}}^{2}\int
d\epsilon\delta(\epsilon-\hat{h}){\partial}_{\epsilon_{d}}\hat{h}\delta(\epsilon-\hat{h}){\partial}_{\epsilon_{d}}\hat{h}{\partial}_{\epsilon}\delta(\epsilon-\hat{h}){\partial}_{\epsilon}f(\epsilon).$

We then proceed to evaluate $N^{(2)}$:

$\displaystyle N^{(2)}=$
$\displaystyle\left<d\right|\sigma^{(2)}\left|d\right>=-\frac{\hbar^{2}}{4\pi}\dot{\epsilon}_{d}^{2}\int
d\epsilon{\partial}_{\epsilon}A(\epsilon)A^{2}(\epsilon){\partial}_{\epsilon}{f}(\epsilon)$
(73) $\displaystyle=$
$\displaystyle-\frac{\hbar^{2}\dot{\epsilon}_{d}^{2}}{12\pi}\int
d\epsilon{\partial}_{\epsilon}A^{3}(\epsilon){\partial}_{\epsilon}{f}(\epsilon)$
$\displaystyle=$
$\displaystyle-\frac{\hbar^{2}\dot{\epsilon}_{d}^{2}}{12\pi}\left[A^{3}(\epsilon){\partial}_{\epsilon}{f}(\epsilon)\right]_{0}^{\infty}+\frac{\hbar^{2}\dot{\epsilon}_{d}^{2}}{12\pi}\int
d\epsilon A^{3}(\epsilon){\partial}_{\epsilon}^{2}{f}(\epsilon)$
$\displaystyle=$
$\displaystyle\frac{\hbar^{2}\dot{\epsilon}_{d}^{2}}{12\pi}\int d\epsilon
A^{3}(\epsilon){\partial}_{\epsilon}^{2}{f}(\epsilon),$

which gives us the results shown in Eq. 35. We have used the fact that
$\partial_{\epsilon_{d}}\hat{h}=|d\rangle\langle d|$.

## References

  * Gemmer  _et al._ (2009) J. Gemmer, M. Michel, and G. Mahler, “Quantum thermodynamics: Emergence of thermodynamic behavior within composite quantum systems, volume 784 of lecture notes in physics,” (2009).
  * Kosloff (2013) R. Kosloff, Entropy 15, 2100 (2013).
  * Esposito  _et al._ (2015a) M. Esposito, M. A. Ochoa, and M. Galperin, Phys. Rev. Lett. 114, 080602 (2015a).
  * Vinjanampathy and Anders (2016) S. Vinjanampathy and J. Anders, Contemporary Physics 57, 545 (2016).
  * Anders and Esposito (2017) J. Anders and M. Esposito, New Journal of Physics 19, 010201 (2017).
  * Carrega  _et al._ (2016) M. Carrega, P. Solinas, M. Sassetti, and U. Weiss, Phys. Rev. Lett. 116, 240403 (2016).
  * Alicki and Kosloff (2018) R. Alicki and R. Kosloff, arXiv preprint arXiv:1801.08314 (2018).
  * Benenti  _et al._ (2017) G. Benenti, G. Casati, K. Saito, and R. S. Whitney, Physics Reports 694, 1 (2017).
  * Campisi  _et al._ (2011) M. Campisi, P. Hänggi, and P. Talkner, Rev. Mod. Phys. 83, 771 (2011).
  * Campisi  _et al._ (2009) M. Campisi, P. Talkner, and P. Hänggi, Phys. Rev. Lett. 102, 210401 (2009).
  * Liu  _et al._ (2021) J. Liu, K. A. Jung, and D. Segal, Physical review letters 127, 200602 (2021).
  * Poot and van der Zant (2012) M. Poot and H. S. van der Zant, Phys. Rep. 511, 273 (2012).
  * Pekola (2015) J. P. Pekola, Nature Physics 11, 118 (2015).
  * Roßnagel  _et al._ (2016) J. Roßnagel, S. T. Dawkins, K. N. Tolazzi, O. Abah, E. Lutz, F. Schmidt-Kaler, and K. Singer, Science 352, 325 (2016).
  * Argun  _et al._ (2017) A. Argun, J. Soni, L. Dabelow, S. Bo, G. Pesce, R. Eichhorn, and G. Volpe, Phys. Rev. E 96, 052106 (2017).
  * Klatzow  _et al._ (2019) J. Klatzow, J. N. Becker, P. M. Ledingham, C. Weinzetl, K. T. Kaczmarek, D. J. Saunders, J. Nunn, I. A. Walmsley, R. Uzdin, and E. Poem, Phys. Rev. Lett. 122, 110601 (2019).
  * Esposito  _et al._ (2015b) M. Esposito, M. A. Ochoa, and M. Galperin, Phys. Rev. B 92, 235440 (2015b).
  * Millen and Xuereb (2016) J. Millen and A. Xuereb, New J. Phys. 18, 011002 (2016).
  * Talkner  _et al._ (2007) P. Talkner, E. Lutz, and P. Hänggi, Phys. Rev. E 75, 050102 (2007).
  * Brandao  _et al._ (2015) F. Brandao, M. Horodecki, N. Ng, J. Oppenheim, and S. Wehner, Proc. Natl. Acad. Sci. U.S.A. 112, 3275 (2015).
  * Jarzynski (2011) C. Jarzynski, Annu. Rev. Condens. Matter Phys. 2, 329 (2011).
  * Solinas  _et al._ (2013) P. Solinas, D. V. Averin, and J. P. Pekola, Physical Review B 87, 060508 (2013).
  * Schmidt  _et al._ (2015) R. Schmidt, M. F. Carusela, J. P. Pekola, S. Suomela, and J. Ankerhold, Physical Review B 91, 224303 (2015).
  * Gogolin and Eisert (2016) C. Gogolin and J. Eisert, Reports on Progress in Physics 79, 056001 (2016).
  * Subaşı  _et al._ (2012) Y. Subaşı, C. Fleming, J. Taylor, and B. L. Hu, Physical Review E 86, 061132 (2012).
  * Ness (2017) H. Ness, Entropy 19, 158 (2017).
  * Strasberg  _et al._ (2017) P. Strasberg, G. Schaller, T. Brandes, and M. Esposito, Phys. Rev. X 7, 021003 (2017).
  * Bruch  _et al._ (2016) A. Bruch, M. Thomas, S. Viola Kusminskiy, F. von Oppen, and A. Nitzan, Phys. Rev. B 93, 115318 (2016).
  * Haughian  _et al._ (2018) P. Haughian, M. Esposito, and T. L. Schmidt, Physical Review B 97, 085435 (2018).
  * Ochoa  _et al._ (2016) M. A. Ochoa, A. Bruch, and A. Nitzan, Phys. Rev. B 94, 035420 (2016).
  * Bruch  _et al._ (2018) A. Bruch, C. Lewenkopf, and F. von Oppen, Physical review letters 120, 107701 (2018).
  * Semenov and Nitzan (2020) A. Semenov and A. Nitzan, The Journal of Chemical Physics 152, 244126 (2020).
  * Bergmann and Galperin (2021) N. Bergmann and M. Galperin, The European Physical Journal Special Topics 230, 859 (2021).
  * Strasberg and Winter (2021) P. Strasberg and A. Winter, PRX Quantum 2, 030202 (2021).
  * Dou  _et al._ (2018) W. Dou, M. A. Ochoa, A. Nitzan, and J. E. Subotnik, Phys. Rev. B 98, 134306 (2018).
  * Dou  _et al._ (2020) W. Dou, J. Bätge, A. Levy, and M. Thoss, Physical Review B 101, 184304 (2020).
  * Schinabeck and Thoss (2020) C. Schinabeck and M. Thoss, Phys. Rev. B 101, 075422 (2020).
  * Jin  _et al._ (2008) J. Jin, X. Zheng, and Y. Yan, The Journal of chemical physics 128, 234703 (2008).
  * Tanimura (2020) Y. Tanimura, The Journal of chemical physics 153, 020901 (2020).
  * Dou  _et al._ (2017) W. Dou, G. Miao, and J. E. Subotnik, Phys. Rev. Lett. 119, 046001 (2017).
  * Oz  _et al._ (2019) A. Oz, O. Hod, and A. Nitzan, Journal of Chemical Theory and Computation (2019).
  * Dou  _et al._ (2015) W. Dou, A. Nitzan, and J. E. Subotnik, J. Chem. Phys. 142, 084110 (2015).
  * Dou and Subotnik (2016) W. Dou and J. E. Subotnik, J. Chem. Phys. 145, 054102 (2016).
  * Hershfield (1993) S. Hershfield, Phys. Rev. Lett. 70, 2134 (1993).
  * Bode  _et al._ (2012) N. Bode, S. V. Kusminskiy, R. Egger, and F. von Oppen, Beilstein J. Nanotechnol 3, 144 (2012).