# Fundamental limitations in Lindblad descriptions of systems weakly coupled
to baths

Devashish Tupkary<EMAIL_ADDRESS>Indian Institute of Science,
Bangalore 560012, India  International Centre for Theoretical Sciences, Tata
Institute of Fundamental Research, Bangalore 560089, India  Institute for
Quantum Computing and Department of Physics and Astronomy, University of
Waterloo, Waterloo, Ontario, Canada, N2L 3G1  Abhishek Dhar
<EMAIL_ADDRESS>International Centre for Theoretical Sciences, Tata
Institute of Fundamental Research, Bangalore 560089, India  Manas Kulkarni
<EMAIL_ADDRESS>International Centre for Theoretical Sciences, Tata
Institute of Fundamental Research, Bangalore 560089, India  Archak Purkayastha
<EMAIL_ADDRESS>Trinity College Dublin, The University of Dublin, College
Green, Dublin, Ireland

###### Abstract

It is very common in the literature to write down a Markovian quantum master
equation in Lindblad form to describe a system with multiple degrees of
freedom and weakly connected to multiple thermal baths which can, in general,
be at different temperatures and chemical potentials. However, the
microscopically derived quantum master equation up to leading order in system-
bath coupling is of the so-called Redfield form which is known to not preserve
complete positivity in most cases. Additional approximations to the Redfield
equation are required to obtain a Lindblad form. We lay down some fundamental
requirements for any further approximations to Redfield equation, which, if
violated, leads to physical inconsistencies like inaccuracies in the leading
order populations and coherences in the energy eigenbasis, violation of
thermalization, violation of local conservation laws at the non-equilibrium
steady state (NESS). We argue that one or more of these conditions will
generically be violated in all the weak system-bath-coupling Lindblad
descriptions existing in literature to our knowledge. As an example, we study
the recently derived Universal Lindblad Equation (ULE) and use these
conditions to show violation of local conservation laws due to inaccurate
coherences but accurate populations in energy eigenbasis. Finally, we
exemplify our analytical results numerically in an interacting open quantum
spin system.

##  I Introduction

A fundamental problem relevant across a wide range of fields including quantum
optics Hur  _et al._ (2016), thermodynamics Binder  _et al._ (2018), chemistry
Cao  _et al._ (2019), engineering Awschalom  _et al._ (2021) and biology Cao
_et al._ (2020) is to describe a quantum system with multiple degrees of
freedom connected to multiple thermal baths. An approach very often taken is
to derive an effective quantum master equation (QME) for the dynamics of the
system assuming weak coupling between the system and the baths. The dynamics
of the quantum system, as described by the QME, is often desired to be
Markovian. It was shown by Gorini, Kossakowski, Sudarshan Gorini  _et al._
(1976), and Lindblad Lindblad (1976) (GKSL) that a QME that preserves all
properties of the density matrix and describes Markovian dynamics has to be of
the form,

$\displaystyle\frac{\partial\hat{\rho}}{\partial
t}=i[\hat{\rho},\hat{H}_{S}]+\hat{\mathcal{L}}(\hat{\rho}),$
$\displaystyle\hat{\mathcal{L}}(\hat{\rho})=i[\hat{\rho},\hat{H}_{LS}]+\sum_{\lambda=1}^{D^{2}-1}\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{)},$
(1)

which is commonly called a Lindblad equation. Here, $\hat{\rho}$ is the
density matrix of the system, $\hat{H}_{S}$ is the system Hamiltonian,
$\hat{H}_{LS}$ is the Hermitian contribution due to the presence of the baths,
commonly called the Lamb shift term, $\hat{L}_{\lambda}$ are the Lindblad
operators and $\gamma_{\lambda}$ are the rates, and $D$ is the dimension of
the Hilbert space. Here we will confine to (effective) finite dimensional
systems. Both $\hat{H}_{LS}$ and $\hat{L}_{\lambda}$ are operators in the
system Hilbert space. This equation preserves Hermiticity and trace of
$\hat{\rho}$, as well as ensures non-negativity of all eigenvalues of
$\hat{\rho}$ at all times for all initial states of the system. The last
condition crucially requires that $\gamma_{\lambda}\geq 0$ Gorini  _et al._
(1976); Lindblad (1976). This result is one of the cornerstones of open
quantum systems. Indeed, a large number of analytical Breuer and Petruccione
(2006); Carmichael (2002); Rivas and Huelga (2012) and numerical Weimer  _et
al._ (2021); Plenio and Knight (1998); Dalibard  _et al._ (1992); Mølmer  _et
al._ (1993) techniques rely on describing experimental systems via Lindblad
equations.

The standard way to obtain a QME to leading order in system-bath coupling is
via the Born-Markov approximation Breuer and Petruccione (2006); Carmichael
(2002); Rivas and Huelga (2012). The QME so obtained, often called the
Redfield equation (RE) Redfield (1965), though reducible to a Lindblad-like
form, is known to generically not satisfy complete positivity, i.e, not
satisfy the requirement $\gamma_{\lambda}\geq 0$ 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 that for certain initial states and at certain times, it
will not give a positive semi-definite density matrix. To rectify this
drawback, typically, further approximations are made to obtain a Lindblad
equation either in the so called local or global forms, which we call local
Lindblad (LLE) and eigenbasis Lindblad (ELE) equations respectively Breuer and
Petruccione (2006); Carmichael (2002); Rivas and Huelga (2012). Several
shortcomings of the LLE and ELE so obtained, in particular their failure to
correctly describe the steady state when coupled to multiple baths, have been
pointed out in the literature, and their regimes of validity discussed Walls
(1970); Wichterich  _et al._ (2007); 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 (2020); Scali  _et al._ (2021); Benatti  _et al._ (2021); Trushechkin
(2021); Davidovic (2021). Despite this, the conventional wisdom is that, in
principle, it should be possible to find a Lindblad equation, different from
both LLE and ELE, which accurately describes the steady state. Indeed, there
has been a number of recent attempts towards developing such new variants of
Lindblad equations Nathan and Rudner (2020); Kleinherbers  _et al._ (2020);
Davidović (2020); Mozgunov and Lidar (2020); McCauley  _et al._ (2020);
Kiršanskas  _et al._ (2018), which are intended to be as accurate as the RE.

Here, we lay down some fundamental requirements on any such attempt to recover
complete positivity, which, if violated, causes physical inconsistencies like
inaccuracies in the leading order populations (diagonal elements of density
matrix) and coherences (off-diagonal elements of density matrix) in the energy
eigenbasis, violation of thermalization and violation of local conservation
laws in the non-equilibrium steady state (NESS). We show that the RE does not
have any of the above physical inconsistencies, despite being generically not
completely positive. On the other hand, we argue that no existing weak system-
bath coupling Lindblad description, to our knowledge, generically satisfies
all the conditions, and therefore has one or more of the physical
inconsistencies despite being completely positive. As an example, we study in
detail the so called Universal Lindblad equation (ULE) Nathan and Rudner
(2020), which has been recently rigorously derived to have accuracy comparable
to that of the RE. We show that the ULE gives inaccurate coherences in the
energy eigenbasis to the leading order and violates local conservation laws.
All the above statements are shown in complete generality for time-independent
Hamiltonians, without writing down any specific system Hamiltonian. This
model-independent discussion sets apart our work from most previous works on
checking accuracy of various QMEs, where the accuracies were studied
numerically in specific models Walls (1970); Wichterich  _et al._ (2007);
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 (2020); Scali  _et al._ (2021);
Benatti  _et al._ (2021). We finally exemplify our discussion numerically in a
three-site XXZ model coupled to two bosonic baths.

The paper is organized as follows. In Sec.II, we discuss the fundamental
requirements. In Sec.III, we discuss how the RE satisfies all the fundamental
requirements except complete positivity. In Sec.IV, we show that the ULE
violates some of the fundamental requirements while restoring complete
positivity and make some general comments about other Lindblad equations. In
Sec.V we numerically exemplify our discussions using the XXZ-model. In Sec.VI,
we summarize and conclude.

##  II Fundamental requirements for an accurate weak-coupling Markovian
description

Figure 1: A typical two-terminal non-equilibrium set-up with not all sites
attached to the baths. The Hamiltonian $\hat{H}_{M}$ describes part of
$\hat{H}_{S}$ that commutes with the system-bath coupling Hamiltonians
$\hat{H}_{SB_{L}}$ and $\hat{H}_{SB_{R}}$, while the Hamiltonian $\hat{H}_{L}$
($\hat{H}_{R}$) describes part of the system Hamiltonian that does not commute
with $\hat{H}_{SB_{L}}$ ($\hat{H}_{SB_{R}}$).

###  II.1 The set-up and assumptions

For simplicity and concreteness, let us consider a typical two-terminal set-up
of the form given in Fig. 1. The Hamiltonian of the full set-up can be written
in the form

$\displaystyle\hat{H}=\hat{H}_{S}+\epsilon\hat{H}_{SB}+\hat{H}_{B}$
$\displaystyle\hat{H}_{SB}=\hat{H}_{SB_{L}}+\hat{H}_{SB_{R}},~{}~{}\hat{H}_{B}=\hat{H}_{B_{L}}+\hat{H}_{B_{R}},$
(2)

where $\hat{H}_{S}$ is the Hamiltonian of the system, $\hat{H}_{B_{L}}$
($\hat{H}_{B_{R}}$) is the Hamiltonian of the left (right) bath,
$\hat{H}_{SB_{L}}$ ($\hat{H}_{SB_{R}}$) is the coupling between the left
(right) bath and the system, and the dimensionless parameter $\epsilon$
controls the strength of the system-bath couplings. We consider the system
bath couplings to be weak, $\epsilon\ll 1$. The system Hamiltonian is further
broken into

$\displaystyle\hat{H}_{S}=\hat{H}_{L}+\hat{H}_{M}+\hat{H}_{R},$ (3)

where $\hat{H}_{M}$ contains the part of the system Hamiltonian that commutes
with the system-bath coupling Hamiltonians $\hat{H}_{SB_{L}}$ and
$\hat{H}_{SB_{R}}$, while, $\hat{H}_{L}$ ($\hat{H}_{R}$) contains the part of
the system Hamiltonian that does not commute with $\hat{H}_{SB_{L}}$
($\hat{H}_{SB_{R}}$). We will further assume for simplicity that the system
Hamiltonian $\hat{H}_{S}$ has no degeneracies. Initially the system is at some
arbitrary state while the baths are at inverse temperatures $\beta_{L}$ and
$\beta_{R}$, and chemical potentials $\mu_{L}$ and $\mu_{R}$. We take the
total number of particles (excitations) in the whole set-up to be conserved.
Without loss of generality, it is possible to assume that ${\rm
Tr}\left(\hat{H}_{SB}\hat{\rho}(0)\otimes\hat{\rho}_{B}\right)=0$ Breuer and
Petruccione (2006), where $\hat{\rho}_{B}$ denotes the composite state of the
two baths given by product of their individual thermal states. The state of
the system at a time $t$ is

$\displaystyle\hat{\rho}(t)=Tr_{\rm
B}\left(e^{-i\hat{H}t}\hat{\rho}(0)\otimes\hat{\rho}_{B}e^{i\hat{H}t}\right),$
(4)

where $Tr_{\rm B}(...)$ implies trace over bath degrees of freedom. We will
assume that in the long-time limit, the system reaches a unique NESS. This
assumption physically necessitates that the system size is finite, while the
baths are in the thermodynamic limit. The NESS density matrix is then defined
as

$\displaystyle\hat{\rho}_{\rm NESS}=\lim_{t\rightarrow\infty}Tr_{\rm
B}\left(e^{-i\hat{H}t}\hat{\rho}(0)\otimes\hat{\rho}_{B}e^{i\hat{H}t}\right).$
(5)

With this concrete, but still fairly general set-up and the assumptions in
mind, we now look at what are the fundamental requirements for a weak-
coupling-QME to accurately describe the NESS of the system.

###  II.2 The fundamental requirements

We will like to have a Markovian QME, written to leading order in system-bath
coupling, that accurately predicts the NESS density matrix. To this end, we
require the following physical conditions.

(a) Preservation of all properties of density matrix — We want a QME of the
form

$\displaystyle\frac{\partial\hat{\rho}}{\partial
t}=i[\hat{\rho},\hat{H}_{S}]+\epsilon^{2}\hat{\mathcal{L}}_{2}(\rho),$ (6)

where we have made it explicit that this equation is written to
$O(\epsilon^{2})$, which can be shown to be leading order Breuer and
Petruccione (2006). To preserve all properties of the density matrix: (i)
Hermiticity, (ii) trace and (iii) complete positivity,
$\hat{\mathcal{L}}_{2}(\rho)$ must be reducible to the form of
$\hat{\mathcal{L}}(\rho)$ given in Eq.(I), with $\gamma_{\lambda}\geq 0$.

