# Quantum thermodynamics with fast driving and strong coupling via the
mesoscopic leads approach

Artur M. Lacerda  Department of Physics, Trinity College Dublin, Dublin 2,
Ireland  Archak Purkayastha<EMAIL_ADDRESS>Department of Physics,
Trinity College Dublin, Dublin 2, Ireland  Center for Complex Quantum Systems,
Department of Physics and Astronomy, Aarhus University, Ny Munkegade 120,
DK-8000 Aarhus C, Denmark.  Michael Kewming  Department of Physics, Trinity
College Dublin, Dublin 2, Ireland  Gabriel T. Landi  Department of Physics,
Trinity College Dublin, Dublin 2, Ireland  Instituto de Física da Universidade
de São Paulo, 05314-970 São Paulo, Brazil  John Goold  Department of Physics,
Trinity College Dublin, Dublin 2, Ireland

###### Abstract

Understanding the thermodynamics of driven quantum systems strongly coupled to
thermal baths is a central focus of quantum thermodynamics and mesoscopic
physics. A variety of different methodological approaches exist in the
literature, all with their own advantages and disadvantages. The mesoscopic
leads approach was recently generalised to steady state thermal machines and
has the ability to replicate Landauer Büttiker theory in the non-interacting
limit. In this approach a set of discretised lead modes, each locally damped,
provide a markovian embedding for the baths. In this work we further
generalise this approach to incorporate an arbitrary time dependence in the
system Hamiltonian. Following a careful discussion of the calculation of
thermodynamic quantities we illustrate the power of our approach by studying
several driven mesoscopic examples coupled to finite temperature fermionic
baths, replicating known results in various limits. In the case of a driven
non interacting quantum dot we show how fast driving can be used to induce
heat rectification.

##  I Introduction

Thermodynamics is a pillar of physical, chemical and biological sciences
primarily due to the unique synergy it offers between the fundamental and
pragmatic.It tells us that heat and work are just two different forms of
energy (the first law), and it sets constraints on the heat generated by
mechanical operations that we are allowed to perform (the second law) Callen
(1998). Importantly, the notion of thermal equilibrium is central for defining
the state functions of conventional thermodynamics. However, when we move
towards the microscopic domain, understanding thermodynamics far from
equilibrium becomes increasingly important De Groot and Mazur (2013);
Jarzynski (2011); Seifert (2012, 2019).

In particular, as we go towards the advent of quantum devices across various
platforms like atoms and ions in optical traps Monroe and Kim (2013),
molecular junctions Aradhya and Venkataraman (2013), quantum dots Kagan  _et
al._ (2021), hybrid semiconductor-superconductor circuits Burkard  _et al._
(2020); Clerk  _et al._ (2020), nitrogen-vacancy centers in diamonds Doherty
_et al._ (2013) etc., it has become necessary to understand the energy and the
entropy cost of such applications Arrachea (2022); Auffèves (2021). These
devices are inherently noisy because complete isolation from surrounding
environment is neither possible at such small length scales, nor desired in
many applications. In fact, in some applications, coupling to multiple baths
is desired, for example, to induce a voltage bias (for example, Liu  _et al._
(2015); Stehlik  _et al._ (2016); Wen  _et al._ (2020)). Moreover, the working
principle of these devices often require control via external driving,
modelled by a time dependent Hamiltonian Koch  _et al._ (2019); James (2021);
Bhandari  _et al._ (2020). To describe and understand the energy and the
entropy cost of such devices, one needs (a) a consistent extension of the
standard notion thermodynamics to externally driven microscopic quantum
systems coupled to multiple macroscopic baths, (b) a prescription to calculate
the required thermodynamic quantities. In this paper, we show that the so-
called mesoscopic leads approach (also known by other names, such as
pseudomodes approach and driven Liouville von-Neumann equation) Imamoglu
(1994); Garraway (1997a, b); Sánchez  _et al._ (2006); Subotnik  _et al._
(2009); Dzhioev and Kosov (2011); Ajisaka  _et al._ (2012); Ajisaka and Barra
(2013); Arrigoni  _et al._ (2013); Dorda  _et al._ (2014); Chen  _et al._
(2014); Zelovich  _et al._ (2014); Dorda  _et al._ (2015); Hod  _et al._
(2016); Gruss  _et al._ (2016); Schwarz  _et al._ (2016); Dorda  _et al._
(2017); Elenewski  _et al._ (2017); Gruss  _et al._ (2017); Zelovich  _et al._
(2017); Tamascelli  _et al._ (2018); Lemmer  _et al._ (2018); Oz  _et al._
(2019); Chen  _et al._ (2019); Zwolak (2020); Wójtowicz  _et al._ (2020);
Chiang and Hsu (2020); Brenes  _et al._ (2020); Lotem  _et al._ (2020); Fugger
_et al._ (2020); Wójtowicz  _et al._ (2021); Elenewski  _et al._ (2021)
achieves these goals, providing a numerically efficient formalism applicable
far beyond regimes accessible by other approaches.

Most standard formalisms to describe such open quantum systems rely on weak
system-bath couplings such that a Born-Markov approximation can be made, and
an effective quantum master equation can be derived Breuer and Petruccione
(2007); Rivas and Huelga (2012). Although usually analytically and numerically
the simplest, various drawbacks of such approaches have been pointed out Walls
(1970); Novotný (2002); Wichterich  _et al._ (2007); Prachar and Novotný
(2010); Rivas  _et al._ (2010); Deçordi and Vidiella-Barranco (2017); Levy and
Kosloff (2014); Purkayastha  _et al._ (2016); Eastham  _et al._ (2016); Hofer
_et al._ (2017); González  _et al._ (2017); Mitchison and Plenio (2018);
Cattaneo  _et al._ (2019); Hartmann and Strunz (2020); Konopik and Lutz
(2022); Scali  _et al._ (2021); Tupkary  _et al._ (2022), and it has been
shown that no such approximation can generically accurately describe the long-
time state of the system in presence of multiple baths Tupkary  _et al._
(2022). Moreover, in many applications, like those involving quantum dots and
molecular junctions, experimentally the system-bath couplings are often not
desired to be weak Aradhya and Venkataraman (2013); Burkard  _et al._ (2020);
Clerk  _et al._ (2020); Kagan  _et al._ (2021). The presence of external
driving makes the situation even more complicated, since driving can introduce
further non-Markovian effects (see for example, Schnell  _et al._ (2020,
2021)).

In the non-Markovian regime, the most widely used formalism are those of non-
equilibrium Green’s functions (NEGF) Haug and Jauho (2008) and Schwinger-
Keldysh path integrals Kamenev (2011). Other formalisms like microscopically
derived quantum Langevin equations Cortés  _et al._ (1985); Ford  _et al._
(1988); Benguria and Kac (1981); Dhar and Sriram Shastry (2003); Dhar and Roy
(2006), and the scattering approach Moskalets (2011); Datta  _et al._ (1987);
Arrachea and Moskalets (2006); Lesovik and Sadovskyy (2011); Brandner (2020);
Potanina  _et al._ (2021) are closely related to these approaches. In absence
of many-body interactions (i.e, higher than quadratic terms in fermionic or
bosonic creation and annihilation operators in Hamiltonian), these approaches
can efficiently give exact results for long-time steady state if there is no
external driving. Treating systems with many-body interaction using these
approaches requires existence of a small parameter, such that a diagrammatic
expansion is possible Haug and Jauho (2008); Kamenev (2011). Further, in
presence of an external driving, often even quadratic Hamiltonians cannot be
treated exactly. Assuming that the drive is periodic, one has to often rely on
perturbative expansions either in frequency or in strength of the drive (for
example, Ludovico  _et al._ (2016a); Bruch  _et al._ (2018); Riera-Campeny
_et al._ (2019); Brandner (2020)). In experimental situations, however, one
may not be just limited to such regimes (for example, Stehlik  _et al._
(2016)). Moreover, even in cases where these approaches are efficiently give
the long-time behavior, obtaining the full behavior from short-times to long-
times is usually quite difficult.

The mesoscopic leads approach is an alternate way to obtain dynamics of the
system in strong-coupling non-Markovian regime Imamoglu (1994); Garraway
(1997a, b); Sánchez  _et al._ (2006); Subotnik  _et al._ (2009); Dzhioev and
Kosov (2011); Ajisaka  _et al._ (2012); Ajisaka and Barra (2013); Arrigoni
_et al._ (2013); Dorda  _et al._ (2014); Chen  _et al._ (2014); Zelovich  _et
al._ (2014); Dorda  _et al._ (2015); Hod  _et al._ (2016); Gruss  _et al._
(2016); Schwarz  _et al._ (2016); Dorda  _et al._ (2017); Elenewski  _et al._
(2017); Gruss  _et al._ (2017); Zelovich  _et al._ (2017); Tamascelli  _et
al._ (2018); Lemmer  _et al._ (2018); Oz  _et al._ (2019); Chen  _et al._
(2019); Zwolak (2020); Wójtowicz  _et al._ (2020); Chiang and Hsu (2020);
Brenes  _et al._ (2020); Lotem  _et al._ (2020); Fugger  _et al._ (2020);
Wójtowicz  _et al._ (2021); Elenewski  _et al._ (2021). In this approach, the
macroscopic baths are systematically approximated by a finite number of damped
modes, which we call leads. Initially introduced for bosonic baths Imamoglu
(1994); Garraway (1997a, b), this approach has been subsequently extensively
used to study quantum transport in fermionic set-ups also Sánchez  _et al._
(2006); Subotnik  _et al._ (2009); Dzhioev and Kosov (2011); Ajisaka  _et al._
(2012); Ajisaka and Barra (2013); Arrigoni  _et al._ (2013); Dorda  _et al._
(2014); Chen  _et al._ (2014); Zelovich  _et al._ (2014); Dorda  _et al._
(2015); Hod  _et al._ (2016); Gruss  _et al._ (2016); Schwarz  _et al._
(2016); Dorda  _et al._ (2017); Elenewski  _et al._ (2017); Gruss  _et al._
(2017); Zelovich  _et al._ (2017); Chen  _et al._ (2019); Oz  _et al._ (2019);
Zwolak (2020); Wójtowicz  _et al._ (2020); Chiang and Hsu (2020); Brenes  _et
al._ (2020); Lotem  _et al._ (2020); Fugger  _et al._ (2020); Wójtowicz  _et
al._ (2021); Elenewski  _et al._ (2021). When combined with tensor network
techniques, it has been recently shown that it is possible to completely non-
perturbatively obtain energy and particle currents in non-equilibrium steady
states (NESS) of interacting quantum many-body systems, and thereby the
thermodynamics at NESS, using this approach Brenes  _et al._ (2020). Related
approaches have been developed to describe impurity models at and beyond Kondo
regimes Arrigoni  _et al._ (2013); Dorda  _et al._ (2014, 2015, 2017); Lotem
_et al._ (2020); Fugger  _et al._ (2020). However, most existing
investigations have been limited to time-independent Hamiltonians, with only
few considering time-dependent system Hamiltonians Chen  _et al._ (2014); Oz
_et al._ (2019).

Formalisms to describe thermodynamics of externally driven systems coupled to
baths have been investigated using various approaches and approximations
Esposito  _et al._ (2015a); Brandner and Seifert (2016); Bruch  _et al._
(2018); Brandner (2020); Potanina  _et al._ (2021); Dann and Kosloff (2021),
each with its own merits and drawbacks, and applied to simple set-ups Esposito
_et al._ (2015b); Ludovico  _et al._ (2016a, b); Bruch  _et al._ (2016); Katz
and Kosloff (2016); Haughian  _et al._ (2018); Riera-Campeny  _et al._ (2019);
Oz  _et al._ (2019). But, to our knowledge, thermodynamics with time dependent
drive in the mesoscopic leads approach has only been previously considered in
the case of a driven resonant level model, where the thermodynamics has been
described only in the regime of slow driving Oz  _et al._ (2019). Here, we
extend the mesoscopic leads approach to arbitrary system Hamiltonians in
presence of arbitrary time dependent driving. We start from a microscopic
description which is common to all standard approaches to open quantum systems
(Sec. II.1). Following existing approaches to thermodynamics of quantum
systems Esposito  _et al._ (2010); Reeb and Wolf (2014); Landi and Paternostro
(2021); Strasberg and Winter (2021), we then give a careful and detailed
discussion of the quantities that are required to be calculated to order to
describe the thermodynamics in presence of an arbitrary external drive and
coupling to multiple baths (Sec. II.2). Next, we discuss how the mesoscopic
leads approach can be systematically derived from the microscopic description
(Sec. III.1), and show that all quantities required for description of
thermodynamics can be obtained via this approach (Sec. III.2). Then, we lay
down an elegant formalism for quadratic Hamiltonians, which holds for
arbitrary driving in the system and arbitrary strengths of system-bath
couplings, and can easily describe dynamics and thermodynamics at all times.
This therefore provides a simple way to access cases beyond both NEGF and
standard weak-coupling quantum master equation descriptions (Secs. III.3,
III.4). We apply the formalism to two examples: a driven resonant level (Sec.
IV), and a two-site fermionic system with one site being externally
periodically driven (Sec. V). We use the former for benchmarking our approach.
In the later, motivated by Riera-Campeny  _et al._ (2019) we investigate
rectification of energy currents. In particular, we obtain the entropy cost of
energy rectification as a function of drive frequency. This would be difficult
to obtain using most other approaches. Finally, we conclude by summarizing and
giving the outlook (Sec. VI).

