# Lyapunov equation in open quantum systems and non-Hermitian physics

Archak Purkayastha<EMAIL_ADDRESS>School of Physics, Trinity College
Dublin, College Green, Dublin 2, Ireland  Centre for complex quantum systems,
Aarhus University, Nordre Ringgade 1, 8000 Aarhus C, Denmark

###### Abstract

The continuous-time differential Lyapunov equation is widely used in linear
optimal control theory, a branch of mathematics and engineering. In quantum
physics, it is known to appear in Markovian descriptions of linear (quadratic
Hamiltonian, linear equations of motion) open quantum systems, typically from
quantum master equations. Despite this, the Lyapunov equation is seldom
considered a fundamental formalism for linear open quantum systems. In this
work we aim to change that. We establish the Lyapunov equation as a
fundamental and efficient formalism for linear open quantum systems that can
go beyond the limitations of various standard quantum master equation
descriptions, while remaining of much less complexity than general exact
formalisms. This also provides valuable insights for non-Hermitian quantum
physics. In particular, we derive the Lyapunov equation for the most general
number conserving linear system in a lattice of arbitrary dimension and
geometry, connected to an arbitrary number of baths at different temperatures
and chemical potentials. Three slightly different forms of the Lyapunov
equation are derived via an equation of motion approach, by making increasing
levels of controlled approximations, without reference to any quantum master
equation. Then we discuss their relation with quantum master equations,
positivity, accuracy and additivity issues, the possibility of describing dark
states, general perturbative solutions in terms of single-particle
eigenvectors and eigenvalues of the system, and quantum regression formulas.
Our derivation gives a clear understanding of the origin of the non-Hermitian
Hamiltonian describing the dynamics and separates it from the effects of
quantum and thermal fluctuations. Many of these results would have been hard
to obtain via standard quantum master equation approaches.

##  I Introduction

General background — As we go towards the age of quantum technology, it has
become crucial to formulate the theory to describe microscopic (i.e,
containing finite number of degrees of freedom) quantum systems coupled to
multiple macroscopic (i.e, containing infinite number of degrees of freedom)
thermal environments (baths) which can all be at different temperatures and
chemical potentials. This is relevant across quantum optics Hur  _et al._
(2016), thermodynamics Binder  _et al._ (2018), chemistry Cao  _et al._
(2019), biologyCao  _et al._ (2020) and engineering Awschalom  _et al._
(2021). However, this is an extremely challenging problem in general because
the dimension of Hilbert space scales exponentially with the number of degrees
of freedom.

In absence of coupling to any macroscopic bath, the situation can be
considerably simplified if the Hamiltonian is quadratic in fermionic or
bosonic creation and annihilation operators. In that case, many properties of
the system can be obtained by calculating the so-called single-particle
density matrix or the correlation matrix, whose elements are the equal time
two-point correlation functions of the system Mahan (2000). The dynamics of
the correlation matrix is governed by the so-called single-particle
Hamiltonian of the system, whose dimension scales linearly with the number of
degrees of freedom in the system. This therefore leads to an exponential
simplification of the problem.

However, in presence of coupling to multiple macroscopic thermal baths,
describing even dynamics of systems governed by quadratic Hamiltonians in
generality becomes quite complicated because of the baths having infinite
degrees of freedom. In this case, standard technique in all open system
formalisms, like non-equilibrium Green’s functions (NEGF) Haug and Jauho
(2008); Wang  _et al._ (2014), Feynman-Vernon influence functional approach
Kamenev (2011), quantum Langevin equations Ford  _et al._ (1965); Benguria and
Kac (1981); Cortés  _et al._ (1985); Ford  _et al._ (1988); Dhar and Sriram
Shastry (2003); Banerjee  _et al._ (2004); Dhar and Sen (2006); Dhar and Roy
(2006) and quantum master equations (QME) Breuer and Petruccione (2007); Rivas
and Huelga (2012), is to derive effective equations for dynamics of the system
after analytically integrating out the baths. For quadratic Hamiltonians, this
can be done in generality, and leads to non-Markovian equations of motion Dhar
and Sen (2006); Dhar and Roy (2006); Tu and Zhang (2008); Jin  _et al._
(2010); Zhang  _et al._ (2012); McCutcheon  _et al._ (2015); Martinez and Paz
(2013); Yang  _et al._ (2013); Ribeiro and Vieira (2015). Due to this non-
Markovianity, for a system with large (but finite) number of degrees of
freedom, obtaining the full dynamics again becomes quite challenging. If the
long-time non-equilibrium steady state (NESS) is unique, and only the NESS
properties are desired, they can be relatively easily found via a Fourier
transform. But obtaining full dynamics requires a Laplace transform that needs
to be finally inverted. Inverting a Laplace transform can become difficult
depending on the number of degrees of freedom and the spectrum of the system.
As a result, to obtain the full dynamics, in many cases it becomes useful to
make further approximations to have a much simpler effective Markovian
description.

Weak-coupling Markovian descriptions and their limitations — The standard
approach to obtain a Markovian description is to assume weak coupling between
the system and the baths and implement the so-called Born-Markov-secular
approximations to obtain a QME in the so-called Lindblad form Breuer and
Petruccione (2007). But several drawbacks of Born-Markov-secular
approximations have been pointed out, in particular, for describing a system
coupled to multiple thermal baths 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); Trushechkin and Volovich (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); Benatti  _et al._ (2020); Konopik
and Lutz (2022); Scali  _et al._ (2021); Benatti  _et al._ (2021); Trushechkin
(2021); Tupkary  _et al._ (2022). It has been shown that not implementing the
secular approximation, thereby obtaining a QME in so-called Redfield form
Redfield (1965), can be more accurate (for example, Purkayastha  _et al._
(2016); Hartmann and Strunz (2020); Tupkary  _et al._ (2022); Davidovic
(2021)). On the other hand, it is known that the Redfield equation violates
complete positivity Davidovic (2021); Purkayastha  _et al._ (2020a); Hartmann
and Strunz (2020); Eastham  _et al._ (2016); Anderloni  _et al._ (2007);
Gaspard and Nagaoka (1999); Kohen  _et al._ (1997); Gnutzmann and Haake
(1996); Suárez  _et al._ (1992). This means it can lead to having unphysical
negative eigenvalues of the system density matrix. But such negative
eigenvalues will be quite close to zero, and can be shown to be below accuracy
level of the Redfield equation Fleming and Cummings (2011); Hartmann and
Strunz (2020); Tupkary  _et al._ (2022); Davidovic (2021). Several more-
refined Lindblad equations have been proposed to circumvent these drawbacks
Nathan and Rudner (2020); Kleinherbers  _et al._ (2020); Davidović (2020);
Mozgunov and Lidar (2020); McCauley  _et al._ (2020); Kiršanskas  _et al._
(2018). Nevertheless, it has been shown that all of these approaches have
inherent limitations even to the leading order in the system-bath coupling
strength when describing a system coupled to multiple thermal baths Tupkary
_et al._ (2022).

The Lyapunov equation from weak-coupling Markovian descriptions — Regardless
of the several limitations, because of their simplicity and the valuable
insight they provide, weak coupling QME descriptions of thermal baths remain
of wide use. Interestingly, for systems governed by time-independent quadratic
Hamiltonians, the weak system-bath coupling QMEs referred to above give an
amazing connection to seemingly unrelated fields in mathematics and
engineering. It turns out that the equation of motion for the correlation
matrix as obtained from such QMEs has the form of a continuous-time Lyapunov
equation (see for example Prosen (2008); Prosen and Žunkovič (2010); Žunkovič
and Prosen (2012); Koga and Yamamoto (2012); Nicacio  _et al._ (2015); Landi
_et al._ (2013); Nicacio  _et al._ (2016); Santos  _et al._ (2017); Malouf
_et al._ (2019); McDonald  _et al._ (2022); Bernal-García  _et al._ (2021);
Roccati  _et al._ (2021)). The Lyapunov equation is an extremely well-studied
equation in mathematics and engineering, where appears extensively in the
field of linear optimal control theory Gajić and Qureshi (1995); Sun and Yong
(2020); Yuan  _et al._ (2021). It is surprising that the same equation appears
in a completely different context, that of linear (i.e, governed by quadratic
Hamiltonians) open quantum systems. Further, efficient numerical methods for
solving Lyapunov equations already exist in all high-level scientific
programming languages, like python, Mathematica, Matlab etc.

Summary of our main results — Although this above fact is already known and
used, Lyapunov equations are not usually discussed as a fundamental formalism
for general linear open quantum systems. This is presumably due to the various
limitations of the underlying QME descriptions. In this manuscript, we attempt
to change that. We establish the continuous-time differential Lyapunov
equation as a rigorously derived efficient description of linear open quantum
systems, that is more fundamental than many of the existing QME descriptions
and can go beyond their limitations.

To this end, we consider a system governed by a number conserving quadratic
time-independent Hamiltonian, bosonic or fermionic, in a lattice of arbitrary
dimension and geometry bilinearly coupled at an arbitrary number of sites to
thermal baths which can all be at different temperatures and chemical
potentials (see Fig. 1). The baths are assumed to be modelled by an infinite
number of bosonic or fermionic modes. By systematically carrying out various
controlled approximations on bath spectral functions and strength of system-
bath couplings, we follow an equation of motion approach to obtain the
continuous-time differential Lyapunov equation without referring to any QME.
In particular, we derive three slightly different Lyapunov equations at three
different increasing levels of approximations. At the first level of
approximation, the Lyapunov equation does not have any positivity problem at
all times, at the second level of approximation, the Lyapunov equation does
not have a positivity problem at the NESS, while at the third level of
approximation, the Lyapunov equation has the same positivity problem as the
Redfield equation. This shows that the Lyapunov equation is more fundamental
for describing such set-ups than the standard QMEs. In fact, the Lyapunov
equations at first two levels suggest corresponding QMEs, which would be
difficult to obtain otherwise. For fermionic systems, at the first two level
of approximations, even the weak system-bath-coupling approximation is not
required, but rather a so-called wide-band-limit approximation on the bath
spectral functions suffices. The controlled microscopic derivation allows us
to specify the validity regime of the Lyapunov equations and the accuracy of
the solutions. A plethora of semi-analytical results follow, some of which
reduces the complexity of the problem of simulating dynamics of this non-
equilibrium open quantum system to the same level as simulating dynamics of
the isolated system. We also give generalized regression formulas for two-time
correlations, which can be easily obtained via our operator equation of motion
approach. These formulas do not correspond to those that would be obtained by
naively applying the standard quantum regression formula Breuer and
Petruccione (2007) at the level of the associated QMEs.

The Lyapunov equation and non-Hermitian quantum physics — The continuous-time
differential Lyapunov equation is a specific form of a linear differential
equation with a homogeneous and an inhomogeneous part. In our microscopic
derivation, it becomes completely clear that homogeneous part is associated
with time-evolution via a non-Hermitian Hamiltonian, while the inhomogeneous
part is associated with quantum and thermal fluctuations due to sources of
loss. The Lyapunov equation therefore gives a natural way of identifying the
non-Hermitian Hamiltonian governing the dynamics and study the effect of
quantum and thermal fluctuations. Non-Hermitian physics is an extremely
rapidly growing field at present (see for example, Bergholtz  _et al._ (2021);
Ashida  _et al._ (2020); El-Ganainy  _et al._ (2018, 2019); Feng  _et al._
(2017)). But majority of works phenomenologically assume a non-Hermitian
Hamiltonian, which is very often expected to govern a classical system, and
quantum and thermal fluctuations are ignored. Only recently, there has been
growing interest in going beyond such descriptions, and realizing dynamics
governed by specific non-Hermitian Hamiltonians in quantum systems (for
example, Quijandría  _et al._ (2018); Purkayastha  _et al._ (2020b); Huber
_et al._ (2020a, b); Jaramillo Ávila  _et al._ (2020); McDonald  _et al._
(2022); Álvaro Gómez-León  _et al._ (2021); Arkhipov  _et al._ (2020, 2021);
Roccati  _et al._ (2021)). On the other hand, classifying all kinds of non-
Hermitian Hamiltonians and their relations to topology has remained an
extremely active direction of research (for example, Wang  _et al._ (2021);
Bergholtz  _et al._ (2021); Kawabata  _et al._ (2019); Torres (2019); Gong
_et al._ (2018); Harari  _et al._ (2018); Bandres  _et al._ (2018)). The
Lyapunov equation then gives the unified way to account for quantum and
thermal fluctuations in all such non-Hermitian systems, as long as there is no
source of gain, i.e, for all passive non-Hermitian systems. (The theory
presented here does not describe sources of gain. See however, Purkayastha
_et al._ (2020b).) If, for some physical reason, the baths can be considered
empty initially, then there are no quantum and thermal fluctuations, perfectly
realizing passive non-Hermitian systems. Interestingly, we find that, at all
the three levels of approximation for the Lyapunov equation, the non-Hermitian
Hamiltonian is the same. The three levels of approximations only change the
inhomogeneous part of the Lyapunov equation. Thus, they specify at what level
of accuracy the quantum and thermal fluctuations are considered.

Plan of the paper — In Sec. II we formally introduce the Lyapunov equation. In
Sec. III we give the set-up and the exact equations of motion. In Sec. IV we
derive the non-Hermitian Hamiltonian and the Lyapunov equations. In Sec. V we
obtain the associated QMEs. In Sec. VI we discuss positivity, accuracy and
additivity issues. In Sec. VII, we discuss the possibility of dark states,
give perturbative solutions and discuss thermalization. In Sec. VIII we derive
the generalized regression formulas. In Sec. IX, we discuss two insightful
examples. In Sec. X, we summarize and discuss the implications of our results.
A summary of our main results is given in Table. 1. Finally, there is an
Appendix, where details of some proofs are given.

##  II The Lyapunov equation

The continuous-time differential Lyapunov equation is an equation of the form

$\displaystyle\frac{d\mathbf{C}}{dt}=-\left(\mathcal{G}\mathbf{C}+\mathbf{C}\mathcal{G}^{\dagger}\right)+\epsilon^{2}\mathbf{Q},$
(1)

where $\mathbf{C}$, $\mathbf{Q}$ and $\mathcal{G}$ are $N\times N$ matrices,
for some given $N$, and $\mathcal{G}^{\dagger}$ is the conjugate transpose of
$\mathcal{G}$. The $\epsilon^{2}$ is used for notational convenience to be
consistent with that used later in the paper. Usually, as will also be in our
case, $\mathbf{C}$ is required to be Hermitian and positive semi-definite at
all times. The necessary and sufficient condition for this $\mathbf{Q}$ being
Hermitian and positive semi-definite. The formal solution of the Lyapunov
equation is

$\displaystyle\mathbf{C}(t)=e^{-\mathcal{G}t}\mathbf{C}(0)e^{-\mathcal{G}^{\dagger}t}+\epsilon^{2}\int_{0}^{t}dt^{\prime}e^{-\mathcal{G}t^{\prime}}\mathbf{Q}e^{-\mathcal{G}^{\dagger}t^{\prime}}.$
(2)

If the real parts of eigenvalues of $\mathcal{G}$ are positive, then are is a
unique steady solution in the long time limit, given by

$\displaystyle\mathbf{C}(\infty)=\epsilon^{2}\int_{0}^{\infty}dt^{\prime}e^{-\mathcal{G}t^{\prime}}\mathbf{Q}e^{-\mathcal{G}^{\dagger}t^{\prime}}.$
(3)

Since this is the steady state, this is the solution of

$\displaystyle\mathcal{G}\mathbf{C}(\infty)+\mathbf{C}(\infty)\mathcal{G}^{\dagger}=\epsilon^{2}\mathbf{Q},$
(4)

which is called the algebraic Lyapunov equation. This greatly simplifies the
problem of finding $\mathbf{C}(\infty)$ numerically. It is complete set of
$N^{2}$ linear equations. Usually, time taken to solve such a system of
equations scales as $N^{6}$. However, for the algebraic Lyapunov equation,
standard efficient algorithms (standard packages in Mathematica, python,
Matlab etc.) are available for which the time taken to solve scales as $N^{3}$
Gajić and Qureshi (1995); Behr  _et al._ (2019). Further, recently, even more
efficient Krylov subspace methods for solving Lyapunov equations have been
investigated Behr  _et al._ (2019). In the following, we derive equations of
the above form for dynamics of a very general linear open quantum system.

Figure 1:  A system on a lattice of arbitrary geometry and dimension coupled
to at an arbitrary number of sites to thermal baths which can all be at
different temperatures (inverse temperatures, say,
$\beta_{\ell},\beta_{m},\beta_{p},\beta_{q},\beta_{r}$) and chemical
potentials (say $\mu_{\ell},\mu_{m},\mu_{p},\mu_{q},\mu_{r}$) .

##  III The set-up and exact equations of motion

We consider the most general number-conserving linear system (quadratic
Hamiltonian) in a lattice of $N$ sites in arbitrary dimension and geometry,

$\displaystyle\hat{\mathcal{H}}_{S}=\sum_{\ell,m=1}^{N}\mathbf{H}_{\ell
m}\hat{c}_{\ell}^{\dagger}\hat{c}_{m}.$ (5)

where $\hat{c}_{\ell}$ is fermionic or bosonic annihilation operator, and
$\mathbf{H}$ is a Hermitian matrix, often called the single-particle
Hamiltonian. We are interested in the case where an arbitrary number of sites
can be attached to baths (see Fig. 1). However, in deriving the theory, for
notational simplicity and generality, we consider that each site of the system
is coupled to a bath, each of which is described by an infinite number of non-
interacting modes,