(b) Correct populations and coherences in energy eigenbasis to the leading
order— Since we are making a weak-coupling approximation at the level of the
QME, we cannot expect to have accurate results at NESS to all orders in
$\epsilon$. But, we will like to have results correct to at least the leading
order. In particular, we will like to have correct results for both
populations and coherences in the energy eigenbasis at least to the leading
order. Coherences in the energy eigenbasis of the system, as well as
population inversions in the energy eigenbasis, have been considered a
resource in quantum thermodynamics and information Nielsen and Chuang (2011);
Bennett and DiVincenzo (2000); Streltsov  _et al._ (2017); Lostaglio  _et al._
(2015); Narasimhachar and Gour (2015); Mitchison  _et al._ (2015);
Allahverdyan  _et al._ (2004); Korzekwa  _et al._ (2016); Kammerlander and
Anders (2016); Santos  _et al._ (2019); Francica  _et al._ (2019). Thus, it is
important to describe populations and coherences correctly. Moreover, as we
will see later, both populations and coherences are crucially linked to having
correct currents at NESS.

(c) Thermalization — Further, having correct populations to leading order is
also linked with the fundamental phenomenon of thermalization of the system
with the baths when they have same temperatures and chemical potentials
$\beta_{L}=\beta_{R}=\beta$, $\mu_{L}=\mu_{R}=\mu$. In that case, on physical
grounds, the steady state density matrix is expected to satisfy the following
condition

$\displaystyle\lim_{\epsilon\rightarrow
0}\Big{(}\lim_{t\rightarrow\infty}\hat{\rho}(t)\Big{)}~{}=~{}\frac{e^{-\beta(\hat{H}_{S}-\mu\hat{N}_{S})}}{{\rm
Tr}\left(e^{-\beta(\hat{H}_{S}-\mu\hat{N}_{S})}\right)}.$ (7)

where $\hat{N}_{S}$ is the total particle or excitation number operator of the
system. In other words, populations of the density matrix in the energy
eigenbasis of the system, to the leading order in system-bath coupling, should
follow a Gibbs distribution specified by the temperature and the chemical
potential of the baths when they are all equal. As we will show later, if the
steady state is unique, this condition is always satisfied in the exact steady
state (Eq.5). Note that, in Eq.(7), the order of limits cannot be
interchanged.

(d) Preservation of local conservation laws — Writing the Heisenberg’s
equation of motion for our set-up and taking trace, it is clear that the rate
of change of expectation value of any operator $\hat{O}$ which commutes with
the system-bath couplings is given by

$\displaystyle\frac{d\langle\hat{O}\rangle}{dt}=-i\langle[\hat{O},\hat{H}_{S}]\rangle,~{}~{}\forall~{}~{}[\hat{O},\hat{H}_{SB}]=0,$
(8)

where $\langle\ldots\rangle={\rm Tr}(\hat{\rho}\ldots)$. This condition is
obviously true irrespective of strength of system-bath couplings. We will like
our effective weak system-bath coupling description to preserve this
condition. Combining Eq.(8) and Eq.(6), we see that this entails

$\displaystyle{\rm
Tr}\left(\hat{O}\hat{\mathcal{L}_{2}}(\hat{\rho})\right)=0,~{}~{}\forall~{}~{}[\hat{O},\hat{H}_{SB}]=0.$
(9)

This condition is not usually discussed in the context of deriving weak-
coupling Markovian descriptions. However, the importance of this condition
becomes clear when it is written for one of the locally conserved quantities,
such as $\hat{H}_{M}$,

$\displaystyle\frac{d\langle\hat{H}_{M}\rangle}{dt}=-i\langle[\hat{H}_{M},\hat{H}_{S}]\rangle=J_{L\rightarrow
M}-J_{M\rightarrow R},$ $\displaystyle J_{L\rightarrow
M}=-i\langle[\hat{H}_{M},\hat{H}_{L}]\rangle,~{}~{}J_{M\rightarrow
R}=i\langle[\hat{H}_{M},\hat{H}_{R}]\rangle,$ (10)

where we have defined the energy current into the central region from the left
as $J_{L\rightarrow M}$, and that from the central region to the right as
$J_{M\rightarrow R}$. These currents are expectation values of some system
operators and do not involve any explicit dependence on system-bath couplings.
At NESS, the rate of change of expectation values of all system operators is
zero. Thus, $\frac{d\langle\hat{H}_{M}\rangle}{dt}=0$, which implies
$J_{L\rightarrow M}=J_{M\rightarrow R}$. So we have,

$\displaystyle\langle[\hat{H}_{M},\hat{H}_{L}]\rangle=-\langle[\hat{H}_{M},\hat{H}_{R}]\rangle\neq
0.$ (11)

This is a fundamental property of the NESS that should exactly hold
irrespective of strength of system bath couplings. The Eq.(II.2), which at
NESS leads to the above condition, is a continuity equation coming from local
conservation of energy. Analogous conditions can be derived for other locally
conserved quantities. The QME written to leading order should satisfy all such
conditions. It is clear that Eq.(9) guarantees this.

It can be seen that satisfying Eq.(9), and giving correct rate of change of
expectation value of $\hat{O}$ in Eq.(8) requires that the coherences in the
energy eigenbasis are given correctly. To see this, we write the system
Hamiltonian in the energy eigenbasis,

$\displaystyle\hat{H}_{S}=\sum_{\alpha}E_{\alpha}|E_{\alpha}\rangle\langle
E_{\alpha}|.$ (12)

Using this basis, Eq.(8) can be written as

$\displaystyle\frac{d\langle\hat{O}\rangle}{dt}=-i\sum_{\alpha,\nu=1}^{D}(E_{\alpha}-E_{\nu})\langle
E_{\alpha}|\hat{\rho}|E_{\nu}\rangle\langle
E_{\nu}|\hat{O}|E_{\alpha}\rangle,$
$\displaystyle\forall~{}~{}[\hat{O},\hat{H}_{SB}]=0.$ (13)

Clearly, the right-hand-side of above equation is governed by the coherences
between in the energy eigenbasis of the system. Thus, satisfying local
conservation laws and obtaining correct currents at NESS requires the
coherences to be correct (at least to leading order).

1. | Correct populations to leading order [$O(\epsilon^{0})$] |  $\left\langle E_{\alpha}\left|\hat{\mathcal{L}}_{2}^{\prime}\left[\sum_{\alpha}p_{\alpha}|E_{\alpha}\rangle\langle E_{\alpha}|\right]\right|E_{\alpha}\right\rangle=\left\langle E_{\alpha}\left|\hat{\mathcal{L}}_{2}\left[\sum_{\alpha}p_{\alpha}|E_{\alpha}\rangle\langle E_{\alpha}|\right]\right|E_{\alpha}\right\rangle$,
---|---|---
2. | Correct coherences to leading order [$O(\epsilon^{2})$] |  $\left\langle E_{\alpha}\left|\hat{\mathcal{L}}_{2}^{\prime}\left[\sum_{\alpha}p_{\alpha}|E_{\alpha}\rangle\langle E_{\alpha}|\right]\right|E_{\nu}\right\rangle=\left\langle E_{\alpha}\left|\hat{\mathcal{L}}_{2}\left[\sum_{\alpha}p_{\alpha}|E_{\alpha}\rangle\langle E_{\alpha}|\right]\right|E_{\nu}\right\rangle~{}\forall~{}\alpha\neq\nu$.
3. | Thermalization [Eq.(7))] |  $\left\langle E_{\alpha}\left|\hat{\mathcal{L}}_{2}^{\prime}\left[\sum_{\alpha}p_{\alpha}|E_{\alpha}\rangle\langle E_{\alpha}|\right]\right|E_{\alpha}\right\rangle=0\Rightarrow p_{\alpha}\propto e^{-\beta(E_{\alpha}-\mu N_{\alpha})}$
when $\beta_{L}=\beta_{R}=\beta$, $\mu_{L}=\mu_{R}=\mu$.
4. | Preservation of local conservation laws |  ${\rm Tr}\left(\hat{O}\hat{\mathcal{L}}_{2}^{\prime}(\hat{\rho})\right)=0,~{}~{}\forall~{}~{}[\hat{O},\hat{H}_{SB}]=0.$
5. | Complete positivity |  $\hat{\mathcal{L}}_{2}^{\prime}[\hat{\rho}]=i[\hat{\rho},\hat{H}_{LS}]+\sum_{\lambda=1}^{D^{2}-1}\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{)}$, $\gamma_{\lambda}\geq 0$.
Table 1: Summary of necessary conditions that a QME written to the leading
order in system-bath coupling must satisfy to accurately describe the steady
state of a set-up of the form in Fig.1. In above, $\hat{\mathcal{L}}_{2}$ is
the $O(\epsilon^{2})$ term obtained by systematically doing time-convolution-
less expansion and taking long-time limit [see Eqs.(15), (16)]. The
$\hat{\mathcal{L}}_{2}[\hat{\rho}]$ coincides with that obtained in Redfield
equation via Born-Markov approximation and is known to generically not satisfy
the complete positivity requirement $\gamma_{\lambda}\geq 0$. The
$\hat{\mathcal{L}}_{2}^{\prime}$ is any modification of
$\hat{\mathcal{L}}_{2}$ which may restore complete positivity.

###  II.3 The equations giving correct populations and coherences to the
leading order

Let us now see what condition in terms of $\hat{\mathcal{L}}_{2}$ must be
satisfied so that the populations are the coherences in the energy eigenbasis
are correct to the leading order. We assume that the NESS density matrix
defined in Eq.(5) and can be expanded in powers of $\epsilon$,

$\hat{\rho}_{\rm NESS}=\sum_{m=0}^{\infty}\epsilon^{2m}\hat{\rho}_{\rm
NESS}^{(2m)}.$ (14)

The QME describing the evolution of $\hat{\rho}(t)$ can also be systematically
expanded in the so-called time-convolution-less form Breuer and Petruccione
(2006),

$\frac{\partial\hat{\rho}}{\partial
t}=\sum_{m=0}^{\infty}\epsilon^{2m}\hat{{\mathcal{L}}}_{2m}(t)[\hat{\rho}(t)],$
(15)

where $\hat{{\mathcal{L}}}_{2m}(t)$ are in general time dependent linear
operators and
$\hat{{\mathcal{L}}}_{0}(t)[\hat{\rho}(t)]=i[\hat{\rho}(t),\hat{H}_{S}]$. By
definition, the steady state $\hat{{\rho}}_{\rm NESS}$ satisfies

$\displaystyle
0=\sum_{m=0}^{\infty}\epsilon^{2m}{\hat{\mathcal{L}}}_{2m}[\hat{\rho}_{\rm
NESS}],$ (16)

where we denote
$\hat{\mathcal{L}}_{2m}\equiv\hat{{\mathcal{L}}}_{2m}(t\to\infty)$.

Using the series expansion of $\hat{{\rho}}_{\rm NESS}$ and matching terms
order-by-order yields, at the lowest order ($\epsilon^{0}$), the condition
$[\hat{{\rho}}_{\rm NESS}^{(0)},\hat{H}_{S}]=0.$ This implies that
$\hat{{\rho}}_{\rm NESS}^{(0)}$ is diagonal in the energy eigenbasis of the
system,

$\displaystyle\langle E_{\alpha}|\hat{{\rho}}_{\rm
NESS}^{(0)}|E_{\nu}\rangle=0~{}~{}\forall~{}\alpha\neq\nu.$ (17)

This shows that while the populations can have $O(\epsilon^{0})$ term, the
leading order contribution to coherences in energy eigenbasis is
$O(\epsilon^{2})$. (Remember that we have assumed no degeneracies in the
system Hamiltonian.) This means $\hat{{\rho}}_{\rm NESS}^{(0)}$ can be written
as

$\displaystyle\hat{\rho}_{\rm
NESS}^{(0)}=\sum_{\alpha}p_{\alpha}|E_{\alpha}\rangle\langle E_{\alpha}|.$
(18)

As a consequence, when the temperatures and the chemical potentials of the
baths are equal, $\beta_{L}=\beta_{R}=\beta$, $\mu_{L}=\mu_{R}=\mu$, the
thermalization condition Eq.(7) requires,

$\displaystyle p_{\alpha}\propto e^{-\beta(E_{\alpha}-\mu N_{\alpha})},$ (19)

where $N_{\alpha}$ is the number of particles or excitations in the energy
eigenstate $E_{\alpha}$.

Let us now look at two higher orders in $\epsilon$. At order $\epsilon^{2}$,
using Eq.(14) in Eq. (16) then leads to the conditions:

$\displaystyle\left\langle
E_{\alpha}\left|\hat{\mathcal{L}}_{2}[\hat{{\rho}}_{\rm
NESS}^{(0)}]\right|E_{\alpha}\right\rangle=0,$ (20) $\displaystyle
i(E_{\alpha}-E_{\nu})\left\langle E_{\alpha}\left|\hat{{\rho}}_{\rm
NESS}^{(2)}\right|E_{\nu}\right\rangle$ $\displaystyle+\left\langle
E_{\alpha}\left|\hat{\mathcal{L}}_{2}[\hat{{\rho}}_{\rm
NESS}^{(0)}]\right|E_{\nu}\right\rangle=0\quad\forall~{}~{}\alpha\neq\nu,$
(21)