##  II Thermodynamics of driven open quantum systems

###  II.1 Microscopic description of open quantum systems

Let us begin by describing a general system configuration; a microscopic
system of interest coupled to multiple macroscopic baths. The former is
controlled by some external protocol, mathematically described by a time-
dependent Hamiltonian $\hat{H}_{\text{S}}(t)$. The Hamiltonian describing the
coupling between the baths and the central system is thus

$\hat{H}_{\text{tot}}=\hat{H}_{\text{S}}(t)+\sum_{\alpha}\quantity(\hat{H}_{\alpha}+\hat{H}_{\text{S}\alpha}),$
(1)

where $\alpha$ labels the different baths, $\hat{H}_{\alpha}$ is the free
Hamiltonian of the baths, and $\hat{H}_{S\alpha}$ describes the coupling
between the baths and the system. We assume that initial state of the system
is in an uncorrelated product state with the baths. Furthermore, we assume
that each bath begins in a thermal state at inverse temperature
$\beta_{\alpha}$ and chemical potential $\mu_{\alpha}$. Thus the joint state
is

$\hat{\rho}_{\text{tot}}(0)=\hat{\rho}_{\text{S}}(0)\prod_{\alpha}\hat{\rho}_{\alpha}^{\text{th}},\mbox{\quad
where\quad}\hat{\rho}_{\alpha}^{\text{th}}=\frac{e^{-\beta_{\alpha}(\hat{H}_{\alpha}-\mu_{\alpha}\hat{N}_{\alpha})}}{Z_{\alpha}}.$
(2)

where $\hat{N}_{\alpha}$ is the number operator of the $\alpha$th bath and
$Z_{\alpha}=\Tr[e^{-\beta_{\alpha}(\hat{H}_{\alpha}-\mu_{\alpha}\hat{N}_{\alpha})}]$
is the corresponding partition function. Additionally, we assume that the bath
Hamiltonian conserves the number of particles, i.e.,
$\commutator{\hat{H}_{\alpha}}{\hat{N}_{\alpha}}=0$. Starting from the initial
state Eq. (2), the entire system is left to evolve unitarily under the
Hamiltonian $\hat{H}_{\text{tot}}$, according to
$\differential\hat{\rho}_{\text{tot}}(t)/\differential
t=-i[\hat{H}_{\text{tot}}(t),\hat{\rho}_{\text{tot}}]$ (we set $\hbar=1$), so
at time $t$ the total density-matrix is

$\hat{\rho}_{\text{tot}}(t)=\hat{U}(t)\,\hat{\rho}_{\text{tot}}(0)\hat{U}^{\dagger}(t),$
(3)

where $\hat{U}(t)=\mathcal{T}e^{-i\int_{0}^{t}\differential
t^{\prime}\hat{H}_{\text{tot}}(t^{\prime})}$ is the time-evolution operator
and $\mathcal{T}$ is the time-ordering operator. The evolution of the system
of interest, $\hat{\rho}_{\text{S}}$, can be obtained by tracing out the the
bath. We will denote $L_{\text{S}}$ and $L_{\text{B}}$ as the number of
degrees of freedom (e.g., the number of sites in a lattice descriptions) in
the system and in each bath, respectively. Since the system is microscopic,
$L_{\text{S}}$ is finite. Each bath, on the other hand, is by assumption
macroscopic, thus $L_{\text{B}}\to\infty$. However, for mathematical
rigorousness and to avoid unwanted divergences, this limit should be taken
only after the bath degrees of freedom have been traced out. That is,

$\hat{\rho}_{\text{S}}(t)=\lim_{L_{\text{B}}\to\infty}\Tr_{\text{B}}[\hat{\rho}_{\text{tot}}(t)],$
(4)

where $\Tr_{\text{B}}[\ldots]$ refers to trace over all bath degrees of
freedom. The order of the operations is crucial, with the limit is taken after
the trace. In order to avoid carrying cumbersome notation, we define the
expected value of some generic operator $\hat{\mathcal{O}}$ as

$\expectationvalue{\hat{\mathcal{O}}(t)}=\lim_{L_{\text{B}}\to\infty}\Tr[\hat{\mathcal{O}}\hat{\rho}_{\text{tot}}(t)].$
(5)

In many cases, we will be interested in the long time behavior of the system
which corresponds to

$\lim_{t\rightarrow\text{large}}\expectationvalue{\hat{\mathcal{O}}(t)}=\lim_{t\rightarrow\text{large}}\left(\lim_{L_{\text{B}}\to\infty}\Tr[\hat{\mathcal{O}}\hat{\rho}_{\text{tot}}(t)]\right).$
(6)

Once again, the order of limits cannot be interchanged.

The above unitary picture of describing thermal baths is the starting point
for microscopic derivations of all standard approaches to open quantum
systems, such as non-equilibrium Green’s functions (NEGF) Haug and Jauho
(2008), Feynman-Vernon influence functional methods Kamenev (2011), quantum
master equations Breuer and Petruccione (2007) and quantum Langevin equations
Cortés  _et al._ (1985); Ford  _et al._ (1988); Benguria and Kac (1981); Dhar
and Sriram Shastry (2003); Dhar and Roy (2006). The crucial point is the
microscopic nature of the system and the macroscopic nature of the baths.
Given it’s microscopic description, the system is assumed to be accessible to
microscopic measurements. This is not true for the macroscopic baths whose
microscopic details are not accessible. Instead, only knowledge of some of
their macroscopic properties such as the temperatures, chemical potentials,
and spectral functions are assumed to be known. Consequently, effective
expressions for system observables and currents from baths are obtained after
integrating out the baths. The various standard techniques of open quantum
systems allow one to compute these, albeit at differing levels of
approximation Haug and Jauho (2008); Kamenev (2011); Breuer and Petruccione
(2007). With this microscopic setting in mind, we will now turn to describing
the necessary thermodynamic quantities and how they can be studied in this
context.

###  II.2 Description of the thermodynamics

We focus on the average thermodynamic currents, such as work, heat and entropy
production. The work performed by the external protocol is defined as the
change in energy of the global set-up,

$W_{\text{ext}}(t)=\expectationvalue{\hat{H}_{\text{tot}}(t)}-\expectationvalue{\hat{H}_{\text{tot}}(0)}.$
(7)

Under this sign convention, work is positive if done on the system and
negative if its done by the system. Since only the system Hamiltonian is time-
dependent, the work can be re-written as follows,

$W_{\text{ext}}(t)=\int_{0}^{t}\differential
t^{\prime}\derivative{\\!\expectationvalue{\hat{H}_{\text{tot}}(t^{\prime})}}{t^{\prime}}=\int_{0}^{t}\differential
t^{\prime}\expectationvalue{\partialderivative{\hat{H}_{\text{S}}(t^{\prime})}{t^{\prime}}}$
(8)

where we used the fact that $\partial\hat{H}_{\text{tot}}/\partial
t=\partial\hat{H}_{\text{S}}/\partial t$.

We define the entropy production as Refs. Esposito  _et al._ (2010); Reeb and
Wolf (2014); Landi and Paternostro (2021); Strasberg and Winter (2021),

$\displaystyle\begin{split}\Sigma(t)&=\lim_{L_{B}\to\infty}D\left(\hat{\rho}_{\text{tot}}(t)\leavevmode\nobreak\
\big{\|}\leavevmode\nobreak\
\hat{\rho}_{\text{S}}(0)\prod_{\alpha}\hat{\rho}_{\alpha}^{\text{th}}\right)\,,\end{split},$
(9)

where $D\left(\hat{\rho}\leavevmode\nobreak\ \middle\|\leavevmode\nobreak\
\hat{\sigma}\right)=\Tr[\hat{\rho}\log\hat{\rho}-\hat{\rho}\log\hat{\sigma}]$
is the quantum relative entropy between the density matrices $\hat{\rho}$ and
$\hat{\sigma}$. Since $D\left(\hat{\rho}\leavevmode\nobreak\
\middle\|\leavevmode\nobreak\ \hat{\sigma}\right)\geq 0$, we have
$\Sigma(t)\geq 0$, the second law of thermodynamics is therefore naturally
satisfied. Given this definition, one can recast the entropy production in
terms of the heat dissipated and the change in the von-Neumann entropy of the
system only

$\Sigma(t)=\delta S_{\text{S}}(t)+\sum_{\alpha}\beta_{\alpha}Q_{\alpha}(t)\,,$
(10)

where $\delta S_{\text{S}}$ is the change in Von-Neumann entropy of the
system, $\delta S_{\text{S}}=S_{\text{S}}(t)-S_{\text{S}}(0)$, with
$S(t)=-{\rm Tr}(\hat{\rho}\log\hat{\rho})$ and $Q_{\alpha}$ is the heat
dissipated into the $\alpha$th bath given by

$Q_{\alpha}(t)=\expectationvalue{\hat{H}_{\alpha}(t)}-\expectationvalue{\hat{H}_{\alpha}(0)}-\mu\quantity(\expectationvalue{\hat{N}_{\alpha}(t)}-\expectationvalue{\hat{N}_{\alpha}(0)})\,.$
(11)

Apart from the work done by the external protocol, there is a chemical work
due to the non-zero chemical potentials of the system. Chemical work done by
the system generates an increase in the number of particles within the baths.

$W_{\text{chem}}(t)=-\sum_{\alpha}\mu_{\alpha}\quantity(\expectationvalue{\hat{N}_{\alpha}(t)}-\expectationvalue{\hat{N}_{\alpha}(0)})\,.$
(12)

The total work is therefore given by the sum of the external work Eq. (7) and
the chemical work

$\displaystyle W(t)=W_{\text{chem}}(t)+W_{\text{ext}}(t).$ (13)

The change in internal energy is defined by

$\delta U(t)=W(t)-\sum_{\alpha}Q_{\alpha}(t)\,,$ (14)

which corresponds to the first law of thermodynamics. By substituting in the
definitions of both the total work, the dissipated heat, and the expression
for $\hat{H}_{\rm tot}$ given by Eq. 1, one finds that the change in internal
energy is

$\delta
U(t)=\expectationvalue{\hat{H}_{\text{S}}(t)+\sum_{\alpha}\hat{H}_{\text{S}\alpha}(t)}-\expectationvalue{\hat{H}_{\text{S}}(0)+\sum_{\alpha}\hat{H}_{\text{S}\alpha}(0)}.$
(15)

Intuitively, the change in the internal energy of the system only depends on
those components of the Hamiltonian that relate to the system.

We also call attention to the limit $L_{B}\to\infty$ in Eq. (9), which permits
definitions of the work, heat and internal energy, in accordance with
expectation value given by Eq. (5). This limit is important for thermodynamic
considerations: Eq. (10) has the form of a standard Clausius’s statement of
second law, if the temperatures and the chemical potentials of the baths can
be considered constant. If the baths are macroscopic, they have an infinite
capacity for heat and particles, and hence their macroscopic properties
temperatures and chemical potentials can be considered constant. So, for
infinitely large baths Eq. (9) is indeed a statement of the second law. On the
other hand, for finite bath sizes, this constancy assumption regarding the
temperatures and chemical potentials does not hold beyond a finite time thus
requires corrections for a consistent thermodynamic description Strasberg  _et
al._ (2021); Strasberg and Winter (2021).

The above definitions of heat and work seem to suggest that one needs to
access the full baths to describe the thermodynamics. However, this is not the
case. The heat and work defined above can be written in terms of currents from
the baths, which can be obtained without having access to the full microscopic
details of the baths. The energy current from the $\alpha$-th bath is

$J^{E}_{\alpha}(t)=-\derivative{\expectationvalue{\hat{H}_{\alpha}}}{t}=-i\expectationvalue{\commutator{\hat{H}_{\text{tot}}(t)}{\hat{H}_{\alpha}}}=i\expectationvalue{\commutator{\hat{H}_{\alpha}}{\hat{H}_{\text{S}\alpha}}},$
(16)

while the particle current is defined

$J^{P}_{\alpha}(t)=-\derivative{\expectationvalue{\hat{N}_{\alpha}}}{t}=-i\expectationvalue{\commutator{\hat{H}_{\text{tot}}(t)}{\hat{N}_{\alpha}}}=i\expectationvalue{\commutator{\hat{N}_{\alpha}}{\hat{H}_{\text{S}\alpha}}}.$
(17)

Finally, the heat current from the $\alpha$-th bath is thus defined as in
terms of the energy and particle currents accordingly

$J^{Q}_{\alpha}(t)=J^{E}_{\alpha}(t)-\mu_{\alpha}J^{P}_{\alpha}(t).$ (18)