$\displaystyle\hat{\mathcal{H}}_{SB}=\epsilon\sum_{\ell=1}^{N}\hat{\mathcal{H}}_{SB_{\ell}},~{}\hat{\mathcal{H}}_{SB_{\ell}}=\epsilon\sum_{r=1}^{\infty}\left(\kappa_{rl}\hat{c}_{\ell}^{\dagger}\hat{B}_{r\ell}+\kappa_{rl}^{*}\hat{B}_{r\ell}^{\dagger}\hat{c}_{\ell}\right)$
$\displaystyle\hat{\mathcal{H}}_{B}=\sum_{\ell=1}^{N}\hat{\mathcal{H}}_{B_{\ell}},~{}\hat{\mathcal{H}}_{B_{\ell}}=\sum_{r=1}^{\infty}\Omega_{rl}\hat{B}_{r\ell}^{\dagger}\hat{B}_{r\ell}.$
(6)

Here $\hat{B}_{r\ell}$ is the fermionic or bosonic annihilation operator of
the $r$th mode of the bath attached to the $\ell$th site of the system. The
factor $\epsilon$ is a dimensionless parameter controlling the strength of
system-bath coupling. The initial state of the whole set-up is taken to be

$\displaystyle\hat{\rho}_{tot}(0)=\hat{\rho}(0)\prod_{\ell=1}^{N}\frac{e^{-\beta_{\ell}\left(\hat{\mathcal{H}}_{B_{\ell}}-\mu_{\ell}\hat{N}_{B_{\ell}}\right)}}{Z_{B_{\ell}}},$
(7)

where
$\hat{N}_{B_{\ell}}=\sum_{r=1}^{\infty}\hat{B}_{r\ell}^{\dagger}\hat{B}_{r\ell}$
is the operator for total number of particles in the bath attached at the
$\ell$th site and $Z_{B_{\ell}}$ is the partition function. So, the system is
initially in an arbitrary state while the baths are all in the thermal states
with their individual temperatures and chemical potentials. The spectral
functions of the baths are defined by

$\displaystyle\mathfrak{J}_{\ell}(\omega)=2\pi\sum_{r=1}^{\infty}|\kappa_{r\ell}|^{2}\delta(\omega-\Omega_{r\ell}).$
(8)

For future reference, we also introduce the notation $f^{H}(\omega)$ for the
Hilbert transform of any arbitrary function $f(\omega)$,

$\displaystyle
f^{H}(\omega)=\frac{1}{\pi}\mathcal{P}\int_{-\infty}^{\infty}d\omega^{\prime}\frac{f(\omega^{\prime})}{\omega-\omega^{\prime}},$
(9)

where $\mathcal{P}$ denotes the principal value.

The dynamics of the system can be described by the an exact quantum Langevin
equation Dhar and Sriram Shastry (2003); Dhar and Sen (2006); Dhar and Roy
(2006). This is derived in two steps. First, the formal solution for
annihilation operators of the baths is written down,

$\displaystyle\hat{B}_{r\ell}(t)=e^{-i\Omega_{r\ell}t}\hat{B}_{r\ell}(0)-i\epsilon\kappa_{r\ell}\int_{0}^{t}dt^{\prime}e^{-\Omega_{r\ell}(t-t^{\prime})}\hat{c}_{\ell}(t^{\prime}).$
(10)

Then this solution is used in the equation of motion for the annihilation
operators of the system

$\displaystyle\frac{d\hat{c}_{\ell}}{dt}$
$\displaystyle=-i\sum_{\ell,m=1}^{N}\mathbf{H}_{\ell
m}\hat{c}_{m}(t)-i\epsilon\sum_{r=1}^{\infty}\kappa_{r\ell}^{*}\hat{B}_{r\ell}(t),$
(11)

to obtain the so-called quantum Langevin equation,

$\displaystyle\frac{d\hat{c}_{\ell}}{dt}$
$\displaystyle=-i\sum_{\ell,m=1}^{N}\mathbf{H}_{\ell
m}\hat{c}_{m}(t)-i\epsilon\hat{\xi}_{\ell}(t)$
$\displaystyle-\epsilon^{2}\int_{0}^{t}dt^{\prime}\int\frac{d\omega}{2\pi}\mathfrak{J}_{\ell}(\omega)e^{-i\omega(t-t^{\prime})}\hat{c}_{\ell}(t^{\prime}),$
(12)

where the integration is over all bath frequencies and
$\hat{\xi}_{\ell}(t)=\sum_{r=1}^{\infty}\kappa_{r\ell}e^{i\Omega_{r\ell}t}\hat{B}_{r\ell}(0)$
is the noise operator. The noise correlation functions are given by

$\displaystyle\langle\hat{\xi}_{\ell}(t)\rangle=0,~{}\langle\hat{\xi}_{\ell}^{\dagger}(t)\hat{\xi}_{m}(t^{\prime})\rangle=\int\frac{d\omega}{2\pi}\mathbf{F}_{\ell
m}(\omega)e^{i\omega(t-t^{\prime})},$ (13)

with

$\displaystyle\mathbf{F}_{\ell
m}(\omega)=\mathfrak{J}_{\ell}(\omega)\mathfrak{n}_{B_{\ell}}(\omega)\delta_{\ell
m}.$ (14)

where

$\displaystyle\mathfrak{n}_{B_{\ell}}(\omega)=[e^{\beta_{\ell}(\omega-\mu_{\ell})}\pm
1]^{-1}$ (15)

is the fermi or bose distribution function, $\delta_{\ell m}$ is Kronecker
delta function and $\langle...\rangle={\rm
Tr}\left(...\hat{\rho}_{tot}(0)\right)$. Note that the quantum Langevin
equation, Eq.(III), is completely exact for our set-up. It does not require on
any further approximation. There is no Markovian approximation, no weak
system-bath coupling approximation, also no approximation on bath spectral
functions, except that they are continuous, and no approximation on the
temperatures and the chemical potentials of the baths. It also holds for both
fermionic and bosonic set-ups. This generality and fundamental nature of such
quantum Langevin equations, although known for a long time Ford  _et al._
(1965); Cortés  _et al._ (1985); Ford  _et al._ (1988); Dhar and Sriram
Shastry (2003); Dhar and Sen (2006); Dhar and Roy (2006), is often not
emphasized in standard textbooks of quantum optics (for example, Walls and
Milburn (2008); Wiseman and Milburn (2010)), where quantum Langevin equations
are instead derived from weak system-bath coupling Lindblad equations.

Another point to mention is that the effect of baths are additive at the level
of the exact quantum Langevin equation. This additivity, which stems from the
fact that the system-bath couplings are bilinear, allows for easy
generalization. We have assumed that each site is coupled to single bath. It
is clear that if some sites are not connected to a bath, it can be easily
incorporated just by setting the system-bath coupling of those sites to zero
in the quantum Langevin equation. This leads to dropping the last two terms
Eq.(III) for the corresponding sites. Further, it can also easily be
generalized to the case one or more sites are coupled to more than one bath.
This can be done just by adding the corresponding terms at the corresponding
sites in the quantum Langevin equation (an example of this will be given later
in Sec. IX.1). Keeping this in mind, for notational simplicity, in derivation
of our results, we stick to the set-up of having each site coupled to a single
bath.

The Eq.(III) can be exactly solved for NESS via a Fourier transform (assuming
the bandwidth of the bath includes all the system modes, a necessary condition
for unique NESS). It yields exactly the same expressions as obtained via a
NEGF approach Dhar and Roy (2006); Dhar and Sen (2006). However, the dynamics
of approach to NESS requires Laplace transform. This can be difficult,
because, inverse Laplace transform is hard. Depending on the number of lattice
sites and the spectral properties of $\mathbf{H}$, this can essentially be
intractable.

On the other hand, in absence of the baths, the dynamics of the system is
numerically tractable up to very large number of sites in the lattice. The
$N\times N$ matrix $\mathbf{H}$ can be diagonalized

$\displaystyle\Phi^{\dagger}\mathbf{H}\Phi=\mathbf{D},$ (16)

where $\mathbf{D}={\rm diag}\\{\omega_{\alpha}\\}$ is a diagonal matrix
containing the eigenvalues of the matrix $\mathbf{H}$, which are the single-
particle eigenvalues, $\Phi$ is a unitary matrix whose columns are the
corresponding single-particle eigenvectors and $\Phi^{\dagger}$ represents the
conjugate transpose. The dynamics of the isolated system is given in terms of
these is given as

$\displaystyle\hat{c}_{\ell}(t)=\sum_{\alpha,m=1}^{N}\Phi_{\ell\alpha}\Phi_{m\alpha}^{*}e^{-\omega_{\alpha}t}\hat{c}_{m}(0).$
(17)

Thus, given the single-particle eigenvectors and eigenvalues, the isolated
system dynamics can be obtained exactly. This is not true for the exact open
system dynamics. In the following, using various levels of Markov
approximations, we simplify the open system time evolution greatly, deriving
the Lyapunov equation.

##  IV Deriving the non-Hermitian Hamiltonian and the Lyapunov equations

The Eq.(III) is our starting point for further Born-Markov-like
approximations. The non-Markovianity of Eq.(III) is encoded in two places, (i)
the last term of Eq.(III), (ii) in the fact that the noise correlation is not
a delta function. Correspondingly, we will do two levels of Markov
approximations, involving two time scales, $\tau_{B_{1}}$ and $\tau_{B_{2}}$.
The corresponding Markov approximations are valid for times larger than these
time scales. As a collorary, we will also identify the non-Hermitian
Hamiltonian governing the time evolution of open quantum system under such
approximations.

###  IV.1 First Markov approximation

####  IV.1.1 Equation of motion and the non-Hermitian Hamiltonian

We call the approximation that deals with the last term in Eq.(III) the first
Markov approximation. The main goal is to make the last term depend only on
$\hat{c}_{\ell}(t)$, and not on its past history. There can be two different
approximations to obtain this. The first is the so-called wide-band limit,
with

$\displaystyle\mathfrak{J}_{\ell}(\omega)=\Gamma_{\ell}.$ (18)

This directly gives

$\displaystyle\int_{0}^{t}dt^{\prime}\int\frac{d\omega}{2\pi}\mathfrak{J}_{\ell}(\omega)e^{-i\omega(t-t^{\prime})}\hat{c}_{\ell}(t^{\prime})=\frac{\Gamma_{\ell}}{2}\hat{c}_{\ell},$
(19)

where the factor of $1/2$ appears because the time integration is from $0$ to
$t$. This approximation does not depend on strength of system-bath coupling,
$\Gamma_{\ell}$ can be arbitrary. While this approximation is often used for
fermionic systems, it is artificial for bosonic systems, where
$\mathfrak{J}_{\ell}(\omega)$ must go to zero at $\omega=0$. So, in the
following, we will focus on the other way of doing the first Markov
approximation, which explicitly depends on weak system-bath coupling. Later,
we will discuss which of the approximations hold in the wide-band limit for
fermionic systems without assuming weak system-bath coupling, and which do
not.

We will write Eq.(III) correct to $O(\epsilon^{2})$. The last term of the
equation is already or $O(\epsilon^{2})$. So, following Eq.(17) we use the
result,

$\displaystyle\hat{c}_{\ell}(t^{\prime})=\sum_{\alpha,m=1}^{N}\Phi_{\ell\alpha}\Phi_{m\alpha}^{*}e^{\omega_{\alpha}(t-t^{\prime})}\hat{c}_{m}(t)+O(\epsilon),$
(20)

and neglect the $O(\epsilon)$ term. This gives

$\displaystyle\epsilon^{2}\int_{0}^{t}dt^{\prime}\int\frac{d\omega}{2\pi}\mathfrak{J}_{\ell}(\omega)e^{-i\omega(t-t^{\prime})}\hat{c}_{\ell}(t^{\prime})$
(21) $\displaystyle\simeq$
$\displaystyle\epsilon^{2}\sum_{\alpha,m=1}^{N}\Phi_{\ell\alpha}\Phi_{m\alpha}^{*}\hat{c}_{m}(t)\int_{0}^{t}dt^{\prime}\int\frac{d\omega}{2\pi}\mathfrak{J}_{\ell}(\omega)e^{-i(\omega-\omega_{\alpha})t^{\prime}}.$

The above is essentially a Born approximation. Next we make a Markov
approximation. Let $\tau_{B_{1}}$ be defined via the condition,

$\displaystyle\left|\int\frac{d\omega}{2\pi}\mathfrak{J}_{\ell}(\omega)e^{-i\omega
t}\right|<\textrm{some tolerance, say $O(\epsilon)$,
}\forall~{}t>\tau_{B_{1}}$ (22)

for all $\ell$. If $t\gg\tau_{B_{1}}$, we can essentially take the upper limit
of the time integral to infinity since this only causes changes in higher
order terms. This, with a little algebra, leads to

$\displaystyle\frac{dc_{vec}}{dt}=-i\mathbf{H}_{\rm
NH}~{}c_{vec}(t)-i\epsilon\xi_{vec}(t),$ (23)

where $c_{vec}(t)$ ($\xi_{vec}(t)$) is a column vector with the $\ell$th
element being $\hat{c}_{\ell}(t)$ ($\hat{\xi}_{\ell}(t)$), and

$\displaystyle\mathbf{H}_{\rm NH}=\mathbf{H}-i\epsilon^{2}\mathbf{v},$ (24)

is the non-Hermitian Hamiltonian governing the open system dynamics. The
elements of the matrix $\mathbf{v}$ are given by,

$\displaystyle\mathbf{v}_{\ell
m}=\frac{1}{2}\sum_{\alpha=1}^{N}\Phi_{\ell\alpha}\Phi_{m\alpha}^{*}\Big{(}\mathfrak{J}_{\ell}(\omega_{\alpha})+i\mathfrak{J}_{\ell}^{H}(\omega_{\alpha})\Big{)},$
(25)

with $\mathfrak{J}_{\ell}^{H}(\omega)$ being the Hilbert transform of
$\mathfrak{J}_{\ell}(\omega)$ (see Eq.(9)). Thus, the first Markov
approximation allows identification of the (single-particle) non-Hermitian
Hamiltonian that governs the dynamics of the open quantum system. All kinds of
non-Hermitian Hamiltonians with losses can be derived by changing the relative
strengths of system-bath couplings and the choosing the spectral functions.
However, baths of this kind are unable to describe sources of gain. So, this
gives a microscopic way of generating all kinds of passive non-Hermitian
systems. In accordance with fluctuation-dissipation theorem, the sources of
loss also give sources of noise, embodied in the term $\xi_{vec}(t)$. Thus,
the effects of quantum and thermal fluctuations on all the classes of passive
non-Hermitian systems, as well as possible transitions between them Bergholtz
_et al._ (2021); Kawabata  _et al._ (2019); Torres (2019); Gong  _et al._
(2018), can be studied in this set-up. It can also be seen that if
temperatures and chemical potentials of the baths are such that they can be
considered empty, i.e, $\mathfrak{n}_{B_{\ell}}(\omega)=0$, then there will be
no noise, thereby perfectly realizing quantum passive non-Hermitian systems.

Note that, the matrix $\mathbf{v}$ in Eq.(25) is not diagonal in the site
basis in general. However, for fermionic wide-band baths, Eq.(18),
$\mathbf{v}$ is diagonal in the site basis,

$\displaystyle\mathbf{v}_{\ell m}=\frac{\Gamma_{\ell}}{2}\delta_{\ell
m}~{}~{}~{}\textrm{for wide band limit}.$ (26)

As mentioned, this case does not require weak system-bath coupling
approximation (i.e, holds with $\epsilon=1$).

Having derived the non-Hermitian Hamiltonian, we now derive the Lyapunov
equations.

####  IV.1.2 The correlation matrix and the Lyapunov-like form

The formal solution of Eq.(23) can be written as

$\displaystyle c_{vec}(t)=e^{-i\mathbf{H}_{\rm
NH}t}c_{vec}(0)-i\epsilon\int_{0}^{t}dt^{\prime}e^{-i\mathbf{H}_{\rm
NH}(t-t^{\prime})}\xi_{vec}(t^{\prime}).$ (27)

We will be interested in the the correlation matrix $\mathbf{C}$ whose
elements are

$\displaystyle\mathbf{C}_{\ell
m}(t)=\langle\hat{c}_{\ell}^{\dagger}(t)\hat{c}_{m}(t)\rangle.$ (28)

This matrix is Hermitian and positive semi-definite by construction. The
expression for $\mathbf{C}(t)$ can be obtained from Eq.(27) after some algebra
(see Appendix A for explicit steps),

$\displaystyle\mathbf{C}(t)=e^{-\mathcal{G}t}\mathbf{C}(0)e^{-\mathcal{G^{\dagger}}t}$
$\displaystyle+\epsilon^{2}\int\frac{d\omega}{2\pi}\left(\frac{1-e^{-(\mathcal{G}+i\omega\mathbb{I})t}}{\mathcal{G}+i\omega\mathbb{I}}\mathbf{F}(\omega)\frac{1-e^{-(\mathcal{G}^{\dagger}-i\omega\mathbb{I})t}}{\mathcal{G}^{\dagger}-i\omega\mathbb{I}}\right),$
(29)

where $\mathbb{I}$ is $N$ dimensional identity matrix,

$\displaystyle\mathcal{G}=-i\mathbf{H}_{\rm NH}^{*},$ (30)