and at order $\epsilon^{4}$ we get:

$\displaystyle\left\langle
E_{\alpha}\left|\hat{\mathcal{L}}_{2}[\hat{{\rho}}_{\rm
NESS}^{(2)}]\right|E_{\alpha}\right\rangle+\left\langle
E_{\alpha}\left|\hat{\mathcal{L}}_{4}[\hat{{\rho}}_{\rm
NESS}^{(0)}]\right|E_{\alpha}\right\rangle=0.$ (22)

The set of equations obtained by writing Eq. (20) for various values of
eigenstate index $\alpha$, gives a complete set of equations which determines
the non-zero elements of $\hat{\rho}_{\rm NESS}^{(0)}$. Given $\hat{\rho}_{\rm
NESS}^{(0)}$, the set of equations obtained by writing Eq. (21) for various
values of eigenstate indices $\alpha$ and $\nu$, determines the off-diagonal
elements of $\hat{\rho}_{\rm NESS}^{(2)}$ in the energy eigenbasis. Finally,
given the off-diagonal elements of $\hat{\rho}_{\rm NESS}^{(2)}$ in the energy
eigenbasis, and the $\hat{\rho}_{\rm NESS}^{(0)}$, the set of equations
obtained by writing Eq. (22) for various values of the eigenstate index
$\alpha$ determines the diagonal elements of $\hat{\rho}_{\rm NESS}^{(2)}$. It
is crucial to note the occurrence of $\hat{\mathcal{L}}_{4}$ in Eq. (22).
Since the leading order populations are $O(\epsilon^{0})$ and the leading
order coherences are $O(\epsilon^{2})$, Eqs.(20) and (22) give the necessary
equations to solve in order to obtain populations and coherences correct to
the leading order. The above discussion correspond to the NESS obtained from
the exact QME, Eq. (15), in the regime of small $\epsilon$.

###  II.4 The need to go beyond standard Lindblad descriptions

In the previous subsection, $\hat{\mathcal{L}}_{2}$ was derived from a time-
convolution-less expansion. This fixes the form of $\hat{\mathcal{L}}_{2}$. If
we truncate the time-convolution-less expansion at $O(\epsilon^{2})$ and take
the long-time limit, we get an equation exactly of our desired form in Eq.(6).
This coincides with Born-Markov approximation and gives the RE Breuer and
Petruccione (2006). However, it is known that the RE does not generically
satisfy the requirement of complete positvity 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). Due to this, very often further approximations are made on
$\hat{\mathcal{L}}_{2}$, modifying it to, say,
$\hat{\mathcal{L}}_{2}^{\prime}$ so as to impose complete positivity. Our
analysis in previous subsections allows us to impose some necessary
restrictions on the approximations by fixing certain components of
$\hat{\mathcal{L}}_{2}^{\prime}$, so that the fundamental requirements are
satisfied. The set of conditions is summarized in Table 1. In particular, we
see that for both populations and coherences to be correct to the leading
order, we require the operation of $\hat{\mathcal{L}}_{2}^{\prime}$ on any
state diagonal in the energy eigenbasis to be same as that of
$\hat{\mathcal{L}}_{2}$.

In most studies of accuracy of results from QMEs, the accuracy is checked by
applying them to some particular system and comparing results with some other
more accurate method Walls (1970); Wichterich  _et al._ (2007); 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 (2020); Scali  _et al._ (2021); Benatti  _et
al._ (2021); Trushechkin (2021). This does not let one clearly comment on
situations where such more accurate methods are unavailable. The set of
conditions in Table 1 allows us to assess the accuracy of
$\hat{\mathcal{L}}_{2}^{\prime}$ by making formal checks, without explicitly
writing it down for any particular system. They therefore allow a completely
general discussion of accuracy of NESS results from QMEs, even in cases where
some more accurate method is unavailable.

The most standard approximations on $\hat{\mathcal{L}}_{2}$ involve reducing
it to either the form of LLE or ELE. As mentioned before, it is already well-
established that the LLE and the ELE do not satisfy all the fundamental
requirements mentioned in Sec. II.2 Walls (1970); Wichterich  _et al._ (2007);
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 (2020); Scali  _et al._ (2021);
Benatti  _et al._ (2021); Trushechkin (2021), and therefore do not satisfy all
the conditions in Table 1. The LLE does not show thermalization, while it
explicitly satisfies the fourth condition in Table 1, thereby preserving local
conservation laws. On the other hand, the ELE explicitly neglects various
coherences in the energy eigenbasis of the system, but satisfies
thermalization Breuer and Petruccione (2006). Moreover, ELE does not also
satisfy the fourth condition in Table 1. Further, the ELE is known to not
satisfy Eq.(11), because it gives zero currents inside the system even in an
out-of-equilibrium condition Wichterich  _et al._ (2007); Purkayastha  _et
al._ (2016). So neither of these standard Lindblad equations satisfy all the
requirements above. This makes it necessary to go beyond these standard
approximations, even at weak system-bath coupling.

In the following section, we consider the RE in more detail and show that even
though the RE does not generically satisfy the requirement of complete
positivity, it satisfies all other of the fundamental requirements.

##  III The Redfield equation

###  III.1 Accuracy of populations and coherences from the Redfield equation
and lack of complete positivity

As mentioned before, the RE is obtained by truncating Eq.(15) at second order
and replacing $\hat{{\mathcal{L}}}_{2}(t)$ by
$\hat{{\mathcal{L}}}_{2}(t\to\infty)\equiv\hat{{\mathcal{L}}}_{2}$ Breuer and
Petruccione (2006),

$\displaystyle\frac{\partial{\hat{{\rho}}}}{\partial
t}=\hat{\mathcal{L}}_{0}[\hat{{\rho}}(t)]+\epsilon^{2}\hat{\mathcal{L}}_{2}[\hat{{\rho}}(t)].$
(23)

The NESS obtained from the RE, which we denote by $\hat{\widetilde{\rho}}_{\rm
NESS}$, is given by

$\displaystyle\hat{\widetilde{\rho}}_{\rm
NESS}=\lim_{t\rightarrow\infty}e^{t(\hat{\mathcal{L}}_{0}+\hat{\mathcal{L}}_{2})}[{\rho}(0)],$
(24)

and satisfies

$\displaystyle 0=\hat{\mathcal{L}}_{0}[\hat{\widetilde{{\rho}}}_{\rm
NESS}]+\epsilon^{2}\hat{\mathcal{L}}_{2}[\hat{\widetilde{{\rho}}}_{\rm
NESS}].$ (25)

The density matrix can also be expanded in powers of $\epsilon$,
$\hat{\widetilde{\rho}}_{\rm
NESS}=\sum_{m=0}^{\infty}\epsilon^{2m}\hat{\widetilde{\rho}}_{\rm
NESS}^{(2m)}$. The question then is, to what order in $\epsilon$ does elements
of $\hat{\widetilde{\rho}}_{\rm NESS}$ and $\hat{{\rho}}_{\rm NESS}$ agree.

Proceeding as before we again get $[\hat{\widetilde{\rho}}^{(0)}_{\rm
NESS},\hat{H}_{S}]=0$ implying a diagonal $\hat{\widetilde{\rho}}^{(0)}_{\rm
NESS}$ in energy eigenbasis. At $O(\epsilon^{2})$ we now obtain:

$\displaystyle\left\langle
E_{\alpha}\left|{\hat{\mathcal{L}}}_{2}[\hat{{\widetilde{\rho}}}_{\rm
NESS}^{(0)}]\right|E_{\alpha}\right\rangle=0,$ (26) $\displaystyle
i(E_{\alpha}-E_{\nu})\left\langle
E_{\alpha}\left|\hat{{\widetilde{\rho}}}_{\rm
NESS}^{(2)}\right|E_{\nu}\right\rangle$ $\displaystyle+\left\langle
E_{\alpha}\left|\hat{\mathcal{L}}_{2}[\hat{{\widetilde{\rho}}}_{\rm
NESS}^{(0)}]\right|E_{\nu}\right\rangle=0\quad\forall~{}~{}\alpha\neq\nu,$
(27)

while at $O(\epsilon^{4})$ we get

$\displaystyle\left\langle
E_{\alpha}\left|\hat{\mathcal{L}}_{2}[\hat{{\widetilde{\rho}}}_{\rm
NESS}^{(2)}]\right|E_{\alpha}\right\rangle=0.$ (28)

Equations (26), (27) and (28) are the analogs of Eqs. (20), (21), and (22)
respectively. We see from Eqs.(20) and (26) that $\hat{{\rho}}^{(0)}_{\rm
NESS}$ and $\hat{\tilde{\rho}}^{(0)}_{\rm NESS}$ satisfies the exact same set
of equations, which implies

$\displaystyle\hat{\widetilde{\rho}}_{\rm NESS}^{(0)}=\hat{\rho}_{\rm
NESS}^{(0)}.$ (29)

Next, from Eqs.(21) and (27), we see that the off-diagonal elements of
$\hat{\rho}_{\rm NESS}^{(2)}$ and $\hat{\widetilde{\rho}}_{\rm NESS}^{(2)}$ in
energy eigenbasis also satisfy exactly same equation therefore implying,

$\displaystyle\left\langle E_{\alpha}\left|\hat{\widetilde{\rho}}_{\rm
NESS}^{(2)}\right|E_{\nu}\right\rangle=\left\langle
E_{\alpha}\left|\hat{\rho}_{\rm
NESS}^{(2)}\right|E_{\nu}\right\rangle,~{}~{}\forall~{}~{}\alpha\neq\nu.$ (30)

Thus, the leading order terms in populations and coherences are given
correctly by the RE. However, Eq.(22) which fixes the diagonal elements of
$\hat{\rho}_{\rm NESS}^{(2)}$ in energy eigenbasis is different from Eq.(28)
which fixes the same for $\hat{\widetilde{\rho}}_{\rm NESS}^{(2)}$. So,

$\displaystyle\left\langle E_{\alpha}\left|\hat{{{\widetilde{\rho}}}}_{\rm
NESS}^{(2)}\right|E_{\alpha}\right\rangle\neq\left\langle
E_{\alpha}\left|\hat{\rho}_{\rm NESS}^{(2)}\right|E_{\alpha}\right\rangle.$
(31)

Specifically, Eq.(22) shows that, to obtain the diagonal elements of NESS
density matrix in the energy eigenbasis correct to $O(\epsilon^{2})$, one
needs the QME up of $O(\epsilon^{4})$. Based on these results, we have,

$\displaystyle||\hat{\rho}_{\rm NESS}-\hat{\widetilde{\rho}}_{\rm NESS}||\sim
O(\epsilon^{2}),$ (32)

where $||\hat{P}||$ is norm of the operator $\hat{P}$. So, the error in
obtaining the NESS from Eq.(23) is $O(\epsilon^{2})$. Despite this, the
leading order coherences in the energy eigenbasis, which are
$O(\epsilon^{2})$, are given correctly. Having $O(\epsilon^{2})$ term in the
coherences accurately while not having the corresponding correction in the
populations can cause at $O(\epsilon^{2})$ violation of positivity of the
density matrix Fleming and Cummings (2011); Purkayastha  _et al._ (2020).
Thus, the lack of complete positivity of the RE stems from the above mismatch
in order of accuracy between populations and coherences.

###  III.2 Thermalization

The explicit form of Eq.(23) is given by (Eq.(9.52) of Ref.Breuer and
Petruccione (2006)),

$\displaystyle\frac{\partial\hat{\rho}}{\partial t}=$ $\displaystyle
i[\hat{\rho}(t),{H}_{S}]$ $\displaystyle+$
$\displaystyle\epsilon^{2}\int_{0}^{\infty}dt^{\prime}[\hat{H}_{SB},[\hat{H}_{SB}(-t^{\prime}),\hat{\rho}(t)\otimes\hat{\rho}_{B}]],$
(33)

where
$\hat{H}_{SB}(t)=e^{i(\hat{H}_{S}+\hat{H}_{B})t}\hat{H}_{SB}e^{-i(\hat{H}_{S}+\hat{H}_{B})t}$.
Let us now consider the following canonical model of thermal baths:

$\displaystyle\hat{H}_{B}$
$\displaystyle=\sum_{\ell}{{}^{{}^{\prime}}}\sum_{r=1}^{\infty}\Omega^{\ell}_{r}\hat{B}^{{(\ell)}\dagger}_{r}\hat{B}^{(\ell)}_{r},$
(34) $\displaystyle\hat{H}_{SB}$
$\displaystyle=\sum_{\ell}{{}^{{}^{\prime}}}\sum_{r=1}^{\infty}(\kappa_{\ell
r}\hat{B}^{(\ell)\dagger}_{r}\hat{S}_{\ell}+\kappa_{\ell
r}^{*}\hat{S}_{\ell}^{\dagger}\hat{B}^{(\ell)}_{r})$ (35)