Using these definitions, the heat dissipated into the $\alpha$th bath and the
chemical work can be recast explicity in terms of these currents as

$\displaystyle
Q_{\alpha}(t)=-\int_{0}^{t}dt^{\prime}J^{Q}_{\alpha}(t^{\prime}),\quad\leavevmode\nobreak\
W_{\text{chem}}(t)=\sum_{\alpha}\mu_{\alpha}\int_{0}^{t}dt^{\prime}J^{P}_{\alpha}(t^{\prime}),$
$\displaystyle W_{\text{ext}}(t)=\delta
U(t)-\int_{0}^{t}dt^{\prime}J^{E}_{\alpha}(t^{\prime}),$ (19)

where the last expression follows from the first law. These expressions, along
with Eqs. (10) and (15), show that the description of average thermodynamics
quantities can be described by knowing, as function of time, the energy and
the particle currents from the baths, the state of the system and the system-
bath coupling energy.

It is important to mention that the above definitions of thermodynamic
quantities, while being completely self-consistent, are not unique, especially
at strong coupling between the system and the baths Talkner and Hänggi (2020).
In particular, there are valid criticisms against the definition of $\delta
U(t)$ including the system-bath coupling energy, and the identification of
von-Neumann entropy as the thermodynamic entropy. However, if the system is
finite-dimensional and is driven (either by an external drive, or via chemical
potential and temperature biases, or both), the work and heat increase
indefinitely, while the change in internal energy $\delta U(t)$ and the change
in entropy of the system $\delta S_{\text{S}}(t)$ remain finite, regardless of
their actual definitions. Consequently, the $\delta U(t)$ and $\delta
S_{\text{S}}(t)$ become negligible in expressions of external work Eq. (II.2)
and entropy production Eq. (10) respectively. As a result, in such cases, the
long-time thermodynamics can be entirely and uniquely described in terms of
the currents from the baths.

So far the formalism we have described holds for arbitrary bath Hamiltonians.
Next we focus specifically on to the case where each bath is described by an
infinite number of fermionic or bosonic modes prepared in a thermal state. In
this case, the Hamiltonian of the $\alpha$th bath is given by

$\hat{H}_{\alpha}=\sum_{m=1}^{\infty}\omega_{m\alpha}\hat{b}_{m\alpha}^{\dagger}\hat{b}_{m\alpha},$
(20)

where $\hat{b}_{m\alpha}$ are fermionic or bosonic annihilation operators for
the $m$th mode of the $\alpha$th bath. We also fix the nature of the system-
bath couplings to

$\hat{H}_{\text{S}\alpha}=\sum_{m=1}^{\infty}\quantity(\lambda_{m\alpha}\hat{S}^{\dagger}_{\alpha}\hat{b}_{m\alpha}+\lambda_{m\alpha}^{*}\hat{b}_{m\alpha}^{\dagger}\hat{S}_{\alpha}),$
(21)

where $\hat{S}_{\alpha}$ is the system operator coupling with the $\alpha$th
bath. With this canonical model of thermal baths, it can be shown that the
dynamics of the system, the currents from the baths, and the system-bath
coupling energy are all uniquely determined by the bath spectral functions,

$\mathcal{J}_{\alpha}(\omega)=2\pi\sum_{m=1}^{\infty}\absolutevalue{\lambda_{m\alpha}}^{2}\delta\quantity(\omega-\omega_{m\alpha}),$
(22)

and the Fermi or the Bose distributions,
$f_{\alpha}(\omega)=(e^{\beta_{\alpha}(\omega-\mu_{\alpha})}\pm 1)^{-1}$.
Crucially, these properties do not depend on the fine microscopic details of
each bath. Any set of baths with the same spectral functions and same initial
temperatures and chemical potentials gives rise to the same dynamics of the
system, the currents from the baths, and the system-bath coupling energy.

While the $L_{B}\to\infty$ limit is crucial for a consistent description of
dynamics and thermodynamics of open quantum systems, it makes calculating the
the above thermodynamic quantities intractable. This is additionally
complicated by the fact that, without any further approximations, the reduced
dynamics of the system is non-Markovian. Only in absence of external driving,
if the system is also quadratic in bosonic or fermionic creation and
annihilation operators, the long-time results can be obtained without making
further approximations using standard analytical techniques like NEGF,
Landauer-Buttiker formalism etc. In presence of external driving, even for
quadratic cases, all approaches often have to rely on perturbative techniques
(for example, Ludovico  _et al._ (2016a); Bruch  _et al._ (2018); Riera-
Campeny  _et al._ (2019); Brandner (2020)).

In absence of exact analytical techniques, arguably, the most direct way to
numerically simulate the dynamics is to consider finite, but large enough,
baths and carry out the full unitary evolution. A systematic way to do this is
given by the chain-mapping of the baths Nazir and Schaller (2018); de Vega
_et al._ (2015); Chin  _et al._ (2010); Strasberg  _et al._ (2018); Prior  _et
al._ (2010); Woods  _et al._ (2015); Woods and Plenio (2016); Mascherpa  _et
al._ (2017); Garg  _et al._ (1985) (see Appendix A). This utilizes the
understanding that we can choose a convenient way to microscopically model the
baths, such that the bath spectral functions are left invariant, without
changing the dynamics or the thermodynamics. However, the size of the baths
required to simulate the dynamics grows with time. If the drive is periodic,
and a periodic steady state is reached within the time possible to simulate
with the largest size baths, the whole dynamics can be inferred. If this is
not the case, for example if the drive is not periodic, the long-time behavior
cannot be obtained by this technique.

In the following section, we introduce the mesoscopic leads approach which
allows us to obtain completely non-perturbative numerically exact results, and
thereby accurately describe the thermodynamics of externally driven open
quantum systems, up to arbitrarily long times. We will benchmark these results
against the brute force simulation with chain-mapping technique.

##  III The mesoscopic leads approach

Figure 1: a) The strong coupling between system $S$ and the bath can be can be
approximated by $L$ leads which are each independently coupled to independent
infinite baths. b) A close up of the one of the leads with energy
$\epsilon_{k\alpha}$ and the system, with coupling strength
$\kappa_{k\alpha}$. The lead is weakly coupled to the infinite bath with decay
rate $\gamma_{k\alpha}$.

###  III.1 Microscopic derivation

We now show how the mesoscopic leads approach Imamoglu (1994); Garraway
(1997a, b); Sánchez  _et al._ (2006); Subotnik  _et al._ (2009); Dzhioev and
Kosov (2011); Ajisaka  _et al._ (2012); Ajisaka and Barra (2013); Arrigoni
_et al._ (2013); Dorda  _et al._ (2014); Chen  _et al._ (2014); Zelovich  _et
al._ (2014); Dorda  _et al._ (2015); Hod  _et al._ (2016); Gruss  _et al._
(2016); Schwarz  _et al._ (2016); Dorda  _et al._ (2017); Elenewski  _et al._
(2017); Gruss  _et al._ (2017); Zelovich  _et al._ (2017); Tamascelli  _et
al._ (2018); Lemmer  _et al._ (2018); Oz  _et al._ (2019); Chen  _et al._
(2019); Zwolak (2020); Wójtowicz  _et al._ (2020); Chiang and Hsu (2020);
Brenes  _et al._ (2020); Lotem  _et al._ (2020); Fugger  _et al._ (2020);
Wójtowicz  _et al._ (2021); Elenewski  _et al._ (2021) can be microscopically
derived in the setting described in previous section. A schematic of how the
each bath is microscopically represented in this approach is given in Fig.
1(a). This can be done for any given bath spectral function. The first step of
this approach is the observation that any continuous bath spectral functions
can be approximated by discretizing it into $L_{\alpha}$ (for simplicity,
equally spaced) points, $\\{\epsilon_{k}\\}$, $k=1$ to $L_{\alpha}$, and
defining

$\mathcal{J}^{\text{eff}}_{\alpha}(\omega)=\sum_{k=1}^{L_{\alpha}}\frac{|\kappa_{k\alpha}|^{2}\gamma_{k}}{(\omega-\epsilon_{k\alpha})^{2}+(\gamma_{k\alpha}/2)^{2}},$
(23)

with

$\displaystyle\kappa_{k\alpha}=\sqrt{\frac{\mathcal{J}_{\alpha}(\epsilon_{k})e_{k\alpha}}{2\pi}},\quad\gamma_{k\alpha}=e_{k\alpha},\quad
e_{k\alpha}=\epsilon_{k+1\leavevmode\nobreak\ \alpha}-\epsilon_{k\alpha}.$
(24)

With this definition, it can be checked that
$\mathcal{J}^{\text{eff}}_{\alpha}(\omega)$ tends to
$\mathcal{J}_{\alpha}(\omega)$ as $L_{\alpha}$ increases. Hence, for a finite
but large enough $L_{\alpha}$ (which corresponds to small enough $e_{k}$)
$\mathcal{J}^{\text{eff}}_{\alpha}(\omega)$ is a controlled approximation to
$\mathcal{J}_{\alpha}(\omega)$. Eq. (23) shows that each bath $\alpha$, with
spectral function $\mathcal{J}_{\alpha}(\omega)$, can be decomposed into
$L_{\alpha}$ different baths, each with a specific Lorentzian spectral
function.

A bath with a Lorentzian spectral density, however, yields the same dynamics
as if the system was coupled with strength $\kappa_{k\alpha}$) to an
additional fermionic mode $\hat{a}_{ka}$, with energy $\epsilon_{k\alpha}$,
which itself is coupled to its own bath, with a flat spectral function
$\gamma_{\alpha}$ Brenes  _et al._ (2020), see Fig. 1(b). This can be viewed
as a single step of the chain mapping of Appendix A, Eq. (A)). We thus arrive
at the mesoscopic lead, consisting of $L_{\alpha}$ modes,

$\displaystyle\hat{H}_{L_{\alpha}}=\sum_{k=1}^{L_{\alpha}}\epsilon_{k\alpha}\hat{a}_{k\alpha}^{\dagger}\hat{a}_{k\alpha},$
(25)

while the system-bath coupling becomes the system-lead coupling,

$\displaystyle\hat{H}_{S\alpha}=\hat{H}_{SL_{\alpha}}=\sum_{k=1}^{L_{\alpha}}\kappa_{k\alpha}\left(\hat{S}_{\alpha}^{\dagger}\hat{a}_{k\alpha}+\hat{a}_{k\alpha}^{\dagger}\hat{S}_{\alpha}\right).$
(26)

Now, we define extended state of the system and the lead modes
$\hat{\rho}_{\text{ext}}$ and the extended Hamiltonian,

$\displaystyle\hat{H}_{\text{ext}}(t)=\hat{H}_{S}(t)+\sum_{\alpha}\left(\hat{H}_{L_{\alpha}}+\hat{H}_{SL_{\alpha}}\right),$
(27)

where $\hat{H}_{L_{\alpha}}$ is the Hamiltonian of the bath attached to the
$\alpha$th lead mode. The final step is to integrate out the residual baths to
obtain an effective equation of motion for the system and the lead modes to
the leading order in $e_{k}$. The crucial point to note in doing so is that
$\kappa_{k\alpha}\propto\sqrt{e_{k\alpha}}$ and $\gamma_{k\alpha}\propto
e_{k\alpha}$. So both $H_{SL_{\alpha}}$ and the coupling between the lead
modes and their own residual baths are small in the limit where
$\mathcal{J}^{\text{eff}}_{\alpha}(\omega)$ is a good approximation to
$\mathcal{J}_{\alpha}(\omega)$. Using standard techniques Brenes  _et al._
(2020), this allows us to obtain the following quantum master equation for
$\hat{\rho}_{\text{ext}}$,

$\derivative{\hat{\rho}_{\text{ext}}}{t}=-i\commutator{\hat{H}_{\text{ext}}(t)}{\hat{\rho}_{\text{ext}}(t)}+\sum_{\alpha}\hat{\mathcal{L}}_{\alpha}\quantity{\hat{\rho}_{\text{ext}}(t)},$
(28)

where

$\displaystyle\begin{split}&\hat{\mathcal{L}}_{\alpha}\quantity{\hat{\rho}_{\text{ext}}(t)}\\\
&=\sum_{k=1}^{L_{\alpha}}\gamma_{k\alpha}e^{\beta_{\alpha}(\epsilon_{k\alpha}-\mu_{\alpha})}f_{k\alpha}\quantity[\hat{a}_{k\alpha}\hat{\rho}_{\text{ext}}(t)\hat{a}_{k\alpha}^{\dagger}-\frac{1}{2}\anticommutator{\hat{a}_{k\alpha}^{\dagger}\hat{a}_{k\alpha}}{\hat{\rho}_{\text{ext}}(t)}]\\\
&+\sum_{k=1}^{L_{\alpha}}\gamma_{k\alpha}f_{k\alpha}\quantity[\hat{a}_{k\alpha}^{\dagger}\hat{\rho}_{\text{ext}}(t)\hat{a}_{k\alpha}-\frac{1}{2}\anticommutator{\hat{a}_{k\alpha}\hat{a}_{k\alpha}^{\dagger}}{\hat{\rho}_{\text{ext}}(t)}],\end{split}$
(29)