with $\mathbf{H}_{\rm NH}^{*}$ representing complex conjugate of
$\mathbf{H}_{\rm NH}$ and $\mathbf{F}(\omega)$ is a $N\times N$ matrix defined
in Eq.(14). Since $\mathbf{F}(\omega)$ is a positive semi-definite diagonal
matrix by construction, the structure of Eq.(IV.1.2) ensures Hermiticity and
positive semi-definiteness of $\mathbf{C}(t)$ at all times. If the NESS is
unique, i.e, if the real parts of eigenvalues of $\mathcal{G}$ are positive,
$\mathbf{C}(\infty)$ is given by

$\displaystyle\mathbf{C}(\infty)=\epsilon^{2}\int\frac{d\omega}{2\pi}\left(\frac{1}{\mathcal{G}+i\omega\mathbb{I}}\mathbf{F}(\omega)\frac{1}{\mathcal{G}^{\dagger}-i\omega\mathbb{I}}\right).$
(31)

This equation could also be obtained by solving Eq.(23) by a fourier transform
and then obtaining the correlation matrix.

While Eq.(IV.1.2) can be directly used to numerically obtain the correlation
matrix, it does not reveal the relation to the Lyapunov equation. To see this,
we write down the equation of motion for the correlation matrix directly from
Eq.(23). After some simplification, this can be written as

$\displaystyle\frac{d\mathbf{C}}{dt}=-(\mathcal{G}\mathbf{C}+\mathbf{C}\mathcal{G}^{\dagger})+\epsilon^{2}\mathbf{Q}(t),$
$\displaystyle\mathbf{Q}(t)=i\left(\mathbf{C}_{\xi}(t)-\mathbf{C}_{\xi}^{\dagger}(t)\right),$
(32)

where elements of the matrix $\mathbf{C}_{\xi}(t)$ are given by
$\mathbf{C}_{\xi_{\ell
m}}(t)=\langle\hat{\xi}_{\ell}^{\dagger}(t)\hat{c}_{m}(t)\rangle$. Using
Eq.(27), $\mathbf{C}_{\xi}(t)$ is obtained as

$\displaystyle\mathbf{C}_{\xi}(t)=-i\int_{0}^{t}dt^{\prime}\int\frac{d\omega}{2\pi}\mathbf{F}(\omega)e^{i\omega
t^{\prime}}e^{-\mathcal{G}^{\dagger}t^{\prime}}.$ (33)

The Eq.(IV.1.2) is exactly of the Lyapunov equation form, except for a time-
dependent inhomogeneous part. It can be checked that Eq.(IV.1.2) is the
solution of this equation (see Appendix B). Therefore, this Lyapunov-like form
preserves the positive semi-definiteness of $\mathbf{C}(t)$ at all times. Also
note that since at the first-Markov approximation level we treat the noise
correlations without further approximations, and the noise correlations are
not delta functions in time, the dynamics of the system is actually explicitly
non-Markovian in this case. Next, we make approximations on the noise
correlations.

###  IV.2 Second Markov approximation

It can be seen that the time-dependence of the inhomogeneous part in
Eq.(IV.1.2) stems from the noise correlation not being a delta function. We do
the second Markov approximation to do away with this time dependence. This can
be done at two levels.

####  IV.2.1 Level-I second Markov approximation

Let $\tau_{B_{2}}$ be defined via the condition

$\displaystyle\left|\int\frac{d\omega}{2\pi}\mathfrak{J}_{\ell}(\omega)\mathfrak{n}_{B_{\ell}}(\omega)e^{i\omega
t}\right|<\textrm{some tolerance, say $O(\epsilon)$, }$
$\displaystyle~{}\forall~{}t>\tau_{B_{2}}$ (34)

for all $\ell$. So, assuming $t\gg\tau_{B_{2}}$, from Eq.(14), we see that we
can essentially extend the upper limit of the time integral in Eq.(33) to
infinity, since this only changes the higher order terms. This gives the
Lyapunov equation,

$\displaystyle\frac{d\mathbf{C}}{dt}=-(\mathcal{G}\mathbf{C}+\mathbf{C}\mathcal{G}^{\dagger})+\epsilon^{2}\mathbf{Q}_{1},$
(35)
$\displaystyle\mathbf{Q}_{1}=\int\frac{d\omega}{2\pi}\left(\mathbf{F}(\omega)\frac{1}{\mathcal{G}^{\dagger}-i\omega\mathbb{I}}+\frac{1}{\mathcal{G}+i\omega\mathbb{I}}\mathbf{F}(\omega)\right).$
(36)

The $\mathbf{Q}_{1}$ above is nothing but $\mathbf{Q}(\infty)$. This level-I
second Markov approximation approximation, by construction, becomes more and
more accurate as $t$ is increased. In the $t\rightarrow\infty$ limit,
corresponding to NESS, the results from Eq.(35) and Eq.(IV.1.2) match. Thus,
$\mathbf{C}(\infty)$ obtained from the algebraic version of above Lyapunov
equation (see Eq.(4)) is Hermitian and positive semi-definite as required. If
the NESS is unique, it is given by Eq.(31) (see Appendix B).

However, $\mathbf{Q}_{1}$ given in the expression in Eq.(35) is not guaranteed
to be positive semi-definite. Thus, unlike Eq.(IV.1.2), Eq.(35) may not
preserve the positive semi-definiteness $\mathbf{C}(t)$ at all times.
Nevertheless, the negative eigenvalues of $\mathbf{C}(t)$, if at all they
exist, will go towards zero on reducing $\epsilon$ and with increase in time.
So, they should be taken to be zero within the accuracy regime of this
equation.

The fact that the NESS obtained from Eq.(35) is positive semi-definite even
though $\mathbf{Q}_{1}$ may not be so is presumably quite interesting on the
mathematical front. This is because, it means that the algebraic Lyapunov
equation, Eq.(4), can provably yield positive semi-definite results even when
its inhomogeneous part is not positive semi-definite.

Next, we make one more level of approximation, which, as we will see later,
makes the description equivalent to that obtained from the Redfield equation.

####  IV.2.2 Level-II second Markov approximation

In this level, a further approximation, using the fact that the inhomogeneous
part is already $O(\epsilon^{2})$, is made to essentially perform the
integration in the definition of $\mathbf{Q}_{1}$ in Eq.(35). To do this, we
define,

$\displaystyle\mathbf{F}^{E}(\omega)=\Phi^{\dagger}\mathbf{F}(\omega)\Phi,~{}\mathbf{v}^{E}=\Phi^{\dagger}\mathbf{v}\Phi,$
$\displaystyle~{}\mathcal{G}^{E}=\Phi^{\dagger}\mathcal{G}\Phi,~{}~{}\mathcal{G}^{E}=-i\mathbf{D}-\epsilon^{2}\mathbf{v}^{E}.$
(37)

From Eq.(33) and above definitions, it follows that

$\displaystyle\mathbf{C}_{\xi}(t)\simeq$
$\displaystyle-i\Phi\left(\int_{0}^{t}dt^{\prime}\int\frac{d\omega}{2\pi}\mathbf{F}^{E}(\omega)e^{i\omega
t^{\prime}}e^{-i\mathbf{D}t^{\prime}}\right)\Phi^{\dagger}$
$\displaystyle+O(\epsilon^{2}).$ (38)

Assuming $t\gg\tau_{B_{2}}$ to extend the upper limit of the time integral to
infinity, performing the integral, and neglecting $O(\epsilon^{2})$ terms, we
obtain,

$\displaystyle\mathbf{C}_{\xi_{\ell
m}}(\infty)\simeq-\frac{i}{2}\sum_{\alpha=1}^{N}\Phi_{\ell\alpha}\Phi_{m\alpha}^{*}\Big{(}\mathbf{F}_{\ell\ell}(\omega_{\alpha})-i\mathbf{F}_{\ell\ell}^{H}(\omega_{\alpha})\Big{)}.$
(39)

The Lyapunov equation is then given by

$\displaystyle\frac{d\mathbf{C}}{dt}=-(\mathcal{G}\mathbf{C}+\mathbf{C}\mathcal{G}^{\dagger})+\epsilon^{2}\mathbf{Q}_{2},$
(40)

with the elements of the matrix $\mathbf{Q}_{2}$ given by

$\displaystyle\mathbf{Q}_{2_{\ell m}}=\frac{1}{2}\sum_{\alpha=1}^{N}$
$\displaystyle\Phi_{\ell\alpha}\Phi_{m\alpha}^{*}\Big{(}\mathbf{F}_{\ell\ell}(\omega_{\alpha})-i\mathbf{F}_{\ell\ell}^{H}(\omega_{\alpha})$
$\displaystyle+\mathbf{F}_{mm}(\omega_{\alpha})+i\mathbf{F}_{mm}^{H}(\omega_{\alpha})\Big{)}.$
(41)

Thus, $\mathbf{Q}_{2}$ is given in terms of single-particle eigenvalues and
eigenvectors of the system. Even though $\mathbf{Q}_{2}$ is Hermitian, like
$\mathbf{Q}_{1}$, it is not guaranteed to be positive semi-definite. Moreover,
unlike Eq.(IV.1.2), even the steady state solution of Eq.(40) may not strictly
be same as that from Eq.(35). So there might be violation of positive semi-
definiteness of $\mathbf{C}(\infty)$. But, the deviations from the results of
Eq.(35) will be small and will decrease on decreasing $\epsilon$.

The level-II second markov approximation explicitly requires weak system-bath
coupling. It requires weak system-bath coupling even for fermionic baths in
the wide-band limit (Eq.(18)) (unless all the baths are either completely full
or completely empty, i.e, the chemical potentials are $\pm\infty$). However,
the first Markov approximation and the level-I second Markov approximation do
not require weak system-bath coupling approximation for fermionic baths at
wide-band limit. The first-Markov approximation gives exact dynamics in this
case, while the level-I second Markov approximation gives approximate dynamics
but exact steady state (even when system-bath coupling is strong, i.e,
$\epsilon=1$).

Having derived the Lyapunov equations, in the next section, we find the QMEs
which can be associated with them.

##  V Relation to quantum master equations

In above, we have derived the Lyapunov equations directly from equations of
motion. Usually, Lyapunov equations arise in open quantum systems when
calculating correlation matrices of linear systems governed by some QME. The
question, then, is what are the QMEs corresponding to the Lyapunov equations
derived in the previous section.

Using standard Born-Markov approximations (without any secular approximation),
the Redfield QME corresponding to our set-up can be derived. This is given by
Purkayastha  _et al._ (2016)

$\displaystyle\frac{\partial\hat{\rho}}{\partial
t}=i[\hat{\rho},\hat{\mathcal{H}}_{S}+\hat{\mathcal{H}}_{LS}]+\epsilon^{2}\sum_{\ell,m=1}^{N}\Big{[}\big{(}\mathbf{v}_{\ell
m}+\mathbf{v}_{m\ell}^{*}\mp\mathbf{Q}_{2_{m\ell}}\big{)}\left(\hat{c}_{m}\rho\hat{c}_{\ell}^{\dagger}-\frac{1}{2}\\{\hat{c}_{\ell}^{\dagger}\hat{c}_{m},\hat{\rho}\\}\right)+\mathbf{Q}_{2_{m\ell}}\left(\hat{c}_{\ell}^{\dagger}\rho\hat{c}_{m}-\frac{1}{2}\\{\hat{c}_{m}\hat{c}_{\ell}^{\dagger},\hat{\rho}\\}\right)\Big{]},$
(42)

where $-$ sign is for fermions and $+$ sign is for bosons,
$\\{\hat{A},\hat{B}\\}=\hat{A}\hat{B}+\hat{B}\hat{A}$ is the anti-commutator,
$\mathbf{Q}_{2}$ and $\mathbf{v}$ are as given in Eqs.(IV.2.2) and (25), and

$\displaystyle\hat{\mathcal{H}}_{LS}=\sum_{\ell,m,\alpha=1}^{N}\Phi_{\ell\alpha}\Phi_{m\alpha}^{*}\left(\frac{\mathfrak{J}_{\ell}^{H}(\omega_{\alpha})+\mathfrak{J}_{m}^{H}(\omega_{\alpha})}{4}\right)\hat{c_{\ell}}^{\dagger}\hat{c}_{m}$
(43)

is the so called Lamb-shift Hamiltonian. It can be checked, directly by
calculating the correlation matrix $\mathbf{C}_{\ell m}(t)={\rm
Tr}\left(\hat{c}_{\ell}^{\dagger}\hat{c}_{m}\hat{\rho}(t)\right)$, that the
corresponding Lyapunov equation is nothing but Eq.(40). This shows equivalence
of the equation of motion approach and the Redfield QME approach for Gaussian
systems. This is of course expected as they describe the same set-up under
same approximations.

From above, it follows by direct inspection that

$\displaystyle\frac{\partial\hat{\rho}}{\partial
t}=i[\hat{\rho},\hat{\mathcal{H}}_{S}+\hat{\mathcal{H}}_{LS}]+\epsilon^{2}\sum_{\ell,m=1}^{N}\Big{[}\big{(}\mathbf{v}_{\ell
m}+\mathbf{v}_{m\ell}^{*}\mp\mathbf{Q}_{1_{m\ell}}\big{)}\left(\hat{c}_{m}\rho\hat{c}_{\ell}^{\dagger}-\frac{1}{2}\\{\hat{c}_{\ell}^{\dagger}\hat{c}_{m},\hat{\rho}\\}\right)+\mathbf{Q}_{1_{m\ell}}\left(\hat{c}_{\ell}^{\dagger}\rho\hat{c}_{m}-\frac{1}{2}\\{\hat{c}_{m}\hat{c}_{\ell}^{\dagger},\hat{\rho}\\}\right)\Big{]}$
(44)

is the QME corresponding to the Lyapunov equation after level-I second Markov
approximation Eq.(35), while,

$\displaystyle\frac{\partial\hat{\rho}}{\partial
t}=i[\rho,\hat{\mathcal{H}}_{S}+\hat{\mathcal{H}}_{LS}]+\epsilon^{2}\sum_{\ell,m=1}^{N}\Big{[}\Big{(}\mathbf{v}_{\ell
m}+\mathbf{v}_{m\ell}^{*}\mp\mathbf{Q}_{{m\ell}}(t)\Big{)}\left(\hat{c}_{m}\rho\hat{c}_{\ell}^{\dagger}-\frac{1}{2}\\{\hat{c}_{\ell}^{\dagger}\hat{c}_{m},\hat{\rho}\\}\right)+\mathbf{Q}_{{m\ell}}(t)\left(\hat{c}_{\ell}^{\dagger}\rho\hat{c}_{m}-\frac{1}{2}\\{\hat{c}_{m}\hat{c}_{\ell}^{\dagger},\hat{\rho}\\}\right)\Big{]}$
(45)

is the QME corresponding to the Lyapunov equation with time-dependent
inhomogeneous part after first Markov approximation, Eq.(IV.1.2). These QMEs
would be hard to derive using the standard approaches. In particular, Eq.(45)
has explicitly time-dependent rates, which means it can potentially describe
non-Markovian dynamics. Further, for fermionic systems in the wide-band limit,
Eq.(45) is essentially exact, i.e, does not require any more approximation
(holds even in strong system-bath coupling, i.e, when $\epsilon=1$). On the
other hand, Eq.(44), although makes some approximation on the dynamics, in
this case, it gives the correct steady state exactly, without any further
approximations. In the next section, we discuss in more detail the positivity,
accuracy and additivity issues of the Lyapunov equations and the associated
QMEs.

##  VI Positivity, accuracy and additivity

The positivity issues of the Lyapunov equation have been mentioned above,
while those in QMEs have already been well studied Lindblad (1976); Gorini
_et al._ (1976, 1978); Breuer and Petruccione (2007); Rivas and Huelga (2012).
In the following, we discuss both kinds of positivity issues, the relation
between them, their origin from the accuracy issue and how they may be
circumvented in practice. We also comment on additivity of the QMEs.

###  VI.1 Positivity

All the QMEs given above are in the Gorini-Kossakowski-Sudarshan-Lindblad
(GKSL) Lindblad (1976); Gorini  _et al._ (1976, 1978) form. For time-
independent rates, Eqs.(42), (44), the condition for complete positivity of
the density matrix $\hat{\rho}$ is given by Breuer and Petruccione (2007);
Rivas and Huelga (2012)

$\displaystyle\mathbf{v}^{*}+\mathbf{v}^{T}\mp\mathbf{Q}_{1,2}\rightarrow\textrm{positive
semi-definite}$ $\displaystyle\mathbf{Q}_{1,2}\rightarrow\textrm{positive
semi-definite}.$ (46)

On the other hand, as we have discussed before, positive semi-definiteness of
the correlation matrix $\mathbf{C}$ requires

$\displaystyle\mathbf{Q}_{1,2}\rightarrow\textrm{positive semi-definite}.$
(47)

Thus, the first condition in Eq.(VI.1) is not required. It follows that, for
Gaussian initial states of the system, only $\mathbf{Q}_{1,2}$ being positive
semi-definite is sufficient for positivity of $\hat{\rho}$. To see this, we
note that QMEs derived above preserve the Gaussianity of the initial state.
This is, of course, consistent, because the actual time evolution, without any
approximation,
$e^{-i\hat{\mathcal{H}}t}\hat{\rho}_{tot}(0)e^{i\hat{\mathcal{H}}t}$ has this
property. Now, for all Gaussian states, the density matrix can be constructed
from the correlation matrix (see, for example,Chung and Peschel (2001); Dhar
_et al._ (2012); Eisler and Peschel (2017)), giving a one-to-one mapping
between the two at all times. So, a valid, positive semi-definite correlation
matrix, guaranteed by Eq.(47), will yield a valid density matrix. If, on the
other hand, the initial state is non-Gaussian, this one-to-one mapping breaks.
The correlation matrix has no information about the non-Gaussianity of the
initial state. To ensure complete positivity of the density matrix at all
times in such cases, both conditions in Eq.(VI.1) are required. It immediately
follows that if the steady state is unique irrespective of the initial
condition (i.e, if the real parts of eigenvalues of $\mathcal{G}$ are
positive), then it must be Gaussian. This means that to ensure positivity in
the steady state, only Eq.(47) is sufficient.