where $\sum_{\ell}^{\prime}$ indicates sum over all sites of the system where
the baths are attached, $\hat{B}^{(\ell)}_{r}$ is bosonic or fermionic
annihilation operator for the $r$-th mode of the bath attached at the
$\ell$-th site and $\hat{S}_{\ell}$ is the system operator coupling to the
bath at site $\ell$. At initial time, the baths are taken to be in their
respective thermal state with inverse temperatures $\beta_{\ell}$ and chemical
potentials $\mu_{\ell}$. The dynamics of the system can be shown to be
governed by the bath spectral functions, defined as

$\displaystyle\mathfrak{J}_{\ell}(\omega)=\sum_{k}2\pi\left|\kappa_{\ell
k}\right|^{2}\delta(\omega-\Omega^{\ell}_{k})$ (36)

and the Fermi or Bose distributions,
$n_{\ell}(\omega)=[e^{\beta_{\ell}(\omega-\mu_{\ell})}\pm 1]^{-1}$,
corresponding to the initial states of the baths. The RE for this set-up is
obtained by simplification of Eq. (III.2), using Eqs. (34),(35):

$\displaystyle\frac{\partial\hat{\rho}}{\partial t}$
$\displaystyle=i[\hat{\rho},\hat{H}_{S}]$
$\displaystyle-\epsilon^{2}\sum_{\ell}{{}^{{}^{\prime}}}\Big{(}\big{[}\hat{S}_{\ell}^{\dagger},\hat{S}_{\ell}^{(1)}\hat{\rho}\big{]}+\big{[}\hat{\rho}\hat{S}_{\ell}^{(2)},\hat{S}_{\ell}^{\dagger}\big{]}+{\rm
h.c.}\Big{)},$ (37)

where

$\displaystyle\hat{S}_{\ell}^{(1)}=\int_{0}^{\infty}dt^{\prime}\int\frac{d\omega}{2\pi}\hat{S}_{\ell}(-t^{\prime})e^{\beta_{\ell}(\omega-\mu_{\ell})}\mathfrak{J}_{\ell}(\omega)n_{\ell}(\omega)e^{-i\omega
t^{\prime}},$
$\displaystyle\hat{S}_{\ell}^{(2)}=\int_{0}^{\infty}dt^{\prime}\int\frac{d\omega}{2\pi}\hat{S}_{\ell}(-t^{\prime})\mathfrak{J}_{\ell}(\omega)n_{\ell}(\omega)e^{-i\omega
t^{\prime}},$ (38)

with $\hat{S}_{\ell}(t)=e^{i\hat{H}_{S}t}\hat{S}_{\ell}e^{-i\hat{H}_{S}t}$,
and h.c. denoting Hermitian conjugate. Note that, in general,
$[\hat{S}_{\ell},\hat{S}_{\ell}^{(1)}]\neq 0$,
$[\hat{S}_{\ell},\hat{S}_{\ell}^{(2)}]\neq 0$. Let us now check if the RE
satisfies the thermalization condition given in Eq.(7). To this end, we need
to find the leading order populations when the temperatures and chemical
potentials of the baths are the same, $\beta_{L}=\beta_{R}=\beta$,
$\mu_{L}=\mu_{R}=\mu$. Defining $E_{\nu\alpha}=E_{\nu}-E_{\alpha}$ and writing
$\hat{\rho}_{\rm NESS}^{(0)}=\sum_{\alpha}p_{\alpha}|E_{\alpha}\rangle\langle
E_{\alpha}|$, it can be checked after some algebra that Eq. (26), for the RE
in Eq. (III.2), takes the form:

$\displaystyle\sum_{\nu}\Big{[}|\langle
E_{\nu}|\hat{S}_{\ell}|E_{\alpha}\rangle|^{2}\left[p_{\alpha}-e^{\beta(E_{\nu\alpha}-\mu)}p_{\nu}\right]\mathfrak{J}_{\ell}(E_{\nu\alpha})n_{\ell}(E_{\nu\alpha})$
$\displaystyle+|\langle
E_{\alpha}|\hat{S}_{\ell}|E_{\nu}\rangle|^{2}\left[p_{\nu}-e^{\beta(E_{\alpha\nu}-\mu)}p_{\alpha}\right]\mathfrak{J}_{\ell}(E_{\alpha\nu})n_{\ell}(E_{\alpha\nu})\Big{]}$
$\displaystyle=0.$ (39)

We now consider the case where the system operator coupling to the bath can
create or annihilate a single particle, i.e $\langle
E_{\nu}|\hat{S}_{\ell}|E_{\alpha}\rangle=0$ $\forall$ $N_{\alpha}-N_{\nu}\neq
1$. This condition ensures that with the system-bath coupling of the form in
Eq.(35), the total number of exictations in the full set-up is conserved. It
can then be checked that the following choice of $p_{\alpha}$ satisfies the
above equation,

$\displaystyle p_{\alpha}\propto e^{-\beta(E_{\alpha}-\mu N_{\alpha})}.$ (40)

Since we have assumed the steady state to be unique, this proves that the
thermalization condition, i.e, the third condition in Table. 1, is satisfied.

Note that, since the leading order populations obtained from RE are exactly
the same as those obtained from the full general time-convolutionless-equation
Eq.(15), the above proof of thermalization is essentially a proof of
thermalization for the full time-convolutionless-equation, Eq.(15). Thus, if
the steady state is unique, the thermalization condition is satisfied in
complete generality. So, any Lindblad equation that does not satisfy this
condition (such as the LLE) cannot describe the steady state of a system
coupled to thermal baths.

###  III.3 Satisfying local conservation laws

Given the RE, Eq.(III.2) one can write down an expression for expectation
value of any system operator ${O}$,

$\displaystyle\frac{d\braket{\hat{O}}}{dt}=$
$\displaystyle-i\braket{[\hat{O},\hat{H}_{S}]}$ $\displaystyle-\epsilon^{2}$
$\displaystyle\sum_{\ell}^{\prime}\Big{(}\left\langle[\hat{O},\hat{S}_{\ell}^{\dagger}]\hat{S}_{\ell}^{(1)}\right\rangle-\left\langle\hat{S}_{\ell}^{(2)}[\hat{O},\hat{S}_{\ell}^{\dagger}]\right\rangle$
$\displaystyle+\left\langle\hat{S}_{\ell}^{(1)^{\dagger}}[\hat{S}_{\ell},\hat{O}]\right\rangle-\left\langle[\hat{S}_{\ell},\hat{O}]\hat{S}_{\ell}^{(2)^{\dagger}}\right\rangle\Big{)}.$
(41)

It is clear from above expression that the fourth condition in Table. 1 is
satisfied. Thus, it is clear that the local conservation laws inside the
system are satisfied. Not only that, as we show below, further conservation
laws relating currents from the baths and currents inside the system are also
satisfied.

For the kind of set-up in Fig.1, the RE can be written as

$\displaystyle\frac{\partial\hat{\rho}}{\partial
t}=i[\hat{\rho},\hat{H}_{S}]+\epsilon^{2}\left(\hat{\mathcal{L}}_{2}^{(L)}[\hat{\rho}(t)]+\hat{\mathcal{L}}_{2}^{(R)}[\hat{\rho}(t)]\right),$
(42)

where $\hat{\mathcal{L}}_{2}^{(L)}[\hat{\rho}(t)]$
($\hat{\mathcal{L}}_{2}^{(R)}[\hat{\rho}(t)]$) encodes the effect of the left
(right) bath, the full superoperator being
$\hat{\mathcal{L}}_{2}[\hat{\rho}(t)]=\hat{\mathcal{L}}_{2}^{(L)}[\hat{\rho}(t)]+\hat{\mathcal{L}}_{2}^{(R)}[\hat{\rho}(t)]$.
Using this, we can write down the equation for rate of change of energy in the
system,

$\displaystyle\frac{d\braket{\hat{H}_{S}}}{dt}=J_{B_{L}}+J_{B_{R}},$
$\displaystyle
J_{B_{L}}=\epsilon^{2}Tr\left(\hat{H}_{S}{\mathcal{L}}_{2}^{(L)}[{\rho}(t)]\right),$
(43) $\displaystyle
J_{B_{R}}=\epsilon^{2}Tr\left(\hat{H}_{S}{\mathcal{L}}_{2}^{(R)}[{\rho}(t)]\right),$
(44)

where $J_{B_{L}}$ ($J_{B_{R}}$) can be interpreted as the energy current from
the left (right) bath into the system.

Next, let us look at the relation between the currents from the baths, and the
currents in the system. Due to local conservation of energy, following the
continuity equations must be satisfied

$\displaystyle\frac{d\braket{\hat{H}_{L}}}{dt}=J_{B_{L}}-J_{L\rightarrow M},$
$\displaystyle\frac{d\braket{\hat{H}_{R}}}{dt}=J_{B_{R}}+J_{M\rightarrow R},$
(45)

where $J_{L\rightarrow M}$ and $J_{M\rightarrow R}$ are the defined in
Eq.(II.2). Satisfying this condition requires, along with Eq.(9), the
following relation is satisfied,

$\displaystyle
Tr\left(\hat{H}_{\ell}\hat{\mathcal{L}}_{2}^{(L)}[{\rho}(t)]\right)+Tr\left(\hat{H}_{\ell}\hat{\mathcal{L}}_{2}^{(R)}[{\rho}(t)]\right)$
$\displaystyle=Tr\left(\hat{H}_{S}\hat{\mathcal{L}}_{2}^{(\ell)}[{\rho}(t)]\right),~{}\ell=L,R.$
(46)

From Eqs.(42) and (44), we can show that this condition is indeed satisfied
and the local continuity equation relating the currents from the bath and the
currents inside the system holds. At NESS, since the rate of change of all the
system operators goes to zero, we have

$\displaystyle J_{B_{L}}=J_{L\rightarrow M}=J_{M\rightarrow R}=-J_{B_{R}}\geq
0.$ (47)

This, once again, is a fundamental property of NESS, which is preserved by the
RE.

The above discussion is with energy currents, but it holds true for any other
currents associated with other local conserved quantities, for example,
particle currents.

###  III.4 Accuracy of currents from RE

Although the local continuity equations are always satisfied by RE, the
currents obtained from RE are only accurate to $O(\epsilon^{2})$. This can be
seen noting that the currents inside the system are given by coherences in the
energy eigenbasis, as shown by Eqs (II.2), (II.2). To see this explicitly, we
write $J_{L\rightarrow M}$ in steady state as

$\displaystyle J_{L\rightarrow
M}=-i\braket{[\hat{H}_{M},\hat{H}_{L}]}=-i\braket{[\hat{H}_{S},\hat{H}_{L}]}$
$\displaystyle=i\sum_{\alpha,\nu=1}^{D}(E_{\alpha}-E_{\nu})\langle
E_{\alpha}|\hat{\widetilde{{\rho}}}_{\rm NESS}|E_{\nu}\rangle\langle
E_{\nu}|\hat{H}_{L}|E_{\alpha}\rangle$
$\displaystyle=i\epsilon^{2}\sum_{\alpha,\nu=1}^{D}(E_{\alpha}-E_{\nu})\langle
E_{\alpha}|\hat{\widetilde{{\rho}}}_{\rm NESS}^{(2)}|E_{\nu}\rangle\langle
E_{\nu}|\hat{H}_{L}|E_{\alpha}\rangle$ $\displaystyle+\textrm{higher orders}.$
(48)

In the last line, we have used the fact that the leading order coherences are
$O(\epsilon^{2})$. Since the coherences obtained from RE are given correctly
only to $O(\epsilon^{2})$, any higher order term obtained from the RE cannot
be trusted. Similar expressions can be written for $J_{M\rightarrow R}$ or any
other currents related to other local conserved quantities.

The fact that currents obtained from RE are correct to $O(\epsilon^{2})$ can
also be seen by calculating the currents from the baths. Carrying out the
trace in Eq.(44) in the energy eigenbasis of the system, the currents from the
baths can be written as

$\displaystyle J_{B_{\ell}}$
$\displaystyle=\epsilon^{2}\sum_{\alpha}E_{\alpha}\left\langle
E_{\alpha}\left|{\mathcal{L}}_{2}^{(\ell)}[\hat{\widetilde{{\rho}}}_{\rm
NESS}^{(0)}]\right|E_{\alpha}\right\rangle$
$\displaystyle+\epsilon^{4}\sum_{\alpha}E_{\alpha}\left\langle
E_{\alpha}\left|\hat{\mathcal{L}}_{2}^{(\ell)}[\hat{\widetilde{{\rho}}}_{\rm
NESS}^{(2)}]\right|E_{\alpha}\right\rangle$ $\displaystyle+\textrm{higher
orders},\textrm{ where}~{}\ell=L,R.$ (49)