with
$f_{k\alpha}=\quantity[e^{\beta_{\alpha}(\epsilon_{k\alpha}-\mu_{\alpha})}\pm
1]^{-1}$ being the Fermi or the Bose distribution function. This is the
central equation of the mesoscopic leads approach. Here, each bath is modelled
by a finite set of modes, each mode being damped via a local Lindblad
operator. The local Lindblad operator is such that it would take the mode to
its thermal state if it was uncoupled from the system. To obtain correct
dynamics at all times, the initial state Eq. (2) will also have to be closely
approximated. For this, we choose

$\displaystyle\hat{\rho}_{\text{ext}}(0)=\hat{\rho}_{\text{S}}(0)\prod_{\alpha}\frac{e^{-\beta_{\alpha}(\hat{H}_{L_{\alpha}}-\mu_{\alpha}\hat{N}_{L_{\alpha}})}}{Z_{L_{\alpha}}},$
(30)

where
$\hat{N}_{L_{\alpha}}=\sum_{k=1}^{L_{\alpha}}\hat{a}_{k\alpha}^{\dagger}\hat{a}_{k\alpha}$
is the total number operator for the lead modes, and $Z_{L_{\alpha}}$ is the
corresponding partition function. Clearly, this initial state becomes a close
approximation to Eq. (2) for large enough lead modes. The above derivation
shows that by increasing then number of lead modes, the dynamics of the system
obtained from the above quantum master equation will converge to the dynamics
of the system in the presence of baths with spectral functions
$\mathcal{J}_{\alpha}(\omega)$.

There are several salient features of the mesoscopic leads approach that makes
it advantageous over most other approaches. Although this approach has been
used in presence of external driving before Chen  _et al._ (2014); Oz  _et
al._ (2019), these important features, to our knowledge, have hardly been
emphasized and appreciated in previous works. The first important point to
note is that the above derivation holds for arbitrary types of driving in the
system. This fact is quite non-trivial because, in microscopic derivations of
quantum master equations, the dissipative part of the quantum master equation
can change depending on the nature of the drive Schnell  _et al._ (2020,
2021). However, in the mesoscopic leads approach, the fact that both
$\kappa_{k\alpha}$ and $\gamma_{k\alpha}$ are small allows us to obtain the
same dissipators to the leading order in $e_{k}$ for all kinds of driving.

The small parameter in the derivation is $e_{k\alpha}$, which reduces on
increasing the number of lead modes. This is not a small parameter of the
physical set-up we want to model, but rather, a small parameter that controls
the error in the numerical simulation. Thus, this approach can give fully non-
perturbative results, with controlled errors which can be reduced by
increasing the number of lead modes. Further, although the individual coupling
between each lead mode and the system is small, overall, the effective
spectral function $\mathcal{J}^{\text{eff}}_{\alpha}(\omega)$ is not small.
Rather, it is a good approximation to $\mathcal{J}_{\alpha}(\omega)$. The
mesoscopic leads approach can thus treat arbitrary strengths of system-bath
coupling. Although the evolution equation for the extended density matrix
$\hat{\rho}_{\rm ext}$ is an additive Lindbald equation, the reduced dynamics
of the system (i.e, not including the leads) can be fully non-Markovian can
have completely non-additive contributions from the leads.

Finally and perhaps most importantly, the extended set-up of the system and
the lead modes is of finite size. Due to the Lindblad damping, it incorporates
the effect of infinite number of degrees of freedom of the baths.
Consequently, unlike the chain-mapping technique described in the previous
section, long time simulation can be done keeping the number of lead modes
finite, provided the spectral function is accurately enough represented. This
is a crucial advantage of this approach over the brute force method using
chain-mapping.

Along with the above salient features, all quantities required for description
of thermodynamics can be obtained from the mescoscopic leads approach, as we
discuss in the next subsection

###  III.2 Thermodynamics from mesoscopic leads approach

The mesoscopic leads approach simulates the evolution of
$\hat{\rho}_{\text{ext}}$. Since this extended state evolves under a Lindblad
equation it is Markovian wheras the orignal system alone still undegoes non-
markovian evolution. This is known as a markovian embedding. As such, it
provides access to expectation values of operators in the extended Hilbert
space of system plus leads, as well as the rate of change of expectation
values of such operators. We need to write all the quantities required for the
description of thermodynamics in terms of such quantities.

As discussed before, we need to know the state of the system, the system-bath
coupling energy, and the particle and the energy currents from the baths as a
function of time. The state of the system is directly obtained from the
mesoscopic leads approach. The system-bath coupling energy is just the
expectation value of the system-lead coupling operator, which is also directly
obtained because it is an operator in the extended Hilbert space. For
obtaining the currents from the baths, we note that the mesoscopic leads
approach consider the following microscopic structure of each bath,

$\displaystyle\hat{H}_{\alpha}=\hat{H}_{L_{\alpha}}+\hat{H}_{LE_{\alpha}}+\hat{H}_{E_{\alpha}},$
(31)

where $\hat{H}_{E_{\alpha}}$ describe the composite Hamiltonian of the
residual baths of all the lead modes, and $\hat{H}_{LE_{\alpha}}$ describe the
corresponding composite coupling between the lead modes and their residual
baths. The particle current from the $\alpha$th bath, defined in Eq. (17), can
now be seen to be given by

$\displaystyle
J^{P}_{\alpha}(t)=i\expectationvalue{\commutator{\hat{N}_{L_{\alpha}}}{\hat{H}_{SL_{\alpha}}}},$
(32)

where
$\hat{N}_{L_{\alpha}}=\sum_{k=1}^{L_{\alpha}}\hat{a}_{k\alpha}^{\dagger}\hat{a}_{k\alpha}$
is the total number operator for the lead modes. This is also the expectation
value of an operator in the extended Hilbert space, so can be directly
obtained from the mesoscopic leads approach. However, for energy current,
defined in Eq. (16), we get

$\displaystyle
J^{E}_{\alpha}(t)=i\expectationvalue{\commutator{\hat{H}_{L_{\alpha}}+\hat{H}_{LE_{\alpha}}}{\hat{H}_{SL_{\alpha}}}}.$
(33)

The fact that $\hat{H}_{SL_{\alpha}}$ and $\hat{H}_{LE_{\alpha}}$ does not
commute does not let us write the energy current as an expectation value in
the extended Hilbert space. However, we note that

$\displaystyle\begin{split}\derivative{\expectationvalue{\hat{H}_{SL_{\alpha}}}}{t}&=-i\expectationvalue{\commutator{\hat{H}_{SL_{\alpha}}}{\hat{H}_{\text{S}}(t)}}-i\expectationvalue{\commutator{\hat{H}_{SL_{\alpha}}}{\hat{H}_{L_{\alpha}}+\hat{H}_{LE_{\alpha}}}}\\\
&=J^{E}_{\alpha}-i\expectationvalue{\commutator{\hat{H}_{SL_{\alpha}}}{\hat{H}_{\text{S}}(t)}},\end{split}.$
(34)

This gives the energy current from the $\alpha$th bath in terms of expectation
values of operators on the extended Hilbert space, and rate of change of such
operators,

$\displaystyle
J^{E}_{\alpha}(t)=\derivative{\expectationvalue{\hat{H}_{SL_{\alpha}}}}{t}+i\expectationvalue{\commutator{\hat{H}_{SL_{\alpha}}}{\hat{H}_{\text{S}}(t)}}.$
(35)

This can therefore be obtained accurately from the mesoscopic lead approach.
Calculating the above expectation values from the quantum master equation, Eq.
(28) finally yields

$J^{E}_{\alpha}(t)=i\expectationvalue{\commutator{\hat{H}_{L_{\alpha}}}{\hat{H}_{SL_{\alpha}}}}+\Tr[\hat{H}_{SL_{\alpha}}\hat{\mathcal{L}}_{\alpha}\quantity{\hat{\rho}_{\text{ext}}}].$
(36)

The second term on the right hand side would naively not be expected.

Figure 2: Schematic representation of a system coupled to mesoscopic leads.
The arrows represent the two set of currents associated with the left bath,
the “external” ($I_{L}^{P}$ and $I_{L}^{E}$) and “internal” ($J_{L}^{P}$ and
$J_{L}^{E}$).

The proper definition of currents is important to describe thermodynamics in
the presence of an external drive. Given the quantum master equation, Eq.
(28), it is tempting to define the particle and the energy currents as

$\displaystyle
I^{P}_{\alpha}(t)=\Tr(\hat{N}_{L_{\alpha}}\hat{\mathcal{L}}_{\alpha}\quantity[\hat{\rho}_{\text{ext}}(t)]),\quad\leavevmode\nobreak\
I^{E}_{\alpha}(t)=\Tr(\hat{H}_{L_{\alpha}}\hat{\mathcal{L}}_{\alpha}\quantity[\hat{\rho}_{\text{ext}}(t)]).$
(37)

However, these describe currents from the residual baths into the mesoscopic
leads. They do not correspond to the current from the baths into the system
defined in Eqs. (17), (16), see Fig. 2. For a driven system,
$I^{P}_{\alpha}(t)\neq J^{P}_{\alpha}(t)$, $I^{E}_{\alpha}(t)\neq
J^{E}_{\alpha}(t)$. In the absence of drive, and at the NESS, due to
continuity equations, these two different definitions will agree. For similar
reasons, with periodic drives, when a periodic NESS (limit cycle) is reached,
the time-period average of these two different currents become the same.
However, instantaneously they will still be different. The currents from the
baths, Eqs. (17), (16), are independent of the microscopic modelling of the
baths. If instead of the mesoscopic leads, a different microscopic modelling
of the baths giving the same spectral functions is used, like the chain-
mapping, the currents would be exactly the same. However, the currents
$I^{P}_{\alpha}(t)$, $I^{E}_{\alpha}(t)$ depend on the particular microscopic
modelling of the baths as mesoscopic leads, and may not be replicated in a
different microscopic model of the baths, like the chain-mapping.

From the above discussion, it is clear that all quantities required for the
description of thermodynamics can be obtained accurately from the mesoscopic
leads approach. A crucial point to note is that, throughout the entire
discussion till now, we have left the system Hamiltonian and the system
operators coupling to the baths completely arbitrary. The entire discussion
therefore holds for arbitrary system Hamiltonians, with arbitrary driving.
However, it is non-trivial to simulate the dynamics of
$\hat{\rho}_{\text{ext}}(t)$. Simulation in the presence of many-body
interactions (i.e, higher than quadratic terms) require tensor network
techniques Wójtowicz  _et al._ (2020); Lotem  _et al._ (2020); Brenes  _et
al._ (2020); Fugger  _et al._ (2020). While this may be numerically expensive,
it is still feasible for a range of parameters where many other techniques
fail. And wherever possible, it gives numerically exact non-perturbative
results. This has been shown for quantum transport and thermodynamics in
absence of any explicit time-dependent drive Brenes  _et al._ (2020).

In this work, however, we will not consider many-body interactions. In absence
of many-body interactions, as will be shown below, numerically exact results
can be obtained quite elegantly by simulating a Lyapunov equation. For driven
systems, such exact results are often not amenable to other numerical or
analytical techniques (except for brute force numerics via chain mapping).

###  III.3 The Lyapunov equation for non-interacting systems

Consider a system Hamiltonian of the form

$\hat{H}_{\text{S}}(t)=\sum_{\ell,m=1}^{L_{S}}[\mathbf{H_{\text{S}}}(t)]_{\ell
m}\hat{c}^{\dagger}_{\ell}\hat{c}_{m},$ (38)

where $\mathbf{H_{\text{S}}}(t)$ is a $L_{S}\times L_{S}$ Hermitian matrix,
sometimes called single-particle Hamiltonian, with time $t$ as a parameter.
This describes a number conserving non-interacting fermionic system of $L_{S}$
sites, in a lattice of arbitrary dimension and geometry with an arbitrary
external type of external drive. Consider this system to be coupled via number
conserving bilinear coupling with baths at various sites. In the mesoscopic
leads approach, the system-lead coupling for the lead describing the bath
attached at the $\alpha$th site is

$\displaystyle\hat{H}_{SL_{\alpha}}=\sum_{k=1}^{L_{\alpha}}\kappa_{k\alpha}\left(\hat{c}_{\alpha}^{\dagger}\hat{a}_{k\alpha}+\hat{a}_{k\alpha}^{\dagger}\hat{c}_{\alpha}\right).$
(39)

The extended set-up of system and leads can now be written in the form

$\displaystyle\hat{H}_{\text{ext}}(t)=\sum_{\ell,m=1}^{L_{S}+\sum_{\alpha}L_{\alpha}}[\mathbf{H_{\text{ext}}}(t)]_{\ell
m}\hat{d}^{\dagger}_{\ell}\hat{d}_{m},$ (40)

where $\hat{d}_{m}$ is a fermionic annihilation operator of the either a
system site or a lead site. For such cases, if the initial state of the system
is Gaussian, the dynamics of the whole set-up remains Gaussian at all times.
All properties of the extended system can be obtained by calculating the
correlation matrix, also called the single-particle density matrix,
$\mathbf{C}(t)$, with elements Landi  _et al._ (2021)