The above discussion leads us to the following two important conclusions. The
first is that, for Gaussian initial states, the dynamics obtained from Eq.(45)
is free from positivity issues at all times. The second is that, the NESS of
Eq.(44) is free from any positivity issues. This is despite the fact that in
general, $\mathbf{Q}_{1}$ may not be positive semi-definite and therefore
Eq.(44) may not be completely positive. In fact, for fermionic wide-band
baths, even in cases where Eq.(44) is not completely positive, it always
yields exact NESS results. This is contrary to the somewhat popular belief
that complete positivity is a necessary requirement for accurately describing
the steady state.

###  VI.2 Accuracy

The positivity issues remain in the Redfield QME Eq.(42), and the
corresponding Lyapunov equation Eq.(40), even for Gaussian initial states and
even in the steady state. But this is related to accuracy of the results. Both
of these are differential equations written correct to $O(\epsilon^{2})$. But,
their solution requires an exponentiation, thereby generating all orders of
$\epsilon$ in the result. Clearly, all orders of $\epsilon$ in the result are
not correct. For the Redfield equation, it can be shown that the diagonal
elements of $\hat{\rho}$ in the eigenbasis of the system Hamiltonian
$\hat{\mathcal{H}}_{S}$ are given correctly to $O(1)$ (leading order), the
error occurring at $O(\epsilon^{2})$, while the off-diagonal elements in that
basis are given correctly to $O(\epsilon^{2})$ (leading order), the error
occurring at $O(\epsilon^{4})$ Fleming and Cummings (2011); Tupkary  _et al._
(2022) (assuming no degeneracy). The correlation matrix in the eigenbasis of
the system Hamiltonian is given by,

$\displaystyle\mathbf{C}^{E}(t)=\Phi^{\dagger}\mathbf{C}(t)\Phi.$ (48)

Following similar arguments, it can be shown that
$\mathbf{C}^{E}_{\alpha\alpha}(t)$ is obtained correct to $O(1)$ while
$\mathbf{C}^{E}_{\alpha\nu}(t)$, $\alpha\neq\nu$ is given correct to
$O(\epsilon^{2})$. It is exactly this mismatch in orders of accuracy between
the diagonal and the off-diagonal elements that leads to positivity violation
both for $\hat{\rho}(t)$ and $\mathbf{C}(t)$. It can be checked that, the
diagonal elements in any basis carry an error of $O(\epsilon^{2})$. So
$O(\epsilon^{2})$ violation of positivity is related with accuracy limits of
the equations.

If positivity is restored by doing ad-hoc approximations at the level of the
differential equations, for example, by neglecting the negative eigenvalues of
$\mathbf{Q}_{1,2}$, or making secular approximations, it does not guarantee
better accuracy. In fact, accuracy usually become worse Tupkary  _et al._
(2022); Davidovic (2021). On the other hand, the results from the Redfield QME
and the corresponding Lyapunov equation can be checked by scaling the required
matrix elements with $\epsilon$ Purkayastha  _et al._ (2020a); Tupkary  _et
al._ (2022). If the scaling observed is higher than that dictated by the
accuracy, then that matrix element is to be taken as zero. This is because, it
would mean that the leading order contribution is zero while the non-zero
values are coming as an artefact of the higher orders present in the solution
of the equation. This gives a completely controlled way of checking and
correcting the results. Further, interesting techniques for obtaining the
$O(\epsilon^{2})$ correction to diagonal elements in the energy eigenbasis at
steady state from the Redfield equation have been developed Thingna  _et al._
(2012, 2013); Xu  _et al._ (2017). These techniques can be used to alleviate
positivity issues at steady state Purkayastha  _et al._ (2020a).

Also, in the first and the second Markov approximations, we have assumed
$t\gg\tau_{B_{1}}$ (Eq.(22)) and $t\gg\tau_{B_{2}}$ (Eq.(IV.2.1))
respectively. So results for times smaller than these time scales will not be
accurate.

###  VI.3 Additivity

A point to note is that in all the forms of Lyapunov equations and the
associated QMEs the contribution from each bath comes as an additive term. It
is sometimes believed that additive QMEs cannot give accurate NESS, especially
at strong system-bath coupling Chan  _et al._ (2014); Mitchison and Plenio
(2018); Giusteri  _et al._ (2017); Kołodyński  _et al._ (2018); Maguire  _et
al._ (2019). However, as mentioned before, Eqs.(45) and (44) give exact NESS
for fermionic case with wide-band baths, irrespective of the strenght of
system-bath coupling. But clearly they are additive. This prompts a deeper
discussion.

Most works exploring inaccuracy of additive QMEs (except Ref. Kołodyński  _et
al._ (2018)), refer to additivity of equations in the diagonal Lindblad form,

$\displaystyle\frac{\partial\hat{\rho}}{\partial
t}=i[\hat{\rho},\mathcal{\hat{H}}_{S}]+\hat{\mathcal{L}}(\hat{\rho}),$
$\displaystyle\hat{\mathcal{L}}(\hat{\rho})=i[\hat{\rho},\mathcal{\hat{H}}_{LS}]+\sum_{\lambda}\gamma_{\lambda}\Big{(}\hat{L}_{\lambda}\hat{\rho}\hat{L}_{\lambda}^{\dagger}-\frac{1}{2}\\{\hat{L}_{\lambda}^{\dagger}\hat{L}_{\lambda},\hat{\rho}\\}\Big{)},$
(49)

where $\hat{L}_{\lambda}$ are called Lindblad operators, and
$\gamma_{\lambda}$ are called rates. Clearly, for a system with more than one
site, none of Eqs.(42), (44) and (45) are explicitly in this form. Rather,
they are in the so-called off-diagonal GKSL form. They can be cast into the
form of Eq.(VI.3) by making a change of basis Lindblad (1976); Gorini  _et
al._ (1976, 1978); Breuer and Petruccione (2007); Rivas and Huelga (2012). The
Lindblad operators and the rates so obtained will no longer have just additive
contribution from each bath, but rather will be a linear combination of
contributions from all the baths. Thus, once converted to form of Eq.(VI.3),
they will be non-additive in this strict sense. Note that, the above
statements also hold true for the Redfield equation Eq.(42), which therefore
would be non-additive in this sense, whenever the system has more than one
site.

A single-site system, then, gives an exception. A single fermionic site
coupled to two wide-band fermionic baths at different temperatures and
chemical potentials is the so-called resonant level model. This is one of the
simplest and extremely well-studied open quantum system. Exact results are
known from various approaches like NEGF Haug and Jauho (2008). As we will see
explicitly later in Sec.IX.1, it turns out, for this system, the QME after
level-I second Markov approximation, i.e, the analog of Eq.(44), is additive,
perfectly of (diagonal) Lindblad form, and yet gives the exact NESS answer.
This simple example therefore shows perfectly additive QMEs, even this the
strict diagonal Lindblad sense, may be able to give exact results in some
cases. This example does not even satisfy the sufficient condition for
microscopically derived additive QMEs found in Ref. Kołodyński  _et al._
(2018).

From above discussion, we see that our equations suggest a non-trivial result:
linear fermionic open systems in the wide-band limit are governed by additive
QMEs (not in the strict diagonal Lindblad sense) at all strengths of system-
bath coupling. Interestingly, the exact QME for linear systems, which can be
derived without wide-band limit approximation, is manifestly non-additive
Zhang  _et al._ (2012). Thus, it seems that, at least for linear systems,
additivity of QMEs depends on the validity of our first Markov approximation.
This can be done exactly for all coupling strengths for wide-band fermionic
baths, which makes the corresponding QME additive. For bosonic baths, it would
require a weak system-bath coupling approximation.

With the positivity, accuracy and additivity issues of the Lyapunov equations
and the associated QMEs clarified, in the next section, we look at the
possibility of having dark states, and also provide perturbative solutions in
the regime of very small system-bath coupling.

##  VII Dark states, the perturbative solution and thermalization

###  VII.1 Dark states

Dark states are eigenstates of the system Hamiltonian which are left invariant
by presence of the baths (for example, see Emary (2007); Buča and Prosen
(2012); Finkelstein-Shapiro  _et al._ (2019); Quach and Munro (2020)).
Depending on the geometry of the lattice, it may so happen that the
eigenfunction of some system mode has nodes at exactly the sites where the
baths are attached. In that case, it can be checked by transforming the exact
quantum Langevin equation (Eq.III)) to the single-particle eigenbasis that the
corresponding mode is completely detached from the baths, and will evolve in
the same way as the isolated system. So, if the system was initially prepared
in that state, it will remain in that state even in presence of the baths,
making it a dark state. This, in turn, means, there is no unique steady state
of the system. By transforming the Lyapunov equations to the single-particle
eigenbasis, it is easy to show that this property is respected by them.

###  VII.2 Perturbative solution: dynamics and NESS

In the limit of very small system-bath coupling, the Level-II second Markov
approximation is good for single-time correlations, the corresponding Lyapunov
equation being Eq.(40). Although this equation can be exactly solved, it gives
correct answers only up to the leading order term in $\epsilon$. So, it is
useful to find analytical expressions for $\mathbf{C}(t)$ up to the leading
order term in $\epsilon$. In order to do so, we transform Eq.(40) to single-
particle eigenbasis

$\displaystyle\frac{d\mathbf{C}^{E}}{dt}=-(\mathcal{G}^{E}\mathbf{C}^{E}+\mathbf{C}{\mathcal{G}^{E}}^{\dagger})+\epsilon^{2}\mathbf{Q}^{E}_{2},$
(50)

where $\mathcal{G}^{E}$, $\mathbf{C}^{E}$ are as defined in Eqs.(IV.2.2),
(48), and $\mathbf{Q}^{E}_{2}=\Phi^{\dagger}\mathbf{Q}_{2}\Phi$ likewise. The
explicit expressions for the elements of $\mathbf{v}^{E}$ and
$\mathbf{Q}_{2}^{E}$ are

$\displaystyle\mathbf{v}^{E}_{\alpha\nu}$
$\displaystyle=\frac{1}{2}\left(\mathbf{f}^{E}_{\alpha\nu}(\omega_{\nu})+i{\mathbf{f}^{E}}^{H}_{\alpha\nu}(\omega_{\nu})\right),$
(51) $\displaystyle\mathbf{Q}^{E}_{2_{\alpha\nu}}$
$\displaystyle=\frac{1}{2}\Big{(}\mathbf{F}_{\alpha\nu}^{E}(\omega_{\nu})-i{\mathbf{F}^{E}}^{H}_{\alpha\nu}(\omega_{\nu})+(\alpha\leftrightarrow\nu)^{*}\Big{)},$
(52)

where $(\alpha\leftrightarrow\nu)^{*}$ notation means that the labels $\alpha$
and $\nu$ are to be interchanged, and the resulting expression is to be
complex-conjuagted. The functions $\mathbf{f}^{E}_{\alpha\nu}(\omega)$ and
$\mathbf{F}^{E}_{\alpha\nu}(\omega)$ are

$\displaystyle\mathbf{f}^{E}_{\alpha\nu}(\omega)=\sum_{\ell=1}^{N}\Phi_{\ell\alpha}^{*}\Phi_{\ell\nu}\mathfrak{J}_{\ell}(\omega),$
$\displaystyle\mathbf{F}^{E}_{\alpha\nu}(\omega)=\sum_{\ell=1}^{N}\Phi_{\ell\alpha}^{*}\Phi_{\ell\nu}\mathfrak{J}_{\ell}(\omega)\mathfrak{n}_{B_{\ell}}(\omega),$
(53)

and ${\mathbf{f}^{E}}^{H}_{\alpha\nu}(\omega)$ and
${\mathbf{F}^{E}}^{H}_{\alpha\nu}(\omega)$ are the corresponding Hilbert
transforms. This form of the Lyapunov equation can be used to find
perturbative solutions. In the following, we will assume that the NESS is
unique, and there is no degeneracy in the system. Results can be easily
generalized to cases without these approximations. In particular, we use the
condition

$\displaystyle|\omega_{\alpha}-\omega_{\nu}|\gg\epsilon^{2}\left|\mathbf{v}^{E}_{\alpha\alpha}+\mathbf{v}^{E*}_{\nu\nu}\right|,~{}~{}~{}\forall~{}~{}\alpha\neq\nu.$
(54)

To this end, first we define

$\displaystyle
w_{\alpha\nu}=i(\omega_{\alpha}-\omega_{\nu})+\epsilon^{2}\left(\mathbf{v}^{E}_{\alpha\alpha}+\mathbf{v}^{E*}_{\nu\nu}\right).$
(55)

The perturbative solutions of Eq.(50) up to leading order in $\epsilon$ are
given by

$\displaystyle\mathbf{C}^{E}_{\alpha\alpha}(t)\simeq\mathbf{C}^{E}_{\alpha\alpha}(0)e^{-2\epsilon^{2}\mathbf{f}^{E}_{\alpha\alpha}(\omega_{\alpha})t}$
$\displaystyle+\frac{\mathbf{F}^{E}_{\alpha\alpha}(\omega_{\alpha})}{\mathbf{f}^{E}_{\alpha\alpha}(\omega_{\alpha})}(1-e^{-2\epsilon^{2}\mathbf{f}^{E}_{\alpha\alpha}(\omega_{\alpha})t}),$
(56)
$\displaystyle\mathbf{C}^{E}_{\alpha\nu}(t)\simeq\mathbf{C}^{E}_{\alpha\nu}(0)e^{-w_{\alpha\nu}t}-\frac{i\epsilon^{2}\mathbf{Q}^{E}_{2_{\alpha\nu}}}{\omega_{\alpha}-\omega_{\nu}}(1-e^{-w_{\alpha\nu}t})$
$\displaystyle+\frac{i\epsilon^{2}}{\omega_{\alpha}-\omega_{\nu}}\Big{[}\mathbf{v}^{E*}_{\nu\alpha}\mathbf{C}^{E}_{\alpha\alpha}(0)(e^{-2\epsilon^{2}\mathbf{f}^{E}_{\alpha\alpha}(\omega_{\alpha})t}-e^{-w_{\alpha\nu}t})$
$\displaystyle+\mathbf{v}^{E}_{\alpha\nu}\mathbf{C}^{E}_{\nu\nu}(0)(e^{-2\epsilon^{2}\mathbf{f}^{E}_{\nu\nu}(\omega_{\alpha})t}-e^{-w_{\alpha\nu}t})\Big{]}$
$\displaystyle+\frac{i\epsilon^{2}}{\omega_{\alpha}-\omega_{\nu}}\Big{[}\mathbf{v}^{E^{*}}_{\nu\alpha}\frac{\mathbf{F}^{E}_{\alpha\alpha}(\omega_{\alpha})}{\mathbf{f}^{E}_{\alpha\alpha}(\omega_{\alpha})}(1-e^{-2\epsilon^{2}\mathbf{f}^{E}_{\alpha\alpha}(\omega_{\alpha})t})$
$\displaystyle+\mathbf{v}^{E}_{\alpha\nu}\frac{\mathbf{F}^{E}_{\nu\nu}(\omega_{\nu})}{\mathbf{f}^{E}_{\nu\nu}(\omega_{\nu})}(1-e^{-2\epsilon^{2}\mathbf{f}^{E}_{\nu\nu}(\omega_{\nu})t})\Big{]}.$
(57)

Real part of $w_{\alpha\nu}$ is
$\mathbf{f}^{E}_{\alpha\alpha}(\omega_{\alpha})+\mathbf{f}^{E}_{\nu\nu}(\omega_{\nu})$,
$\mathbf{f}^{E}_{\alpha\alpha}(\omega_{\alpha})>0$ by construction. Thus, for
time
$t\gg\textrm{max}\\{[\epsilon^{2}\mathbf{f}^{E}_{\alpha\alpha}(\omega_{\alpha})]^{-1}\\}$,
the steady state is reached. The perturbative results for the steady state are
given by

$\displaystyle\mathbf{C}^{E}_{\alpha\alpha}(\infty)=\frac{\mathbf{F}^{E}_{\alpha\alpha}(\omega_{\alpha})}{\mathbf{f}^{E}_{\alpha\alpha}(\omega_{\alpha})}=\frac{\sum_{\ell}|\Phi_{\ell\alpha}|^{2}\mathfrak{J}_{\ell}(\omega_{\alpha})\mathfrak{n}_{\ell}(\omega_{\alpha})}{\sum_{\ell}|\Phi_{\ell\alpha}|^{2}\mathfrak{J}_{\ell}(\omega_{\alpha})}$
(58)
$\displaystyle\mathbf{C}^{E}_{\alpha\nu}(\infty)=\frac{i\epsilon^{2}}{\omega_{\alpha}-\omega_{\nu}}\Big{[}\frac{\mathbf{v}^{E*}_{\nu\alpha}\mathbf{F}^{E}_{\alpha\alpha}(\omega_{\alpha})}{\mathbf{f}^{E}_{\alpha\alpha}(\omega_{\alpha})}$
$\displaystyle\qquad\qquad\qquad\qquad+\frac{\mathbf{v}^{E}_{\alpha\nu}\mathbf{F}^{E}_{\nu\nu}(\omega_{\nu})}{\mathbf{f}^{E}_{\nu\nu}(\omega_{\nu})}-\mathbf{Q}^{E}_{2_{\alpha\nu}}\Big{]}.$
(59)