This shows that the leading order currents from the baths are
$O(\epsilon^{2})$ and depend on $\hat{\widetilde{{\rho}}}_{\rm NESS}^{(0)}$
which, as shown before, is diagonal in the energy eigenbasis and is given
correctly by the RE. In other words, the leading order currents from the baths
are determined by the leading order populations in the energy eigenbasis. On
the other hand, the $O(\epsilon^{4})$ term contains contributions from both
the diagonal and the off-diagonal elements of $\hat{\widetilde{{\rho}}}_{\rm
NESS}^{(2)}$ in the energy eigenbasis of the system. Since the diagonal
elements of $\hat{\widetilde{{\rho}}}_{\rm NESS}^{(2)}$ in the energy
eigenbasis of the system are not given correctly by RE, currents from RE will
carry $O(\epsilon^{4})$ error.

A particularly critical condition arises when, on physical grounds, the
currents are zero, for example, when the temperatures and the chemical
potentials of the baths are same, $\beta_{L}=\beta_{R}=\beta$,
$\mu_{L}=\mu_{R}=\mu$. In this case, the $O(\epsilon^{2})$ contribution to
currents obtained from RE will be zero, since this is given correctly. But,
the $O(\epsilon^{4})$ term, which can contain an error, may not be zero. This
will give unphysical $O(\epsilon^{4})$ currents even when the current is
expected to be zero on physical grounds. Knowing this, it is however possible
to clearly identify such spurious unphysical currents stemming from
inaccuracies of RE by checking the scaling of the currents with $\epsilon$. If
the scaling is $O(\epsilon^{2})$, the currents obtained from RE can be
trusted. If the scaling in $O(\epsilon^{4})$, it shows that the
$O(\epsilon^{2})$ contribution to currents is zero. The currents must be taken
as zero within the accuracy of RE in that case. This was explicitly used in
Ref.Purkayastha  _et al._ (2020) (see the supplemental material of
Ref.Purkayastha  _et al._ (2020)).

Another interesting point to note from Eqs.(III.4),(III.4) and (47) is that
the $O(\epsilon^{2})$ coherences, the $O(\epsilon^{0})$ populations in energy
eigenbasis and the preservation of local conservation laws are tightly related
to one another. If either one of the populations and the coherences are given
accurately, while the other is not, it will lead to violation of local
continuity equations.

It is clear from the discussion in this section that the only drawback of the
RE is the lack of complete positivity. In the next section, we look in detail
at one recent attempt to rigorously rectify this drawback, and show that it
violates some of the other required conditions in Table.1. We also make some
general comments on all existing weak system-bath coupling Lindblad
descriptions.

##  IV The Universal Lindblad Equation

As mentioned before, we need to go beyond the standard Lindblad equations in
local and eigenbasis forms, LLE and ELE, to satisfy the required conditions in
Table.1. Since the RE is microscopically derived and satisfies all other
requirements except for the complete positivity, it is intuitive that any QME
satisfying the required conditions should give results close to RE. Recently,
a Lindblad equation, called the Universal Lindblad Equation (ULE) Nathan and
Rudner (2020) has been rigorously derived such that it is of Lindblad form,
while maintaining

$\displaystyle||\hat{{\rho}}_{\rm NESS}^{\rm ULE}-\hat{\rho}_{\rm
NESS}||\sim||\hat{\widetilde{{\rho}}}_{\rm NESS}-\hat{\rho}_{\rm NESS}||\sim
O(\epsilon^{2}),$ (50)

where $\hat{{\rho}}_{\rm NESS}^{\rm ULE}$ is the NESS density matrix obtained
from ULE. The ULE reduces to the LLE and the ELE with further approximations
in appropriate limits Nathan and Rudner (2020). Moreover, as discussed in
Nathan and Rudner (2020); Davidović (2020) it is closely related to other
recently derived Lindblad equations to go beyond ELE and LLE. This makes the
study of ULE in the light of the fundamental requirements representative of
existing Lindblad equation approaches.

###  IV.1 The general form of ULE

The ULE approach requires the system-bath coupling to be written in terms of
Hermitian operators. So we write the system-bath coupling as,
$\hat{H}_{SB}=\sum_{\ell}^{\prime}\sum_{k=1,2}\hat{X}_{(\ell,k)}\hat{\mathcal{B}}_{(\ell,k)}$,
where

$\displaystyle\hat{X}_{(\ell,1)}=\hat{S}_{\ell}^{\dagger}+\hat{S}_{\ell},~{}\hat{X}_{(\ell,2)}=i(\hat{S}_{\ell}-\hat{S}_{\ell}^{\dagger}),$
$\displaystyle\hat{\mathcal{B}}_{(\ell,1)}=\sum_{r=1}^{\infty}\frac{\kappa_{r\ell}\hat{B}_{r}^{(\ell)\dagger}+\kappa_{r\ell}^{*}\hat{B}_{r}^{(\ell)}}{2},$
$\displaystyle\hat{\mathcal{B}}_{(\ell,2)}=i\sum_{r=1}^{\infty}\frac{\kappa_{r\ell}^{*}\hat{B}_{r}^{(\ell)}-\kappa_{r\ell}\hat{B}_{r}^{(\ell)\dagger}}{2}.$
(51)

The system-bath coupling is then given in the form

$\displaystyle\hat{H}_{SB}=\sum_{\lambda}\hat{X}_{\lambda}\hat{\mathcal{B}}_{\lambda},$
(52)

where $\lambda=(\ell,k)$ is the combined index. The ULE is of the form of

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

with the Lindblad operators and the Lamb-shift Hamiltonian being given by

$\displaystyle\hat{L}_{\lambda}=\sum_{\lambda^{\prime}}\int_{-\infty}^{\infty}ds\widetilde{\bm{g}}_{\lambda,\lambda^{\prime}}(s)\hat{X}_{\lambda^{\prime}}(-s)$
(54) $\displaystyle\hat{H}_{\rm
LS}=\frac{1}{2i}\int_{-\infty}^{\infty}ds\int_{-\infty}^{\infty}ds^{\prime}\sum_{\lambda,\lambda^{\prime}}\hat{X}_{\lambda}(s)\hat{X}_{\lambda^{\prime}}(s^{\prime})\bm{\phi}_{\lambda\lambda^{\prime}}(s,s^{\prime})$

where
$\hat{X}_{\lambda}(t)=e^{i\hat{H}_{S}t}\hat{X}_{\lambda}e^{-i\hat{H}_{S}t}$
denotes the interaction picture operator. The matrices $\bm{\phi}(t,s)$ and
$\widetilde{\bm{g}}(t)$ are,

$\displaystyle\widetilde{\bm{g}}(t)$
$\displaystyle=\int_{-\infty}^{\infty}d\omega\bm{g}(\omega)e^{-i\omega
t},~{}\bm{\phi}(t,s)=\widetilde{\bm{g}}(t)\widetilde{\bm{g}}(-s)\text{sgn}(t-s),$
$\displaystyle\bm{g}(\omega)$
$\displaystyle=\sqrt{\frac{\bm{G}(\omega)}{2\pi}},~{}\bm{G}_{\lambda,\lambda^{\prime}}(\omega)=\int_{-\infty}^{\infty}\frac{dt}{2\pi}\braket{\hat{\mathcal{B}}_{\lambda}(t)\hat{\mathcal{B}}_{\lambda^{\prime}}(0)}_{B}e^{i\omega
t}$ (55)

with $\text{sgn}(x)$ being the sign function, and $\braket{...}_{B}$ denoting
the expectation value taken over the bath initial state. Given this general
form the ULE we now investigate if the ULE satisfies the fundamental
requirements in Table 1.

###  IV.2 Accuracy of ULE

####  IV.2.1 Accurate populations but inaccurate coherences

As mentioned before, the ULE has been derived so that the steady state it
predicts differs from the exact steady state in $O(\epsilon^{2})$. This means
that the $O(\epsilon^{0})$ term of the NESS density matrix, must be given
correctly by the ULE, and hence would match with that from RE. Here we explain
how to see this explicitly in general.

For this, we write $\hat{\widetilde{\rho}}_{\rm NESS}^{(0)}$ in the energy
eigenbasis, as in Eq.(18) and explicitly evaluate the term required by the
first condition in Table 1,

$\displaystyle\left\langle
E_{\alpha}\left|\hat{\mathcal{L}}_{2}^{\prime}[\hat{\widetilde{{\rho}}}_{\rm
NESS}^{(0)}]\right|E_{\alpha}\right\rangle=$
$\displaystyle\sum_{\nu,\lambda,\lambda^{\prime}}2\pi
p_{\nu}\Big{\\{}\bra{E_{\alpha}}\hat{X}_{\lambda^{\prime}}\ket{E_{\nu}}\bra{E_{\nu}}\hat{X}_{\lambda}\ket{E_{\alpha}}\bm{G}_{\lambda,\lambda^{\prime}}(E_{\nu}-E_{\alpha})$
$\displaystyle-\sum_{\gamma}\bra{E_{\alpha}}\hat{X}_{\lambda}\ket{E_{\gamma}}\bra{E_{\gamma}}\hat{X}_{\lambda^{\prime}}\ket{E_{\alpha}}\bm{G}_{\lambda,\lambda^{\prime}}(E_{\alpha}-E_{\gamma})\delta_{\alpha\nu}\Big{\\}}.$
(56)

We need to compare this result with that obtained from RE to check the first
condition in Table 1. To do this, it is easier to write the RE also in the
form where the system-bath coupling terms are Hermitian operators, as in
Eq.(52). The RE for this form of coupling is given by

$\displaystyle\frac{\partial\hat{\rho}}{\partial
t}=i[\hat{\rho},\hat{H}_{S}]-\epsilon^{2}\sum_{\lambda,\lambda^{\prime}}\Big{(}\big{[}\hat{X}_{\lambda},\hat{X}_{\lambda^{\prime}\lambda}\hat{\rho}\big{]}+{\rm
h.c.}\Big{)},$ (57)

with

$\displaystyle\hat{X}_{\lambda^{\prime}\lambda}=\int_{0}^{\infty}dt^{\prime}\langle\hat{\mathcal{B}}_{\lambda}(t^{\prime})\hat{\mathcal{B}}_{\lambda^{\prime}}(0)\rangle\hat{X}_{\lambda^{\prime}}(-t^{\prime}),$
$\displaystyle=\int_{0}^{\infty}dt^{\prime}\int_{-\infty}^{\infty}\bm{G}_{\lambda,\lambda^{\prime}}(\omega)e^{-i\omega
t^{\prime}}\hat{X}_{\lambda^{\prime}}(-t^{\prime}).$ (58)

Using Eq.(IV.1), it can be checked, that Eq.(57) can be reduced to Eq.(III.2).
With Eq.(57), direct evaluation of $\left\langle
E_{\alpha}\left|\hat{\mathcal{L}}_{2}[\hat{\widetilde{{\rho}}}_{\rm
NESS}^{(0)}]\right|E_{\alpha}\right\rangle$ gives exactly the same expression
as Eq.(IV.2.1). Thus the ULE explicitly satisfies the condition first
condition in Table 1. It is important to note that we do not explicitly
require the values of the populations $p_{\nu}$ to show that the condition is
satisfied.

Using exactly similar steps, it can be checked that the ULE does not satisfy
the second condition in Table 1. Thus, the leading order coherences from ULE
are not given correctly. Since the leading order coherences in energy
eigenbasis are $O(\epsilon^{2})$, this is consistent with the ULE having an
error of $O(\epsilon^{2})$ (see Eq.(50)).

Figure 2:  The top row [(a), (b), (c)] is for the equilibrium case
($\beta_{L}=1$, $\beta_{R}=1$) and the bottom row [(d), (e), (f)] is for the
non-equilibrium case ($\beta_{L}=5$, $\beta_{R}=0.5$). Apart from (c), in all
other cases, $\omega_{0}^{(1)}=\omega_{0}^{(2)}=\omega_{0}^{(3)}=1$. (a) The
trace distance between the expected thermal state $\hat{\rho}_{\rm th}$ and
the equilibrium steady-state $\hat{\rho}_{\rm eq}$ as obtained by the LLE, ULE
and RE as a function of $g$. (b) The local magnetization in equilibrium steady
state as obtained by LLE, ULE, RE, ELE and by the expected thermal state as a
function of $g$. (c) The scaling of spin currents in equilibrium as given by
LLE, ULE, RE as a function of system-bath coupling strength $\epsilon$, for
the case where $\omega_{0}^{(1)}=1,\omega_{0}^{(2)}=1.5,\omega_{0}^{(3)}=2$.
Here $I_{1}^{ULE}$, $I_{2}^{ULE}$ and $I_{B}^{ULE}$ refer to the two bond-
currents and the boundary current as obtained by ULE. (d) Spin currents in
NESS as obtained by LLE, ULE, RE and ELE. For ULE and ELE, only the boundary
currents $I_{B}^{ULE}$ and $I_{B}^{ELE}$ are plotted. (For ELE, the bond
currents are exactly zero.) (e) The local magnetization in NESS as obtained by
LLE, ULE, RE, ELE as a function of $g$. (f) The scaling of difference between
the two bond currents at NESS from ULE as a function of system-bath coupling
strength $\epsilon$. Other parameters: $\mu_{L}=\mu_{R}=-0.5$,
$\omega_{c}=10$. All energy parameters are in units of $\omega_{0}^{(1)}$.