$\displaystyle\mathbf{C}_{pq}(t)=\Tr[\hat{d}_{q}^{\dagger}\hat{d}_{p}\hat{\rho}_{\text{ext}}(t)].$
(41)

Since $\hat{\rho}_{\text{ext}}(t)$ is Gaussian at all times, it can be
obtained exactly from $\mathbf{C}_{pq}(t)$. The crucial simplification is that
one can directly write down the equation of motion for $\mathbf{C}_{pq}(t)$
from Eq. (28), in the form of a Lyapunov equation, as we discuss below. For
simplicity and relevance to later examples, in the following, we consider two
baths, but the discussion can be straightforwardly generalized to more than
two baths.

For two baths, the matrices $\mathbf{H_{\text{ext}}}(t)$ and $\mathbf{C}(t)$
can be written in the block form

$\displaystyle\mathbf{H_{\text{ext}}}(t)=\left[\begin{array}[]{ccc}\mathbf{H}_{\text{S}}(t)&\mathbf{H}_{SL_{1}}&\mathbf{H}_{SL_{2}}\\\
\mathbf{H}_{SL_{1}}^{\dagger}&\mathbf{H}_{L_{1}}&\mathbf{0}\\\
\mathbf{H}_{SL_{2}}^{\dagger}&\mathbf{0}&\mathbf{H}_{L_{2}}\end{array}\right],$
(45)
$\displaystyle\mathbf{C}(t)=\left[\begin{array}[]{ccc}\mathbf{C}_{\text{S}}(t)&\mathbf{C}_{SL_{1}}(t)&\mathbf{C}_{SL_{2}}(t)\\\
\mathbf{C}_{SL_{1}}^{\dagger}(t)&\mathbf{C}_{L_{1}}(t)&\mathbf{C}_{L_{1}L_{2}}(t)\\\
\mathbf{C}_{SL_{2}}^{\dagger}(t)&\mathbf{C}_{L_{1}L_{2}}^{\dagger}(t)&\mathbf{C}_{L_{2}}(t)\end{array}\right].$
(49)

Here, the $L_{S}\times L_{\alpha}$ matrix $\mathbf{H}_{SL_{\alpha}}$ gives the
coupling between the system and the lead representing the $\alpha$th bath, and
$\mathbf{H}_{B_{\alpha}}=\text{diag}\\{\epsilon_{k\alpha}\\}$ gives the
Hamiltonian of the $\alpha$th lead. Similarly, the $L_{S}\times L_{S}$
dimensional matrix $\mathbf{C}_{S}(t)$ gives the system correlation matrix,
the $L_{S}\times L_{\alpha}$ dimensional matrix $\mathbf{C}_{SL_{\alpha}}(t)$
gives the correlations between the system and the $\alpha$th lead,
$\mathbf{C}_{L_{\alpha}}(t)$ gives the correlation matrix corresponding to the
$\alpha$th lead, and $\mathbf{C}_{L_{1}L_{2}}(t)$ gives the correlations
between the two leads. With the correlation matrix written in this form, its
equation of motion, obtained from Eq. (28), is given by the following
continuous-time differential Lyapunov equation

$\derivative{\mathbf{C}}{t}=-\left(\mathbf{W}(t)\mathbf{C}(t)+\mathbf{C}(t)\mathbf{W}^{\dagger}(t)\right)+\mathbf{F},$
(50)

with $\mathbf{W}(t)=i\mathbf{H_{\text{ext}}}(t)+\mathbf{\Upsilon}/2$, where

$\displaystyle\mathbf{\Upsilon}=\left[\begin{array}[]{ccc}\mathbf{0}&\mathbf{0}&\mathbf{0}\\\
\mathbf{0}&\mathbf{\Upsilon}_{1}&\mathbf{0}\\\
\mathbf{0}&\mathbf{0}&\mathbf{\Upsilon}_{2}\end{array}\right],\leavevmode\nobreak\
\leavevmode\nobreak\
\mathbf{F}=\left[\begin{array}[]{ccc}\mathbf{0}&\mathbf{0}&\mathbf{0}\\\
\mathbf{0}&\mathbf{F}_{1}&\mathbf{0}\\\
\mathbf{0}&\mathbf{0}&\mathbf{F}_{2}\end{array}\right]$ (57)

and $\mathbf{\Upsilon}_{\alpha}=\text{diag}\\{\gamma_{k\alpha}\\}$,
$\mathbf{F}_{\alpha}=\text{diag}\\{\gamma_{k\alpha}f_{k\alpha}\\}$. The
dynamics can be obtained by integrating this matrix equation numerically, say
via Runge-Kutta methods. The size of the matrices required scales linearly
with total number of sites in the system and the leads, and not exponentially,
as it would have been if we were to directly consider the evolution of
$\hat{\rho}_{\text{ext}}(t)$. All quantities required for thermodynamics can
also be obtained by knowing $\mathbf{C}(t)$. Moreover, this is possible for
arbitrary types of driving, which is beyond the reach of most other analytical
or numerical techniques. Further simplification is possible if the drive is
periodic, as we show in the next subsection.

###  III.4 Floquet solution of the periodically-driven Lyapunov equation

The dynamics of Eq. (50), in the presence of a periodic drive, will be
characterized by a transient, followed by a limit cycle where $\mathbf{C}(t)$
becomes time-periodic. If one is interested only in the limit cycle, naively
this would require solving Eq. (50) over many periods, which may be expensive.
Instead, one can employ the following method For simplicity, we assume the
drive takes the form $\mathbf{W}(t)=\mathbf{W}_{0}+\mathbf{W}_{1}\cos(\omega
t)$. The generalisation to multiple harmonics is straightforward. We then
attempt a solution of the form
$\mathbf{C}(t)=\sum_{n=-\infty}^{\infty}\mathbf{C}_{n}e^{in\omega t}$.
Hermiticity of $\mathbf{C}$ implies that
$\mathbf{C}_{n}^{\dagger}=\mathbf{C}_{-n}$. Plugging this in Eq. (50) yields a
set of recursive algebraic equations

$\displaystyle
in\omega\mathbf{C}_{n}+\mathbf{W}_{0}\mathbf{C}_{n}+\mathbf{C}_{n}\mathbf{W}_{0}^{\dagger}=\mathbf{F}\delta_{n,0}$
$\displaystyle-\frac{1}{2}\mathbf{W}_{1}(\mathbf{C}_{n-1}+\mathbf{C}_{n+1})$
(58)
$\displaystyle-\frac{1}{2}(\mathbf{C}_{n-1}+\mathbf{C}_{n+1})\mathbf{W}_{1}^{\dagger}.$

This can now be solved iteratively, as in the Gauss-Siedel method Golub and
Loan (1996). First, one takes $\mathbf{C}_{n\neq 0}=0$ and solve Eq. (58) for
$\mathbf{C}_{0}$. Then the result is used to set up 3 new equations for
$\mathbf{C}_{-1},\mathbf{C}_{0},\mathbf{C}_{1}$, which in turn is used to set
up 5 equations for
$\mathbf{C}_{-2},\mathbf{C}_{-1},\mathbf{C}_{0},\mathbf{C}_{1},\mathbf{C}_{2}$,
and so on.

The solution is greatly facilitated when $\mathbf{W}_{0}$ is diagonalizable
(which is almost always the case); that is,
$\mathbf{W}_{0}=\mathbf{S}\mathbf{\Lambda}\mathbf{S}^{-1}$, where
$\mathbf{\Lambda}$ is a diagonal matrix containing the eigenvalues (which are
generally complex). Define
$\tilde{\mathbf{C}}_{n}=\mathbf{S}^{-1}\mathbf{C}_{n}(\mathbf{S}^{-1})^{\dagger}$
and similarly for $\tilde{\mathbf{F}}$ and $\tilde{\mathbf{W}}_{1}$. Then Eq.
(58) can be rewritten element-wise as

$\displaystyle[\tilde{\mathbf{C}}_{n}]_{i,j}=\frac{1}{in\omega+\Lambda_{i}+\Lambda_{j}^{*}}\Big{[}\delta_{n,0}\tilde{\mathbf{F}}$
$\displaystyle-\frac{1}{2}\tilde{\mathbf{W}}_{1}(\tilde{\mathbf{C}}_{n-1}+\tilde{\mathbf{C}}_{n+1})$
(59)
$\displaystyle-\frac{1}{2}(\tilde{\mathbf{C}}_{n-1}+\tilde{\mathbf{C}}_{n+1})\tilde{\mathbf{W}}_{1}\Big{]}_{i,j}.$

The reason why this is advantageous is because the largest overhead in Eq.
(58) is the linear system of equations associated with the fact that
$\mathbf{W}_{0}$ is not diagonal. In Eq. (59), the only computationally
expensive part is the diagonalization of $\mathbf{W}_{0}$, which only has to
be performed once. Afterwards, all operations involve only simple matrix
multiplications. After convergence, one recovers
$\mathbf{C}_{n}=\mathbf{S}\tilde{\mathbf{C}}_{n}\mathbf{S}^{\dagger}$ and use
this to construct
$\mathbf{C}(t)=\sum_{n=-\infty}^{\infty}\mathbf{C}_{n}e^{in\omega t}$. This
yields the solution within the limit cycle; that is, valid through an entire
period of oscillation. We now test out our technique on a number of models
that highlight the power and versatility of the methodology to explore
thermodynamics beyond the reach of most other approaches.

Figure 3: (a) Particle current as a function of time in the driven-resonant
model coupled to two baths for increasing number of modes $L$ in each lead,
where both leads have the same number of modes. The red curve shows the
results obtained the chain-mapping (CM) simulation, which is exact up to a
finite time, indicated with a dashed line. (b) Same as (a), but for the energy
current. The inset shows a zoomed vision of the short-time scale. (c) Same as
(a) and (b), but for the entropy production rate. In these calculations we
used $\Gamma=1$, $W=8$, $W^{*}=4$, $T_{L}=T_{R}=1$, $\mu_{L}=-\mu_{R}=0.5$,
$A=1$, $\omega=0.25$, $p_{0}=0.5$ and $L_{\text{lin}}/L_{\text{log}}=0.1$.

Figure 4: (a) Comparison of the particle current from the left bath in the
two-terminal driven resonant level with the currents from the residual baths
into the mesoscopic leads. (b) Same, but for the energy currents. In these
calculations, both leads had $L=100$ modes and the other parameters were
$\Gamma=1$, $W=8$, $W^{*}=4$, $T_{L}=T_{R}=1$, $\mu_{L}=-\mu_{R}=0.5$, $A=1$,
$\omega=0.25$, $p_{0}=0.5$ and $L_{\text{lin}}/L_{\text{log}}=0.1$.

Figure 5: Comparison of the entropy of the dot in the driven resonant-level
model with two terminals using the mesoscopic lead approach (blue curve)
against the behaviour predicted by the Pauli master equation [Eq (63)] (orange
curve) for increasing values of temperature, where the temperatures of both
baths are set to the same value. The following temperatures were used: (a)
$T=0.1$ (b) $T=1$ (c) $T=5$ (d) $T=10$. The amplitude of the driving is $A=1$
and the driving frequency is $\omega=25$. In the mesoscopic lead simulations
the other parameters were $W=8$, $W^{*}=4$, $L=100$,
$L_{\text{log}}/L_{\text{lin}}=0.1$, $\Gamma=1$ and $p_{0}=0.5$.

##  IV Thermodynamics of a driven resonant-level

As a first example, we apply the formalism to study the thermodynamics of the
driven resonant-level model Esposito  _et al._ (2015a, b); Ludovico  _et al._
(2016a, b); Bruch  _et al._ (2018); Haughian  _et al._ (2018); Oz  _et al._
(2019). This provides a simple model in which we can benchmark our results
against other methodologies employed in the literature. The system is composed
of a single dot with an externally controlled, time-dependent energy
$\epsilon(t)$. The Hamiltonian is

$\hat{H}_{\text{S}}(t)=\epsilon(t)\hat{c}^{\dagger}\hat{c},$ (60)

where $\hat{c}^{\dagger}$ and $\hat{c}$ are the fermionic creator and
annihilator operators of the system. We focus on the two-terminal case, in
which the dot is coupled to two baths, labelled by L and R. Our framework can,
in principle, handle the case in which each bath is described by an arbitrary,
structured spectral density. For simplicity, however, we consider the case in
which both baths are described by the same flat spectral density, given by

$\mathcal{J}(\omega)=\begin{cases}\Gamma,&\omega\in\quantity[-W,W]\\\
0,&\text{otherwise},\end{cases}$ (61)

where $\Gamma$ is the system-bath coupling strength and $W$ is a hard cut-off.