We clearly see that the leading order diagonal elements in the single-particle
eigenbasis of the system are $O(1)$, while the leading-order off-diagonal
elements are $O(\epsilon^{2})$, as already mentioned before. Given the single-
particle eigenvalues and eigenvectors of the system, the equations in this
subsection are almost closed-form solutions for dynamics and NESS of our very
general set-up under Born-Markov approximation and Eq.(54). Thus, we see that,
in this regime, all results for both the transient dynamics and the NESS can
be written in terms of the single-particle eigenvalues and eigenvectors of the
system. This simplifies the extremely complicated problem of the system in
arbitrary dimension connected to an arbitrary number of baths at different
temperatures and chemical potentials, to the same level as obtaining the
dynamics of the isolated system in absence of the baths. This is a very
significant simplification. Moreover, as we will see later in Sec. IX.2, if
the system is one-dimensional with nearest neighbour hopping, further
simplification is possible and a very general and insightful expression for
current at NESS in terms of single-particle eigenvalues and eigenvectors can
be obtained. Below, we use Eqs.(58), (59) to discuss thermalization in the
arbitrary dimensions and geometry setting.

###  VII.3 Thermalization

The Eqs.(58), (59) reveal the very important physics of thermalization. In
equilibrium all baths have same temperatures and chemical potentials, i.e, the
Bose or Fermi distributions of all the baths are exactly the same,
$\mathfrak{n}_{\ell}(\omega)\rightarrow\mathfrak{n}(\omega)$. Thus, from
Eq.(58),
$\mathbf{C}^{E}_{\alpha\alpha}(\infty)=\mathfrak{n}(\omega_{\alpha})$. So we
get the non-trivial and physically important result that

$\displaystyle\textrm{in equilibrium},~{}~{}~{}\lim_{\epsilon\rightarrow
0}\left(\lim_{t\rightarrow\infty}\hat{\rho}(t)\right)=\frac{e^{-\beta(\hat{\mathcal{H}}_{S}-\mu\hat{N}_{S})}}{Tr\left(e^{-\beta(\hat{\mathcal{H}}_{S}-\mu\hat{N}_{S})}\right)},$
(60)

where the order of limits cannot be changed, and
$\hat{N}_{S}=\sum_{\ell=1}^{N}\hat{c}_{\ell}^{\dagger}\hat{c}_{\ell}$. Note
that the $\epsilon\rightarrow 0$ limit is consistent with Eq.(54). Although
Eq.(58) is written for the case where each site is attached to the bath, as
mentioned before, it is possible to send an arbitrary number of system-bath
couplings to zero to obtain an arbitrary distribution for location of baths.
Even if a bath at one site is kept connected, while all others are
disconnected, Eq.(58) shows that Eq.(60) holds. To appreciate the non-
triviality of the result, remember that we are working with an extremely
general system in arbitrary lattice and geometry with arbitrary number of site
attached to baths. What we showed above is that even if a thermal bath is
attached to one site of such a system, all the modes attached to it will reach
thermal equilibrium, and if the NESS is unique, the system will thermalize in
the sense of Eq.(60), irrespective of any further details of the system.

Away from equilibrium, when the temperatures and chemical potentials of the
baths are different, there will be non-zero current in NESS. For systems with
time-reversal symmetry, $\mathbf{H}$ will be real symmetric and consequently,
$\Phi$ can be chosen to be real orthogonal. In this case, the current in NESS
depends directly on the imaginary part of the off-diagonal elements
$\mathbf{C}^{E}_{\alpha\nu}$. After some algebra, the imaginary part of
$\mathbf{C}^{E}_{\alpha\nu}$ can be explicitly written as

$\displaystyle\textrm{Im}\left(\mathbf{C}^{E}_{\alpha\nu}(\infty)\right)=\frac{\epsilon^{2}}{\omega_{\alpha}-\omega_{\nu}}\Big{[}$
$\displaystyle\frac{\sum_{\ell,m}\Phi_{m\alpha}^{2}\Phi_{\ell\alpha}\Phi_{\ell\nu}\mathfrak{J}_{\ell}(\omega_{\alpha})\mathfrak{J}_{m}(\omega_{\alpha})\left(\mathfrak{n}_{m}(\omega_{\alpha})-\mathfrak{n}_{\ell}(\omega_{\alpha})\right)}{\sum_{\ell}^{\prime}\Phi_{\ell\alpha}^{2}\mathfrak{J}_{\ell}(\omega_{\alpha})}$
$\displaystyle+(\alpha\leftrightarrow\nu)\Big{]}.$ (61)

In equilibrium,
$\mathfrak{n}_{\ell}(\omega_{\alpha})=\mathfrak{n}_{m}(\omega_{\alpha})=\mathfrak{n}(\omega_{\alpha})$,
so Im$\left(\mathbf{C}^{E}_{\alpha\nu}(\infty)\right)=0$, which is consistent
with the fact that there is no steady state current in equilibrium for systems
with time-reversal symmetry.

Until now, we have discussed only equal time correlations of the system. The
understanding in terms of equation of motion gives us a natural way to discuss
two-time correlation functions and regression formulas, which we discuss in
the next section.

Approximation level | First Markov, | Level-I second Markov, | Level-II second Markov,
---|---|---|---
| Eqs.(IV.1.2), (33) | Eqs.(35),(36) | Eqs.(40), (IV.2.2)
Assumption | $t\gg\tau_{B_{1}}$ | $t\gg\tau_{B_{1}},\tau_{B_{2}}$ | $t\gg\tau_{B_{1}},\tau_{B_{2}}$
| see Eq.(22) for $\tau_{B_{1}}$ | see Eq.(IV.2.1) for $\tau_{B_{2}}$ | \+ weak system-bath coupling
Weak system-bath coupling | Holds | Holds | Holds
Strong coupling wide-band limit | Holds and exact | Holds, exact at NESS | Does not hold
fermionic case [Eq.(18)] |  |  |
Positivity of Lyapunov equation | Preserved at all times | Preserved at NESS | Violated below accuracy level
| Eq.(45) | Eq.(44) | Eq.(42)
Associated QME | time-dependent rates | time-independent rates, | Redfield equation
|  | may not be completely positive, |
|  | but always gives positive NESS |
Additivity | Additive | Additive | Additive
(both Lyapunov equation and QME) |  |  |
Generalized regression formula | Can be derived, | Can be derived, | Cannot be derived
from equations of motion | Eqs.(65), (66) | Eqs.(70), (71) |
Table 1:  Summary of the three level of approximations for the three slightly
different continuous-time Lyapunov equations. The homogeneous part of the
Lyapunov equations, which is associated with the non-Hermitian Hamiltonian
governing the dynamics (see Eqs.(23), (24), (25), (30)), is the same in all
cases. The approximations change the inhomogeneous part of the Lyapunov
equations which embody the quantum and thermal fluctuations. The level-II
second Markov approximation is identical to the Redfield equation and has the
same, controlled, accuracy and positivity issues. It can be used to obtain
semi-analytical results (Eqs.(56), (57), (58), (59)) in the regime of very
small system-bath coupling (Eq.(54)), which give open system dynamics in terms
of single-particle eigenvalues and eigenvectors of the system.

##  VIII Two-time correlations and regression formulas

The matrix $\mathbf{C}(t)$ gives equal time correlations. In this section, we
move to calculating two-time correlations

$\displaystyle\mathbf{C}_{\ell
m}(t,t^{\prime})=\langle\hat{c}_{\ell}^{\dagger}(t)\hat{c}_{m}(t^{\prime})\rangle.$
(62)

These quantities are directly related to NEGF Haug and Jauho (2008). For
calculating two-time correlations from QMEs one has to resort to the
corresponding quantum regression formulas, which often relies on assumptions
with their own set of issues Talkner (1986); Ford and O’Connell (1996);
Guarnieri  _et al._ (2014). On the other hand, having an operator equation of
motion description makes obtaining two-time correlations completely
straightforward.

###  VIII.1 First Markov approximation

After first Markov approximation, calculating two-time correlation functions
from Eq.(27) gives (see Appendix A)

$\displaystyle\mathbf{C}(t,t^{\prime})=e^{-\mathcal{G}t}\mathbf{C}(0)e^{-\mathcal{G^{\dagger}}t^{\prime}}$
$\displaystyle+\epsilon^{2}\int\frac{d\omega}{2\pi}\left(\frac{1-e^{-(\mathcal{G}+i\omega\mathbb{I})t}}{\mathcal{G}+i\omega\mathbb{I}}\mathbf{F}(\omega)\frac{1-e^{-(\mathcal{G}^{\dagger}-i\omega\mathbb{I})t^{\prime}}}{\mathcal{G}^{\dagger}-i\omega\mathbb{I}}\right),$
(63)

which is exactly same as Eq.(IV.1.2) except with the two time arguments being
different. The regression formula deals with evolution of
$\mathbf{C}(t+\tau,t)$ as a function of $\tau$, $\tau>0$. If we naively used
the quantum regression relation on the QMEs Breuer and Petruccione (2007), we
would get,

$\displaystyle\frac{d\mathbf{C}(t+\tau,t)}{d\tau}=-\mathcal{G}\mathbf{C}(t+\tau,t).$
(64)

This would be same for all three QMEs Eq.(45), (44), (42). However, we can
derive the differential equation for evolution of $\mathbf{C}(t+\tau,t)$ as a
function of $\tau$ directly from equations of motion, exactly as Eq.(IV.1.2)
was obtained, without further approximations. This gives,

$\displaystyle\frac{d\mathbf{C}(t+\tau,t)}{d\tau}=-\mathcal{G}\mathbf{C}(t+\tau,t)+i\epsilon^{2}\mathbf{C}_{\xi}(t+\tau,t).$
(65)

Here, $\mathbf{C}_{\xi_{\ell
m}}(t+\tau,t)=\langle\hat{\xi}_{\ell}^{\dagger}(t+\tau)\hat{c}_{m}(t)\rangle$,
which can be calculated from Eq.(27) as

$\displaystyle\mathbf{C}_{\xi}(t+\tau,t)=-i\int_{0}^{t}dt^{\prime}\int\frac{d\omega}{2\pi}\mathbf{F}(\omega)e^{i\omega(t^{\prime}+\tau)}e^{-\mathcal{G}^{\dagger}t^{\prime}}.$
(66)

Thus, naively using quantum regression at the QMEs would not have given the
inhomogeneous part of Eq.(65). Such deviations from the quantum regression
formula at the level of QMEs has been associated with non-Markovian behavior
Guarnieri  _et al._ (2014). We call equations of the form of Eq.(65)
generalized regression formulas. The Eq.(65) is under the first Markov
approximation, and requires $t+\tau\gg\tau_{B_{1}}$ with $\tau_{B_{1}}$
defined in Eq.(22). So, if $\tau\gg\tau_{B_{1}}$, only then $t=0$ is allowed,
otherwise not. If we put $t=0$, the inhomogeneous part becomes zero. Thus,
only in this case the equation becomes same as the one expected by using
quantum regression from the QME.

The formal solution of Eq.(65) can be written as

$\displaystyle\mathbf{C}(t+\tau,t)=e^{-\mathcal{G}\tau}\mathbf{C}(t)$
$\displaystyle+\epsilon^{2}\int\frac{d\omega}{2\pi}\left(\frac{e^{i\omega\tau}-e^{-\mathcal{G}\tau}}{\mathcal{G}+i\omega\mathbb{I}}\mathbf{F}(\omega)\frac{1-e^{-(\mathcal{G}^{\dagger}-i\omega\mathbb{I})t}}{\mathcal{G}^{\dagger}-i\omega\mathbb{I}}\right).$
(67)

If real parts of eigenvalues of $\mathcal{G}$ are positive, i.e, if there is
unique NESS, then noting that $\mathbf{C}(\infty)$ is given by Eq.(31), we can
write the above equation in the following suggestive form,

$\displaystyle\mathbf{C}(t+\tau,t)=e^{-\mathcal{G}\tau}\left[\mathbf{C}(t)-\mathbf{C}(\infty)\right]$
$\displaystyle+\epsilon^{2}\int\frac{d\omega}{2\pi}e^{i\omega\tau}\left(\frac{1}{\mathcal{G}+i\omega\mathbb{I}}\mathbf{F}(\omega)\frac{1}{\mathcal{G}^{\dagger}-i\omega\mathbb{I}}\right)$
(68)
$\displaystyle-\epsilon^{2}\int\frac{d\omega}{2\pi}\left(\frac{e^{i\omega\tau}-e^{-\mathcal{G}\tau}}{\mathcal{G}+i\omega\mathbb{I}}\mathbf{F}(\omega)\frac{e^{-(\mathcal{G}^{\dagger}-i\omega\mathbb{I})t}}{\mathcal{G}^{\dagger}-i\omega\mathbb{I}}\right).$

So the contribution from the term that would be obtained via naive application
of quantum regression at the level of the QME actually goes to zero at NESS.
We get the following result for two-time correlations at NESS,

$\displaystyle\lim_{t\rightarrow\infty}\mathbf{C}(t+\tau,t)=\epsilon^{2}\int\frac{d\omega}{2\pi}e^{i\omega\tau}\left(\frac{1}{\mathcal{G}+i\omega\mathbb{I}}\mathbf{F}(\omega)\frac{1}{\mathcal{G}^{\dagger}-i\omega\mathbb{I}}\right).$
(69)

This is exactly the same as would be obtained by solving Eq.(23) by a Fourier
transform and then obtaining the two-time correlation functions. Clearly,
quantum regression at the level of the QME would not give this result.

###  VIII.2 Second Markov approximation

The level-I second Markov approximation would assume $t+\tau\gg\tau_{B_{2}}$,
with $\tau_{B_{2}}$ defined in Eq.(IV.2.1), and take $t\rightarrow\infty$ in
Eq.(66). This yields,

$\displaystyle\frac{d\mathbf{C}(t+\tau,t)}{d\tau}=-\mathcal{G}\mathbf{C}(t+\tau,t)+i\epsilon^{2}\widetilde{\mathbf{C}}_{\xi}(\tau),$
(70)

where

$\displaystyle\widetilde{\mathbf{C}}_{\xi}(\tau)=-i\int\frac{d\omega}{2\pi}e^{i\omega\tau}\mathbf{F}(\omega)\frac{1}{\mathcal{G}^{\dagger}-i\omega\mathbb{I}}.$
(71)

The formal solution at for the above equations is

$\displaystyle\mathbf{C}(t+\tau,t)=e^{-\mathcal{G}\tau}\mathbf{C}(t)$
$\displaystyle+\epsilon^{2}\int\frac{d\omega}{2\pi}\left(\frac{e^{i\omega\tau}-e^{-\mathcal{G}\tau}}{\mathcal{G}+i\omega\mathbb{I}}\mathbf{F}(\omega)\frac{1}{\mathcal{G}^{\dagger}-i\omega\mathbb{I}}\right).$
(72)

If real parts of eigenvalues of $\mathcal{G}$ are positive, the above
expression can be written in the same form as Eq.(VIII.1) except the last line
being set to zero. So, in the long time limit, the expression again reduces to
Eq.(69).

A level-II second Markov approximation would further approximate Eq.(71) and
carry out an analogous procedure of Eq.(IV.2.2). However, carrying out this
approximation gives inconsistent results for two-time correlations. This is
because, $\widetilde{\mathbf{C}}_{\xi}(\tau)$ typically would decay with
$\tau$. Setting $\epsilon=0$ in the expression with finite $\tau$ would not
give this behavior and therefore would be inconsistent. As a result, the
level-II second Markov approximation cannot be performed here. So, to obtain
two-time correlations corresponding to the Redfield equation, Eq.(42), one has
to resort to the regression relation at the level of QMEs Eq.(64), which is
likely to impose additional restrictions Guarnieri  _et al._ (2014). In
particular, we see from Eq.(VIII.2) and Eq.(31) that only if, due to some
additional approximations over and above the level-I second Markov
approximation, the following relation approximately holds

$\displaystyle\int\frac{d\omega}{2\pi}e^{i\omega\tau}\left(\frac{1}{\mathcal{G}+i\omega\mathbb{I}}\mathbf{F}(\omega)\frac{1}{\mathcal{G}^{\dagger}-i\omega\mathbb{I}}\right)\approx
e^{-\mathcal{G}\tau}\mathbf{C}(\infty),$ (73)

will Eq.(64) be satisfied at all times. Whether these additional
approximations are same as those required for deriving Redfield equation is
not clear and requires further work Khan  _et al._ (2021).

It is usually believed that if the NESS is given by a QME of Lindblad form,
the quantum regression formula must be valid. However, as we will see below in
Sec. IX.1, the simple example of the resonant level model with wide-band baths
shows that it is not so. Although its NESS can be obtained from a Lindblad
equation for arbitrary strengths of system-bath coupling, the exact two-time
correlations at NESS are given by the analog of Eq.(69), which is different
from what is obtained via analog of Eq.(64). Only if further approximations,
like weak system-bath coupling, or high temperatures, are made, can the two
results be reduced to the same form. So, clearly, regression relations require
further set of approximations, which can sometimes be different from those
required to obtain Lindblad descriptions.

The formulas given in this section allow calculation of the
$\mathbf{C}(t+\tau,t)$, $\tau>0$, knowing $\mathbf{C}(t)$. If
$\mathbf{C}(t,t+\tau)$ is desired instead, it can be obtained by simply noting
that $\mathbf{C}^{\dagger}(t+\tau,t)=\mathbf{C}(t,t+\tau)$.