####  IV.2.2 Violation of local conservation laws

As we have already discussed, having correct populations to leading order
while having incorrect coherences leads to violation of local conservation
laws. So, the ULE does not satisfy the fourth requirement in Table. 1. This
can be explicitly seen by writing the evolution equation for any system
operator,

$\displaystyle\frac{d\braket{\hat{O}}}{dt}$
$\displaystyle=-i\braket{[\hat{O},\hat{H}_{S}+\hat{H}_{\rm LS}]}$
$\displaystyle+\sum_{\lambda}\frac{1}{2}\left(\braket{[\hat{L}_{\lambda}^{\dagger},\hat{O}]\hat{L}_{\lambda}}+\braket{\hat{L}_{\lambda}^{\dagger}[\hat{O},\hat{L}_{\lambda}]}\right).$
(59)

From the form of the operators $\hat{L}_{\lambda}$ in Eq.(54), we see that
even if $[\hat{O},\hat{H}_{SB}]=0$, ${\rm
Tr}\left(\hat{O}\hat{\mathcal{L}}_{2}^{\prime}(\hat{\rho})\right)\neq 0$.

####  IV.2.3 Accuracy of currents from ULE

Despite this violation of local conservation laws, the ULE gives the correct
currents from the baths to the leading order. This can be seen by writing the
ULE in the form of Eq.(42) and defining the corresponding currents as in
Eq.(44). The currents can then be written in energy eigenbasis, as in
Eq.(III.4), to see that the leading order contributions are given by the
$O(\epsilon^{0})$ populations in the energy eigenbasis. Since these terms are
given accurately by the ULE, it follows that the currents from the baths are
given accurately to $O(\epsilon^{2})$ and agree with the same currents
obtained from RE. The error in the currents from the baths as obtained from
ULE is therefore $O(\epsilon^{4})$, which is same as that from RE.

However, since the coherences in energy eigenbasis contain $O(\epsilon^{2})$
error in ULE, the currents inside the system, defined in Eq.(II.2), will have
an error in the leading order. This, once again, shows that local conservation
laws will be violated.

###  IV.3 General comments regarding Lindblad equations

We end this section with some general comments regarding Lindblad equations.
We have already mentioned that the ULE is closely related to some of the other
Lindblad equations devised to go beyond the limitations of ELE and LLE. It is
worth mentioning that, to our knowledge, all existing Lindblad equations
Nathan and Rudner (2020); Kleinherbers  _et al._ (2020); Davidović (2020);
Mozgunov and Lidar (2020); McCauley  _et al._ (2020); Kiršanskas  _et al._
(2018) except the LLE violate the fourth condition in Table. 1, thereby
violating local conservation laws. The LLE on the other hand is known to not
satisfy the thermalization requirement, the third condition in Table 1.
Further, although explicitly satisfying the local conservation laws, it does
not always give the correct currents in NESS Walls (1970); Wichterich  _et
al._ (2007); 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 (2020);
Scali  _et al._ (2021); Benatti  _et al._ (2021); Trushechkin (2021). This is
contrary to the ULE, which despite violating local conservation laws, gives
correct currents from the baths to the leading order. The RE, on the other
hand, except for generically violating the requirement of complete positivity,
satisfies all other the conditions in Table 1. It therefore seems that a weak
system-bath coupling QME respecting all the conditions in Table 1 is
generically impossible. In particular, since currents inside the system are
related to coherences and currents from the baths are related to populations,
it seems that, no weak system-bath coupling QME can generically satisfy the
requirement of complete positivity and simultaneously give correct populations
and coherences to the leading order. Thus all weak system-bath coupling
Lindblad equations seem fundamentally limited. While this is true for the
existing QMEs to our knowledge, a general proof of the fact is missing till
now. We leave the general proof of this to future work. Although all the
results above have been derived keeping a set-up of the form of Fig. 1 in
mind, these can be easily carried forward to cases where there are more than
two baths.

Our discussion till now has been general. We have not yet written down any
particular system Hamiltonian, or chosen any particular bath spectral
function. This completely general discussion of accuracy of QMEs sets apart
our work from most previous worksWalls (1970); Wichterich  _et al._ (2007);
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 (2020); Scali  _et al._ (2021);
Benatti  _et al._ (2021); Trushechkin (2021), where the accuracy is discussed
referring to a particular model. In the following, we numerically check our
general discussion in the three site Heisenberg model coupled to bosonic
baths.

Figure 3:  The top row [(a), (b), (c)] is for the equilibrium case
($\beta_{L}=1$, $\beta_{R}=1$) and the bottom row [(d), (e), (f)] is for the
non-equilibrium case ($\beta_{L}=5$, $\beta_{R}=0.5$). Apart from (c),
everywhere, $\omega_{0}^{(1)}=1,\omega_{0}^{(2)}=1.5,\omega_{0}^{(3)}=2$. (a)
The trace distance between the expected thermal state $\hat{\rho}_{th}$ and
the equilibrium steady-state $\hat{\rho}_{eq}$ as obtained by the local
Lindblad (LLE), Universal Lindblad (ULE) and and Redfield (RE) equations as a
function of $g$. (b) The local magnetization in equilibrium steady state as
obtained by LLE, ULE, RE, ELE and by the expected thermal state as a function
of $g$. (c) The scaling of bond spin currents in equilibrium as obtained from
ULE, $I_{1}^{ULE}$, $I_{2}^{ULE}$, for the case where
$\omega_{0}^{(1)}=\omega_{0}^{(2)}=\omega_{0}^{(3)}=1$. The LLE, RE and ULE
boundary currents are zero correct to numerical precision of $10^{-16}$. (d)
Spin currents in NESS as obtained by LLE, ULE, RE and ELE. For ULE and ELE,
the boundary currents $I_{B}^{ULE}$ and $I_{B}^{ELE}$ are plotted. (e) The
local magnetization in NESS as obtained by LLE, ULE, RE, ELE as a function of
$g$. (f) The scaling of difference between the two bond currents at NESS from
ULE as a function of system-bath coupling strength $\epsilon$. Other
parameters: $\mu_{1}=\mu_{2}=-0.5$, $\omega_{c}=10$. All energy parameters are
in units of $\omega_{0}^{(1)}$.

##  V Numerical Results

We now numerically exemplify the above discussion by using a XXZ spin-chain in
the presence of a magnetic field, with the first and the last sites attached
to baths modelled by infinite number of bosonic modes,

$\displaystyle\hat{H}_{S}=\sum_{\ell=1}^{N}\frac{\omega^{(\ell)}_{0}}{2}\hat{\sigma}_{z}^{\ell}-\sum_{\ell=1}^{N-1}g(\hat{\sigma}_{x}^{\ell}\hat{\sigma}_{x}^{\ell+1}+\hat{\sigma}_{y}^{\ell}\hat{\sigma}_{y}^{\ell+1}+\Delta\hat{\sigma}_{z}^{\ell}\hat{\sigma}_{z}^{\ell+1}),$
$\displaystyle\hat{H}_{SB}=\sum_{\ell=1,N}\sum_{r=1}^{\infty}(\kappa_{\ell
r}\hat{B}^{{(\ell)}\dagger}_{r}\hat{\sigma}^{\ell}_{-}+\kappa_{lr}^{*}\hat{B}^{(\ell)}_{r}\hat{\sigma}^{\ell}_{+}),$
(60)

where $\hat{\sigma}^{\ell}_{x,y,z}$ denotes the Pauli matrices acting on the
$\ell^{\text{th}}$ spin,
$\hat{\sigma}^{\ell}_{+}=(\hat{\sigma}^{\ell}_{x}+i\hat{\sigma}^{\ell}_{y})/2$,
$\hat{\sigma}^{\ell}_{-}=(\hat{\sigma}^{\ell}_{x}-i\hat{\sigma}^{\ell}_{y})/2$,
$\hat{B}^{(\ell)}_{r}$ is bosonic annihilation operator for the
$r^{\text{th}}$ mode of the bath attached at the $\ell^{\text{th}}$ site.
Here, $\omega^{(\ell)}_{0}$, $g$, and $g\Delta$ represent the magnetic field,
the overall spin-spin coupling strength and the anisotropy respectively. The
total number of excitations in the system is given by
$\hat{N}_{S}=\sum_{\ell=1}^{N}\hat{\sigma}^{\ell}_{+}\hat{\sigma}^{\ell}_{-}$
and it satisfies $[\hat{N}_{S},\hat{H}_{S}]=0$. We consider bosonic baths
described by Ohmic spectral functions with Gaussian cut-offs,
$\mathfrak{J}_{\ell}(\omega)=\sum_{k}2\pi\left|\kappa_{\ell
k}\right|^{2}\delta(\omega-\Omega^{\ell}_{k})=\omega
e^{-(\omega/\omega_{c})^{2}}\Theta(\omega),$ where, $\Theta(\omega)$ is
Heaviside step function, and $\omega_{c}$ is the cut-off frequency. The set-up
is exactly of the form of Fig. 1, with the initial inverse temperatures and
chemical potentials of the baths as shown. We look at the equilibrium and non-
equilibrium steady states as obtained by RE, LLE, ELE, and ULE. The explicit
forms of all these equations for the XXZ spin-chain are given in the Appendix
(Appedices. A, B, C and D respectively). For numerical simplicity, we consider
the $N=3$ case, which is sufficient to demonstrate the fundamental limitations
of all the three Lindblad equations. All numerical results below are obtained
using QuTiP Johansson  _et al._ (2012, 2013).

| LLE | ELE | ULE | RE
---|---|---|---|---
Diagonal (Populations) | O(1) wrong | O(1) wrong in NESS | O(1) correct | O(1) correct
Off-diagonal (Coherences) | O($\epsilon^{2}$) wrong | O($\epsilon^{2}$) wrong | O($\epsilon^{2}$) wrong | O($\epsilon^{2}$) correct
Conservation Laws | Conserved | Violated | Violated | Conserved
Thermalization | Violated | Preserved | Preserved | Preserved
Currents from baths | valid for $g\ll\epsilon$ | valid for $|E_{\alpha}-E_{\gamma}|\gg\epsilon$ | valid for all $g$ | valid for all $g$
Currents in system | valid for $g\ll\epsilon$ | zero | O($\epsilon^{2}$) wrong | O($\epsilon^{2}$) correct
Complete positivity | Preserved | Preserved | Preserved | Violated
Table 2: Accuracy of various QMEs in various settings.

First, we look at the equilibrium case (top row of Fig. 2),
$\beta_{L}=\beta_{R}=\beta$, $\mu_{L}=\mu_{R}=\mu$, and calculate the trace
distance Nielsen and Chuang (2011), $D(\hat{\rho}_{\rm eq},\hat{\rho}_{\rm
th})={\rm Tr}[\sqrt{(\hat{\rho}_{\rm eq}-\hat{\rho}_{\rm th})^{2}}]/2$,
between the steady state $\hat{\rho}_{\rm eq}$ given by the three QMEs and the
thermal state $\hat{\rho}_{\rm
th}~{}=~{}\frac{e^{-\beta(\hat{H}_{S}-\mu\hat{N}_{S})}}{{\rm
Tr}\left(e^{-\beta(\hat{H}_{S}-\mu\hat{N}_{S})}\right)}$. The trace distance
$D(\hat{\rho}_{\rm eq},\hat{\rho}_{\rm th})$ is plotted as a function of $g$
for $\epsilon=0.1$ in Fig. 2(a). We see that the trace distance is quite small
and of the same order for RE and ULE, and does not vary much with $g$. On the
other hand, the trace distance for LLE is larger and grows with $g$ while that
for ELE is identically zero. Note that generically any system is expected to
have steady-state coherences in energy eigenbasis in equilibrium Guarnieri
_et al._ (2018); Purkayastha  _et al._ (2020) for any finite $\epsilon$ and
convergence to $\hat{\rho}_{\rm th}$ is physically only expected for
$\epsilon\rightarrow 0$. Nevertheless, as shown in Fig 2(b), observables like
local magnetizations, $\braket{\hat{\sigma}^{\ell}_{z}}$, obtained from ULE
and RE show very small difference from the thermal expectation values, since
they depend on populations in leading order.