Within the mesoscopic leads approach, each bath is described by a damped lead
with $L$ modes. Fixing a set of energies $\\{\epsilon_{k}\\}$ that
appropriately sample the spectral density, the couplings
$\\{\kappa_{k\alpha}\\}$ and the damping rates $\\{\gamma_{k\alpha}\\}$ are
fixed by Eq. (24). Here, we follow Refs. Schwarz  _et al._ (2016);
Weichselbaum  _et al._ (2009) and choose the energies $\\{\epsilon_{k}\\}$ via
linear-logarithmic discretization. In this procedure, $L_{\text{lin}}$
energies are chosen to be linearly spaced in some frequency window
$[-W^{*},W^{*}]$, then $L_{\text{log}}$ energies are chosen to be
logarithmically spaced in each of the intervals $[-W,W^{*}]$ and $[W^{*},W]$,
so that the total number of modes if $L_{\text{lin}}+2L_{\text{log}}$.
Throughout all the examples in this section, we fix $\Gamma=1$ and use this to
set the energy scale. Furthermore, we fix $W=8$ and $W^{*}=4$.

We take the energy of the dot to be driven by a sine-protocol of the form

$\epsilon(t)=A\sin(\omega t).$ (62)

We study the full time-evolution of the extended system by numerically
integrating the Lyapunov equation (50). In Fig. 3 (a) and (b), we show the
particle and energy current, respectively, computed with $A=1$ and
$\omega=0.25$. We set the baths to have equal temperature
$T_{\text{L}}=T_{\text{R}}=1$, but opposite chemical potentials
$\mu_{\text{L}}=-\mu_{\text{R}}=0.5$, and the initial population of the dot to
be $p_{0}=0.5$. We take both leads to have the same number of modes $L$, and
show the results for increasing values of $L$. For comparison, we also present
results obtained via brute-force simulation with chain mapping (Appendix A),
which is numerically exact only up to a finite time, determined by the size of
the chain taken, indicated with a vertical dashed lines. It is completely
clear that with increase in number of lead modes, the results converge to that
obtained by chain-mapping method up to the time it remains valid. In Fig. 3
(b), the inset shows a zoom of the short time-scale. Remarkably, even the
short-time dynamics is well described.

It is worth re-stating that in the mesoscopic-leads formalism one has full
access to the state of the system at all times and for all driving
frequencies. This is in contrast to most other formalisms, like NEGF,
Landauer-Büttiker formalism, where it is difficult to calculate dynamics at
all times and for all driving frequencies. This means, for instance, that we
can compute the entropy production rate at all times, and not only the average
over a cycle. In the case of the resonant-level, since the correlation matrix
reduces to only a number, the entropy of the system is given by
$S_{\text{S}}(t)=-p_{t}\log p_{t}-(1-p_{t})\log(1-p_{t})$, where
$p_{t}=\langle c^{\dagger}c\rangle$ is the population of the dot at time $t$.
Therefore, by simply differentiating this function with respect to time, we
can compute the instantaneous entropy production rate [c.f. Eq. (10)] as
$\dot{\Sigma}(t)=\dot{S}_{\text{S}}(t)-\sum_{\alpha}\beta_{\alpha}J_{\alpha}^{Q}$.
This is shown in Fig. 3 (c), where convergence to chain-mapping results up to
a finite time is also shown.

The convergence of results obtained from mesoscopic leads to those obtained by
chain-mapping also demonstrates the fact these results are independent of the
particular microscopic modelling of the baths. All microscopic models of the
baths which give the same spectral functions lead to same dynamics and
thermodynamics at all times. Using the proper definition of currents, given in
Eqs.(32) and (36), is crucial for this. If, instead, we used Eq. (37), which
is tempting to use given the Lindblad description of mesoscopic leads, the
converged result from mesoscopic leads approach would be different, as shown
in Fig. 4. Therefore, they do not agree with results from the chain-mapping
approach. The expressions in Eq. (37) give currents from the residual baths to
the lead modes (see Fig. 2) which depend on the particular microscopic
modelling on the baths.

Next, as an important check, we compare the quantities computed for the driven
resonant-level model against their classical counterpart in the high-
temperature. We consider the same setup as before, with exactly the same
parameters except for the temperature of the two baths. In the limit
$T\gg\Gamma$, the population of the dot is well described by a classical Pauli
master equation of the form

$\derivative{p_{t}}{t}=\sum_{\alpha=\text{L,R}}\big{[}-\Gamma(1-f_{\alpha}(t))p_{t}+\Gamma
f_{\alpha}(t)(1-p_{t})\big{]},$ (63)

with a time-dependent Fermi distribution
$f_{\alpha}(t)=[e^{(\epsilon(t)-\mu_{\alpha})/T}+1]^{-1}$. In Fig. 5, we show
the comparison of the entropy of the dot computed using our approach against
the results obtained by numerically integrating the classical master equation
[Eq. (63)] for increasing values of $T$. The plots clearly show the
convergence of the results with increasing temperature, providing another
benchmark for the consistency of our approach.

Having benchmarked our approach in a prototype example, in the next section,
we proceed to a highly non-trivial case.

Figure 6: Schematic of the rectification model. Two double dots are coupled to
one another with coupling strength $\lambda$. Here the left dot is driven by a
time dependent, oscillatory force at a frequency $\omega$. We consider the
difference in the particle current when the hot and cold baths are swapped.

##  V Entropy cost of energy rectification in a driven non-interacting system

Motivated by Ref. Riera-Campeny  _et al._ (2019), we study heat rectification
in a periodically driven non-interacting system. There, the authors studied
the problem for bosonic degrees of freedom. Here, we consider instead a chain
of two fermionic modes, in which only the first mode is driven. A similar
problem was also considered in Ref. Kohler  _et al._ (2004); Chen  _et al._
(2014), but focusing only on particle currents. With our framework, we can now
also extend this to energy current, and hence heat rectification Li  _et al._
(2012).

The Hamiltonian of the system is given by

$\hat{H}_{\text{S}}(t)=\epsilon_{1}(t)\hat{c}^{\dagger}_{1}\hat{c}_{1}+\epsilon_{2}\hat{c}^{\dagger}_{2}\hat{c}_{2}+\lambda\quantity(\hat{c}^{\dagger}_{1}\hat{c}_{2}+\hat{c}^{\dagger}_{2}\hat{c}_{1}),$
(64)

where $\lambda$ is the coupling strength between the two modes, $\epsilon_{2}$
is the energy of the second and $\epsilon_{1}(t)$ is the energy of the first
dot, which is controlled by periodic protocol of the form
$\epsilon_{1}(t)=\epsilon_{1}(0)+A\sin(\omega t)$. In our simulations, we
always fixed $\epsilon_{1}(0)=\epsilon_{2}=0$. Additionally, each mode is
coupled to its own bath. By simplicity, we take both baths to be described by
the same flat spectral density, given by Eq. (61). Again, we fix $\Gamma=1$
for both baths, and use this parameter to set the energy scales. We also
choose the chemical potentials of both leads to be zero, so that the heat
current matches the energy current.

In order to study heat rectification, two configurations of the setup are
considered. In the forward configuration, we set $T_{\text{L}}=T_{\text{hot}}$
and $T_{\text{R}}=T_{\text{cold}}$, where $T_{\text{hot}}>T_{\text{cold}}$. In
the reverse setup, the temperatures are swapped, so
$T_{\text{L}}=T_{\text{cold}}$ and $T_{\text{R}}=T_{\text{hot}}$ (see Fig. 6).
Since only the first dot is driven, the system is not symmetric with respect
to left-right inversion. As a consequence, the currents in the two scenarios
will generally be different.

Figure 7: Energy current as a function of time in the forward configuration
(blue curve) and opposite of the energy current as a function of time in the
reverse configuration (orange curve) for different values of driving
frequency: (a) $\omega=3$ (b) $\omega=4$ (c) $\omega=5$ (d) $\omega=6$. The
dashed lines indicate the average current in the limit cycle. The coupling
strength between the two modes if fixed to $\lambda=3$ and the other
parameters are $T_{\text{hot}}=2$, $T_{\text{cold}}=1$, $A=1$, $\Gamma=1$,
$\mu_{\text{L}}=\mu_{\text{R}}=0$, $W=8$, $W^{*}=4$ and $L=100$ with
$L_{\text{lin}}/L_{\text{log}}=0.1$ Figure 8: (a) Rectification coefficient as
a function of the driving frequency divided by the coupling strength $\lambda$
between the two modes $\lambda$, for increasing values of $\lambda$. (b)
Cycle-averaged entropy production rate as a function of the driving frequency
divided by the coupling strength between the two modes $\lambda$, for
increasing values of $\lambda$. The parameters used are $T_{\text{hot}}=2$,
$T_{\text{cold}}=1$, $A=1$, $\Gamma=1$, $\mu_{\text{L}}=\mu_{\text{R}}=0$,
$W=8$, $W^{*}=4$ and $L=100$ with $L_{\text{lin}}/L_{\text{log}}=0.1$

Since the system is time-dependent, one must be careful with the notion of
rectification, since the heat flowing through the system is ill-defined. Due
to the presence of the driving, the magnitude of the time period averaged
energy currents associated with the left and right bath are not the same, even
in the long-time limit. In fact, they differ precisely by the time period
averaged power associated with the driving. For this reason, we shall focus
only on the energy current flowing towards the left mode, which is the one
being driven. As a particular example, let us choose $\lambda=3$. We use the
method described in Sec. III.3 to plot the energy currents in the forward and
reversed configuration as a function of time for increasing values of
$\omega$, which is shown in Fig. 7. In our simulations, the amplitude of the
driving is $A$ and the temperatures are $T_{\text{hot}}=2$ and
$T_{\text{cold}}=1$. We fixed the number of modes in both leads to be $L=100$,
which was checked to be enough to ensure proper convergence of the
simulations. We denote by $J^{E}_{L}$ the energy current from the left bath in
the forward configuration, while by $-\tilde{J}^{E}_{L}$ energy current into
the left bath in the backward configuration (see Fig. 6). Both currents
eventually go into a limit cycle show periodic oscillations with the same
frequency as the drive. The horizontal dashed lines correspond to the cycle-
average values in the limit cycle, computed using the method of Sec. III.4.
The plots clearly show that the forward and backward currents are not the same
in the limit cycle. Instead, in the the gap them depends on the driving
frequency. We have found that analogous rectification does not occur in the
particle current for our chosen set-up.

Another interesting point to note is that, in Fig. 7(d), by our sign
conventions, the cycle-averaged heat flows into the left bath in both forward
and backward configurations. In forward configuration the left bath is hot
(see Fig. 6). So, it shows that, due to the relatively high frequency of the
drive, the hot bath is getting heated in this configuration.

In order to better study the frequency dependence of rectification, we define
the rectification coefficient as

$R=\frac{\absolutevalue{\overline{J^{E}_{L}}+\overline{\tilde{J}^{E}_{L}}}}{\absolutevalue{\overline{J^{E}_{L}}-\overline{\tilde{J}^{E}_{L}}}},$
(65)

where $\overline{J^{E}_{L}}$ ($\overline{\tilde{J}^{E}_{L}}$) is the cycle-
averaged energy current in the forward (backward) configuration in the long-
time limit. In the case of a perfect diode, i.e., where there is some current
flowing in one configuration but none in the other, we have that $R=1$. If
there is no rectification, i.e., both currents have the same magnitude but
opposite signs, then $R=0$. Interestingly, this quantity, which is analogous
to that used to characterize rectification in absence of external driving, is
not upper bounded by $1$. Due to external driving, it can happen that the
direction of cycle-averaged heat current is the same in both configurations,
as we saw in Fig. 7(d). In such cases, we will have $R>1$.

Using the method described in Sec. III.4, we now compute the limit-cycle
solution for a wide range of couplings $\lambda$ and driving frequencies
$\omega$, and look at both the rectification coefficient and the entropy
production. In Fig. 8 (a) we show the rectification coefficient as a function
of $\omega/\lambda$, for increasing values of $\lambda$. As can be seen, in
the low and high frequency regimes the rectification coefficient goes to zero,
but for intermediate values it peaks. For $\lambda>1$, this peak occurs in the
same value of $\omega/\lambda$.

In the limit cycle, the cycle-averaged entropy production rate is given simply
by $\overline{\sigma}=\sum_{\alpha}\beta_{\alpha}\overline{J}^{Q}_{\alpha}$,
where $\overline{J}^{Q}_{\alpha}$ denotes the cycle-averaged heat current
associated with each bath. Interestingly, we have found that, even though
there is rectification of energy current flowing into the left bath, the
cycle-averaged entropy production rate, i.e., $\overline{\sigma}$ is the same
in both the forward and reversed configurations. In Fig. 8 (b) we plot
$\overline{\sigma}$ as a function of $\omega/\lambda$. We find that there is a
strong peak in the entropy production rate at the same frequency in which the
maximum rectification coefficient occurs. This clearly shows that an increase
in rectification of energy current in our set-up comes at the cost of an
increase in the entropy production rate.

We will like to mention that plots in Fig. 8, which give the entropy cost of
energy rectification as a function of driving frequency in our set-up, would
be difficult to obtain via most other techniques. This is far beyond Markovian
regime of system dynamics and additionally does not have any small parameter
in the Hamiltonian to allow perturbation techniques. The frequency of the
drive is also a free parameter, which has been varied from low to a high
value, ruling out any possibility of perturbation in frequency also. This
clearly highlights the power of our approach.