This concludes our main general results. A summary of our results is given in
Table 1. In the next section, we explicitly discuss the two examples we have
referred to before, viz., the resonant level model, and a one-dimensional
system with nearest neighbour hopping.

##  IX Insightful examples

###  IX.1 Resonant level model

In this section, we work out the extremely simple example of a resonant level
model, which is a single fermionic site, coupled to two wide-band fermionic
reservoirs at different temperatures and chemical potentials, see Fig. 2. The
Hamiltonian is given by

$\displaystyle\hat{\mathcal{H}}_{SB}=\varepsilon\hat{c}^{\dagger}\hat{c},$
$\displaystyle\hat{\mathcal{H}}_{SB}=\sum_{\ell=L,R}\hat{\mathcal{H}}_{SB_{\ell}},~{}\hat{\mathcal{H}}_{SB_{\ell}}=\sum_{r=1}^{\infty}\left(\kappa_{rl}\hat{c}^{\dagger}\hat{B}_{r\ell}+\kappa_{rl}^{*}\hat{B}_{r\ell}^{\dagger}\hat{c}\right)$
$\displaystyle\hat{\mathcal{H}}_{B}=\sum_{\ell=L,R}\hat{\mathcal{H}}_{B_{\ell}},~{}\hat{\mathcal{H}}_{B_{\ell}}=\sum_{r=1}^{\infty}\Omega_{rl}\hat{B}_{r\ell}^{\dagger}\hat{B}_{r\ell}.$
(74)

In above, $\ell=L,R$ labels the left and the right baths attached to the
system site. The system-bath coupling and the baths are of the same form as
Eq.(III) except the following differences. Since we will consider the wide-
band limit fermionic case,

$\displaystyle\mathfrak{J}_{\ell}(\omega)=\Gamma_{\ell},~{}~{}\ell=L,R,$ (75)

we have set $\epsilon=1$. Further, unlike in previous sections, there are two
baths attached to the same single site. As mentioned before, since all the
equations are additive this hardly complicates anything. We will specifically
focus on the level-I second Markov approximation. For each bath, we have an
analog of $\mathbf{v}$ and $\mathbf{Q}_{1}$ matrices, and they are just summed
over.

Figure 2: Schematic of the resonant level model.

$\displaystyle\mathbf{v}=\mathbf{v}^{(L)}+\mathbf{v}^{(R)},~{}~{}\mathbf{Q}_{1}=\mathbf{Q}_{1}^{(L)}+\mathbf{Q}_{1}^{(R)},$
(76)

where $\mathbf{v}^{(L,R)}$ and $\mathbf{Q}_{1}^{(L,R)}$ are calculated by
using Eq.(26) and Eq.(36). For a single site, they are just scalars,
$\mathbf{v}^{(L,R)}$ being given by

$\displaystyle\mathbf{v}^{(\ell)}=\frac{\Gamma_{\ell}}{2},~{}~{}\ell=L,R,$
(77)

see Eq.(26). The non-Hermitian Hamiltonian (Eq.(24)) and $\mathcal{G}$
(Eq.(30)) are then

$\displaystyle\mathbf{H}_{\rm
NH}=\varepsilon-i\frac{\Gamma_{L}+\Gamma_{R}}{2},~{}\mathcal{G}=-i\varepsilon+\frac{\Gamma_{L}+\Gamma_{R}}{2}.$
(78)

Using Eq.(36), we get the following expressions for $\mathbf{Q}_{2}^{(L,R)}$,

$\displaystyle\mathbf{Q}_{1}^{(\ell)}=\int\frac{d\omega}{2\pi}\frac{(\Gamma_{L}+\Gamma_{R})\Gamma_{\ell}\mathfrak{n}_{B_{\ell}}(\omega)}{(\omega-\varepsilon)^{2}+\left(\frac{\Gamma_{L}+\Gamma_{R}}{2}\right)^{2}}~{}~{}\ell=L,R,$
(79)

where $\mathfrak{n}_{B_{L}}(\omega)$ ($\mathfrak{n}_{B_{R}}(\omega)$) is the
Fermi distribution corresponding to the left (right) bath. For a single-site,
the correlation matrix is also a scalar,
$\langle\hat{n}\rangle=\langle\hat{c}^{\dagger}\hat{c}\rangle$ being the only
element. The Lyapunov equation then becomes

$\displaystyle\frac{d\langle\hat{n}\rangle}{dt}=-(\Gamma_{L}+\Gamma_{R})\langle\hat{n}\rangle+\mathbf{Q}_{1}^{(L)}+\mathbf{Q}_{1}^{(R)}.$
(80)

The corresponding QME is

$\displaystyle\frac{\partial\hat{\rho}}{\partial t}$
$\displaystyle=i[\hat{\rho},\varepsilon\hat{n}]+\sum_{\ell=L,R}\Big{[}\big{(}\Gamma_{\ell}-\mathbf{Q}_{1}^{(\ell)}\big{)}\left(\hat{c}\rho\hat{c}^{\dagger}-\frac{1}{2}\\{\hat{c}^{\dagger}\hat{c},\hat{\rho}\\}\right)$
$\displaystyle+\mathbf{Q}_{1}^{(\ell)}\left(\hat{c}^{\dagger}\rho\hat{c}-\frac{1}{2}\\{\hat{c}\hat{c}^{\dagger},\hat{\rho}\\}\right)\Big{]}.$
(81)

This is clearly additive and of Lindblad form. It can be checked explicitly
that above QME gives Eq.(80).

The NESS is obtained by setting the left-hand-side of Eq.(80) to zero, which,
gives

$\displaystyle\langle\hat{n}\rangle=\frac{\mathbf{Q}_{1}^{(L)}+\mathbf{Q}_{1}^{(R)}}{\Gamma_{L}+\Gamma_{R}}.$
(82)

This, along with Eq.(79), gives,

$\displaystyle\langle\hat{n}\rangle=\int\frac{d\omega}{2\pi}\frac{\Gamma_{L}\mathfrak{n}_{B_{L}}(\omega)+\Gamma_{R}\mathfrak{n}_{B_{R}}(\omega)}{(\omega-\varepsilon)^{2}+\left(\frac{\Gamma_{L}+\Gamma_{R}}{2}\right)^{2}}.$
(83)

This is the well-known correct expression for occupation at NESS for the
resonant level model Haug and Jauho (2008).

Let us now calculate the current at NESS. The Eq.(80) is essentially a
continuity equation for occupation, and currents from the left and the right
baths, $I_{L}$ and $I_{R}$, can be identified as
$I_{\ell}=-\Gamma_{\ell}\langle\hat{n}\rangle+\mathbf{Q}_{1}^{(\ell)},~{}~{}\ell=L,R.$
At NESS, since the left-hand-side of Eq.(80) is zero, we have
$I_{L}=-I_{R}=I$. Using Eq.(82), the expression for current $I$ is

$\displaystyle
I=\frac{\Gamma_{R}\mathbf{Q}_{1}^{(L)}-\Gamma_{L}\mathbf{Q}_{1}^{(R)}}{\Gamma_{L}+\Gamma_{R}}.$
(84)

Substituting the expression for $\mathbf{Q}_{1}^{(\ell)}$ from Eq.(79) gives

$\displaystyle
I=\int\frac{d\omega}{2\pi}\frac{\Gamma_{L}\Gamma_{R}[\mathfrak{n}_{B_{L}}(\omega)-\mathfrak{n}_{B_{R}}(\omega)]}{(\omega-\varepsilon)^{2}+\left(\frac{\Gamma_{L}+\Gamma_{R}}{2}\right)^{2}},$
(85)

which is the well-known correct expression for current in NESS resonant level
model with wide-band baths Haug and Jauho (2008). Note that these results are
true irrespective of strength of system-bath coupling. Therefore the perfectly
additive Lindblad equation, Eq.(IX.1) gives completely exact results for NESS
of the resonant level model (Refer to discussion in Sec. VI.3). Note that, the
system-bath coupling Hamiltonians neither commute with each other, nor with
the system Hamiltonian. So it does not satisfy the sufficient condition for
accurate additive dynamics given in Ref. Kołodyński  _et al._ (2018).

It is interesting to note that even though the NESS is given exactly by a QME
of perfectly Lindblad form, Eq.(IX.1), the two-time correlations at NESS
cannot be obtained by using quantum regression relation at the level of the
QME (Refer to discussions in Sec. VIII). For simplicity, here we will assume
$\Gamma_{L}=\Gamma_{R}=\Gamma$. Let us define the notation

$\displaystyle\langle\hat{c}^{\dagger}(\tau)\hat{c}\rangle=\lim_{t\rightarrow\infty}\langle\hat{c}^{\dagger}(t+\tau)\hat{c}(t)\rangle.$
(86)

Using quantum regression at the level of the QME gives

$\displaystyle\langle\hat{c}^{\dagger}(\tau)\hat{c}\rangle=e^{(i\varepsilon-\Gamma)\tau}\langle\hat{n}\rangle$
$\displaystyle=e^{(i\varepsilon-\Gamma)\tau}\int\frac{d\omega}{2\pi}\frac{\Gamma\left[\mathfrak{n}_{B_{L}}(\omega)+\mathfrak{n}_{B_{R}}(\omega)\right]}{(\omega-\varepsilon)^{2}+\Gamma^{2}}.$
(87)

On the other hand, using Eq.(69) gives

$\displaystyle\langle\hat{c}^{\dagger}(\tau)\hat{c}\rangle=\int\frac{d\omega}{2\pi}e^{i\omega\tau}\frac{\Gamma\left[\mathfrak{n}_{B_{L}}(\omega)+\mathfrak{n}_{B_{R}}(\omega)\right]}{(\omega-\varepsilon)^{2}+\Gamma^{2}},$
(88)

which is the correct expression obtained from NEGF Haug and Jauho (2008).
Clearly this is not the same as Eq.(IX.1). So, despite the NESS being exactly
given by a Lindblad equation, the regression relation does not hold.

The Eq.(88) has the form of Fourier transform of sum of Fermi distributions
weighted by a Lorentzian function. This can approximately reduce to Eq.(IX.1)
if the Fermi distributions are reasonably flat within the width of the
Lorentzian. Since Fermi distributions vary in a scale of $1/\beta$, this can
be satisfied if

$\displaystyle\beta_{\ell}\Gamma\ll 1,~{}~{}\ell=L,R.$ (89)

The above condition can either be satisfied at weak coupling, or at high
temperatures. In either case, it will be justified to approximately set
$\mathfrak{n}_{L,R}(\omega)\approx\mathfrak{n}_{L,R}(\varepsilon)$ in the
integrations in Eqs.(IX.1) and (88). Due to properties of Lorentzian
functions, both integrations will now yield

$\displaystyle\langle\hat{c}^{\dagger}(\tau)\hat{c}\rangle\approx
e^{(i\varepsilon-\Gamma)\tau}\frac{\mathfrak{n}_{B_{L}}(\varepsilon)+\mathfrak{n}_{B_{R}}(\varepsilon)}{2}.$
(90)

This same result could be obtained by using quantum regression on the Redfield
equation (which will also be of Lindblad form in this case), deriving which
will explicitly require weak system-bath coupling approximation. Clearly,
validity of quantum regression requires additional approximations.

###  IX.2 Simple expression for current and dimensionless conductance in one-
dimensional nearest neighbour systems

Figure 3: (a) General linear one-dimensional nearest neighbour system coupled
to multiple thermal baths. (b) The same system but with baths attached at only
the two ends.

Here we consider a one-dimensional chain with nearest-neighbour hopping and
derive a simple and insightful expression for the particle current using the
result in Eq.(59). The system Hamiltonian we will be considering now is given
by

$\displaystyle\hat{\mathcal{H}}_{S}=\sum_{\ell=1}^{N}\varepsilon_{\ell}\hat{c}_{\ell}^{\dagger}\hat{c}_{\ell}+\sum_{\ell=1}^{N-1}g_{\ell}\left(\hat{c}_{\ell}^{\dagger}\hat{c}_{\ell+1}+\hat{c}_{\ell+1}^{\dagger}\hat{c}_{\ell}\right).$
(91)

In other words, $\mathbf{H}$ in Eq.(5) is now the tridiagonal matrix

$\displaystyle\mathbf{H}=\begin{bmatrix}\varepsilon_{1}&g_{1}&0&\ldots&\ldots\\\
g_{1}&\varepsilon_{2}&g_{2}&0&\ldots\\\ 0&g_{2}&\varepsilon_{3}&g_{3}&0\\\
\vdots&\ddots&\ddots&\ddots&\ddots\\\
\ldots&\ldots&0&g_{N-1}&\varepsilon_{N}\end{bmatrix}.$ (92)

This system is coupled at an arbitrary number of sites to different baths
which are all at different temperatures and chemical potentials, see Fig.
3(a). We will find an expression for particle current at a system bond in
terms of the single-particle eigenvalues and eigenfunctions of the system. The
particle current $I_{p}$ at the $p$th bond of the one-dimensional system is
given by

$\displaystyle
I_{p}=g_{p}\textrm{Im}(\langle\hat{c}_{p}^{\dagger}\hat{c}_{p+1}\rangle)=g_{p}\sum_{\alpha,\nu=1}^{N}\Phi_{p\alpha}\Phi_{p\nu}\textrm{Im}(\mathbf{C}^{E}_{\alpha\nu}),$
(93)

where ${\rm Im}(A)$ refers to imaginary part of $A$. Using Eq.(59) and
simplifying utilizing some properties of a tridiagonal matrix, we obtain the
following expression for current at the $p$th bond (see Appendix C),

$\displaystyle
I_{p}=\epsilon^{2}\sum_{\alpha=1}^{N}\frac{1}{\sum_{\ell}^{\prime}\Phi_{\ell\alpha}^{2}\mathfrak{J}_{\ell}(\omega_{\alpha})}\Big{[}$
$\displaystyle\sum_{\ell,m}^{\prime}\sum_{k=1}^{p}\delta_{\ell
k}\Phi_{m\alpha}^{2}\Phi_{\ell\alpha}^{2}\mathfrak{J}_{\ell}(\omega_{\alpha})\mathfrak{J}_{m}(\omega_{\alpha})\left[\mathfrak{n}_{\ell}(\omega_{\alpha})-\mathfrak{n}_{m}(\omega_{\alpha})\right]\Big{]}.$
(94)

Here $\sum^{\prime}$ refers to sum over the points where the baths are
attached. So, we see that current in the limit of very small system-bath
couplings (Eq.(54)) is governed by the amplitudes of the single-particle
eigenfunctions at the sites where the baths are attached.

Now let us further specialize to the case where the only two baths are
attached, which are at the first and last sites of the system, see Fig. 3(b).
In this case, the current should be independent of the bond where it is
calculated, $I_{p}=I$. Indeed it is so, as is seen from carrying out the sum
over $k$ in Eq.(IX.2), which gives

$\displaystyle
I=\epsilon^{2}\sum_{\alpha=1}^{N}\frac{\Phi_{1\alpha}^{2}\Phi_{N\alpha}^{2}\mathfrak{J}_{1}(\omega_{\alpha})\mathfrak{J}_{N}(\omega_{\alpha})\left[\mathfrak{n}_{1}(\omega_{\alpha})-\mathfrak{n}_{N}(\omega_{\alpha})\right]}{\Phi_{1\alpha}^{2}\mathfrak{J}_{1}(\omega_{\alpha})+\Phi_{N\alpha}^{2}\mathfrak{J}_{N}(\omega_{\alpha})}.$
(95)

The above expressions for current are valid for both bosonic and fermionic
cases. We can now obtain an insightful expression for particle conductance of
a fermionic system. Let set-up be fermionic with the two baths having same
temperature, $\beta_{1}=\beta_{N}=\beta$, but different chemical potentials,
$\mu_{1}=\mu+\Delta\mu,\mu_{N}=\mu$ and being described by wide-band baths
coupled at same strength,
$\mathfrak{J}_{1}(\omega_{\alpha})=\mathfrak{J}_{N}(\omega_{\alpha})=\Gamma$.
Then particle conductance is given by,

$\displaystyle G=$
$\displaystyle\lim_{\Delta\mu_{\rightarrow}0}\frac{dI}{d\mu}$
$\displaystyle=\epsilon^{2}\beta\Gamma\sum_{\alpha=1}^{N}\left[\frac{\Phi_{1\alpha}^{2}\Phi_{N\alpha}^{2}}{\Phi_{1\alpha}^{2}+\Phi_{N\alpha}^{2}}\mathfrak{n}^{2}(\omega_{\alpha})e^{\beta(\omega_{\alpha}-\mu)}\right]$
(96)

Note, here we have already imposed very small system-bath coupling condition
before, so we are not in the regime of Eq.(35), even though we are using wide-
band fermionic baths now. In the high temperature limit,
$\mathfrak{n}^{2}(\omega_{\alpha})e^{\beta(\omega_{\alpha}-\mu)}\simeq 1/4$.
So particle conductance of the fermionic system in wide-band and high
temperature limit is given by

$\displaystyle G=\epsilon^{2}\frac{\Gamma\beta}{4}W(1,N),~{}{\rm
where}~{}W(r,s)=\sum_{\alpha=1}^{N}\left[\frac{\Phi_{r\alpha}^{2}\Phi_{s\alpha}^{2}}{\Phi_{r\alpha}^{2}+\Phi_{s\alpha}^{2}}\right].$
(97)

Thus, in this limit, $W(1,N)$ is proportional to conductance. It can be termed
a dimensionless conductance. It is interesting to note that $W(1,N)$ depends
only on system eigenfunctions and is independent of the baths. It is
essentially an isolated system quantity, but to derive Eq.(97), we required to
consider an open system. This expression for particle conductance was used in
Ref. Purkayastha  _et al._ (2018), without proof, to explain the origin of
sub-diffusive scaling of conductance with system size at the critical point of
the Aubry-André-Harper model. The results were also checked against exact
calculations.