The same is not true for the spin currents. The bond spin currents in the
system $I_{j}$ are defined from the continuity equation
$\frac{d\braket{\hat{\sigma}^{j}_{z}}}{dt}=-i\braket{[\hat{\sigma}^{j}_{z},\hat{H}_{S}]}=I_{j}-I_{j-1}$,
which gives
$I_{j}=4ig(\braket{\hat{\sigma}^{j}_{+}\hat{\sigma}^{j+1}_{-}}-\braket{\hat{\sigma}^{j}_{-}\hat{\sigma}^{j+1}_{+}}).$
On the other hand, boundary spin-currents are defined by the continuity
equation, $\frac{d\braket{\hat{M}_{z}}}{dt}=I_{B_{L}}+I_{B_{R}}$, where
$\hat{M}_{z}=\sum_{\ell=1}^{N}\hat{\sigma}^{\ell}_{z}$. For equations of the
form Eq.(42),
$I_{B_{\ell}}=\epsilon^{2}Tr\left(\hat{M}_{z}\hat{\mathcal{L}}_{2}^{(\ell)}[\hat{\rho}(t)]\right)$,
$\ell=L,R$ where $\hat{\mathcal{L}}_{2}^{(\ell)}$ changes depending on the QME
used. In steady state, local conservation laws require
$I_{1}=-I_{B_{L}}=I_{2}=I_{B_{R}}$. Further, in equilibrium, currents should
be zero Levy and Kosloff (2014). We find that, for
$\omega_{0}^{(1)}\neq\omega_{0}^{(2)}\neq\omega_{0}^{(3)}$, all QMEs except
ELE gives unphysical currents in even though we have
$\beta_{L}=\beta_{R}=\beta$, $\mu_{L}=\mu_{R}=\mu$. However, it is paramount
to note that, as shown in Fig 2(c), the currents from RE and the boundary
currents from ULE scale as $\epsilon^{4}$, clearly showing that the leading
order, $O(\epsilon^{2})$, term is zero. This is completely consistent with our
discussion in Sec. III.4 and Sec. IV.2.3. On the other hand, the bond currents
from ULE are different, and scale as $\epsilon^{2}$, thereby showing error in
leading order, and violating local conservation laws. The LLE does not violate
local conservation laws, but still gives equilibrium currents scaling as
$\epsilon^{2}$ Levy and Kosloff (2014), thereby showing leading order
inaccuracies in both populations and coherences.

Now we discuss the results for the non-equilibrium set-up (bottom row of Fig.
2), $\beta_{L}\neq\beta_{R}$. The currents in NESS as a function of $g$ are
shown in Fig. 2(d). The current from RE matches the boundary current from ULE
for all $g$ and shows a non-monotonic behavior with $g$. On the other hand,
LLE fail to capture the non-monotonicity of current as a function of $g$, and
matches only at very small values of $g$. The ELE identically gives zero
currents inside the system, even in NESS Wichterich  _et al._ (2007). However,
the boundary currents from ELE seem to match reasonably well with RE for
sufficiently high $g$. Thus, the ELE also violates the local conservation laws
in NESS. The local magnetizations in NESS also match from RE and ULE and are
significantly different from the LLE and ELE results at high and low $g$
respectively, as shown in Fig. 2(e). Finally, in Fig. 2(f), we demonstrate
that ULE bond currents are different in NESS with the difference scaling as
$O(\epsilon^{2})$, which highlights a clear violation of local conservation
laws in ULE.

All plots in Fig. 2 are made with
$\omega_{0}^{(1)}=\omega_{0}^{(2)}=\omega_{0}^{(3)}=1$, except for panel (c)
which was plotted for
$\omega_{0}^{(1)}=1,\omega_{0}^{(2)}=1.5,\omega_{0}^{(3)}=2$. In Fig. 3, we
give the complimentary plot, where all numerical results execpt panel (c) are
for $\omega_{0}^{(1)}=\omega_{0}^{(2)}=\omega_{0}^{(3)}=1$, while panel (c) is
for $\omega_{0}^{(1)}=1,\omega_{0}^{(2)}=1.5,\omega_{0}^{(3)}=2$. Two
important things are worth noticing. First, when
$\omega_{0}^{(1)}=\omega_{0}^{(2)}=\omega_{0}^{(3)}$, none of LLE, RE and ULE
boundary currents give any unphysical non-zero result (correct to a numerical
precision of $10^{-16}$) if $\beta_{L}=\beta_{R}=\beta$,
$\mu_{L}=\mu_{R}=\mu$. However, the ULE bond currents still gives non-zero
$O(\epsilon^{2})$ result and thereby violate local conservation laws. Second,
for $\omega_{0}^{(1)}=1,\omega_{0}^{(2)}=1.5,\omega_{0}^{(3)}=2$, because of
already large separation between the on-site magnetic fields, the secular
approximation performs much better, and as a result, ELE gives reasonably
accurate results both in equilibrium and out-of-equilibrium steady states, for
a wide range of $g$. However, ELE, by construction, gives zero currents inside
the system Wichterich  _et al._ (2007). Table 2 summarizes the accuracies and
validity regimes of the various QMEs. Although our numerical demonstration is
for $N=3$, our analytical understanding in previous sections shows that these
issues persist for larger $N$, as well as generically for other systems.

##  VI Summary and discussions

In this work, we have investigated fundamental limitations of QMEs obtained in
weak system-bath coupling approach to describe a quantum system coupled to
multiple macroscopic baths which can be at different temperatures and chemical
potentials, in absence of any external driving (i.e, the when the Hamiltonian
is time-independent). Microscopically obtaining the QME up to leading order in
system-bath coupling leads to the so-called Redfield form, which is known to
generically violate the requirement of complete positivity. Very often further
approximations are done on the Redfield form to reduce it to a Lindblad form,
so that complete positivity is restored. We have laid down some fundamental
requirements (Table. 1) that any such approximation must satisfy. If they are
violated there occurs physical inconsistencies like inaccuracies in leading
order populations and coherences in the energy eigenbasis, violation of
thermalization when all baths have same temperatures and chemical potentials,
violation of local conservation laws at NESS. This has allowed us to check for
occurrence of these physical inconsistencies in weak-system-bath coupling QMEs
in general, without writing them down for any particular model. This model
independent general discussion distinguishes our work from most previous works
investigating accuracies of various QMEs Walls (1970); Wichterich  _et al._
(2007); 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 (2020); Scali  _et al._
(2021); Benatti  _et al._ (2021).

We have found that the Redfield equation, despite being generically not
completely positive, does not violate any other of the fundamental
requirements. It therefore gives correct populations and coherences in energy
eigenbasis to the leading order, shows thermalization, and preserves local
conservation laws. On the other hand, weak system-bath coupling Lindblad
descriptions violate one or more the fundamental requirements and thereby have
one or more of the physical inconsistencies mentioned above, despite being
completely positive. In particular, it can be argued that no existing weak
system-bath coupling Lindblad description, to our knowledge, can generically
give both correct populations and correct coherences to the leading order. As
an example, we have explicitly shown these violations in generality for the
so-called universal Lindblad equation which has been recently derived. We have
also numerically demonstrated these statements in a three site XXZ chain
coupled to two bosonic baths. From our results, it seems that there is no
consistent way to correct the violation of complete positivity of the Redfield
equation, without going to higher orders in system-bath coupling.

These results are extremely significant because weak system-bath coupling
Lindblad equations remain the most widely-used descriptions of open quantum
systems, and accurate description of populations, coherences and currents is
important for quantum information and thermodynamics Nielsen and Chuang
(2011); Bennett and DiVincenzo (2000); Streltsov  _et al._ (2017); Lostaglio
_et al._ (2015); Narasimhachar and Gour (2015); Mitchison  _et al._ (2015);
Allahverdyan  _et al._ (2004); Korzekwa  _et al._ (2016); Kammerlander and
Anders (2016); Santos  _et al._ (2019); Francica  _et al._ (2019). It seems
the Redfield equation, despite not being completely positive, should be the
description of choice. The various errors in the Redfield equations are
controlled. Whether a result obtained from the Redfield equation can be
trusted or not, can be determined by checking its scaling with system-bath
coupling strength. Further, there have been techniques suggested to infer the
next order corrections to populations in steady-state from Redfield equation
without actually writing down the next higher order QME Thingna  _et al._
(2012, 2013); Xu  _et al._ (2017). These techniques may be able to correct the
positivity issues in the steady-state obtained from the Redfield equation
Purkayastha  _et al._ (2020).

It is paramount to mention that the fundamental limitations of Lindblad
descriptions discussed in this paper pertain to the case of a system weakly
coupled to multiple thermal baths in absence of any external driving. In
presence of particular kinds of driving, Lindblad equations can be
microscopically derived even for strong system-bath coupling Cattaneo  _et
al._ (2021). Our results do not pertain to those descriptions. Further, our
results do not pertain to cases where the Lindblad dissipators do not act on
the system, but rather on the bath degrees of freedom. Such approaches can be
shown to accurately describe the steady state of the system for arbitrary
strength of system-bath coupling Brenes  _et al._ (2020). Our results,
however, do show that it is imperative to go beyond weak system-bath coupling
QMEs for physically consistent description of open quantum systems.

## Acknowledgements

MK would like to acknowledge support from the project 6004-1 of the Indo-
French Centre for the Promotion of Advanced Research (IFCPAR), Ramanujan
Fellowship (SB/S2/RJN-114/2016), SERB Early Career Research Award
(ECR/2018/002085) and SERB Matrics Grant (MTR/2019/001101) from the Science
and Engineering Research Board (SERB), Department of Science and Technology,
Government of India. AD and MK acknowledge support of the Department of Atomic
Energy, Government of India, under Project No. RTI4001. AP acknowledges
funding from the European Research Council (ERC) under the European Unions
Horizon 2020 research and innovation program (Grant Agreement No. 758403).

## Appendix

The example we looked at in the main text is the three-site Heisenberg model,
coupled to two different bosonic baths at the first and the third sites. Here,
we will give the QMEs used for this set-up, and write them for Heisenberg
spin-chain of length $N$.

##  Appendix A Redfield equation

The Redfield equation for our set-up is given by

$\displaystyle\frac{\partial\hat{\rho}}{\partial t}$
$\displaystyle=i[\hat{\rho},\hat{H}_{S}]$ (61) $\displaystyle-\epsilon^{2}$
$\displaystyle\sum_{\ell=1,N}\sum_{\alpha,\gamma=1}^{2^{N}}\Big{\\{}\left[\hat{\rho}\ket{E_{\alpha}}\bra{E_{\alpha}}\hat{\sigma}^{\ell}_{-}\ket{E_{\gamma}}\bra{E_{\gamma}},\hat{\sigma}^{l}_{+}\right]C_{\ell}(\alpha,\gamma)$
$\displaystyle+$
$\displaystyle\left[\hat{\sigma}^{\ell}_{+},\ket{E_{\alpha}}\bra{E_{\alpha}}\hat{\sigma}^{l}_{-}\ket{E_{\gamma}}\bra{E_{\gamma}}\hat{\rho}\right]D_{\ell}(\alpha,\gamma)+\text{H.c}\Big{\\}}$

with

$\displaystyle C_{\ell}(\alpha,\gamma)$
$\displaystyle=\frac{\mathfrak{J}_{\ell}(E_{\gamma\alpha})n_{\ell}(E_{\gamma\alpha})}{2}-i\mathcal{P}\int_{0}^{\infty}d\omega\frac{\mathfrak{J}_{\ell}(\omega)n_{\ell}(\omega)}{\omega-
E_{\gamma\alpha}},$ (62) $\displaystyle D_{\ell}(\alpha,\gamma)$
$\displaystyle=\frac{e^{\beta_{\ell}(E_{\gamma\alpha}-\mu_{\ell})}\mathfrak{J}_{\ell}(E_{\gamma\alpha})n_{\ell}(E_{\gamma\alpha})}{2}$
$\displaystyle-i\mathcal{P}\int_{0}^{\infty}d\omega\frac{e^{\beta_{\ell}(\omega-\mu_{\ell})}\mathfrak{J}_{\ell}(\omega)n_{\ell}(\omega)}{\omega-
E_{\gamma\alpha}},$

where $\mathcal{P}$ denotes the Cauchy Principal value. Here
$\ket{E_{\alpha}}$ and $\ket{E_{\gamma}}$ are simultaneous eigenkets of the
system hamiltonian and the magnetization operator, with eigenenergies
$E_{\alpha}$ and $E_{\gamma}$, where $E_{\gamma\alpha}=E_{\gamma}-E_{\alpha}$.

##  Appendix B Local Lindblad equation

The local-Lindblad equation, which can be derived microscopically for
$g\ll\epsilon$, from the Redfield equation, is given by,

$\displaystyle\frac{\partial\hat{\rho}}{\partial t}=$ $\displaystyle
i[\hat{\rho},\hat{H}+\hat{H}_{LS}]$ (63) $\displaystyle+\epsilon^{2}$
$\displaystyle\sum_{\ell=1,N}\mathfrak{J}_{\ell}(\omega^{(\ell)}_{0})(n_{\ell}(\omega^{(\ell)}_{0})+1)\left(\hat{\sigma}^{\ell}_{-}\hat{\rho}\hat{\sigma}^{\ell}_{+}-\frac{1}{2}\\{\hat{\sigma}^{\ell}_{+}\hat{\sigma}^{\ell}_{-},\hat{\rho}\\}\right)$
$\displaystyle+$
$\displaystyle\mathfrak{J}_{\ell}(\omega^{(\ell)}_{0})n_{\ell}(\omega^{(\ell)}_{0})\left(\hat{\sigma}^{\ell}_{+}\hat{\rho}\hat{\sigma}^{\ell}_{-}-\frac{1}{2}\\{\hat{\sigma}^{\ell}_{-}\hat{\sigma}^{\ell}_{+},\hat{\rho}\\}\right)$