Our results in this section open many interesting questions, even for this
rather simple set-up. The various interesting features, like the dependence of
the position of the peak on various system parameters, the temperature
dependence of rectification, the effect of having additional chemical
potential bias, all can now be explored completely non-perturbatively. These
directions deserve thorough separate investigations, which we delegate to
future works.

##  VI Conclusion

In this work we have performed a detailed analysis of the mescoscopic leads in
the context of quantum thermodynamics. Our main contribution is that the
methodology allows for the inclusion of an arbitrary time dependence on the
system Hamiltonian. This is a formidable challenge in quantum thermodynamics
which we have demonstrated this technique can overcome. Following an overview
of the thermodynamics of driven open quantum systems we reviewed the
mesoscopic leads formalism and show how the definitions of the energy and
particle current behave in the presence of driving. We emphasis that this
technique has the formidable feature of being able to cope with both strong
coupling to the leads and fast driving. We have demonstrated this focusing in
particular on quadratic systems however an extension to interacting central
systems by means of tensor networks is possible. We believe that the power of
the methodology will allow for extensive explorations of thermodynamics of
quantum systems in regimes that have so far been inaccessible with other
approaches. To demonstrate the power of the methodology we apply our method to
both the driven resonant level model and driven tunnel coupled quantum dot
models respectively. In the case of the driven resonant level we were able to
show that our results replicate known results in the literature in the high
temperature limit. In the case of the driven quantum dot we showed how driving
of one of the dot’s energy can be used to induce a powerful heat rectification
effect in a parameter regime which would be inaccessible to conventional
techniques. In future work we plan to extend this technique further to access
higher moments of the currents both in the presence and absence of driving as
well as incorporating non quadratic interactions in the central system.

Acknowledgements. The authors acknowledge Mark T. Mitchison for fruitful
discussions. This work was funded by the European Research Council Starting
Grant ODYSSEY (Grant Agreement No. 758403) and the EPSRC-SFI joint project
QuamNESS. J. G. is supported by a SFI-Royal Society University Research
Fellowship. and GTL acknowledges the financial support of the São Paulo
Funding Agency FAPESP (Grant No. 2019/14072-0.), and the Brazilian funding
agency CNPq (Grant No. INCT-IQ 246569/2014-0). A.P acknowledges funding from
the European Union’s Horizon 2020 research and innovation programme under the
Marie Sklodowska-Curie Grant Agreement No. 890884. A.P also acknowledges
funding from the Danish National Research Foundation through the Center of
Excellence “CCQ” (Grant agreement no.: DNRF156). J.G. would aslo like to thank
Rosario Fazio for suggestions regarding rectification in the presence of
driving and M. T. Mitchison for general discussions on the project.

## Appendix

##  Appendix A Brute force numerics with chain-mapping

The systematic way to implement the brute force simulation of unitary dynamics
of the system and baths using finite, but large enough, baths is offered by
the chain mapping of the baths. Any bath spectral function, with finite upper
and lower cut-offs, can be mapped into a semi-infinite one-dimensional nearest
neighbour tight-binding chain with only the first site of the chain coupled
with the system Nazir and Schaller (2018); de Vega  _et al._ (2015); Chin  _et
al._ (2010); Strasberg  _et al._ (2018); Prior  _et al._ (2010); Woods  _et
al._ (2015); Woods and Plenio (2016); Mascherpa  _et al._ (2017); Garg  _et
al._ (1985),

$\displaystyle\hat{H}_{\alpha}=\sum_{p=1}^{L_{B}}\varepsilon_{p,\alpha}\hat{a}_{p,\alpha}^{\dagger}\hat{a}_{p,\alpha}+g_{p,\alpha}(\hat{a}_{p,\alpha}^{\dagger}\hat{a}_{p+1,\alpha}+\hat{a}_{p+1,\alpha}^{\dagger}\hat{a}_{p\alpha}),$
$\displaystyle\hat{H}_{S\alpha}=g_{0,\alpha}(S_{\alpha}^{\dagger}\hat{a}_{1,\alpha}+\hat{a}_{1,\alpha}^{\dagger}S_{\alpha}),$
(66)

with $L_{B}\to\infty$. Such a chain gives the spectral function
$\mathcal{J}_{\alpha}(\omega)$ if the on-site potentials
$\varepsilon_{p,\alpha}$ and the hoppings $g_{p,\alpha}$ are obtained from the
following set of recursion relations

$\displaystyle g_{p,\alpha}^{2}=\frac{1}{2\pi}\int
d\omega\mathfrak{J}_{p,\alpha}(\omega),$
$\displaystyle\varepsilon_{p+1,\alpha}=\frac{1}{2\pi g_{p,\alpha}^{2}}\int
d\omega\leavevmode\nobreak\ \omega\mathfrak{J}_{p,\alpha}(\omega),$
$\displaystyle\mathfrak{J}_{p+1,\alpha}(\omega)=\frac{4g_{p,\alpha}^{2}\mathfrak{J}_{p,\alpha}(\omega)}{\left[\mathfrak{J}_{p,\alpha}^{H}(\omega)\right]^{2}+\left[\mathfrak{J}_{p,\alpha}(\omega)\right]^{2}},\leavevmode\nobreak\
\leavevmode\nobreak\ $ (67)

with $p$ going from $0$ to $L_{B}$,
$\mathfrak{J}_{0,\alpha}(\omega)=\mathfrak{J}_{\alpha}(\omega)$ and
$\mathfrak{J}_{p,\alpha}^{H}(\omega)$ being the Hilbert transform of
$\mathfrak{J}_{p,\alpha}(\omega)$,

$\displaystyle\mathfrak{J}_{p,\alpha}^{H}(\omega)=\frac{1}{\pi}\mathcal{P}\int_{-\infty}^{\infty}d\omega^{\prime}\frac{\mathfrak{J}_{p,\alpha}(\omega^{\prime})}{\omega-\omega^{\prime}},$
(68)

where $\mathcal{P}$ denotes the principal value. With finite high and low
frequency cut-offs, the parameters $\varepsilon_{p,\alpha}$ and $g_{p,\alpha}$
quickly tend to constants for increasing $p$. Let these constants be
$\varepsilon_{B_{\alpha}}$ and $g_{B_{\alpha}}$. With baths modelled as such
chains, the process to be simulated involves switching on the system-bath
couplings at initial time. Due to Lieb-Robinson bounds, the information about
this spreads at a finite speed proportional to $g_{B_{\alpha}}$. Consequently,
to simulate up to a time $t$, baths of size

$\displaystyle L_{B}\sim g_{B,\alpha}t$ (69)

suffices to accurately mimic the limit $L_{B}\to\infty$. We see that to
simulate accurately up to a longer time, a larger bath size is required.