##  X Summary and outlook

Summary and distinguishing features from previous works — In this work, we
have derived the Lyapunov equation for describing the dynamics of number
conserving linear systems (quadratic Hamiltonians) in a lattice of arbitrary
dimension and geometry, coupled to an arbitrary number of macroscopic thermal
baths which can all be at different temperatures and chemical potentials.
Three slightly different forms of the Lyapunov equation are derived. Table 1
gives the summary of our results. The following points distinguish our work
from previous works involving Lyapunov equations. We have given a detailed,
systematic, controlled derivation starting from a fully microscopic
Hamiltonian (Hermitian) model of the system coupled with the baths. Unlike
most previous works Prosen (2008); Prosen and Žunkovič (2010); Žunkovič and
Prosen (2012); Koga and Yamamoto (2012); Nicacio  _et al._ (2015); Landi  _et
al._ (2013); Nicacio  _et al._ (2016); Santos  _et al._ (2017); Malouf  _et
al._ (2019); McDonald  _et al._ (2022); Bernal-García  _et al._ (2021)), we
have done so without referring to any QME. Exact approaches require inverse
Laplace transformation Tu and Zhang (2008); Jin  _et al._ (2010); Zhang  _et
al._ (2012); Martinez and Paz (2013); Yang  _et al._ (2013), inverting which
can become challenging depending on number of sites in the system and the
spectrum of the system. The Lyapunov equations greatly simplify the problem of
obtaining dynamics of such open quantum systems via bypassing the need for
Laplace transform. In particular, for fermionic reservoirs in the so-called
wide-band limit, our results do not require a weak system-bath coupling
approximation. Our microscopic derivation makes the validity regimes of the
Lyapunov equations clear. We have then found the associated QMEs which yields
the corresponding Lyapunov equations. For Gaussian initial states of the
system, two of the associated QMEs also allow us to resolve the positivity
problem of the Redfield equation. These two QMEs would have been difficult to
obtain otherwise. On the other hand, the third QME is the Redfield equation,
which shows the equivalence of QME and equation of motion approaches. In the
limit of very small coupling to the baths, we have found semi-analytical
results for the system which are written only in terms of the single-particle
eigenvalues and eigenvectors of the system. This therefore reduces the problem
of obtaining open system dynamics to the same difficulty level as solving the
isolated system dynamics. Finally, we have given the generalized regression
formulas from the correlation matrix, which allows calculation of all two-time
correlation functions for a Gaussian system. These formulas would not be
possible to derive by naively applying quantum regression at the level of the
QMEs. Finally, we have worked out two insightful examples which highlight
several features of our results.

Implications for non-Hermitian quantum physics — One of most immediate
implications of our work is that our microscopic derivation explicitly
separates the non-Hermitian Hamiltonian governing the dynamics of the system
from the quantum and thermal fluctuations due to the presence of the baths.
Therefore, it shows how dynamics governed by effective non-Hermitian
Hamiltonians can microscopically arise out of Hermitian quantum mechanics, and
gives an unified way to treat quantum and thermal fluctuations on such
systems. Since our formulation is very general, all kinds of non-Hermitian
Hamiltonians can be microscopically designed in this way, as long as there are
no sources of gain. The homogeneous part of the Lyapunov equation is
associated with the non-Hermitian Hamiltonian, while the inhomogeneous part of
the Lyapunov equation is associated with the quantum and thermal fluctuations.
Interestingly, the three levels of approximation only change the inhomogeneous
part. They therefore give various accuracy levels of treating the quantum and
thermal fluctuations. While Lyapunov equations from Lindblad QMEs have been
used to explore non-Hermitian physics in few works McDonald  _et al._ (2022);
Arkhipov  _et al._ (2021); Roccati  _et al._ (2021), our formulation goes
beyond the validity regimes of such Lindblad QMEs deriving which usually
require further approximations over the Redfield equation Breuer and
Petruccione (2007); Rivas and Huelga (2012); Tupkary  _et al._ (2022).

Sources of gain typically require a non-linear coupling with bath, which is
beyond the scope of the present paper. However, Lyapunov equations can be
obtained from microscopic considerations, completely out of Hermitian quantum
mechanics, in such cases also, see Purkayastha  _et al._ (2020b). But, such
Lyapunov equations, which give an effective linearized description, may be
unstable and may not be valid up to long times. Nevertheless, at least up to
some finite time (which can be estimated), our present work, in combination
with Purkayastha  _et al._ (2020b), suggests that systems described by all
kinds of effective non-Hermitian Hamiltonians can be obtained from standard
quantum mechanics, as long as quantum and thermal fluctuations are properly
accounted for via a Lyapunov equation. However, in presence of gain, the non-
linear coupling can make the state non-Gaussian. So, the Lyapunov equation in
presence of gain may not describe the full state of the system.

Possible implications for dissipative quantum many-body systems — Systems
governed by quadratic Hamiltonians form the starting point for much of our
understanding of physics, especially, in higher than one dimension. Much of
analytical techniques in physics deal with formulating sophisticated methods
to obtain corrections above such quadratic descriptions Mahan (2000); Haug and
Jauho (2008); Kamenev (2011). These techniques may be combined with our
general Lyapunov equations and regression formulas to treat many-body
interactions (i.e higher than quadratic terms). The simplest of these is the
Hartree-Fock mean field approximation, which is very often made for realistic
three-dimensional systems. At such mean-field level, the dynamics and NESS can
be obtained by solving the Lyapunov equation self-consistently. The Lyapunov
equations may therefore be used as a natural and simple way to include
dissipation into the mean-field description of realistic quantum many-body
systems. This also gives a microscopic meaning, consistent with standard
quantum mechanics, to band structures of non-Hermitian Hamiltonians, which
have been explored in various works (for example, Kawabata  _et al._ (2019);
Gong  _et al._ (2018); Wang  _et al._ (2021)).

Relevance to quantum thermodynamics of Gaussian systems — Lyapunov equations
have been recently used to introduce and describe concepts related to quantum
thermodynamics in Gaussian systems like Wigner entropy production Santos  _et
al._ (2017); Malouf  _et al._ (2019), and non-linear Onsager relations Salazar
and Landi (2020). Once again, this was done from the viewpoint of Markovian
Lindblad equations. Our results therefore suggest possible generalizations of
those results beyond their present validity regimes.

Further generalizations and future works — Here we have considered a number
conserving quadratic Hamiltonian. But most of our results can be readily
generalized to the case of number non-conserving systems, such as
superconductors and squeezed harmonic oscillators, using a slightly different
definition for the correlation matrix Ribeiro and Vieira (2015). They may also
be possible generalize to the case of athermal Gaussian baths, instead of
thermal baths. Another interesting direction is generalization of the Lyapunov
equation to the case where the temperatures and the chemical potentials of the
baths are time-dependent. At the level of the Redfield equation and the
associated Lyapunov equation, this has already been achieved Purkayastha and
Dubi (2017). It will be interesting to see if a generalization beyond the
validity regime of Redfield equation would be possible. Moreover, the entire
rich mathematical understanding of Lyapunov equations Gajić and Qureshi
(1995); Sun and Yong (2020) can now be carried over to open quantum systems
and non-Hermitian physics. Investigations in these directions, as well as in
the direction of quantum thermodynamics, will be taken up in future works.
Overall, we find it fascinating that an equation used in daily life for
control of macroscopic objects Gajić and Qureshi (1995) can be shown to play
such a fundamental role in describing physics of microscopic quantum systems.

Acknowledgements — The author acknowledges funding from the European Union’s
Horizon 2020 research and innovation programme under the Marie Sklodowska-
Curie Grant Agreement No. 890884. The author also acknowledges funding from
the Centre for complex quantum systems, Aarhus University, Denmark.

## Appendix

##  Appendix A Derivation of Eq.(IV.1.2) from Eq.(27)

In order to derive Eq.(IV.1.2) from Eq.(27), first we take the transpose of
Eq.(27),

$\displaystyle c_{vec}^{T}(t)=c_{vec}^{T}(0)e^{-i\mathbf{H}_{\rm
NH}^{T}t}+i\epsilon\int_{0}^{t}dt^{\prime}\xi_{vec}^{T}(t^{\prime})e^{-i\mathbf{H}_{\rm
NH}^{T}(t-t^{\prime})}.$ (98)

Note that with the transpose, the column vectors $c_{vec}(t)$ and
$\xi_{vec}(t)$ with elements $\\{\hat{c}_{\ell}(t)\\}$ and
$\\{\hat{\xi}_{\ell}(t)\\}$ have now become row vectors. Next, we take the
Hermitian conjugate of the above equation. This gives,

$\displaystyle c^{\dagger}_{vec}(t)=e^{i\mathbf{H}_{\rm
NH}^{*}t}c_{vec}^{\dagger}(0)-i\epsilon\int_{0}^{t}dt^{\prime}e^{i\mathbf{H}_{\rm
NH}^{*}(t-t^{\prime})}\xi_{vec}^{\dagger}(t^{\prime}),$ (99)

where $c^{\dagger}_{vec}(t)$ and $\xi_{vec}^{\dagger}(t)$ are column vectors
with elements $\\{\hat{c}_{\ell}^{\dagger}(t)\\}$ and
$\\{\hat{\xi}_{\ell}^{\dagger}(t)\\}$. Multiplying Eqs.(99) and (98) and
taking expectation values with respect to the initial state Eq.(7), we get

$\displaystyle\mathbf{C}(t)=e^{-\mathcal{G}t}\mathbf{C}(0)e^{-\mathcal{G^{\dagger}}t}$
$\displaystyle+\epsilon^{2}\int_{0}^{t}dt_{1}\int_{0}^{t}dt_{2}e^{-i(t-t_{1})\mathcal{G}}\langle\xi_{vec}^{\dagger}(t_{1})\xi_{vec}^{T}(t_{2})\rangle
e^{-i(t-t_{2})\mathcal{G}^{\dagger}},$ (100)

where $\mathcal{G}=-i\mathbf{H}_{\rm NH}^{*}$, as defined in Eq.(30), and
$\langle\xi_{vec}^{\dagger}(t_{1})\xi_{vec}^{T}(t_{2})\rangle$ is a matrix
with elements
$\langle\hat{\xi}_{\ell}^{\dagger}(t_{1})\hat{\xi}_{m}(t_{2})\rangle$. Using
the expression for
$\langle\hat{\xi}_{\ell}^{\dagger}(t_{1})\hat{\xi}_{m}(t_{2})\rangle$ from
Eqs.(13) and (14), and carrying out the integrations over time, we obtain
Eq.(IV.1.2).

Note that Eq.(VIII.1) can also be obtained as above, simply by making the time
arguments in Eqs.(98) and (99) different.

##  Appendix B Some explicit proofs

In the main text, it was mentioned that Eq.(IV.1.2) is the solution for
Eq.(IV.1.2). Further, it was mentioned that although Eq.(IV.1.2) is not the
solution of Eq.(35), the long-time solution of Eq.(35) is given by Eq.(31),
and is hence positive semi-definite. Here we explicitly show the steps for
proving these statements.

Figure 4: The domain of the integrations in Eq.(B).

The formal solution of Eq.(IV.1.2) is

$\displaystyle\mathbf{C}(t)=e^{-\mathcal{G}t}\mathbf{C}(0)e^{-\mathcal{G}^{\dagger}t}+\epsilon^{2}\mathbf{R}(t),$
(101)

where

$\displaystyle\mathbf{R}(t)=\int_{0}^{t}dt^{\prime}e^{-(t-t^{\prime})\mathcal{G}}\mathbf{Q}(t^{\prime})e^{-(t-t^{\prime})\mathcal{G}^{\dagger}}.$
(102)

Then, we have to show that $\mathbf{R}(t)$ simplifies to the second term in
Eq.(IV.1.2), when $\mathbf{Q}(t)$ is as defined in Eqs.(IV.1.2), (33). To show
this, first, we use the expression for $\mathbf{Q}(t)$ in $\mathbf{R}$, and
after some algebra, it can be written in the form

$\displaystyle\mathbf{R}(t)=$
$\displaystyle\int_{0}^{t}dt_{1}\int_{0}^{t_{1}}dt_{2}~{}\mathcal{I}(t_{1},t_{2};t)$
$\displaystyle+\int_{0}^{t}dt_{2}\int_{0}^{t_{2}}dt_{1}~{}\mathcal{I}(t_{1},t_{2};t),$
(103)

with

$\displaystyle\mathcal{I}(t_{1},t_{2};t)=\int\frac{d\omega}{2\pi}e^{-(t-t_{1})\mathcal{G}}e^{i(t_{1}-t_{2})\omega}\mathbf{F}(\omega)e^{-(t-t_{1})\mathcal{G}^{\dagger}}.$
(104)

At this point, visualizing the domain of the integrations, shown in Fig.(4)
becomes useful. Since the two terms in Eq.(B) correspond to exactly the same
function getting integrated over the lower triangular and the upper triangular
parts of square shown in Fig.(4), their sum exactly corresponds to integration
of the function over the entire square. Thus, we have

$\displaystyle\mathbf{R}(t)=\int_{0}^{t}dt_{1}\int_{0}^{t}dt_{2}~{}\mathcal{I}(t_{1},t_{2};t).$
(105)

Now, explicitly carrying out the integrations over time gives exactly second
term in Eq.(IV.1.2). This therefore proves that Eq.(IV.1.2) is indeed the
solution for Eq.(IV.1.2).

The formal solution of Eq.(35) is

$\displaystyle\mathbf{C}(t)=e^{-\mathcal{G}t}\mathbf{C}(0)e^{-\mathcal{G}^{\dagger}t}+\epsilon^{2}\mathbf{R}_{1}(t),$
(106)

where

$\displaystyle\mathbf{R}_{1}(t)=\int_{0}^{t}dt^{\prime}e^{-\mathcal{G}t^{\prime}}\mathbf{Q}_{1}e^{-\mathcal{G}^{\dagger}t^{\prime}}.$
(107)

Noting that $\mathbf{Q}_{1}$ (Eq.(36)) is nothing but $\mathbf{Q}(\infty)$,
similar algebraic simplifications as done to obtain Eq.(B) gives

$\displaystyle\mathbf{R}_{1}(t)=$
$\displaystyle\int_{0}^{t}dt_{1}\int_{0}^{\infty}dt_{2}~{}\mathcal{I}(t_{1},t_{2};t)$
$\displaystyle+\int_{0}^{t}dt_{2}\int_{0}^{\infty}dt_{1}~{}\mathcal{I}(t_{1},t_{2};t).$
(108)

Clearly, using the geometric visualization of domain to combine the two terms
in the sum into one integration is no longer possible for finite $t$. But for
$t\rightarrow\infty$, this can be done and it gives exactly same result as
$\mathbf{R}(\infty)$. By carrying out the time integration, it can be checked
that the result is exactly Eq.(31), which is manifestly positive semi-
definite.

##  Appendix C Derivation of Eq.(IX.2)

In this section, we give the derivation of Eq.(IX.2). In going to the single-
particle eigenbasis we have to diagonalise $\mathbf{H}$ as in Eq.(16). Writing
out Eq.(16) explicitly for Eq.(92), we have the following set of equations

$\displaystyle(\omega_{\alpha}-\varepsilon_{\ell})\Phi_{\ell\alpha}=g_{\ell}\Phi_{\ell+1\alpha}+g_{\ell-1}\Phi_{\ell-1\alpha}~{}~{}~{}\forall~{}~{}~{}\ell\neq
1,N;$
$\displaystyle(\omega_{\alpha}-\varepsilon_{1})\Phi_{1\alpha}=g_{1}\Phi_{2\alpha},~{}~{}(\omega_{\alpha}-\varepsilon_{N})\Phi_{N\alpha}=g_{N-1}\Phi_{N-1\alpha}.$
(109)

Using Eq.(59), and simplifying, the expression for NESS current becomes,

$\displaystyle
I=\epsilon^{2}g_{p}\sum_{\alpha,\nu=1}^{N}\frac{\Phi_{p\alpha}\Phi_{p+1\nu}-\Phi_{p\nu}\Phi_{p+1\alpha}}{\omega_{\alpha}-\omega_{\nu}}\Big{[}$
$\displaystyle\frac{\sum_{\ell,m}^{\prime}\Phi_{m\alpha}^{2}\Phi_{\ell\alpha}\Phi_{\ell\nu}\mathfrak{J}_{\ell}(\omega_{\alpha})\mathfrak{J}_{m}(\omega_{\alpha})\left[\mathfrak{n}_{m}(\omega_{\alpha})-\mathfrak{n}_{\ell}(\omega_{\alpha})\right]}{\sum_{\ell}^{\prime}\Phi_{\ell\alpha}^{2}\mathfrak{J}_{\ell}(\omega_{\alpha})}\Big{]},$
(110)

where $\sum^{\prime}$ refers to sum over the points where the baths are
attached. Remember that since we have used Eq.(VII.3), above expression for
current is only valid when Eq.(54) is satisfied. This expression can be
further simplified using an interesting result which we find for tridiagonal
matrices using Eq.(C),

$\displaystyle\Phi_{p\alpha}\Phi_{p+1\nu}-\Phi_{p\nu}\Phi_{p+1\alpha}=\frac{\omega_{\nu}-\omega_{\alpha}}{g_{p}}\sum_{k=1}^{p}\Phi_{k\alpha}\Phi_{k\nu}.$
(111)

Also, the eigenvectors are orthogonal, so