where
$\hat{H}_{LS}=\sum_{\ell=1,N}(\overline{\Delta}_{\ell}+2\overline{\Delta}^{\prime}_{\ell})\frac{\omega^{(\ell)}_{0}}{2}$,
and $\overline{\Delta}_{\ell}$ and $\overline{\Delta}_{\ell}^{\prime}$ are
given by,

$\displaystyle\overline{\Delta}_{\ell}$
$\displaystyle=\mathcal{P}\frac{\epsilon^{2}}{2\pi}\int_{0}^{\infty}d\omega\frac{\mathfrak{J}_{\ell}(\omega)}{\omega^{(\ell)}_{0}-\omega},$
(64) $\displaystyle\overline{\Delta}_{i}^{\prime}$
$\displaystyle=\mathcal{P}\frac{\epsilon^{2}}{2\pi}\int_{0}^{\infty}d\omega\frac{\mathfrak{J}_{\ell}(\omega)n_{\ell}(\omega)}{\omega^{(\ell)}_{0}-\omega}.$

##  Appendix C Eigenbasis Lindblad equation

The eigenbasis Lindblad equation is obtained from the Redfield equation by the
doing the so-called secular approximation Breuer and Petruccione (2006). This
amounts to an equation of the form

$\displaystyle\frac{\partial\hat{\rho}}{\partial t}$
$\displaystyle=i[\hat{\rho},\hat{H}_{S}]-\epsilon^{2}\sum_{\ell=1,N}\sum_{\alpha,\gamma=1}^{2^{N}}\sum_{\begin{subarray}{c}\alpha^{\prime},\gamma^{\prime}=1\\\
E_{\gamma^{\prime}\alpha^{\prime}}=E_{\gamma\alpha}\end{subarray}}^{2^{N}}$
(65) $\displaystyle\Big{\\{}\big{[}$
$\displaystyle\hat{\rho}\ket{E_{\alpha}}\bra{E_{\alpha}}\hat{\sigma}^{\ell}_{-}\ket{E_{\gamma}}\bra{E_{\gamma}},\ket{E_{\gamma^{\prime}}}\bra{E_{\gamma^{\prime}}}\hat{\sigma}^{\ell}_{+}\ket{E_{\alpha^{\prime}}}\bra{E_{\alpha^{\prime}}}\big{]}\times$
$\displaystyle C_{\ell}$
$\displaystyle(\alpha,\gamma)+\big{[}\ket{E_{\gamma^{\prime}}}\bra{E_{\gamma^{\prime}}}\hat{\sigma}^{\ell}_{+}\ket{E_{\alpha^{\prime}}}\bra{E_{\alpha^{\prime}}},$
$\displaystyle\ket{E_{\alpha}}\bra{E_{\alpha}}\hat{\sigma}^{\ell}_{-}\ket{E_{\gamma}}\bra{E_{\gamma}}\hat{\rho}\big{]}D_{\ell}(\alpha,\gamma)+\text{H.c}\Big{\\}},$

where the $C_{l}$ and $D_{l}$ constants are the same as in Redfield, given by
Eq.(62).

##  Appendix D Universal Lindblad equation

The universal lindblad equation for our setup is derived from the steps
described in Nathan and Rudner (2020). The ULE equation is then given by

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

The lamb shift $\hat{H}_{LS}$ and Lindblad operators $\hat{L}_{(\ell,k)}$ are
given by,

$\displaystyle\hat{L}_{(\ell,k)}=$
$\displaystyle{2\pi\epsilon}\sum_{\alpha,\gamma,\ell^{\prime},k^{\prime}}\bm{g}_{(\ell,\ell^{\prime},k,k^{\prime})}(E_{\gamma\alpha})\bm{X}^{(\ell^{\prime},k^{\prime})}_{\alpha\gamma}\ket{E_{\alpha}}\bra{E_{\gamma}}$
(67) $\displaystyle\hat{H}_{LS}=$ $\displaystyle\sum_{\alpha\gamma\eta;\ell
k\ell^{\prime}k^{\prime}}\bm{X}^{(\ell,k)}_{\alpha\eta}\bm{X}^{(\ell^{\prime},k^{\prime})}_{\eta\gamma}\bm{f}_{\ell,\ell^{\prime},k,k^{\prime}}(E_{\alpha\eta},E_{\gamma\eta})\ket{E_{\alpha}}\bra{E_{\gamma}}$

where
$\bm{X}^{(\ell,k)}_{\alpha\gamma}=\bra{E_{\alpha}}\hat{X}_{(\ell,k)}\ket{E_{\gamma}}$,
and

$\bm{f}(p,q)=-2\pi\epsilon^{2}\mathcal{P}\int_{-\infty}^{\infty}d\omega\frac{\bm{g}(\omega-p)\bm{g}(\omega+q)}{\omega}$
(68)

Note that $\bm{g}_{(\ell,\ell^{\prime},k,k^{\prime})}$ should be treated as a
matrix with row index $(\ell,k)$ and column index
$(\ell^{\prime},k^{\prime})$. Thus,
$[\bm{g}(\omega-p)\bm{g}(\omega+q)]_{(\ell,\ell^{\prime\prime},k,k^{\prime\prime})}=\sum_{(\ell^{\prime},k^{\prime})}\bm{g}(\omega-p)_{(\ell,\ell^{\prime},k,k^{\prime})}\bm{g}(\omega+q)_{(\ell^{\prime},\ell^{\prime\prime},k^{\prime},k^{\prime\prime})}$.
The $\bm{g}(\omega)$ matrix captures the effect of the bath on the system, and
can be evaluated to be

$\bm{g}_{(\ell,\ell^{\prime},k,k^{\prime})}(\omega)=\delta_{\ell,\ell^{\prime}}\begin{cases*}\frac{\sqrt{\mathfrak{J}_{\ell}(-\omega)n_{\ell}(-\omega)}}{4\sqrt{2}\pi}\begin{pmatrix}1&i\\\
-i&1\end{pmatrix}_{k,k^{\prime}}& if $\omega<0$\\\
\frac{\sqrt{\mathfrak{J}_{\ell}(\omega)(1+n_{\ell}(\omega))}}{4\sqrt{2}\pi}\begin{pmatrix}1&-i\\\
i&1\end{pmatrix}_{k,k^{\prime}}& if $\omega>0$\end{cases*}$ (69)

As in the Redfield case, $\ket{E_{\alpha}}$, $\ket{E_{\gamma}}$ are eigenkets
of the system Hamiltonian with eigenenergies $E_{\alpha}$ and $E_{\gamma}$,
and $E_{\gamma\alpha}=E_{\gamma}-E_{\alpha}$.

## 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).
  * 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).
  * 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).
  * Gorini  _et al._ (1976) V. Gorini, A. Kossakowski, and E. C. G. Sudarshan, Journal of Mathematical Physics 17, 821 (1976).
  * Lindblad (1976) G. Lindblad, Communications in Mathematical Physics 48, 119 (1976).
  * Breuer and Petruccione (2006) H.-P. Breuer and F. Petruccione,  _The Theory of Open Quantum Systems_ (Oxford University Press, Oxford, 2006).
  * Carmichael (2002) H. Carmichael,  _Statistical Methods in Quantum Optics 1. Master Equations and Fokker–Planck Equations_ (Springer-Verlag Berlin Heidelberg, 2002).
  * Rivas and Huelga (2012) Á. Rivas and S. F. Huelga, SpringerBriefs in Physics (2012).
  * Weimer  _et al._ (2021) H. Weimer, A. Kshetrimayum, and R. Orús, Rev. Mod. Phys. 93, 015008 (2021).
  * Plenio and Knight (1998) M. B. Plenio and P. L. Knight, Rev. Mod. Phys. 70, 101 (1998).
  * Dalibard  _et al._ (1992) J. Dalibard, Y. Castin, and K. Mølmer, Phys. Rev. Lett. 68, 580 (1992).
  * Mølmer  _et al._ (1993) K. Mølmer, Y. Castin, and J. Dalibard, Journal of the Optical Society of America B 10, 524 (1993).
  * Redfield (1965) A. Redfield, in  _Advances in Magnetic Resonance_ (Elsevier, 1965) pp. 1–32.
  * Hartmann and Strunz (2020) R. Hartmann and W. T. Strunz, Phys. Rev. A 101, 012103 (2020).
  * Eastham  _et al._ (2016) P. R. Eastham, P. Kirton, H. M. Cammack, B. W. Lovett, and J. Keeling, Phys. Rev. A 94, 012110 (2016).
  * 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).
  * Walls (1970) D. F. Walls, Zeitschrift für Physik A Hadrons and nuclei 234, 231 (1970).
  * Wichterich  _et al._ (2007) H. Wichterich, M. J. Henrich, H.-P. Breuer, J. Gemmer, and M. Michel, Phys. Rev. E 76, 031115 (2007).
  * 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).
  * 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).
  * Benatti  _et al._ (2020) F. Benatti, R. Floreanini, and L. Memarzadeh, Phys. Rev. A 102, 042219 (2020).
  * Konopik and Lutz (2020) M. Konopik and E. Lutz, (2020), arXiv:2012.09907 [quant-ph] .
  * 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, (2021), arXiv:2102.10036 [quant-ph] .
  * Trushechkin (2021) A. Trushechkin, Phys. Rev. A 103, 062226 (2021).
  * Davidovic (2021) D. Davidovic,  (2021), arXiv:2112.07863 [quant-ph] .
  * 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).
  * Nielsen and Chuang (2011) M. A. Nielsen and I. L. Chuang,  _Quantum Computation and Quantum Information: 10th Anniversary Edition_ , 10th ed. (Cambridge University Press, USA, 2011).
  * Bennett and DiVincenzo (2000) C. H. Bennett and D. P. DiVincenzo, Nature 404, 247 (2000).
  * Streltsov  _et al._ (2017) A. Streltsov, G. Adesso, and M. B. Plenio, Rev. Mod. Phys. 89, 041003 (2017).
  * Lostaglio  _et al._ (2015) M. Lostaglio, K. Korzekwa, D. Jennings, and T. Rudolph, Phys. Rev. X 5, 021001 (2015).
  * Narasimhachar and Gour (2015) V. Narasimhachar and G. Gour, Nature Communications 6 (2015).
  * Mitchison  _et al._ (2015) M. T. Mitchison, M. P. Woods, J. Prior, and M. Huber, New Journal of Physics 17, 115013 (2015).
  * Allahverdyan  _et al._ (2004) A. E. Allahverdyan, R. Balian, and T. M. Nieuwenhuizen, Europhysics Letters (EPL) 67, 565–571 (2004).
  * Korzekwa  _et al._ (2016) K. Korzekwa, M. Lostaglio, J. Oppenheim, and D. Jennings, New Journal of Physics 18, 023045 (2016).
  * Kammerlander and Anders (2016) P. Kammerlander and J. Anders, Scientific Reports 6 (2016).
  * Santos  _et al._ (2019) J. P. Santos, L. C. Céleri, G. T. Landi, and M. Paternostro, npj Quantum Information 5 (2019).
  * Francica  _et al._ (2019) G. Francica, J. Goold, and F. Plastina, Phys. Rev. E 99, 042105 (2019).
  * Fleming and Cummings (2011) C. H. Fleming and N. I. Cummings, Phys. Rev. E 83, 031117 (2011).
  * Purkayastha  _et al._ (2020) A. Purkayastha, G. Guarnieri, M. T. Mitchison, R. Filip, and J. Goold, npj Quantum Information 6 (2020).
  * Johansson  _et al._ (2012) J. Johansson, P. Nation, and F. Nori, Computer Physics Communications 183, 1760 (2012).
  * Johansson  _et al._ (2013) J. Johansson, P. Nation, and F. Nori, Computer Physics Communications 184, 1234 (2013).
  * Guarnieri  _et al._ (2018) G. Guarnieri, M. Kolář, and R. Filip, Phys. Rev. Lett. 121, 070401 (2018).
  * 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).
  * Cattaneo  _et al._ (2021) M. Cattaneo, G. De Chiara, S. Maniscalco, R. Zambrini, and G. L. Giorgi, Phys. Rev. Lett. 126, 130403 (2021).
  * Brenes  _et al._ (2020) M. Brenes, J. J. Mendoza-Arenas, A. Purkayastha, M. T. Mitchison, S. R. Clark, and J. Goold, Phys. Rev. X 10, 031040 (2020).