## References

  * Callen (1998) H. B. Callen, “Thermodynamics and an introduction to thermostatistics,” (1998).
  * De Groot and Mazur (2013) S. R. De Groot and P. Mazur,  _Non-equilibrium thermodynamics_ (Courier Corporation, 2013).
  * Jarzynski (2011) C. Jarzynski, Annual Review of Condensed Matter Physics 2, 329 (2011).
  * Seifert (2012) U. Seifert, Reports on Progress in Physics 75, 126001 (2012).
  * Seifert (2019) U. Seifert, Annual Review of Condensed Matter Physics 10, 171 (2019).
  * Monroe and Kim (2013) C. Monroe and J. Kim, Science 339, 1164 (2013).
  * Aradhya and Venkataraman (2013) S. V. Aradhya and L. Venkataraman, Nature Nanotechnology 8, 399 (2013).
  * Kagan  _et al._ (2021) C. R. Kagan, L. C. Bassett, C. B. Murray, and S. M. Thompson, Chemical Reviews 121, 3186 (2021).
  * Burkard  _et al._ (2020) G. Burkard, M. J. Gullans, X. Mi, and J. R. Petta, Nature Reviews Physics 2, 129 (2020).
  * Clerk  _et al._ (2020) A. A. Clerk, K. W. Lehnert, P. Bertet, J. R. Petta, and Y. Nakamura, Nature Physics 16, 257 (2020).
  * Doherty  _et al._ (2013) M. W. Doherty, N. B. Manson, P. Delaney, F. Jelezko, J. Wrachtrup, and L. C. Hollenberg, Physics Reports 528, 1 (2013), the nitrogen-vacancy colour centre in diamond.
  * Arrachea (2022) L. Arrachea, arXiv e-prints , arXiv:2205.14200 (2022), arXiv:2205.14200 [quant-ph] .
  * Auffèves (2021) A. Auffèves, arXiv e-prints , arXiv:2111.09241 (2021), arXiv:2111.09241 [quant-ph] .
  * Liu  _et al._ (2015) Y.-Y. Liu, J. Stehlik, C. Eichler, M. J. Gullans, J. M. Taylor, and J. R. Petta, Science 347, 285 (2015).
  * Stehlik  _et al._ (2016) J. Stehlik, Y.-Y. Liu, C. Eichler, T. R. Hartke, X. Mi, M. J. Gullans, J. M. Taylor, and J. R. Petta, Phys. Rev. X 6, 041027 (2016).
  * Wen  _et al._ (2020) Y. Wen, N. Ares, F. J. Schupp, T. Pei, G. A. D. Briggs, and E. A. Laird, Nature Physics 16, 75 (2020).
  * Koch  _et al._ (2019) C. P. Koch, M. Lemeshko, and D. Sugny, Rev. Mod. Phys. 91, 035005 (2019).
  * James (2021) M. James, Annual Review of Control, Robotics, and Autonomous Systems 4, 343 (2021).
  * Bhandari  _et al._ (2020) B. Bhandari, P. T. Alonso, F. Taddei, F. von Oppen, R. Fazio, and L. Arrachea, Phys. Rev. B 102, 155407 (2020).
  * Imamoglu (1994) A. Imamoglu, Physical Review A 50, 3650 (1994).
  * Garraway (1997a) B. M. Garraway, Physical Review A 55, 2290 (1997a).
  * Garraway (1997b) B. M. Garraway, Physical Review A 55, 4636 (1997b).
  * Sánchez  _et al._ (2006) C. G. Sánchez, M. Stamenova, S. Sanvito, D. R. Bowler, A. P. Horsfield, and T. N. Todorov, The Journal of Chemical Physics 124, 214708 (2006).
  * Subotnik  _et al._ (2009) J. E. Subotnik, T. Hansen, M. A. Ratner, and A. Nitzan, The Journal of Chemical Physics 130, 144105 (2009).
  * Dzhioev and Kosov (2011) A. A. Dzhioev and D. S. Kosov, The Journal of Chemical Physics 134, 044121 (2011).
  * Ajisaka  _et al._ (2012) S. Ajisaka, F. Barra, C. Mejía-Monasterio, and T. Prosen, Physical Review B 86 (2012), 10.1103/physrevb.86.125111.
  * Ajisaka and Barra (2013) S. Ajisaka and F. Barra, Physical Review B 87 (2013), 10.1103/physrevb.87.195114.
  * Arrigoni  _et al._ (2013) E. Arrigoni, M. Knap, and W. von der Linden, Phys. Rev. Lett. 110, 086403 (2013).
  * Dorda  _et al._ (2014) A. Dorda, M. Nuss, W. von der Linden, and E. Arrigoni, Phys. Rev. B 89, 165105 (2014).
  * Chen  _et al._ (2014) L. Chen, T. Hansen, and I. Franco, The Journal of Physical Chemistry C 118, 20009 (2014).
  * Zelovich  _et al._ (2014) T. Zelovich, L. Kronik, and O. Hod, Journal of Chemical Theory and Computation 10, 2927 (2014), pMID: 26588268.
  * Dorda  _et al._ (2015) A. Dorda, M. Ganahl, H. G. Evertz, W. von der Linden, and E. Arrigoni, Phys. Rev. B 92, 125145 (2015).
  * Hod  _et al._ (2016) O. Hod, C. A. Rodríguez-Rosario, T. Zelovich, and T. Frauenheim, The Journal of Physical Chemistry A 120, 3278 (2016), pMID: 26807992.
  * Gruss  _et al._ (2016) D. Gruss, K. A. Velizhanin, and M. Zwolak, Scientific Reports 6 (2016), 10.1038/srep24514.
  * Schwarz  _et al._ (2016) F. Schwarz, M. Goldstein, A. Dorda, E. Arrigoni, A. Weichselbaum, and J. von Delft, Phys. Rev. B 94, 155142 (2016).
  * Dorda  _et al._ (2017) A. Dorda, M. Sorantin, W. von der Linden, and E. Arrigoni, New Journal of Physics 19, 063005 (2017).
  * Elenewski  _et al._ (2017) J. E. Elenewski, D. Gruss, and M. Zwolak, The Journal of Chemical Physics 147, 151101 (2017).
  * Gruss  _et al._ (2017) D. Gruss, A. Smolyanitsky, and M. Zwolak, The Journal of Chemical Physics 147, 141102 (2017).
  * Zelovich  _et al._ (2017) T. Zelovich, T. Hansen, Z.-F. Liu, J. B. Neaton, L. Kronik, and O. Hod, The Journal of Chemical Physics 146, 092331 (2017).
  * Tamascelli  _et al._ (2018) D. Tamascelli, A. Smirne, S. Huelga, and M. Plenio, Physical Review Letters 120 (2018), 10.1103/physrevlett.120.030402.
  * Lemmer  _et al._ (2018) A. Lemmer, C. Cormick, D. Tamascelli, T. Schaetz, S. F. Huelga, and M. B. Plenio, New Journal of Physics 20, 073002 (2018).
  * Oz  _et al._ (2019) A. Oz, O. Hod, and A. Nitzan, Journal of Chemical Theory and Computation 16, 1232 (2019).
  * Chen  _et al._ (2019) F. Chen, E. Arrigoni, and M. Galperin, New Journal of Physics 21, 123035 (2019).
  * Zwolak (2020) M. Zwolak, The Journal of Chemical Physics 153, 224107 (2020).
  * Wójtowicz  _et al._ (2020) G. Wójtowicz, J. E. Elenewski, M. M. Rams, and M. Zwolak, Physical Review A 101 (2020), 10.1103/physreva.101.050301.
  * Chiang and Hsu (2020) T.-M. Chiang and L.-Y. Hsu, The Journal of Chemical Physics 153, 044103 (2020).
  * Brenes  _et al._ (2020) M. Brenes, J. J. Mendoza-Arenas, A. Purkayastha, M. T. Mitchison, S. R. Clark, and J. Goold, Phys. Rev. X 10, 031040 (2020).
  * Lotem  _et al._ (2020) M. Lotem, A. Weichselbaum, J. von Delft, and M. Goldstein, Phys. Rev. Research 2, 043052 (2020).
  * Fugger  _et al._ (2020) D. M. Fugger, D. Bauernfeind, M. E. Sorantin, and E. Arrigoni, Phys. Rev. B 101, 165132 (2020).
  * Wójtowicz  _et al._ (2021) G. Wójtowicz, J. E. Elenewski, M. M. Rams, and M. Zwolak, Phys. Rev. B 104, 165131 (2021).
  * Elenewski  _et al._ (2021) J. E. Elenewski, G. Wójtowicz, M. M. Rams, and M. Zwolak, The Journal of Chemical Physics 155, 124117 (2021).
  * Breuer and Petruccione (2007) H.-P. Breuer and F. Petruccione,  _The Theory of Open Quantum Systems_ (Oxford University Press, Oxford, 2007).
  * Rivas and Huelga (2012) A. Rivas and S. F. Huelga,  _Open Quantum Systems_ (Springer Berlin Heidelberg, 2012).
  * Walls (1970) D. F. Walls, Zeitschrift für Physik A Hadrons and nuclei 234, 231 (1970).
  * Novotný (2002) T. Novotný, Europhysics Letters (EPL) 59, 648 (2002).
  * Wichterich  _et al._ (2007) H. Wichterich, M. J. Henrich, H.-P. Breuer, J. Gemmer, and M. Michel, Phys. Rev. E 76, 031115 (2007).
  * Prachar and Novotný (2010) J. Prachar and T. Novotný, Physica E: Low-dimensional Systems and Nanostructures 42, 565 (2010).
  * Rivas  _et al._ (2010) Á. Rivas, A. D. K. Plato, S. F. Huelga, and M. B. Plenio, New Journal of Physics 12, 113032 (2010).
  * Deçordi and Vidiella-Barranco (2017) G. Deçordi and A. Vidiella-Barranco, Optics Communications 387, 366 (2017).
  * Levy and Kosloff (2014) A. Levy and R. Kosloff, EPL (Europhysics Letters) 107, 20004 (2014).
  * Purkayastha  _et al._ (2016) A. Purkayastha, A. Dhar, and M. Kulkarni, Phys. Rev. A 93, 062114 (2016).
  * Eastham  _et al._ (2016) P. R. Eastham, P. Kirton, H. M. Cammack, B. W. Lovett, and J. Keeling, Physical Review A 94 (2016).
  * Hofer  _et al._ (2017) P. P. Hofer, M. Perarnau-Llobet, L. D. M. Miranda, G. Haack, R. Silva, J. B. Brask, and N. Brunner, New Journal of Physics 19, 123037 (2017).
  * González  _et al._ (2017) J. O. González, L. A. Correa, G. Nocerino, J. P. Palao, D. Alonso, and G. Adesso, Open Systems & Information Dynamics 24, 1740010 (2017).
  * Mitchison and Plenio (2018) M. T. Mitchison and M. B. Plenio, New Journal of Physics 20, 033005 (2018).
  * Cattaneo  _et al._ (2019) M. Cattaneo, G. L. Giorgi, S. Maniscalco, and R. Zambrini, New Journal of Physics 21, 113045 (2019).
  * Hartmann and Strunz (2020) R. Hartmann and W. T. Strunz, Phys. Rev. A 101, 012103 (2020).
  * Konopik and Lutz (2022) M. Konopik and E. Lutz, Phys. Rev. Research 4, 013171 (2022).
  * Scali  _et al._ (2021) S. Scali, J. Anders, and L. A. Correa, Quantum 5, 451 (2021).
  * Tupkary  _et al._ (2022) D. Tupkary, A. Dhar, M. Kulkarni, and A. Purkayastha, Phys. Rev. A 105, 032208 (2022).
  * Schnell  _et al._ (2020) A. Schnell, A. Eckardt, and S. Denisov, Phys. Rev. B 101, 100301 (2020).
  * Schnell  _et al._ (2021) A. Schnell, S. Denisov, and A. Eckardt, Phys. Rev. B 104, 165414 (2021).
  * Haug and Jauho (2008) H. Haug and A.-P. Jauho,  _Quantum Kinetics in Transport and Optics of Semiconductors_ (Springer-Verlag Berlin Heidelberg, 2008).
  * Kamenev (2011) A. Kamenev,  _Field Theory of Non-Equilibrium Systems_ (Cambridge University Press, Cambridge, 2011).
  * Cortés  _et al._ (1985) E. Cortés, B. J. West, and K. Lindenberg, The Journal of Chemical Physics 82, 2708 (1985).
  * Ford  _et al._ (1988) G. W. Ford, J. T. Lewis, and R. F. O’Connell, Phys. Rev. A 37, 4419 (1988).
  * Benguria and Kac (1981) R. Benguria and M. Kac, Phys. Rev. Lett. 46, 1 (1981).
  * Dhar and Sriram Shastry (2003) A. Dhar and B. Sriram Shastry, Phys. Rev. B 67, 195405 (2003).
  * Dhar and Roy (2006) A. Dhar and D. Roy, Journal of Statistical Physics 125, 801 (2006).
  * Moskalets (2011) M. V. Moskalets,  _Scattering Matrix Approach to Non-Stationary Quantum Transport_ (IMPERIAL COLLEGE PRESS, 2011).
  * Datta  _et al._ (1987) S. Datta, M. Cahay, and M. McLennan, Phys. Rev. B 36, 5655 (1987).
  * Arrachea and Moskalets (2006) L. Arrachea and M. Moskalets, Phys. Rev. B 74, 245322 (2006).
  * Lesovik and Sadovskyy (2011) G. B. Lesovik and I. A. Sadovskyy, Physics-Uspekhi 54, 1007 (2011).
  * Brandner (2020) K. Brandner, Zeitschrift für Naturforschung A 75, 483 (2020).
  * Potanina  _et al._ (2021) E. Potanina, C. Flindt, M. Moskalets, and K. Brandner, Phys. Rev. X 11, 021013 (2021).
  * Ludovico  _et al._ (2016a) M. F. Ludovico, M. Moskalets, D. Sánchez, and L. Arrachea, Physical Review B 94 (2016a), 10.1103/physrevb.94.035436.
  * Bruch  _et al._ (2018) A. Bruch, C. Lewenkopf, and F. von Oppen, Physical Review Letters 120 (2018), 10.1103/physrevlett.120.107701.
  * Riera-Campeny  _et al._ (2019) A. Riera-Campeny, M. Mehboudi, M. Pons, and A. Sanpera, Phys. Rev. E 99, 032126 (2019).
  * Esposito  _et al._ (2015a) M. Esposito, M. A. Ochoa, and M. Galperin, Physical Review Letters 114 (2015a), 10.1103/physrevlett.114.080602.
  * Brandner and Seifert (2016) K. Brandner and U. Seifert, Phys. Rev. E 93, 062134 (2016).
  * Dann and Kosloff (2021) R. Dann and R. Kosloff, Quantum 5, 590 (2021).
  * Esposito  _et al._ (2015b) M. Esposito, M. A. Ochoa, and M. Galperin, Phys. Rev. B 92, 235440 (2015b).
  * Ludovico  _et al._ (2016b) M. F. Ludovico, L. Arrachea, M. Moskalets, and D. Sánchez, Entropy 2016, Vol. 18, Page 419 18, 419 (2016b), 1610.05143 .
  * Bruch  _et al._ (2016) A. Bruch, M. Thomas, S. Viola Kusminskiy, F. von Oppen, and A. Nitzan, Phys. Rev. B 93, 115318 (2016).
  * Katz and Kosloff (2016) G. Katz and R. Kosloff, Entropy 18, 186 (2016).
  * Haughian  _et al._ (2018) P. Haughian, M. Esposito, and T. L. Schmidt, Physical Review B 97 (2018), 10.1103/physrevb.97.085435.
  * Esposito  _et al._ (2010) M. Esposito, K. Lindenberg, and C. V. den Broeck, New Journal of Physics 12, 013013 (2010).
  * Reeb and Wolf (2014) D. Reeb and M. M. Wolf, New Journal of Physics 16, 103011 (2014).
  * Landi and Paternostro (2021) G. T. Landi and M. Paternostro, Reviews of Modern Physics 93 (2021), 10.1103/revmodphys.93.035008.
  * Strasberg and Winter (2021) P. Strasberg and A. Winter, PRX Quantum 2, 030202 (2021).
  * Strasberg  _et al._ (2021) P. Strasberg, M. G. Díaz, and A. Riera-Campeny, Phys. Rev. E 104, L022103 (2021).
  * Talkner and Hänggi (2020) P. Talkner and P. Hänggi, Rev. Mod. Phys. 92, 041002 (2020).
  * Nazir and Schaller (2018) A. Nazir and G. Schaller, “The Reaction Coordinate Mapping in Quantum Thermodynamics,” in  _Thermodynamics in the Quantum Regime: Fundamental Aspects and New Directions_, Vol. 195, edited by F. Binder, L. A. Correa, C. Gogolin, J. Anders, and G. Adesso (2018) p. 551.
  * de Vega  _et al._ (2015) I. de Vega, U. Schollwöck, and F. A. Wolf, Phys. Rev. B 92, 155126 (2015).
  * Chin  _et al._ (2010) A. W. Chin, A. Rivas, S. F. Huelga, and M. B. Plenio, Journal of Mathematical Physics 51, 092109 (2010).
  * Strasberg  _et al._ (2018) P. Strasberg, G. Schaller, T. L. Schmidt, and M. Esposito, Phys. Rev. B 97, 205405 (2018).
  * Prior  _et al._ (2010) J. Prior, A. W. Chin, S. F. Huelga, and M. B. Plenio, Phys. Rev. Lett. 105, 050404 (2010).
  * Woods  _et al._ (2015) M. P. Woods, M. Cramer, and M. B. Plenio, Phys. Rev. Lett. 115, 130401 (2015).
  * Woods and Plenio (2016) M. P. Woods and M. B. Plenio, Journal of Mathematical Physics 57, 022105 (2016), https://doi.org/10.1063/1.4940436 .
  * Mascherpa  _et al._ (2017) F. Mascherpa, A. Smirne, S. F. Huelga, and M. B. Plenio, Phys. Rev. Lett. 118, 100401 (2017).
  * Garg  _et al._ (1985) A. Garg, J. N. Onuchic, and V. Ambegaokar, The Journal of Chemical Physics 83, 4491 (1985).
  * Landi  _et al._ (2021) G. T. Landi, D. Poletti, and G. Schaller, (2021), arXiv:2104.14350v3 [quant-ph] .
  * Golub and Loan (1996) G. H. Golub and C. F. V. Loan,  _Matrix Computations_ , 3rd ed. (Johns Hopkins University Press, 1996).
  * Weichselbaum  _et al._ (2009) A. Weichselbaum, F. Verstraete, U. Schollwöck, J. I. Cirac, and J. von Delft, Phys. Rev. B 80, 165117 (2009).
  * Kohler  _et al._ (2004) S. Kohler, S. Camalet, M. Strass, J. Lehmann, G.-L. Ingold, and P. Hänggi, Chemical Physics 296, 243 (2004).
  * Li  _et al._ (2012) N. Li, J. Ren, L. Wang, G. Zhang, P. Hänggi, and B. Li, Reviews of Modern Physics 84, 1045 (2012).