$\displaystyle\sum_{\nu=1}^{N}\Phi_{\ell\nu}\Phi_{m\nu}=\delta_{\ell m}.$
(112)

Using Eqs.(111), (112) in Eq.(C), we obtain Eq.(IX.2).

## References

  * Hur  _et al._ (2016) K. L. Hur, L. Henriet, A. Petrescu, K. Plekhanov, G. Roux, and M. Schiró, Comptes Rendus Physique 17, 808 (2016).
  * Binder  _et al._ (2018) F. Binder, L. A. Correa, C. Gogolin, J. Anders, and G. Adesso, eds.,  _Thermodynamics in the Quantum Regime_ (Springer, Cham, 2018).
  * Cao  _et al._ (2019) Y. Cao, J. Romero, J. P. Olson, M. Degroote, P. D. Johnson, M. Kieferová, I. D. Kivlichan, T. Menke, B. Peropadre, N. P. D. Sawaya, and et al., Chemical Reviews 119, 10856–10915 (2019).
  * Cao  _et al._ (2020) J. Cao, R. J. Cogdell, D. F. Coker, H.-G. Duan, J. Hauer, U. Kleinekathöfer, T. L. C. Jansen, T. Mančal, R. J. D. Miller, J. P. Ogilvie, V. I. Prokhorenko, T. Renger, H.-S. Tan, R. Tempelaar, M. Thorwart, E. Thyrhaug, S. Westenhoff, and D. Zigmantas, Science Advances 6 (2020).
  * Awschalom  _et al._ (2021) D. D. Awschalom, C. H. R. Du, R. He, J. Heremans, A. Hoffmann, J. Hou, H. Kurebayashi, Y. Li, L. Liu, V. Novosad, J. Sklenar, S. Sullivan, D. Sun, H. Tang, V. Tyberkevych, C. Trevillian, A. W. Tsen, L. Weiss, W. Zhang, X. Zhang, L. Zhao, and C. W. Zollitsch, IEEE Transactions on Quantum Engineering , 1 (2021).
  * Mahan (2000) G. D. Mahan,  _Many-particle physics_ (Springer, Boston, MA, 2000).
  * Haug and Jauho (2008) H. Haug and A.-P. Jauho,  _Quantum Kinetics in Transport and Optics of Semiconductors_ (Springer-Verlag Berlin Heidelberg, 2008).
  * Wang  _et al._ (2014) J.-S. Wang, B. K. Agarwalla, H. Li, and J. Thingna, Frontiers of Physics 9, 673 (2014).
  * Kamenev (2011) A. Kamenev,  _Field Theory of Non-Equilibrium Systems_ (Cambridge University Press, Cambridge, 2011).
  * Ford  _et al._ (1965) G. W. Ford, M. Kac, and P. Mazur, Journal of Mathematical Physics 6, 504 (1965).
  * Benguria and Kac (1981) R. Benguria and M. Kac, Phys. Rev. Lett. 46, 1 (1981).
  * 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).
  * Dhar and Sriram Shastry (2003) A. Dhar and B. Sriram Shastry, Phys. Rev. B 67, 195405 (2003).
  * Banerjee  _et al._ (2004) D. Banerjee, B. C. Bag, S. K. Banik, and D. S. Ray, The Journal of Chemical Physics 120, 8960 (2004).
  * Dhar and Sen (2006) A. Dhar and D. Sen, Phys. Rev. B 73, 085119 (2006).
  * Dhar and Roy (2006) A. Dhar and D. Roy, Journal of Statistical Physics 125, 801 (2006).
  * 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).
  * Tu and Zhang (2008) M. W. Y. Tu and W.-M. Zhang, Phys. Rev. B 78, 235311 (2008).
  * Jin  _et al._ (2010) J. Jin, M. W.-Y. Tu, W.-M. Zhang, and Y. Yan, New Journal of Physics 12, 083013 (2010).
  * Zhang  _et al._ (2012) W.-M. Zhang, P.-Y. Lo, H.-N. Xiong, M. W.-Y. Tu, and F. Nori, Phys. Rev. Lett. 109, 170402 (2012).
  * McCutcheon  _et al._ (2015) D. P. S. McCutcheon, J. P. Paz, and A. J. Roncaglia, Phys. Rev. Lett. 115, 168901 (2015).
  * Martinez and Paz (2013) E. A. Martinez and J. P. Paz, Phys. Rev. Lett. 110, 130406 (2013).
  * Yang  _et al._ (2013) L.-P. Yang, C. Y. Cai, D. Z. Xu, W.-M. Zhang, and C. P. Sun, Phys. Rev. A 87, 012110 (2013).
  * Ribeiro and Vieira (2015) P. Ribeiro and V. R. Vieira, Phys. Rev. B 92, 100302 (2015).
  * 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).
  * Trushechkin and Volovich (2016) A. S. Trushechkin and I. V. Volovich, EPL (Europhysics Letters) 113, 30005 (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).
  * Benatti  _et al._ (2020) F. Benatti, R. Floreanini, and L. Memarzadeh, Phys. Rev. A 102, 042219 (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).
  * Benatti  _et al._ (2021) F. Benatti, R. Floreanini, and L. Memarzadeh, PRX Quantum 2, 030344 (2021).
  * Trushechkin (2021) A. Trushechkin, Phys. Rev. A 103, 062226 (2021).
  * Tupkary  _et al._ (2022) D. Tupkary, A. Dhar, M. Kulkarni, and A. Purkayastha, Phys. Rev. A 105, 032208 (2022).
  * Redfield (1965) A. Redfield, in  _Advances in Magnetic Resonance_ (Elsevier, 1965) pp. 1–32.
  * Davidovic (2021) D. Davidovic, (2021), arXiv:2112.07863 [quant-ph] .
  * Purkayastha  _et al._ (2020a) A. Purkayastha, G. Guarnieri, M. T. Mitchison, R. Filip, and J. Goold, npj Quantum Information 6 (2020a).
  * Anderloni  _et al._ (2007) S. Anderloni, F. Benatti, and R. Floreanini, Journal of Physics A: Mathematical and Theoretical 40, 1625 (2007).
  * Gaspard and Nagaoka (1999) P. Gaspard and M. Nagaoka, The Journal of Chemical Physics 111, 5668 (1999).
  * Kohen  _et al._ (1997) D. Kohen, C. C. Marston, and D. J. Tannor, The Journal of Chemical Physics 107, 5236 (1997).
  * Gnutzmann and Haake (1996) S. Gnutzmann and F. Haake, Zeitschrift für Physik B Condensed Matter 101, 263 (1996).
  * Suárez  _et al._ (1992) A. Suárez, R. Silbey, and I. Oppenheim, The Journal of Chemical Physics 97, 5101 (1992).
  * Fleming and Cummings (2011) C. H. Fleming and N. I. Cummings, Phys. Rev. E 83, 031117 (2011).
  * Nathan and Rudner (2020) F. Nathan and M. S. Rudner, Phys. Rev. B 102, 115109 (2020).
  * Kleinherbers  _et al._ (2020) E. Kleinherbers, N. Szpak, J. König, and R. Schützhold, Phys. Rev. B 101, 125131 (2020).
  * Davidović (2020) D. Davidović, Quantum 4, 326 (2020).
  * Mozgunov and Lidar (2020) E. Mozgunov and D. Lidar, Quantum 4, 227 (2020).
  * McCauley  _et al._ (2020) G. McCauley, B. Cruikshank, D. I. Bondar, and K. Jacobs, npj Quantum Information 6 (2020).
  * Kiršanskas  _et al._ (2018) G. Kiršanskas, M. Franckié, and A. Wacker, Phys. Rev. B 97, 035432 (2018).
  * Prosen (2008) T. Prosen, New Journal of Physics 10, 043026 (2008).
  * Prosen and Žunkovič (2010) T. Prosen and B. Žunkovič, New Journal of Physics 12, 025016 (2010).
  * Žunkovič and Prosen (2012) B. Žunkovič and T. Prosen, AIP Conference Proceedings 1468, 350 (2012).
  * Koga and Yamamoto (2012) K. Koga and N. Yamamoto, Phys. Rev. A 85, 022103 (2012).
  * Nicacio  _et al._ (2015) F. Nicacio, A. Ferraro, A. Imparato, M. Paternostro, and F. L. Semião, Phys. Rev. E 91, 042116 (2015).
  * Landi  _et al._ (2013) G. T. Landi, T. Tomé, and M. J. de Oliveira, Journal of Physics A: Mathematical and Theoretical 46, 395001 (2013).
  * Nicacio  _et al._ (2016) F. Nicacio, M. Paternostro, and A. Ferraro, Phys. Rev. A 94, 052129 (2016).
  * Santos  _et al._ (2017) J. P. Santos, G. T. Landi, and M. Paternostro, Phys. Rev. Lett. 118, 220601 (2017).
  * Malouf  _et al._ (2019) W. T. B. Malouf, J. P. Santos, L. A. Correa, M. Paternostro, and G. T. Landi, Phys. Rev. A 99, 052104 (2019).
  * McDonald  _et al._ (2022) A. McDonald, R. Hanai, and A. A. Clerk, Phys. Rev. B 105, 064302 (2022).
  * Bernal-García  _et al._ (2021) D. N. Bernal-García, L. Huang, A. E. Miroshnichenko, and M. J. Woolley, (2021), arXiv:2111.04041 [quant-ph] .
  * Roccati  _et al._ (2021) F. Roccati, S. Lorenzo, G. M. Palma, G. T. Landi, M. Brunelli, and F. Ciccarello, Quantum Science and Technology 6, 025005 (2021).
  * Gajić and Qureshi (1995) Z. Gajić and M. T. J. Qureshi,  _Lyapunov matrix equation in system stability and control_ (Academic Press, Inc. San Diego, California, 1995).
  * Sun and Yong (2020) J. Sun and J. Yong,  _Stochastic linear-quadratic optimal control theory: Open-loop and closed-loop solutions_ (Springer Nature, Switzerland, 2020).
  * Yuan  _et al._ (2021) S. Yuan, M. Lv, S. Baldi, and L. Zhang, IEEE Transactions on Automatic Control 66, 2250 (2021).
  * Bergholtz  _et al._ (2021) E. J. Bergholtz, J. C. Budich, and F. K. Kunst, Rev. Mod. Phys. 93, 015005 (2021).
  * Ashida  _et al._ (2020) Y. Ashida, Z. Gong, and M. Ueda, Advances in Physics 69, 249 (2020).
  * El-Ganainy  _et al._ (2018) R. El-Ganainy, K. G. Makris, M. Khajavikhan, Z. H. Musslimani, S. Rotter, and D. N. Christodoulides, Nature Physics 14, 11 (2018).
  * El-Ganainy  _et al._ (2019) R. El-Ganainy, M. Khajavikhan, D. N. Christodoulides, and S. K. Ozdemir, Communications Physics 2, 37 (2019).
  * Feng  _et al._ (2017) L. Feng, R. El-Ganainy, and L. Ge, Nature Photonics 11, 752 (2017).
  * Quijandría  _et al._ (2018) F. Quijandría, U. Naether, S. K. Özdemir, F. Nori, and D. Zueco, Phys. Rev. A 97, 053846 (2018).
  * Purkayastha  _et al._ (2020b) A. Purkayastha, M. Kulkarni, and Y. N. Joglekar, Phys. Rev. Research 2, 043075 (2020b).
  * Huber  _et al._ (2020a) J. Huber, P. Kirton, S. Rotter, and P. Rabl, SciPost Physics 9 (2020a), 10.21468/scipostphys.9.4.052.
  * Huber  _et al._ (2020b) J. Huber, P. Kirton, and P. Rabl, Phys. Rev. A 102, 012219 (2020b).
  * Jaramillo Ávila  _et al._ (2020) B. Jaramillo Ávila, C. Ventura-Velázquez, R. d. J. León-Montiel, Y. N. Joglekar, and B. M. Rodríguez-Lara, Scientific Reports 10, 1761 (2020).
  * Álvaro Gómez-León  _et al._ (2021) Álvaro Gómez-León, T. Ramos, A. González-Tudela, and D. Porras, (2021), arXiv:2109.10930 [quant-ph] .
  * Arkhipov  _et al._ (2020) I. I. Arkhipov, A. Miranowicz, F. Minganti, and F. Nori, Phys. Rev. A 102, 033715 (2020).
  * Arkhipov  _et al._ (2021) I. I. Arkhipov, F. Minganti, A. Miranowicz, and F. Nori, Phys. Rev. A 104, 012205 (2021).
  * Wang  _et al._ (2021) K. Wang, A. Dutt, C. C. Wojcik, and S. Fan, Nature 598, 59 (2021).
  * Kawabata  _et al._ (2019) K. Kawabata, K. Shiozaki, M. Ueda, and M. Sato, Phys. Rev. X 9, 041015 (2019).
  * Torres (2019) L. E. F. F. Torres, Journal of Physics: Materials 3, 014002 (2019).
  * Gong  _et al._ (2018) Z. Gong, Y. Ashida, K. Kawabata, K. Takasan, S. Higashikawa, and M. Ueda, Phys. Rev. X 8, 031079 (2018).
  * Harari  _et al._ (2018) G. Harari, M. A. Bandres, Y. Lumer, M. C. Rechtsman, Y. D. Chong, M. Khajavikhan, D. N. Christodoulides, and M. Segev, Science 359, eaar4003 (2018).
  * Bandres  _et al._ (2018) M. A. Bandres, S. Wittek, G. Harari, M. Parto, J. Ren, M. Segev, D. N. Christodoulides, and M. Khajavikhan, Science 359, eaar4005 (2018).
  * Behr  _et al._ (2019) M. Behr, P. Benner, and J. Heiland, Calcolo 56, 51 (2019).
  * Walls and Milburn (2008) D. Walls and G. J. Milburn,  _Quantum pptics_ (Springer-Verlag Berlin Heidelberg, 2008).
  * Wiseman and Milburn (2010) H. M. Wiseman and G. J. Milburn,  _Quantum measurement and control_ (Cambridge Universtiy Press, New York, 2010).
  * Lindblad (1976) G. Lindblad, Communications in Mathematical Physics 48, 119 (1976).
  * Gorini  _et al._ (1976) V. Gorini, A. Kossakowski, and E. C. G. Sudarshan, Journal of Mathematical Physics 17, 821 (1976).
  * Gorini  _et al._ (1978) V. Gorini, A. Frigerio, M. Verri, A. Kossakowski, and E. Sudarshan, Reports on Mathematical Physics 13, 149 (1978).
  * Chung and Peschel (2001) M.-C. Chung and I. Peschel, Phys. Rev. B 64, 064412 (2001).
  * Dhar  _et al._ (2012) A. Dhar, K. Saito, and P. Hänggi, Phys. Rev. E 85, 011126 (2012).
  * Eisler and Peschel (2017) V. Eisler and I. Peschel, Journal of Physics A: Mathematical and Theoretical 50, 284003 (2017).
  * Thingna  _et al._ (2012) J. Thingna, J.-S. Wang, and P. Hänggi, The Journal of Chemical Physics 136, 194110 (2012).
  * Thingna  _et al._ (2013) J. Thingna, J.-S. Wang, and P. Hänggi, Phys. Rev. E 88, 052127 (2013).
  * Xu  _et al._ (2017) X. Xu, J. Thingna, and J.-S. Wang, Phys. Rev. B 95, 035428 (2017).
  * Chan  _et al._ (2014) C.-K. Chan, G.-D. Lin, S. F. Yelin, and M. D. Lukin, Phys. Rev. A 89, 042117 (2014).
  * Giusteri  _et al._ (2017) G. G. Giusteri, F. Recrosi, G. Schaller, and G. L. Celardo, Phys. Rev. E 96, 012113 (2017).
  * Kołodyński  _et al._ (2018) J. Kołodyński, J. B. Brask, M. Perarnau-Llobet, and B. Bylicka, Phys. Rev. A 97, 062124 (2018).
  * Maguire  _et al._ (2019) H. Maguire, J. Iles-Smith, and A. Nazir, Phys. Rev. Lett. 123, 093601 (2019).
  * Emary (2007) C. Emary, Phys. Rev. B 76, 245319 (2007).
  * Buča and Prosen (2012) B. Buča and T. Prosen, New Journal of Physics 14, 073007 (2012).
  * Finkelstein-Shapiro  _et al._ (2019) D. Finkelstein-Shapiro, S. Felicetti, T. Hansen, T. o. Pullerits, and A. Keller, Phys. Rev. A 99, 053829 (2019).
  * Quach and Munro (2020) J. Q. Quach and W. J. Munro, Phys. Rev. Applied 14, 024092 (2020).
  * Talkner (1986) P. Talkner, Annals of Physics 167, 390 (1986).
  * Ford and O’Connell (1996) G. W. Ford and R. F. O’Connell, Phys. Rev. Lett. 77, 798 (1996).
  * Guarnieri  _et al._ (2014) G. Guarnieri, A. Smirne, and B. Vacchini, Phys. Rev. A 90, 022110 (2014).
  * Khan  _et al._ (2021) S. Khan, B. K. Agarwalla, and S. Jain, (2021), arXiv:2111.14879 [quant-ph] .
  * Purkayastha  _et al._ (2018) A. Purkayastha, S. Sanyal, A. Dhar, and M. Kulkarni, Phys. Rev. B 97, 174206 (2018).
  * Salazar and Landi (2020) D. S. P. Salazar and G. T. Landi, Phys. Rev. Research 2, 033090 (2020).
  * Purkayastha and Dubi (2017) A. Purkayastha and Y. Dubi, Phys. Rev. B 96, 085425 (2017).